python - 我可以在 fortran 中使用像 numpy 和 scipy 这样的库吗？

Question

我试图通过将我的一堆嵌套循环移植到 Fortran 并将它们作为子例程调用来加速我的 python 代码。

但是我的很多循环调用 numpy，以及来自 scipy 的特殊函数，比如贝塞尔函数。

在我尝试使用 fortran 之前，我想知道是否可以将 scipy 和 numpy 导入到我的 fortran 子例程中并为 bessel 函数调用模块？

否则我是否必须在 Fortran 中创建贝塞尔函数才能使用它？

理想情况下，我会创建某种子程序来优化下面的代码。这只是我整个项目的一个片段，目的是让您了解我要完成的工作。

我知道我应该实现其他做法来提高速度，但现在我正在调查在我的主 python 程序中调用 fortran 子例程的好处。

    for m in range(self.MaxNum_Eigen):
        #looping throught the eigenvalues for the given maximum number of eigenvalues allotted
        bm = self.beta[m]


        #not sure
        #*note: rprime = r. BUT tprime ~= t.

        #K is a list of 31 elements for this particular case

        K = (bm / math.sqrt( (self.H2**2) + (bm**2) ))*(math.sqrt(2) / self.b)*((scipy.special.jv(0, bm * self.r))/ (scipy.special.jv(0, bm * self.b))) # Kernel, K0(bm, r).

        #initial condition
        F = [37] * (self.n1)

        # Integral transform of the initial condition
        #Fbar = (np.trapz(self.r,self.r*K*F))

        '''
            matlab syntax trapz(X,Y), x ethier spacing or vector
            matlab:     trapz(r,r.*K.*F)                trapz(X,Y)
            python:     np.trapz(self.r*K*F, self.r)    trapz(Y,X)

        '''

        #*(np.trapz(self.r,self.r*K*F))
        Fbar = np.ones((self.n1,self.n2))*(np.trapz(self.r*K*F, self.r))

        #steady state condition: integral is in steady state
        SS = np.zeros((sz[0],sz[1]))

        coeff = 5000000*math.exp(-(10**3)) #defining value outside of loop with higher precision

        for i in range(sz[0]):
            for j in range(sz[1]):

                '''
                    matlab reshape(Array, size1, size2) takes multiple arguments the item its resizeing and the new desired shape

                    create self variables and so we are not re-initializing them over and over agaian?

                    using generators? How to use generators

                '''
                s = np.reshape(tau[i,j,:],(1,n3))


                # will be used for rprime and tprime in Ozisik solution.
                [RR,TT] = np.meshgrid(self.r,s)



                '''
                    ##### ERROR DUE TO ROUNDING OF HEAT SOURCE ####

                    error in rounding  5000000*math.exp(-(10**3)) becomes zero

                    #log10(e−10000)=−10000∗(0.4342944819)=−4342.944819

                    #e−1000=10−4342.944819=10−4343100.05518=1.13548386531×10−4343

                '''

                #g = 5000000*math.exp(-(10**3)) #*(RR - self.c*TT)**2) #[W / m^2] heat source.
                g = coeff * (RR - self.c*TT)**2


                K = (bm/math.sqrt(self.H2**2 + bm**2))*(math.sqrt(2)/self.b)*((scipy.special.jv(0,bm*RR))/(scipy.special.jv(0,bm*self.b)))

                #integral transform of heat source
                gbar = np.trapz(RR*K*g, self.r, 2) #trapz(Y,X,dx (spacing) )
                gbar = gbar.transpose()



                #boundary condition. BE SURE TO WRITE IN TERMS OF s!!!
                f2 = self.h2 * 37

                A = (self.alpha/self.k)*gbar + ((self.alpha*self.b)/self.k2)*((bm/math.sqrt(self.H2**2 + bm**2))*(math.sqrt(2)/self.b)*((scipy.special.jv(0,bm*self.b))/(scipy.special.jv(0,bm*self.b))))*f2
                #A is essentially a constant is this correct all the time?
                #What does A represent

                SS[i, j] = np.trapz(np.exp( (-self.alpha*bm**2)*(T[i,j] - s) )*A, s)

        #INSIDE M loop
        K = (bm / math.sqrt((self.H2 ** 2) + (bm ** 2)))*(math.sqrt(2) /self.b)*((scipy.special.jv(0, bm * R))/ (scipy.special.jv(0, bm * self.b)))


        U[:,:, m] = np.exp(-self.alpha * bm ** 2 * T)* K* Fbar + K* SS

        #print(['Eigenvalue ' num2str(m) ', found at time ' num2str(toc) ' seconds'])

Answer 1

评论中给出的答案汇编

特定于我的代码的答案: 正如 vorticity 所提到的，我的代码本身并没有最大程度地使用 numpy 和 scipy 包。

关于 Bessel，函数“royvib”提到使用 scipy 中的 .jo 而不是 .jv。调用特殊贝塞尔函数 jv。在计算上要昂贵得多，特别是因为我知道我会为我的许多声明使用零阶贝塞尔函数，从 jv -> j0 解决的微小变化加速了这个过程。

此外，我在循环外声明了变量，以防止为搜索我的适当函数而进行的昂贵调用。示例如下。

之前

for i in range(SomeLength):
    some_var = scipy.special.jv(1,coeff)

之后

Bessel = scipy.special.jv
for i in range(SomeLength):
    some_var = Bessel(1,coeff)

通过不使用点 ('.') 命令在每个循环中查看库，存储函数节省了时间。但是请记住，这确实会降低 python 的可读性，这是我选择在 python 中执行此项目的主要原因。我不知道这个步骤从我的过程中减少了多少时间。

特定于 Fortran: 因为我能够改进我的 python 代码，所以我没有走这条路，因为缺乏细节，但是“高性能标记”所说的一般答案是，是的，已经有一些库可以处理 Fortran 中的贝塞尔函数。

如果我将代码移植到 Fortran 或使用 f2py 混合 Fortran 和 python，我将相应地更新此答案。