F2PY¶
F2py是一个非常好的包，它可以自动包装fortran代码，并使其可以从Python调用。斐波那契的例子取自f2py网页http://cens.ioc.ee/projects/f2py2e/。
在笔记本上，神奇的%fortran会自动编译单元格中的任何fortran代码，所有的子例程都将成为可调用函数（尽管名称将转换为小写）。例如，将以下内容粘贴到单元格中。间距是正确的，这一点很重要，因为默认情况下，代码被视为固定格式的fortran，如果列中的内容不正确，编译器会抱怨。为了避免这种情况，您可以编写fortran 90代码，方法是编写第一行代码！90层。稍后会有一个例子。
%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
现在评估一下。它将被自动编译并导入Sage（虽然导入函数的名称将是小写的）。现在我们要尝试调用它，我们需要以某种方式传递给它一个数组 \(A\) ，以及数组的长度 \(N\) . 其工作方式是numpy数组将自动转换为fortran数组，Python标量将自动转换为fortran标量。所以，要叫fib，我们要做以下事情。
import numpy
m=numpy.array([0]*10,dtype=float)
print(m)
fib(m,10)
print(m)
注意fortran是一个可以在任何字符串上调用的函数。因此，如果你在myprog.f文件中有一个fortran程序，那么你可以做以下的事情
f=open('my_prog.f','r')
s=f.read()
fortran(s)
现在我的prog.f中的所有函数都是可调用的。
可以在fortran代码中调用外部库。你只需要告诉f2py把它们连接起来。例如，假设我们想用lapack（一个线性代数库）编写一个求解线性方程的程序。我们要使用的函数名为dgesv，它具有以下签名。
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.
我们可以做以下事情。注意，库在列表中的顺序实际上很重要，因为它是传递给gcc的顺序。也fortran.库只是链接到中的库的名称列表。你可以直接设置这个列表。以便fortran.库= ['lapack'，'blas'] 相当于以下内容。
fortran.add_library('lapack')
fortran.add_library('blas')
现在
%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
有几件事需要注意。如前所述，如果代码的第一行是！f90，则它将被视为fortran 90代码，不需要采用固定格式。使用上述方法
a=numpy.random.randn(10,10)
b=numpy.array(range(10),dtype=float)
x=b.copy()
linearequations(a,x,10)
numpy.dot(a,x)
这将求解线性系统ax=b并将结果存储在b中。如果库不在Sage的local/lib或路径中，则可以使用将其添加到搜索路径中
fortran.add_library_path('path').
也可以直接设置fortran.库指定路径。它应该是要传递给gcc的路径（字符串）列表。为了让您了解f2py可以做的更多事情，请注意，使用intent语句，您可以控制生成的Python函数的行为方式。例如，考虑对原始fibonacci代码的以下修改。
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
注意注释和意向语句。这告诉f2py \(n\) 是一个输入参数，并且 \(a\) 是输出。这叫做
a=fib(10)
通常，您将把所有声明的intent（In）传递给fortran函数，所有声明的intent（out）都将以元组形式返回。请注意，声明intent（in）意味着您只关心它的值，而不是在调用函数之后。在上面的n告诉我们要计算多少个fibonci数，我们需要将它指定为输入，但是我们不需要返回n，因为它不包含任何新的内容。类似地，A是intent（out），所以我们不需要A预先有一个特定的值，我们只关心后面的内容。F2py生成一个Python函数，因此您只传递声明的intent（in），并为其余参数提供空的工作空间，它只返回intent（out）的那些。默认情况下，所有参数都是intent（in）。
现在考虑以下内容
%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
你可能期望重缩放（a，n）来重缩放numpy矩阵a。唉，这行不通。你传递的任何信息都不会改变。请注意，在上面的fibonacci示例中，一维数组被fortran代码更改，类似地，在我们调用lapack的示例中，一维向量b被它的解决方案替换，而矩阵A没有改变，即使当时dgesv说它修改了输入矩阵。为什么二维数组不会发生这种情况呢。要理解这一点，就需要了解fortran和C存储数组的方式之间的区别。Fortran使用列排序存储矩阵，而C使用行排序存储矩阵。这就是矩阵
存储为
\((0\, 1\, 2\, 3\, 4\, 5\,) \,\,\,\, \text{ in C}\)
\((0\, 3\,1\, 4\, 2\, 5) \,\,\,\, \text{ in Fortran}\)
一维数组在C和Fortran中是相同的。因为f2py允许fortran代码对一维向量进行操作，所以fortran代码将更改传递给它的一维numpy数组。但是，由于二维数组在默认情况下是不同的，f2py将numpy数组（以C格式存储）复制到第二个fortran格式的数组中（即，接受转置），这就是传递给fortran函数的内容。我们稍后会看到一种方法来避免这种复制。首先让我们指出一种编写重缩放函数的方法。
%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
请注意，要调用此函数，您将使用
b=rescale(a,2.0).
请注意，我不是路过这里 \(n\) 它是 \(a\) . 通常f2py可以解决这个问题。这是一个很好的时机，f2py会自动为函数的Python版本生成一些文档，这样您就可以检查需要传递给它的内容以及它将返回什么。使用这个尝试
rescale?
intent（in，out）指令告诉f2py获取 \(a\) 并以numpy数组的形式返回它们。这可能仍然不是你想要的。原来的 \(a\) 你传入的是未修改的。如果你想修改原稿 \(a\) 您在use intent（inout）中传递的。这基本上让fortran代码在适当的地方处理数据。
%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
如果希望fortran代码在适当的地方使用numpy数组，那么应该确保numpy数组是以fortran格式存储的。您可以通过在创建数组时使用order='FORTRAN'关键字来确保这一点，如下所示。
a=numpy.array([[1,2],[3,4]],dtype=float,order='FORTRAN')
rescale(a,2.0)
执行此操作后，a将拥有其自身的重新缩放版本。有一个最终版本结合了前两个版本。
%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
（进，出，过火）意图说如果 \(a\) 在FORTRAN命令中，我们在适当的地方工作，但是如果不是，我们复制它并在之后返回内容。这是两个世界中最好的。请注意，如果您反复向fortran代码传递大型numy数组，那么避免使用（inout）或（in，out，overwrite）复制数组是非常重要的。请记住，numpy数组必须使用Fortran排序来避免复制。
有关F2py的更多示例和更高级的用法，请参阅F2py网页http://cens.ioc.ee/projects/f2py2e/。f2py文档中提到的命令行f2py工具可以从sageshell调用，方法是使用
!f2py