f2py：如何将二维列表传递给Fortran 77

Question

我在把二维数组传给Fortran时遇到了麻烦。我想把一些不重叠的光谱合并在一起。首先，我选择x轴上的一些点，然后把所有数据插值到这个新的公共网格上。我把光谱存储在Python中的一个二维列表里。

在Python 2.7中，这个方法可以用，但速度很慢：

for i in range(len(wlp)):
  print wlp[i],
  for a in range(len(datax)):
    inrange = 0
    if datax[a][0] >= wlp[i] or datax[a][-1] <= wlp[i]:
      for b in range(len(datax[a])-1):
        if float(datax[a][b]) <= wlp[i] and float(datax[a][b+1]) >= wlp[i]:
          sp = float(datax[a][b]); ep = float(datax[a][b+1])
          delx = ep-sp; dely = float(data[a][b+1])-float(data[a][b])
          ji =  (dely/delx)*(wlp[i]-sp)+float(data[a][b])
          inrange = 1
    if inrange == 0: ji = '?0'
    else: ji = ji * weights[a]
    print ji,
  print

公共的x网格在第一列打印，所有插值后的光谱在后面的列中打印。如果有一些光谱的长度较短，超出了范围，就会打印“?0”。这有助于后面为每个数据点设置合适的权重。

为了加快速度，我最终写了这个Fortran子程序，并用f2py来调用：

c wlp = x axis points (wavelength)
c lenwlp = length of list wlp, len(wlp)
c datay = 2D python list with flux
c datax = 2D python list with wavelength
c lendatax = number of spectra, len(datax)
c datax_pl = list of the lengths of all spectra
c weights = list of optional weights
c maxi = length of the longest spectrum   

C============================================================================80
       SUBROUTINE DOIT(wlp,lenwlp,datay,datax,lendatax,datax_pl,
     .                 weights,maxi)
C============================================================================80
       INTEGER  I,a,b,lenwlp,inrange,datax_pl(*),maxi,lendatax
       DOUBLE PRECISION WLP(*),SP,EP,DELY,DELX,ji
       DOUBLE PRECISION WEIGHTS(*)
       DOUBLE PRECISION DATAY(lendatax,maxi)
       DOUBLE PRECISION DATAX(lendatax,maxi)   

 2     FORMAT (E20.12,  1X, $)
 3     FORMAT (A, $)
 4     FORMAT (1X)
       I = 1
       DO WHILE (I.LE.lenwlp)
         WRITE(*,2) WLP(I)
         DO a=1,lendatax
           inrange = 0
           ji = 0.0
           IF (datax(a,1).ge.WLP(I) .or.
     .         datax(a,datax_pl(a)).le.WLP(I)) THEN
             DO b=1,datax_pl(a)-1
               IF (DATAX(a,b).LE.WLP(I) .and.
     .           DATAX(a,b+1).GE.WLP(I)) THEN
                   SP = DATAX(a,b); EP = DATAX(a,b+1)
                   DELX = EP - SP; DELY = datay(a,b+1)-datay(a,b)
                   if (delx.eq.0.0) then
                     ji = datay(a,b)
                   else
                     ji = (DELY/DELX)*(WLP(I)-SP)+datay(a,b)
                   end if
                   inrange = 1
               END IF
             END DO
           END IF
           IF (inrange.eq.0) THEN
             WRITE(*,3) ' ?0'
           ELSE
             WRITE(*,2) ji*WEIGHTS(a)
           END IF
         END DO
         I = I + 1
         write(*,4)
       END DO
       END

这个程序在gfortran 4.8中编译得很好。然后我在Python代码中导入它，设置好列表并运行这个子程序：

import subroutines
wlp = [...]
data = [[...],[...],[...]]
datax = [[...],[...],[...]]
datax_pl = [...]
weights = [...]
maxi = max(datax_pl)

subroutines.doit(wlp,len(wlp),data,datax,len(datax),datax_pl,weights,maxi)

它返回了：

ValueError: setting an array element with a sequence.

我传递了列表和最长光谱的长度（maxi），这应该定义Fortran中的最大维度（？）。我不需要返回值，所有的输出都在标准输出上打印。

问题可能出在数组声明的最开始。我对这方面没有经验……任何建议都很受欢迎。

性能优化二维数组 fortran f2py 数据插值 gfortran 数组声明光谱合并

f2py：如何将二维列表传递给Fortran 77

1 个回答

撰写回答