调用f2py封装的radar5后Python代码停止执行问题排查

已通过f2py将Fortran的RADAR5延迟微分方程积分器封装为Python模块，调用后积分结果与官方示例一致，但调用后的所有Python代码均无法执行，交互式环境下直接导致内核崩溃。未调用RADAR5时代码运行完全正常。

相关代码与编译命令 #

Python调用代码 #

import numpy
import radar5_f2pytest
print(radar5_f2pytest.__doc__)
i=[2]
x=[0]
rpar=numpy.array([1.34, 1.6E9, 8.0E3, 4.0E7, 1, 1, 6E-2, 6E-2, 15E-2],order='F')
ND=2
NRDENS=1
NGRID=1
NLAGS=1
NJACL=2
MXST=4000
LWORK=11
LIWORK=16
IJAC=1
MLJAC=ND
IMAS=0
IOUT=1
X=0
Y=numpy.array([1E-10,1E-5],order='F')
TAU=rpar[8]
XEND=100.5
ITOL=0
RTOL=1E-9
ATOL=RTOL*1E-9
H=1E-6
IWORK=numpy.zeros((16), order='F', dtype=numpy.int64)

WORK=numpy.array([1E-16, 0.9, 0.001, min(0.03, RTOL**0.5), 1, 1.2, XEND-X, 0.2, 8, 0, 5], order='F')
IWORK[1]=1000000
IWORK[4]=ND
IWORK[2]=7
IWORK[8]=1
IWORK[10]=2
IWORK[11]=MXST
IWORK[12]=NGRID
IWORK[13]=1
IWORK[14]=NRDENS
IPAST=numpy.zeros((NRDENS+1), order='F')
IPAST[0]=2

GRID=numpy.array([TAU], order='F')

ipar=numpy.zeros((1))
print("WHAT")
radar5_f2pytest.radar5_f2py(X, Y, XEND, H, 
                                  RTOL, ATOL, ITOL, IJAC, 
                                  MLJAC, NLAGS, NJACL,
                                  IMAS,ND,0, IOUT, WORK, 
                                  IWORK,GRID, IPAST, rpar, ipar)

print("WHAT")  # 该行代码从未执行

Fortran封装子例程 #

! -*- f90 -*-

SUBROUTINE radar5_f2py(N, X, Y, XEND, H, RTOL, ATOL, ITOL,&
&                          IJAC, MLJAC, MUJAC, NLAGS, NJACL, &
&                          IMAS, MLMAS, MUMAS, IOUT, WORK, IWORK, &
&                          GRID, IPAST, RPAR, IPAR, IDID, RESULT)
   
    IMPLICIT NONE
              
    INTEGER, PARAMETER :: DP=kind(1D0) 
    INTEGER, intent(in) :: N, NLAGS, NJACL
    REAL(kind=8), intent(inout) :: X
    REAL(kind=8), intent(in) :: XEND
    REAL(kind=8), intent(in) :: H
    LOGICAL, intent(in) :: ITOL, IJAC, IMAS, IOUT
    INTEGER, intent(in) :: MLJAC, MLMAS, MUMAS

    INTEGER, intent(in) :: MUJAC
    !f2py integer, optional, intent(in) :: MUJAC


    REAL(kind=8), dimension(N), intent(inout) :: Y 
    REAL(kind=8), intent(in), dimension(11) :: WORK 
    REAL(kind=8), intent(in) ::  ATOL,RTOL 
    INTEGER, intent(in), dimension(16) :: IWORK 
    REAL(kind=8), intent(inout) :: GRID 
    INTEGER, intent(inout) :: IPAST 
    REAL(kind=8), intent(in), dimension(9) :: RPAR 
    INTEGER, intent(in) :: IPAR
    INTEGER, intent(out) :: IDID
    REAL*8, intent(out), dimension(1,N) :: RESULT
    REAL(kind=8), EXTERNAL :: PHI
    REAL(kind=8), EXTERNAL :: ARGLAG
    EXTERNAL :: FCN, JFCN, JACLAG, SOLOUT, DUMMY

    CALL RADAR5(N,FCN,PHI,ARGLAG,X,Y,XEND,H,&
     &                  RTOL,ATOL,ITOL, &
     &                  JFCN,IJAC,MLJAC,MUJAC, &
     &                  JACLAG,NLAGS,NJACL, &
     &                  IMAS,SOLOUT,IOUT, &
     &                  WORK,IWORK,RPAR,IPAR,IDID,&
     &                  GRID,IPAST,DUMMY,MLMAS,MUMAS,RESULT)
    
     RESULT(1,1)=1
     RETURN   

END SUBROUTINE

编译命令 #

flang -c radar5.f contr5.f dc_decdel.f decsol.f dontr5.f 
================================================================
start f2py...
f2py -c -m radar5_f2pytest radar5_template.f90 radar5.o contr5.o dc_decdel.o decsol.o dontr5.o only: radar5_f2py :

问题排查与解决方案 #

1. 修正参数Intent声明 #

RADAR5积分器会修改WORK和IWORK数组的内容，但当前Fortran封装中这两个参数被声明为intent(in)，导致内存写入越界：

! 修改前
REAL(kind=8), intent(in), dimension(11) :: WORK 
INTEGER, intent(in), dimension(16) :: IWORK 

! 修改后
REAL(kind=8), intent(inout), dimension(11) :: WORK 
INTEGER, intent(inout), dimension(16) :: IWORK 

2. 统一参数类型匹配 #

Python传入的ITOL、IJAC等是整数0/1，但Fortran中声明为LOGICAL，存在类型转换风险，改为整数类型更稳妥：

! 修改前
LOGICAL, intent(in) :: ITOL, IJAC, IMAS, IOUT

! 修改后
INTEGER, intent(in) :: ITOL, IJAC, IMAS, IOUT

同时，Python中ipar是长度为1的数组，但Fortran中是标量，调整Fortran声明为数组：

! 修改前
INTEGER, intent(in) :: IPAR

! 修改后
INTEGER, intent(in), dimension(1) :: IPAR

3. 接收RADAR5的输出参数 #

当前Python调用未接收IDID（积分状态码）和RESULT，无法确认积分是否正常退出，修改调用代码：

idid, result = radar5_f2pytest.radar5_f2py(X, Y, XEND, H, 
                                  RTOL, ATOL, ITOL, IJAC, 
                                  MLJAC, NLAGS, NJACL,
                                  IMAS,ND,0, IOUT, WORK, 
                                  IWORK,GRID, IPAST, rpar, ipar)
print(f"积分状态码IDID: {idid}")

根据RADAR5文档，IDID=1表示成功完成积分，负数则对应错误类型，可快速定位问题。

4. 解决栈溢出问题 #

RADAR5可能占用较大栈空间，导致交互式环境内核崩溃，编译时增加栈大小参数：

flang -c -Wl,-stack_size,0x10000000 radar5.f contr5.f dc_decdel.f decsol.f dontr5.f 

5. 调试崩溃位置 #

若问题仍存在，用gdb定位崩溃点：

gdb python
run your_script.py
bt  # 打印崩溃调用栈

内容的提问来源于stack exchange，提问作者Karsten