python - Fortran 版本的 python 循环给出了错误的答案

Question

我一直在研究国际贸易模型，但该模型变得非常慢(有时需要几周时间才能完成)。大多数情况下，有一个大的 for 循环减慢了进程，如下所示(在 Python 中):

for s in xrange(self.S-1):
                c_matrix[:,s+1,s+1:self.T+s+1] = ((self.beta * (1-self.MortalityRates[:,s,s:self.T+s]) * (1 + r_path[s+1:self.T+s+1] - self.delta)\
                                                * psi[:,s+1,s+1:self.T+s+1])/psi[:,s,s:self.T+s])**(1/self.sigma) * c_matrix[:,s,s:self.T+s]*np.exp(-self.g_A)

                a_matrix[:,s+1,s+1:self.T+s+1] = (  (we[:,s,s:self.T+s] + (1 + r_path[s:self.T+s] - self.delta)*a_matrix[:,s,s:self.T+s] + bqvec_path[:,s,s:self.T+s])\
                                                -c_matrix[:,s,s:self.T+s]*(1+we[:,s,s:self.T+s]*(self.chi/we[:,s,s:self.T+s])**self.rho)  )*np.exp(-self.g_A)

我不会详细介绍每个部分，但我决定制作一个 Fortran 模块，理论上它应该做完全相同的事情，但运行速度应该明显更快。我使用 Cython 和 distutils 来包装 Fortran 模块。这是 .f90 文件(特别关注上面 for 循环的 a_matrix 部分)供您引用:

MODULE FILLA

USE ISO_C_BINDING, ONLY: C_DOUBLE, C_INT

IMPLICIT NONE

CONTAINS

SUBROUTINE      

FILLASSETS(A_MATRIX,WE,R_PATH,BQ_VECPATH,C_MATRIX,I,S,T,DELTA,CHI,RHO,G_A) BIND(C)

  INTEGER(C_INT), INTENT(IN) :: I,S,T
  INTEGER :: X, E, L
  REAL(C_DOUBLE), INTENT(IN) :: DELTA, CHI, RHO, G_A
  REAL(C_DOUBLE), DIMENSION(0:I-1,0:S-1,0:T+S-1),INTENT(IN) :: C_MATRIX, WE, BQ_VECPATH
  REAL(C_DOUBLE), DIMENSION(0:T+S-1), INTENT(IN) :: R_PATH
  REAL(C_DOUBLE), DIMENSION(0:I-1,0:S,0:T+S-1), INTENT(INOUT) :: A_MATRIX
  REAL(C_DOUBLE), DIMENSION(0:I-1,0:S-1,0:T+S-1):: R_PATHM

  DO E = 0, I-1
    DO L = 0, S-1
        R_PATHM(E,L,:)=R_PATH(:)
    ENDDO
  ENDDO

  DO X = 0, S-2
        A_MATRIX(:,X+1,X+1:T+X+1) = ((WE(:,X,X:T+X)+(1.0D0+R_PATHM(:,X,X:T+X)-DELTA)&
                                    &*A_MATRIX(:,X,X:T+X)+BQ_VECPATH(:,X,X:T+X))&
                                    &-C_MATRIX(:,X,X:T+X)*(1.0D0+WE(:,X,X:T+X)*(CHI/WE(:,X,X:T+X))&
                                    &**RHO))*EXP(-G_A)
  ENDDO
END SUBROUTINE FILLASSETS

END MODULE

这个 .pyx 文件是如何将 Python 中的数组传递到 Fortran 中的:

from numpy cimport ndarray as ar
cdef extern from "filla.h":
        void fillassets(double* A_Matrix, double* we, double* r_path, double* bq_vecpath, double* c_matrix, int* I, int* S, int* T, double* delta, double* chi, double* rho, double* g_A)

cpdef f_filla(ar[double, ndim=3] A_Mat, ar[double, ndim=3] w_e, ar[double, ndim=1] r_path, ar[double, ndim=3] bq_vecpath, ar[double,ndim=3] c_mat, double delta, double chi, double rho, double g_A):
        cdef int I,S,T,Q
        if A_Mat.flags['C_CONTIGUOUS'] and w_e.flags['C_CONTIGUOUS'] and r_path.flags['C_CONTIGUOUS'] and bq_vecpath.flags['C_CONTIGUOUS'] and c_mat.flags['C_CONTIGUOUS']:
                I=c_mat.shape[0]
                S=c_mat.shape[1]
                Q=c_mat.shape[2]
                T=Q-S
                fillassets(&A_Mat[0,0,0],&w_e[0,0,0],&r_path[0],&bq_vecpath[0,0,0],&c_mat[0,0,0],&I,&S,&T,&delta,&chi,&rho,&g_A)
        else:
                raise ValueError("Input array U is not C-contiguous")

所以我的设置没有任何错误，但是当我并排运行它们时，Fortran 模块由于某种原因在 Assets 矩阵中返回不正确的值。我认为这与 Fortran 和 Python 如何处理数组有关，但我还没有找到简单的解释，希望我能在这里找到一些帮助。我对 Fortran 非常不熟悉，所以我可能缺少一些简单的东西。我尽量保持简短，所以如果我遗漏了任何细节，请告诉我。

编辑:我正在使用 Mac OS X El Capitan 和 gfortran 编译器。

Answer 1

Fortran 使用相反的数组元素布局:如果有一个维度为 (m,n) 的数组，则元素将逐列布局，即第二个维度变化最慢。 C(也许还有Python？)则相反，逐行布局，因此第一个维度变化最慢。因此，如果您有一个大小为 (m,n) 的 C 矩阵并将其传递给 Fortran，则需要告诉它维度为 (n,m)。