f2py#
F2py is a very nice package that automatically wraps fortran code and makes it callable from Python. The Fibonacci examples are taken from the f2py webpage http://cens.ioc.ee/projects/f2py2e/.
From the notebook the magic %fortran will automatically compile any fortran code in a cell and all the subroutines will become callable functions (though the names will be converted to lowercase.) As an example paste the following into a cell. It is important that the spacing is correct as by default the code is treated as fixed format fortran and the compiler will complain if things are not in the correct column. To avoid this, you can write fortran 90 code instead by making your first line !f90. There will be an example of this later.
%fortran
C FILE: FIB1.F
SUBROUTINE FIB(A,N)
C
C CALCULATE FIRST N FIBONACCI NUMBERS
C
INTEGER N
REAL*8 A(N)
DO I=1,N
IF (I.EQ.1) THEN
A(I) = 0.0D0
ELSEIF (I.EQ.2) THEN
A(I) = 1.0D0
ELSE
A(I) = A(I-1) + A(I-2)
ENDIF
ENDDO
END
C END FILE FIB1.F
Now evaluate it. It will be automatically compiled and imported into Sage (though the name of imported function will be lowercase). Now we want to try to call it, we need to somehow pass it an array \(A\), and the length of the array \(N\). The way it works is that numpy arrays will be automatically converted to fortran arrays, and Python scalars converted to fortran scalars. So to call fib we do the following.
import numpy
m=numpy.array([0]*10,dtype=float)
print(m)
fib(m,10)
print(m)
Note that fortran is a function that can be called on any string. So if you have a fortran program in a file my prog.f. Then you could do the following
f=open('my_prog.f','r')
s=f.read()
fortran(s)
Now all the functions in my prog.f are callable.
It is possible to call external libraries in your fortran code. You simply need to tell f2py to link them in. For example suppose we wish to write a program to solve a linear equation using lapack (a linear algebra library). The function we want to use is called dgesv and it has the following signature.
SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
* N (input) INTEGER
* The number of linear equations, i.e., the order of the
* matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
* On entry, the N-by-N coefficient matrix A.
* On exit, the factors L and U from the factorization
* A = P*L*U; the unit diagonal elements of L are not stored.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* IPIV (output) INTEGER array, dimension (N)
* The pivot indices that define the permutation matrix P;
* row i of the matrix was interchanged with row IPIV(i).
*
* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
* On entry, the N-by-NRHS matrix of right hand side matrix B.
* On exit, if INFO = 0, the N-by-NRHS solution matrix X.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, U(i,i) is exactly zero. The factorization
* has been completed, but the factor U is exactly
* singular, so the solution could not be computed.
we could do the following. Note that the order that library are in the list actually matters as it is the order in which they are passed to gcc. Also fortran.libraries is simply a list of names of libraries that are linked in. You can just directly set this list. So that fortran.libraries=[‘lapack’,’blas’]is equivalent to the following.
fortran.add_library('lapack')
fortran.add_library('blas')
Now
%fortran
!f90
Subroutine LinearEquations(A,b,n)
Integer n
Real*8 A(n,n), b(n)
Integer i, j, pivot(n), ok
call DGESV(n, 1, A, n, pivot, b, n, ok)
end
There are a couple things to note about this. As we remarked earlier, if the first line of the code is !f90, then it will be treated as fortran 90 code and does not need to be in fixed format. To use the above try
a=numpy.random.randn(10,10)
b=numpy.array(range(10),dtype=float)
x=b.copy()
linearequations(a,x,10)
numpy.dot(a,x)
This will solve the linear system ax=b and store the result in b. If your library is not in Sage’s local/lib or in your path you can add it to the search path using
fortran.add_library_path('path').
You can also directly set fortran.library paths by assignment. It should be a list of paths (strings) to be passed to gcc. To give you an idea of some more things you can do with f2py, note that using intent statements you can control the way the resulting Python function behaves a bit bitter. For example consider the following modification of our original fibonacci code.
C FILE: FIB3.F
SUBROUTINE FIB(A,N)
C
C CALCULATE FIRST N FIBONACCI NUMBERS
C
INTEGER N
REAL*8 A(N)
Cf2py intent(in) n
Cf2py intent(out) a
Cf2py depend(n) a
DO I=1,N
IF (I.EQ.1) THEN
A(I) = 0.0D0
ELSEIF (I.EQ.2) THEN
A(I) = 1.0D0
ELSE
A(I) = A(I-1) + A(I-2)
ENDIF
ENDDO
END
C END FILE FIB3.F
Note the comments with the intent statements. This tells f2py that \(n\) is an input parameter and \(a\) is the output. This is called as
a=fib(10)
In general you will pass everything declared intent(in) to the fortran function and everything declared intent(out) will be returned in a tuple. Note that declaring something intent(in) means you only care about its value before the function is called not afterwards. So in the above n tells us how many fiboncci numbers to compute we need to specify this as an input, however we don’t need to get n back as it doesn’t contain anything new. Similarly A is intent(out) so we don’t need A to have an specific value beforehand, we just care about the contents afterwards. F2py generates a Python function so you only pass those declared intent(in) and supplies empty workspaces for the remaining arguments and it only returns those that are intent(out). All arguments are intent(in) by default.
Consider now the following
%fortran
Subroutine Rescale(a,b,n)
Implicit none
Integer n,i,j
Real*8 a(n,n), b
do i = 1,n
do j=1,n
a(i,j)=b*a(i,j)
end do
end do
end
You might be expecting Rescale(a,n) to rescale a numpy matrix a. Alas this doesn’t work. Anything you pass in is unchanged afterwards. Note that in the fibonacci example above, the one dimensional array was changed by the fortran code, similarly the one dimensional vector b was replaced by its solution in the example where we called lapack while the matrix A was not changed even then dgesv says it modifies the input matrix. Why does this not happen with the two dimensional array. Understanding this requires that you are aware of the difference between how fortran and C store arrays. Fortran stores a matrices using column ordering while C stores them using row ordering. That is the matrix
is stored as
\((0\, 1\, 2\, 3\, 4\, 5\,) \,\,\,\, \text{ in C}\)
\((0\, 3\,1\, 4\, 2\, 5) \,\,\,\, \text{ in Fortran}\)
One dimensional arrays are stored the same in C and Fortran. Because of this f2py allows the fortran code to operate on one dimensional vectors in place, so your fortran code will change one dimensional numpy arrays passed to it. However, since two dimensional arrays are different by default f2py copies the numpy array (which is stored in C format) into a second array that is in the fortran format (i.e. takes the transpose) and that is what is passed to the fortran function. We will see a way to get around this copying later. First let us point one way of writing the rescale function.
%fortran
Subroutine Rescale(a,b,n)
Implicit none
Integer n,i,j
Real*8 a(n,n), b
Cf2py intent(in,out) a
do i = 1,n
do j=1,n
a(i,j)=b*a(i,j)
end do
end do
end
Note that to call this you would use
b=rescale(a,2.0).
Note here I am not passing in \(n\) which is the dimension of \(a\). Often f2py can figure this out. This is a good time to mention that f2py automatically generates some documentation for the Python version of the function so you can check what you need to pass to it and what it will return. To use this try
rescale?
The intent(in,out) directives tells f2py to take the contents of \(a\) at the end of the subroutine and return them in a numpy array. This still may not be what you want. The original \(a\) that you pass in is unmodified. If you want to modify the original \(a\) that you passed in use intent(inout). This essentially lets your fortran code work with the data inplace.
%fortran
Subroutine Rescale(a,b,n)
Implicit none
Integer n,i,j
Real*8 a(n,n), b
Cf2py intent(inout) a
do i = 1,n
do j=1,n
a(i,j)=b*a(i,j)
end do
end do
end
If you wish to have fortran code work with numpy arrays in place, you should make sure that your numpy arrays are stored in fortran’s format. You can ensure this by using the order=’FORTRAN’ keyword when creating the arrays, as follows.
a=numpy.array([[1,2],[3,4]],dtype=float,order='FORTRAN')
rescale(a,2.0)
After this executes, a will have the rescaled version of itself. There is one final version which combines the previous two.
%fortran
Subroutine Rescale(a,b,n)
Implicit none
Integer n,i,j
Real*8 a(n,n), b
Cf2py intent(in,out,overwrite) a
do i = 1,n
do j=1,n
a(i,j)=b*a(i,j)
end do
end do
end
The (in,out,overwrite) intent says that if \(a\) is in FORTRAN ordering we work in place, however if its not we copy it and return the contents afterwards. This is sort of the best of both worlds. Note that if you are repeatedly passing large numpy arrays to fortran code, it is very important to avoiding copying the array by using (inout) or (in,out,overwrite). Remember though that your numpy array must use Fortran ordering to avoid the copying.
For more examples and more advanced usage of F2py you should refer to the f2py webpage http://cens.ioc.ee/projects/f2py2e/. The command line f2py tool which is referred to in the f2py documentation can be called from the Sage shell using
!f2py