Fortran子程序在Python中返回None但R中正常运行的问题求助

问题场景 #

我在实现二元正态CDF的概率计算，R环境下使用pbivnorm包的.Fortran调用能正常返回结果：

library(pbivnorm)
.Fortran("PBIVNORM", as.double(0), c(0,0), as.double(-0.1), as.double(-0.2), as.integer(c(0,0)), as.double(0), as.integer(1), PACKAGE="pbivnorm")

返回结果：

[[1]]
[1] 0.193613

[[2]]
[1] 0 0

[[3]]
[1] -0.1

[[4]]
[1] -0.2

[[5]]
[1] 0 0

[[6]]
[1] 0

[[7]]
[1] 1

提取该包底层Fortran代码，用F2py编译为fortran.so后，Python调用代码如下：

import fortran
fortran.pbivnorm(float(0), [float(0), float(0)], float(-0.1), float(-0.2), [int(0), int(0)], float(0),int(1))

但调用返回None，已知该Fortran子程序会通过引用修改传入值，尝试使用命名常量后值仍未变化，怀疑是输入类型或传递方式的问题。

可行解决方案 #

1. 用Numpy数组包装可修改参数 #

Fortran子程序依赖引用传递修改变量，Python原生标量（如float(0)）是不可变类型，无法被Fortran修改。需用Numpy数组包装所有参数，确保传递可修改的内存缓冲区：

import numpy as np
import fortran

# 严格匹配Fortran的变量类型：double precision对应np.float64，integer对应np.int32
p = np.array(0.0, dtype=np.float64)
xy = np.array([0.0, 0.0], dtype=np.float64)
a = np.array(-0.1, dtype=np.float64)
b = np.array(-0.2, dtype=np.float64)
infin = np.array([0, 0], dtype=np.int32)
maxpts = np.array(0.0, dtype=np.float64)
inform = np.array(1, dtype=np.int32)

# 调用子程序
fortran.pbivnorm(p, xy, a, b, infin, maxpts, inform)

# 查看计算结果（对应R中返回的第一个值）
print(p)  # 预期输出: 0.193613

2. 明确F2py编译的参数意图 #

F2py自动生成的接口可能未正确识别参数的输入/输出属性，需手动指定参数意图。可以在Fortran代码中添加注释，或修改.pyf接口文件：
在Fortran代码开头添加：

subroutine PBIVNORM(p, xy, a, b, infin, maxpts, inform)
    double precision p, xy(2), a, b, maxpts
    integer infin(2), inform
    !f2py intent(inout) p, xy, infin, maxpts, inform
    !f2py intent(in) a, b
end subroutine PBIVNORM

然后重新编译：

f2py -c -m fortran your_fortran_code.f90

这样F2py会明确哪些参数是输入、哪些是输入输出，确保Python端能正确获取修改后的值。

3. 严格匹配变量类型 #

Fortran的double precision必须对应Numpy的np.float64，integer对应np.int32（或np.int64，取决于Fortran编译器默认整数长度），类型不兼容会导致内存错误或参数无法被修改。

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