" segmentation fault " in fortran

Question

I have questions : Is it possible to define final-value of do loop(ngal) as a free parameter with out exact value or not when I don't know the exact size of my catalog either I don't know the size of nm and nz? ( the whole code as you wished )

integer i,iz,im,ngal,n_thresh,nmax,nz,nm,ok,total,ii,j
integer, dimension(:,:),allocatable :: ngalb
double precision ramin,ramax,decmin,decmax,zmin,zmax,mag_min,mag_max,step,bin,S,dix,fi,Vmax,zbin,magbin
double precision, dimension (:),allocatable :: ra,dec,mag_g,mag_r,mag_i,redshift,zup
double precision, dimension (:,:,:),allocatable :: magb,zb
character plotname*100,zbin_str*4

open(unit=1,file='Lagos12.i.asc')
open(2,file='bin.asc')
open(3,file='mag_vol.asc')
open(4,file='luminosity_func.asc')



ramin=0.0
ramax=2.0
decmin=-2.0
    decmax=2.0
zmin=0.0
zmax=3.2
mag_min=13
mag_max=26
step=0.1
bin=0.5
n_thresh=20

 allocate (ra(ngal),stat=ok)
 allocate (dec(ngal),stat=ok)
 allocate (mag_g(ngal),stat=ok)
 allocate (mag_r(ngal),stat=ok)
 allocate (mag_i(ngal),stat=ok)
 allocate (redshift(ngal),stat=ok)
 allocate(ngalb(nz,nm),stat=ok)




     ngalb=0

     do i=1,ngal
     read(1,*)  ra(i),dec(i),mag_g(i),mag_r(i),mag_i(i),redshift(i)

     iz=int((redshift(i)-zmin)/step)+1
     im=int((mag_i(i)-mag_min)/bin)+1
     ngalb(iz,im)=ngalb(iz,im)+1
     enddo

 nmax=0
     do iz=1,nz
        do im=1,nm
           nmax=max(nmax,ngalb(iz,im))
        enddo
     enddo


     allocate(magb(nz,nm,nmax),stat=ok)
     allocate(zb(nz,nm,nmax),stat=ok)
   allocate(zup(nz),stat=ok)


     ngalb=0

     do i=1,ngal

     read(1,*)ra(i),dec(i),mag_g(i),mag_r(i),mag_i(i),redshift(i)

     iz=int((redshift(i)-zmin)/step)+1
     im=int((mag_i(i)-mag_min)/bin)+1


    ngalb(iz,im)=ngalb(iz,im)+1
     j=ngalb(iz,im)
     magb(iz,im,j)=mag_i(i)
     zb(iz,im,j)=redshift(i)
 write(2,'(2x,3i5,2x,i5,2x,f10.7,2x,f10.5)')iz,im,j,ngalb(im,iz),magb(im,iz,j),zb(im,iz,j)
     enddo





deallocate(magb,zb)
deallocate(ngalb)
deallocate (ra,dec,mag_g,mag_r,mag_i,redshift,zup)

the question is in itself not understandable, but from the code I would assume OP's problem is that the size of most arrays should equal the number of lines read (unknown when allocatin), and the size of ngalb depends on the shifts in the data (also unkown when allocating). So, the solution is to think a little bit. — steabert
– steabert, Commented Sep 23, 2013 at 15:29
@steabert: I agree. I see two possibilities to do that: Use fixed-size arrays and fill them until there are no more datasets in the file (using iostat or end), or reading in the file first using dummy-arguments to determine the number of lines (records?), then doing the allocations and finally read in the file for real. — Alexander Vogt
– Alexander Vogt, Commented Sep 23, 2013 at 15:32
Or even reallocate the arrays, when they are full and continue reading. But for a beginner I would choose 1 pass for number of lines, rewind and read the lines. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Sep 23, 2013 at 16:08

Alexander Vogt · Accepted Answer · 2013-09-23 17:40:59Z

4

To sum up the comments by High Performance Mark, steabert, Vladimir F, and myself:

! First pass: Determine sizes
nm = 0 ; nz = 0 ; ngal = 0
do
  read(1,*,iostat=ierr)  dummy_ra, dummy_dec, dummy_mag_g, dummy_mag_r, &
                         dummy_mag_i, dummy_redshift
  if ( ierr < 0 ) exit ! End of file
  if ( ierr > 0 ) stop 'An error occured reading from unit 1'
  nz = max(nz,int((dummy_redshift-zmin)/step)+1)
  nm = max(nm,int((dummy_mag_i-mag_min)/bin)+1)
  ngal = ngal + 1
enddo

! Allocate arrays
allocate (ra(ngal), dec(ngal), mag_g(ngal), mag_r(ngal), mag_i(ngal), &
          redshift(ngal), ngalb(nm,nz), stat=ierr) 
if ( ierr /= 0 ) stop 'Cannot allocate memory'
ngalb=0

! Second pass: Read in
rewind(1)
do i=1,ngal
  read(1,*) ra(i), dec(i), mag_g(i), mag_r(i), mag_i(i), redshift(i)
  iz = int((redshift(i)-zmin)/step)+1
  im = int((mag_i(i)-mag_min)/bin)+1
  ngalb(im,iz) = ngalb(im,iz)+1
enddo

edited Sep 23, 2013 at 17:40

answered Sep 23, 2013 at 16:23

Alexander Vogt

18.2k13 gold badges56 silver badges72 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

" segmentation fault " in fortran

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related