F2py¶
F2py是一个非常好的包，它自动包装Fortran代码，并使其可以从Python调用。斐波纳契数例取自f2py网页http://cens.ioc.ee/projects/f2py2e/.
从笔记本中，神奇的%Fortran将自动编译单元格中的任何Fortran代码，所有子例程将成为可调用的函数(尽管名称将转换为小写)。例如，将以下内容粘贴到一个单元格中。空格是正确的，这一点很重要，因为默认情况下，代码被视为固定格式的Fortran，如果不在正确的列中，编译器将发出警告。为了避免这种情况，您可以编写Fortran 90代码，方法是将第一行！f90。稍后将会有一个这样的例子。
%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是一个可以在任何字符串上调用的函数。所以，如果你有一个Fortran程序在我的程序文件中。然后，您可以执行以下操作
f=open('my_prog.f','r')
s=f.read()
fortran(s)
现在，我的pro.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的顺序。Fortr.Librters也只是一个链接的库的名称列表。您可以直接设置此列表。因此fortr.Librarys= ['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中。
fortran.add_library_path('path').
您也可以通过赋值直接设置fortr.库路径。它应该是要传递给GCC的路径(字符串)列表。为了让您了解可以使用f2py做的更多事情，请注意，使用Intent语句可以控制生成的Python函数的行为方式。例如，考虑对我们的原始斐波纳契代码进行以下修改。
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)都将在一个元组中返回。请注意，声明某个意图(在中)意味着您只关心它在函数被调用之前的值，而不是在调用后。所以在上面的n告诉我们要计算多少fiboncci数，我们需要将其指定为输入，但是我们不需要返回n，因为它不包含任何新的东西。类似地，A是意图(Out)，所以我们不需要A在事前有一个特定值，我们只关心事后的内容。F2py生成一个Python函数，因此您只传递那些声明的Intent(In)，并为其余参数提供空的工作区，并且它只返回Intent的那些(Out)。默认情况下，所有参数都是有意的(在中)。
现在考虑以下几点
%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
你可能期待着renale(a，n)重缩放一个数量级的矩阵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\) 在子例程的末尾，并在Numy数组中返回它们。这可能仍不是你想要的。原版 \(a\) 您传入的不会被修改。如果要修改原始的 \(a\) 您传入了使用意图(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
(In，Out，Overwrite)意图是如果 \(a\) 是FORTRAN排序的，我们工作在适当的位置，如果不是，我们复制它，然后返回内容。这在某种程度上是两个世界的最佳选择。请注意，如果要反复向Fortran代码传递大的数值数组，避免使用(INOUT)或(IN，OUT，OVERWRITE)复制数组非常重要。但是请记住，您的NumPy数组必须使用Fortran排序以避免复制。
有关F2py的更多示例和更高级用法，请参考f2py网页http://cens.ioc.ee/projects/f2py2e/.可以使用以下命令从SageShell调用f2py文档中引用的命令行f2py工具
!f2py