FB040510
General
We are using Fortran 90/95 subroutines, wrapped to be called from python using the f2py package. This is because f90 code executes much faster than python scripts. There are some subtelties to pay attention to when wrapping fortran code, else you will add large overheads from the py-f90 interface, as arrays are copied and reindexed. For an introduction to F90/95 (only minor differences between the two), I recommend the compact and rather comprehensive (and free!) ``Fortran 90 course notes'' by AC Marshall from the University of Liverpool. It contains all you probably need to know. I wrote a simple fortran method in BoA/fortran/BoaTest1.f90 to illustrate some basic features and give you a chance to test the wrapper without BoA. Look at its header for details. For an online F90/95 language reference the best I found is at the NCSA resources page, describing IBM's XL Fortran for AIX 8.1 - which is close to the Intel compiler.
F90 in BoA
For BoA our general idea is to have one f90.so extension module, which includes all the f90 methods (called subroutines and functions in fortran). This is necessitated by that the f90.data module, which contains much of a scans data, is connected (through an ``use data'') to the other f90 program modules, and therefore they all need to be linked together.
The f90 methods may be split into different modules (classes) for convenience. We now have the first operational modules BoaF1.f90, BoaChannelAnalyser.f90, BoaBaseLine.f90, and the data module BoaData.f90. Each module may include any number of subroutines or functions. The data module BoaData.f90 is like a common block that contains all the data which does not change during data reduction. All data which does change is passed to the fortran subroutines as call arguments.
The BoaData.f90 (f90.data from python) module is filled in BoaDataEntity.FillF90. It must be refilled if you change data object, else the fortran methods will work on a different scan. This re-filling must be implemented still. Currently the f90.data is only filled upon read of a new data file.
The CVS directory BoA/fortran contains the fortran source code. You will need to wrap/compile the BoA modules on your local system (see below), since it links to local libraries that have no standard address. This will create the extension module f90.so which you import to BoA. From the CVS directory BoA start BoA, then
>>> from fortran import f90This is how to import any module from a subdirectory, which for this needs to include an empty file __init__.py
The python script fortran/ftest.py contains a series of calls to the fortran subroutines. To run it:
>>> read() # read in some scan
>>> op() # open plot device
[enter]
>>> execfile('ftest.py') # start the script
which is followed with lots of output.
To illustrate the use of new python methods that use fortran, you find
BoA/TestFB.py, which you run like ftest.py. It goes through a number
of data reduction steps and plots the data.
Wrapping F90 code with f2py
To wrap the f90 modules to produce f90.so:
ifc -c -w svd.f90
f2py -c -m f90 BoaData.f90 BoaF1.f90 BoaChannelAnalyser.f90 BoaBaseLine.f90 svd.o
or on some installations alternatively:
f2py -c --fcompiler=intel -m f90 BoaData.f90 BoaF1.f90 BoaChannelAnalyser.f90 BoaBaseLine.f90 svd.o
The first command recompiles the svd.o.
On the f2py line there are some diagnostic options you may add if you
debug your code:
-DF2PY_REPORT_ATEXIT : gives time statistics upon exit from python. -DF2PY_REPORT_ON_ARRAY_COPY=1000 : reports when the f2py interface copies an array. -DNUMARRAY : must be used for numarray support. Default is Numeric.
If the wrapping fails, one of the following may be wrong:
1. You have not initiated the ifc compiler properly. In your shell initialization file (e.g. .cshrc for tcsh) you need
if (-e /opt/intel/compiler60/ia32/bin/ifcvars.sh) then
source /opt/intel/compiler60/ia32/bin/ifcvars.csh
endif
or something equivalent.
2. Your python path does not include the intel fortran compiler:
setenv PYTHONPATH ".:/opt/intel/compiler60/ia32/lib/:
/usr/local/lib/python2.3:
/usr/local/lib/python2.3/site-packages:
/home/bertoldi/bin:
/opt:
/usr/lib"
3. You use an old version of f2py.
<fortran> f2py -version 2.39.235_1644
Once you have successfully imported f90 in BoA, you can inquire about the use of a given method by typing
print f90.f1.NAME.__doc__
Fortran attributes are called f90.data.name_of_attribute. To inquire which ones are available:
boa> print f90.data.__doc__ el - 'f'-array(218) track_el - 'f'-array(218) ffcf_gain - 'f'-array(120) subscan_time - 'f'-array(4) az_p - 'f'-array(109,3) lst - 'f'-array(218) lon_p - 'f'-array(109,3) track_az - 'f'-array(218) lat - 'f'-array(218) az - 'f'-array(218) lat_p - 'f'-array(109,3) lst_p - 'f'-array(109,3) array_gain - 'f'-array(120) lon - 'f'-array(218) ffcf_cn - 'f'-array(120) ut_p - 'f'-array(109,3) nodding_sta - 'i'-array(218) subscan_index - 'i'-array(4) subscan_num - 'i'-array(4) weights - 'f'-array(0), not allocated el_p - 'f'-array(109,3) ut - 'f'-array(218) wobbler_pos - 'f'-array(218)They are filled in in BoaBusiness.py: BoaB.FillF90
Use f90 methods in BoA
To call a fortran method, here an example:
compressed_array,nmax = f90.f1.compress(array,flag_array,0)Two objects are returned as a tuple, an array and an integer. They both are not in the call argument list, they are hidden to python, but are listed in the f90 code call argument list - have a look at the source code.
Limitations
This particular example illustrates one of the limitations of wrapping f90 code: you cannot return an array with a length that is determined upon execution. The wrapper needs to specify the size of an array somehow. It does not have to be fixed, but specified through the size of an input attribute at least. In this example we try to return an array that is a compression of the input array, determined by the condition that the corresponding flag is 0. The trick to still do this here is to return a comressed_array with the same size as array, plus an integer telling the size of the compressed array, so that the final answer is compressed_array[0:nmax].
Fortran vs. C-contiguous
If a Numeric array is proper-contiguous and has a proper type then it is directly passed to the wrapped Fortran function. Otherwise, an element-wise copy of an input array is made and the copy, being proper-contiguous and with proper type, is used as an array argument. There are two types of proper-contiguous Numeric arrays: Fortran-contiguous arrays when data is stored column-wise, i.e. indexing of data as stored in memory starts from the lowest dimension; C-contiguous when data is stored row-wise, i.e. indexing of data as stored in memory starts from the highest dimension. For one-dimensional arrays these notions coincide. To transform input arrays to column major storage order before passing them to Fortran routines, one may use the function as_column_major_storage(<array>) that is provided by all F2PY generated extension modules, such as the BoA f90. If you call a fortran method repeatedly with the same input array, you should convert the array first to avoid conversion by the wrapper interface on each call - which could dominate the execution time here. If you add the option -DF2PY_REPORT_ON_ARRAY_COPY=1000 when wrapping, you will be informed on each copy that the wrapper interface performs. The option -DF2PY_REPORT_ATEXIT gives an execution time summary upon exit that splits up the time used in fortran and in the interface. If the interface time is large or comparable to the fortran execution time, your code is not efficient because it copies arrays too often. Look at examples in BoaBaseLine.py, e.g.:
Data = f90.as_column_major_storage(self.Data.Data_Red_p)
Flag = f90.as_column_major_storage(self.Data.Data_Flag_p)
...
for i_ch in ch_range: # loop over channels and phases
for i_ph in ph_range:
Data = f90.baseline.addpoly(Data,Poly,Mean,Rms,i_ph,i_ch)
The input arrays are copied once into fortran-contiguous arrays before the loop,
so in the loop there is no overhead from copying. Note also the general scheme of
calling a fortran method here: Data is in- and output argument.