keflavich / pyradex

Python interface to RADEX
BSD 3-Clause "New" or "Revised" License
19 stars 12 forks source link

Segmentation faults May 2021 #34

Closed keflavich closed 3 years ago

keflavich commented 3 years ago

I've tracked down segmentation faults happening during the readdata phase to this line:

            read(11,*) dummy,lcu(icoll),lcl(icoll),
     $           (coll(ipart,icoll,itemp),itemp=1,ntemp)

which appears to be caused by accessing the array coll in any way; attempting to read or write to/from that array, even within acceptable dimensions, triggers a segfault. I can reproduce it with something as simple as:

         coll(1,1,1)=1e5
         write(*,*) coll

I don't have a solution right now. This clearly isn't an actual fortran error, since it didn't occur on older fortran/f2py versions. Maybe real*8 isn't supported any more?

keflavich commented 3 years ago

The code in question isn't in the repository, so I reproduce it here:

c     readdata.f
c
      SUBROUTINE readdata
      implicit none
      include 'radex.inc'

c This file is part of the RADEX software package
c to calculate molecular excitation and radiative
c transfer in a homogeneous medium.
c
c Documentation for the program is posted at
c https://sron.rug.nl/~vdtak/radex/index.shtml  
c
c Although this program has been thoroughly tested, the
c authors do not claim that it is free of errors and
c gives correct results in all situations.
c
c Publications using this program should make a reference
c to our paper: A&A 468, 627 (2007).
c
c     ---------------------------------------------------------
c     
c     Reads molecular data files (2003 format)

      integer ilev,jlev   ! to loop over energy levels
      integer iline       ! to loop over lines
      integer ipart,jpart ! to loop over collision partners
      integer itemp       ! to loop over collision temperatures
      integer icoll       ! to loop over collisional transitions
      integer dummy       ! to skip part of the file

      integer id(maxpart)      ! to identify collision partners

c     upper/lower levels of collisional transition 
      integer lcu(maxcoll),lcl(maxcoll) 
      real*8 coll(maxpart,maxcoll,maxtemp)
      real*8 colld(maxpart,maxlev,maxlev)
      real*8 ediff
      real*8 temp(maxtemp) ! collision temperatures

      character*120 collref ! text about source of collisional data
      character*120 blah

c     to interpolate rate coeffs
      integer iup,ilo,nint
      real*8 tupp,tlow,fint

c     to verify matching of densities and rates
      logical found

c     to calculate thermal o/p ratio for H2
      real*8 opr

c     Executable part begins here.
      write(*,*) 'readdata'
      write(*,*) 'molfile: "',molfile,'"'

      open(unit=11,file=molfile,status='old',err=99)

      write(*,*) "molfile successfully opened"
      debug=1

c     in the header, every second line is a comment
  101 format(a)
      read(11,*) 
      read(11,101) specref
      read(11,*) 
      read(11,*) amass
      read(11,*) 
      read(11,*) nlev

      if (nlev.lt.1) stop 'error: too few energy levels defined' 
      if (nlev.gt.maxlev) stop 'error: too many energy levels defined' 
      if (debug) write(*,*) 'readdata: basics'

c     Term energies and statistical weights
      read(11,*)
      do ilev=1,nlev
         read(11,*) dummy,eterm(ilev),gstat(ilev),qnum(ilev)
         if ((dummy.lt.1).or.(dummy.gt.nlev))
     $        stop 'error:illegal level number'
      enddo

      if (debug) write(*,*) 'readdata: levels'!,(eterm(ilev),ilev=1,nlev)

c     Radiative upper & lower levels and Einstein coefficients
      read(11,*) 
      read(11,*) nline
      read(11,*)

      if (nline.lt.1) stop 'error: too few spectral lines defined' 
      if (nline.gt.maxline) stop
     $     'error: too many spectral lines defined'

      do iline=1,nline
         read(11,*) dummy,iupp(iline),ilow(iline),aeinst(iline)
     $      ,spfreq(iline),eup(iline)
         if ((dummy.lt.1).or.(dummy.gt.nline))
     $        stop 'error:illegal line number'
         xnu(iline)=(eterm(iupp(iline))-eterm(ilow(iline)))
         if ((xnu(iline).lt.eps))
     $        stop 'error:illegal line frequency'
      enddo

      if (debug) write(*,*) 'readdata: lines'!,(xnu(iline),iline=1,nline)

c     Number of collision partners
      read(11,*)
      read(11,*) npart
      if (npart.lt.1) stop 'error: too few collision partners defined' 
      if (npart.gt.maxpart) stop 'error: too many collision partners'

      write(*,*) 'readdata: npart=',npart
      write(*,*) "maxtemp=",maxtemp
      write(*,*) "temp=",temp
      write(*,*) "tkin=",tkin

 102  format(i1,a)
      do ipart=1,1
         read(11,*)
         read(11,102) id(ipart),collref 
         read(11,*)
         read(11,*) ncoll
         write(*,*) 'ncoll',ncoll
      if (ncoll.lt.1) stop 'error: too few collision rates defined' 
      if (ncoll.gt.maxcoll) stop 'error: too many collision rates'
         read(11,*) blah
         write(*,*) 'blah',blah
         read(11,*) ntemp
         write(*,*) 'ntemp',ntemp
         if (ntemp.lt.0) stop
     $        'error: too few collision temperatures defined'
         if (ntemp.gt.maxtemp) stop
     $        'error: too many collision temperatures'
         read(11,*)
         read(11,*) (temp(itemp),itemp=1,ntemp)
         write(*,*) 'temp',temp
         write(*,*) 'itemp',itemp
         read(11,*)

         if (debug) write(*,*) 'ready to read ',ncoll,' rates for '
     $        ,ntemp,' temperatures for partner ',ipart

c        these next two lines trigger segfault         

         coll(1,1,1)=1e5
         write(*,*) coll

c the text commented out below here is the functional "real" code that
c fails
c        do icoll=1,ncoll
c           read(11,*) dummy,lcu(icoll),lcl(icoll),
c    $           (coll(ipart,icoll,itemp),itemp=1,ntemp)
c           write(*,*) dummy,lcu(icoll),lcl(icoll)
c        if ((dummy.lt.1).or.(dummy.gt.ncoll))
c    $        stop 'error:illegal collision number'
c        enddo
      enddo
      stop
keflavich commented 3 years ago

SOLVED by https://github.com/numpy/numpy/issues/19112#issuecomment-850039068

keflavich commented 3 years ago

and closed by 7c1035d8f920dbe841522936456de70383d6285d