f2py有更多有趣的例子¶
现在让我们看看使用f2py的一些更有趣的示例。我们将实现一个简单的迭代方法来求解拉普拉斯方程的平方。实际上，这个实现是从http://www.scipy.org/PerformancePython？突出显示scipy网站上的%28performance%29页面。它有很多关于在python中实现数值算法的信息。
下面的fortran代码实现了求解正方形拉普拉斯方程的松弛方法的一次迭代。
%fortran
subroutine timestep(u,n,m,dx,dy,error)
double precision u(n,m)
double precision dx,dy,dx2,dy2,dnr_inv,tmp,diff
integer n,m,i,j
cf2py intent(in) :: dx,dy
cf2py intent(in,out) :: u
cf2py intent(out) :: error
cf2py intent(hide) :: n,m
dx2 = dx*dx
dy2 = dy*dy
dnr_inv = 0.5d0 / (dx2+dy2)
error = 0d0
do 200,j=2,m-1
do 100,i=2,n-1
tmp = u(i,j)
u(i,j) = ((u(i-1,j) + u(i+1,j))*dy2+
& (u(i,j-1) + u(i,j+1))*dx2)*dnr_inv
diff = u(i,j) - tmp
error = error + diff*diff
100 continue
200 continue
error = sqrt(error)
end
如果你这么做的话
timestep?
你发现你需要传递timestep一个numpy数组u，网格间距dx，dy，它将返回更新后的u，以及一个错误估计。要将此应用于实际解决问题，请使用以下驱动程序代码
import numpy
j=numpy.complex(0,1)
num_points=50
u=numpy.zeros((num_points,num_points),dtype=float)
pi_c=float(pi)
x=numpy.r_[0.0:pi_c:num_points*j]
u[0,:]=numpy.sin(x)
u[num_points-1,:]=numpy.sin(x)
def solve_laplace(u,dx,dy):
iter =0
err = 2
while(iter <10000 and err>1e-6):
(u,err)=timestep(u,dx,dy)
iter+=1
return (u,err,iter)
现在调用例程
(sol,err,iter)=solve_laplace(u,pi_c/(num_points-1),pi_c/(num_points-1))
该方法求解的边界条件为：正方形的左、右边缘为正弦函数的半振荡。对于我们正在使用的51x51网格，我发现需要大约0.2秒来解决这个问题，需要2750次迭代。如果您已经安装了gnuplot包（使用可选的packages（）来查找它的名称，并使用install package来安装它），那么可以使用
import Gnuplot
g=Gnuplot.Gnuplot(persist=1)
g('set parametric')
g('set data style lines')
g('set hidden')
g('set contour base')
g('set zrange [-.2:1.2]')
data=Gnuplot.GridData(sol,x,x,binary=0)
g.splot(data)
为了了解使用f2py所取得的成果，让我们比较纯python中的相同算法和使用numpy数组的矢量化版本。
def slowTimeStep(u,dx,dy):
"""Takes a time step using straight forward Python loops."""
nx, ny = u.shape
dx2, dy2 = dx**2, dy**2
dnr_inv = 0.5/(dx2 + dy2)
err = 0.0
for i in range(1, nx-1):
for j in range(1, ny-1):
tmp = u[i,j]
u[i,j] = ((u[i-1, j] + u[i+1, j])*dy2 +
(u[i, j-1] + u[i, j+1])*dx2)*dnr_inv
diff = u[i,j] - tmp
err += diff*diff
return u,numpy.sqrt(err)
def numpyTimeStep(u,dx,dy):
dx2, dy2 = dx**2, dy**2
dnr_inv = 0.5/(dx2 + dy2)
u_old=u.copy()
# The actual iteration
u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
(u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv
v = (u - u_old).flat
return u,numpy.sqrt(numpy.dot(v,v))
您可以通过更改驱动程序例程中使用的timestep函数来尝试这些方法。即使在50x50网格上，python版本也很慢。在3000次迭代中求解系统需要70秒。在大约5000次迭代中，numpy例程需要2秒才能达到误差容限。相比之下，使用3000次迭代，f2py例程大约需要0.2秒才能达到误差容限。我要指出的是，numpy例程不是完全相同的算法，因为它是jacobi迭代，而f2py例程是gauss seidel。这就是为什么numpy版本需要更多的迭代。即便如此，你也可以看到f2py版本的速度似乎是numpy版本的5倍左右。实际上，如果你在一个500x500的网格上尝试这个，我发现numpy例程需要30秒来完成500次迭代，而f2py只需要2秒。所以f2py版本的速度要快15倍。在较小的网格上，每个实际的迭代都相对便宜，因此调用f2py的开销更为明显，在迭代开销较大的示例中，f2py的优势显而易见。即使是在小例子上，它仍然非常快。请注意，使用python编写的500x500网格需要大约半个小时来进行500次迭代。
据我所知，在matlab中实现这个算法最快的方法就是将它矢量化，就像我们的numpy例程一样。matlab和numpy中的向量加法具有可比性。所以除非有我不知道的诀窍，否则使用f2py你可以以交互方式编写比用matlab编写的任何代码快15倍的速度（如果我错了，请纠正我）。实际上，通过使用intent（in、out、overwrite）并使用order='FORTRAN'创建初始numpy数组，可以使f2py版本快一点。这就消除了内存中不必要的复制。