使用Python调用Fortran程序 1 前言 网上关于Python调用Fortran程序的方法通常分为三种：1）基于f2py；2）生成动态链接库；3）生成可执行文件。 其中第1种方法在涉及到“祖传”代码时，通常会出现各种报错；第3种方法在进行数据传递时基本只能通过操作文件的方式，很不方便。 由于涉及Fortran程序时，一般都逃不开“祖传”代码，因此本文将介绍最为稳定可靠的第2种方法。 2 方法详情 2.1 示例代码 新建test01.f90文件，创建子例程sub_test01以及函数func_test01，详细内容如下 subroutine sub_test01(x,y,z) bind(C,name="sub_test01") use iso_c_binding real(c_double), intent(in), value :: x,y real(c_double), intent(out) :: z(2) z(1) = x + x z(2) = y*y end subroutine sub_test01 function func_test01(x,y) result(z) bind(c,name="func_test01") use iso_c_binding real(c_double), intent(in), value :: x,y real(c_double) :: z z = x + y end function func_test01 bind：用于声明外部调用时子例程/函数名称 iso_c_binding：Fortran自带的模组，必须引用 intent：声明变量属性，输入为in，输出为out，即是输入也是输出为inout c_double：变量类型，real对应c_double，integer对应c_int value：输入变量为单个值时，需添加此标记 2.2 生成动态链接库 与正常编译相比增加-shared，生成后缀为.so的文件，如下 gfortran -shared test01.f90 -o test01.so 如果此步骤报错recompile with -fPIC，则在-shared后加上-fPIC。 2.3 使用Python调用 调用sub_test01子例程，需要引用ctypes和numpy，示例代码如下 import ctypes as ct import numpy as np # 加载动态链接库 fortlib = ct.CDLL('test01.so') # 引用sub_test01子例程 f_sub = fortlib.sub_test01 # 声明变量类型 f_sub.argtypes = [ct.c_double, ct.c_double, ct.POINTER(ct.c_double)] # 输入变量赋值 x = ct.c_double(3) y = ct.c_double(4) # 输出变量初始化 z = np.ones(2) z_p = z.ctypes.data_as(ct.POINTER(ct.c_double)) # 调用sub_test01子例程 f_sub(x,y,z_p) print(z) 使用ctypes.CDLL(<so name>)加载动态链接库，其中<so name>为上一节生成的动态链接库名称 子例程的引用名称为上一节bind中name定义的名称 argtypes用于声明变量类型，其中ctypes.POINTER表示指针。当变量为单个值（Fortran代码中value）时，声明为相应类型；当变量为数组（或输出变量，即intent(out)）时，声明为指针 ctypes.c_double(<value>)，其中<value>为变量的值 打印结果为[6. 16.] 调用func_test01函数，与子例程调用方式基本相同，示例代码如下 import ctypes as ct # 加载动态链接库 fortlib = ct.CDLL('test01.so') # 引用sub_test01子例程 f_sub = fortlib.func_test01 # 声明变量类型 f_sub.argtypes = [ct.c_double, ct.c_double] # 声明结果类型 f_sub.restype = ct.c_double # 输入变量赋值 x = ct.c_double(3) y = ct.c_double(4) # 调用sub_test01子例程 z = f_sub(x,y) print(z) restype声明返回值的类型 打印结果为7.0 参考 python调用fortran的3种形式【f2py，动态链接库，os命令】 How to Call Fortran from Python Using Python as glue

非线性最小二乘法拟合函数-3 在前面两个博客中，推导得到了Levenberg-Marquardt算法（简称LM算法）的迭代公式，这里将讲述如何使用Fortran编写一个简单的LM算法。 此处介绍待定系数只有一个的情况。 1 程序设计 1.1 Levenberg_Maquardt_Fit 简介：LM算法子例程，输入离散数据，迭代拟合待定系数 传递参数： n 离散数据的长度 x 离散数据的自变量 y 离散数据的因变量 a 待定系数 1.2 myfunc 简介：待拟合函数的数值计算子例程。为了具有普适性，假定待拟合函数的解析式未知，函数结果只能由此子例程数值计算得到。 传递参数： n 离散数据的长度 x 离散数据的自变量 a 待定系数的当前值 fx 待拟合函数的数值计算结果 1.3 Calculate_Jacobian 简介：计算待拟合函数的Jacobian矩阵的子例程。为了具有普适性，使用差分法计算一阶导数。 传递参数： n 离散数据的长度 x 离散数据的自变量 y 待拟合函数的数值计算结果 a 待定系数的当前值 J 待拟合函数的Jacobian矩阵 2 源代码 2.1 Levenberg_Maquardt_Fit 此处迭代次数上限设置为30次，可根据需求自行更改 对于自变量有多个的情况，可将输入x改为x1,x2,...，满足自变量个数即可 SUBROUTINE Levenberg_Maquardt_Fit(n, x, y, a) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 ! REAL, DIMENSION(3), INTENT(INOUT) :: a ! 待求的拟合系数 REAL, INTENT(INOUT) :: a ! 待求的拟合系数 ! 定义常量 INTEGER, PARAMETER :: m = 1 ! 待求的系数个数 REAL, PARAMETER :: eps = 1.0E-6 ! 收敛阈值 ! 定义变量 REAL :: da(m), r(n), J(n,m), H(m,m) ! 拟合系数、残差、雅可比矩阵、Hessian矩阵 REAL :: lambda, alpha(m) ! 调节因子、步长 INTEGER :: i, iter ! 循环计数器 REAL :: fx(n) ! 初始化调节因子 lambda = 0.001 ! 开始迭代 iter = 0 DO WHILE(iter < 30) ! 迭代次数上限为10000 iter = iter + 1 ! 计算被拟合函数值 CALL myfunc(n,x,a,fx) ! 计算残差向量 r = y - fx ! 计算雅可比矩阵 CALL Calculate_Jacobian(n, x, fx, a, J) ! 计算Hessian矩阵 H = MATMUL(TRANSPOSE(J), J) ! 计算梯度向量 da = MATMUL(TRANSPOSE(J), r) ! 计算搜索方向 DO i = 1, m H(i,i) = H(i,i) + lambda END DO ! CALL SOLVE_LINEAR_SYSTEM(m, H, da, alpha) alpha = da(1)/H(1,1) ! 计算新的拟合系数 a = a + alpha(1) ! 更新调节因子 IF (NORM2(r) < eps) THEN EXIT ELSE IF (NORM2(r) < NORM2(r - MATMUL(J, alpha))/2) THEN lambda = lambda/10. ELSE lambda = lambda*10. END IF END DO END SUBROUTINE 2.2 myfunc 此处使用了一个非常简单的一次函数作为示例，使用时需要将fx = a*x修改为实际的待拟合函数的形式。若待拟合函数是一个子例程，则在此处使用CALL调用 若在Levenberg_Maquardt_Fit子例程中修改了自变量个数，此处需要对应修改 SUBROUTINE myfunc(n,x,a,fx) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n) ! 自变量 REAL, INTENT(IN) :: a ! 拟合系数 REAL, INTENT(OUT) :: fx(n) fx = a*x END SUBROUTINE 2.3 Calculate_Jacobian 本子例程使用了差分法计算一阶导数，因此几乎适用于所有的待拟合函数 若在前面修改了自变量的个数，此处需对应修改 SUBROUTINE Calculate_Jacobian(n, x, y, a, J) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 REAL, INTENT(IN) :: a ! 拟合系数 ! 输出参数 REAL, DIMENSION(n,1), INTENT(OUT) :: J ! 雅可比矩阵 ! 定义变量 INTEGER :: i ! 循环计数器 REAL :: dela, y1(n) ! 步进 dela = 0.01 ! 计算 a+dela 对应的 y1 CALL myfunc(n,x,a+dela,y1) ! 计算雅可比 DO i = 1, n J(i,1) = (y1(i) - y(i))/dela END DO END SUBROUTINE 3 测试 编写程序主体，设定y约为x的两倍，调用LM子例程计算待定系数a program name implicit none INTEGER n real :: x(5),y(5),a x = (/1.,2.,3.,4.,5./) y = (/2.1,4.1,5.9,8.1,9.9/) ! y = 2*x**2 + n = 5 a = 1.0 call Levenberg_Maquardt_Fit(n, x, y, a) print *,a end program name 使用gfortran编译生成可执行文件，运行可执行文件后，在终端打印出 1.99818182 离散数据的因变量约为自变量的两倍，因此上述计算结果符合预期