I am trying to wrap an old and somewhat messy f77 code with f2py. I get everything to compile with 2 different compilers (gfortran and ifort), nevertheless the code segfaults when executed from python, while the original f77 works. The reason seems to be unallocated memory for an array. To be more precise, in the main subroutine the variable are declared as follows:
implicit double precision (a-h,o-z)
include 'gmm01f.par.f'
include 'ampld.par.f'
parameter (nangmax=1801)
parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
parameter (ni0=np*(np+1)*(2*np+1)/3+np*np)
parameter (ng0=np*(2*np**3+10*np**2+19*np+5)/6)
parameter (nrc=4*np*(np+1)*(np+2)/3+np)
integer u,v,u0,nmax(nLp),uvmax(nLp),ind(nLp),idshp(nLp)
double precision k,r0(9,nLp),x(nLp),xc(nLp),r00(3,nLp),
+ besj(0:2*np+1),besy(0:2*np+1),smue(4,4),i11(nangmax),
+ i21(nangmax),i22(nangmax),i12(nangmax),inat(nangmax),
+ pol(nangmax),dang(nangmax),
+ cscaxi(nLp),cscayi(nLp),cextxi(nLp),cextyi(nLp),
+ cabsxi(nLp),cabsyi(nLp),cexti(nLp),cabsi(nLp),cscai(nLp),
+ assymi(nLp),assymxi(nLp),assymyi(nLp),cprxi(nLp),
+ cpryi(nLp),cpri(nLp),drot(nrc),c0i(nLp),c1i(nLp)
double precision shp(3,nLp),
& RT11(NPN6,NPN4,NPN4),RT12(NPN6,NPN4,NPN4),
& RT21(NPN6,NPN4,NPN4),RT22(NPN6,NPN4,NPN4),
& IT11(NPN6,NPN4,NPN4),IT12(NPN6,NPN4,NPN4),
& IT21(NPN6,NPN4,NPN4),IT22(NPN6,NPN4,NPN4)
complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ci,cin,
+ s2x,s4x,s3y,s1y,s2xj,s4xj,s3yj,s1yj,
+ ref(nLp),refc(nLp),atr(2,np,nmp),at(nmp),bt(nmp),ek(np),
+ p0(nLp,nmp),q0(nLp,nmp),an(np),bn(np),atj(nmp),btj(nmp),
+ as00(nLp,nmp),bs00(nLp,nmp),as0(nLp,nmp),bs0(nLp,nmp),
+ as(nLp,nmp),bs(nLp,nmp),asc(nLp,nmp),bsc(nLp,nmp),
+ as1(nLp,nmp),bs1(nLp,nmp),ast(nLp,nmp),bst(nLp,nmp),
+ asp(nLp,nmp),bsp(nLp,nmp),asv(nLp,nmp),bsv(nLp,nmp),
+ asx(nLp,nmp),bsx(nLp,nmp),aty(nmp),bty(nmp),atyj(nmp),
+ btyj(nmp),Ay(2),By(2),Ayj(2),Byj(2),
+ tbar(nLp,2,2,nmp,nmp),tbar0(2,2,nmp,nmp),bar(2,2),
+ ekt(-2*np:2*np)
CHARACTER FLNAME*20,fileout*20,fileout1*19,flout*22
COMMON/MIESUB/ twopi,pih
common/rot/bcof(0:np+2),dc(-np:np,0:nmp)
common/fnr/fnr(0:2*(np+2))
common/pitau/pi(nmp0),tau(nmp0)
common/tran/atr
common/ig0/iga0(ni0)
common/g0/ga0(ng0)
common/cofmnv0/cof0(ni0)
common/crot/cofsr(nmp)
COMMON /TMAT/ RT11,RT12,RT21,RT22,IT11,IT12,IT21,IT22
where gmm01f.par.f is as follows:
parameter (nLp=10,np=30)
and ampld.par.f is:
Parameter (NPN1=100, NPNG1=500, NPNG2=2*NPNG1, NPN2=2*NPN1,
& NPL=NPN2+1, NPN3=NPN1+1,
& NPN4=NPN1, NPN5=2*NPN4, NPN6=NPN4+1)
What happens is that the memory for the cofsr(nmp) array seems, for some reason, not to be allocated, and anytime I try to assign a value to the array the program segfaults. This only happens for cofsr(nmp) and not for the other arrays, that are initialized to 0 at the beginning of the program. I can't really see any reason for this behavior, the only exception being the fact that in a subroutine cofsr(nmp) is initialized as
subroutine cofsrd(nmax)
include 'gmm01f.par.f'
parameter (nmp=np*(np+2))
double precision cofsr(nmp),lnfacd,c
common/crot/cofsr
while otherwise we only have the
implicit double precision (a-h,o-z)
statement.
As additional informations the beginning of the subroutine is the following:
SUBROUTINE gmm03f(idMie,
& idAlpha,idBeta,iseed,
& small,MXINT,
& NADD,
& fint,
& sang,pang,
& w,irat,
& nL,
& idshp,shp,r0,
& cextxi,cextyi,cscaxi,
& cscayi,cabsxi,cabsyi)
Cf2py intent(in) idMie
Cf2py intent(in) idAlpha,idBeta,iseed
Cf2py intent(in) small,MXINT
Cf2py intent(in) NADD
Cf2py intent(in) fint
Cf2py intent(in) sang,pang
Cf2py intent(in) w,irat
Cf2py intent(in) nL
Cf2py intent(in) idshp,shp,r0
Cf2py intent(out) cextxi,cextyi,cscaxi,cscayi,cabsxi,cabsyi
and the signature file is:
python module stm ! in
interface ! in :stm
subroutine gmm03f(idmie,idalpha,idbeta,iseed,small,mxint,nadd,fint,sang,pang,w,irat,nl,idshp,shp,r0,cextxi,cextyi,cscaxi,cscayi,cabsxi,cabsyi) ! in :stm:stm.f
integer intent(in) :: idmie
integer intent(in) :: idalpha
integer intent(in) :: idbeta
integer intent(in) :: iseed
double precision intent(in) :: small
integer intent(in) :: mxint
integer intent(in) :: nadd
double precision intent(in) :: fint
double precision intent(in) :: sang
double precision intent(in) :: pang
double precision intent(in) :: w
integer intent(in) :: irat
integer intent(in) :: nl
integer dimension(10),intent(in) :: idshp
double precision dimension(3,10),intent(in) :: shp
double precision dimension(9,10),intent(in) :: r0
double precision dimension(10),intent(out) :: cextxi
double precision dimension(10),intent(out) :: cextyi
double precision dimension(10),intent(out) :: cscaxi
double precision dimension(10),intent(out) :: cscayi
double precision dimension(10),intent(out) :: cabsxi
double precision dimension(10),intent(out) :: cabsyi
double precision dimension(527) :: pi
double precision dimension(527) :: tau
double precision dimension(65) :: fnr
double precision dimension(101,100,100) :: rt11
double precision dimension(101,100,100) :: rt12
double precision dimension(101,100,100) :: rt21
double precision dimension(101,100,100) :: rt22
double precision dimension(101,100,100) :: it11
double precision dimension(101,100,100) :: it12
double precision dimension(101,100,100) :: it21
double precision dimension(101,100,100) :: it22
integer dimension(19810) :: iga0
double precision :: twopi
double precision :: pih
double precision dimension(960) :: cofsr
complex*16 dimension(2,30,960) :: atr
double precision dimension(19810) :: cof0
double precision dimension(317875) :: ga0
double precision dimension(33) :: bcof
double precision dimension(np + 31,961) :: dc
common /pitau/ pi,tau
common /fnr/ fnr
common /tmat/ rt11,rt12,rt21,rt22,it11,it12,it21,it22
common /ig0/ iga0
common /miesub/ twopi,pih
common /crot/ cofsr
common /tran/ atr
common /cofmnv0/ cof0
common /g0/ ga0
common /rot/ bcof,dc
include 'gmm01f.par.f'
include 'ampld.par.f'
end subroutine gmm03f
end interface
end python module stm
Thanks a Lot
I found an answer, even if I can't really say that I understand the rationale behind it. What causes the problem is the name of the common block:
common/crot/cofsr(nmp)
when switched with a different name such as
common/testok/cofsr(nmp)
everything works as expected, therefore my problem is solved.