19990826-4.f   [plain text]

* Culled from 970528-1.f in Burley's g77 test suite.  Copyright
* status not clear.  Feel free to chop down if the bug is still
* reproducible (see end of test case for how bug shows up in gdb
* run of f771).  No particular reason it should be a noncompile
* case, other than that I didn't want to spend time "fixing" it
* to compile cleanly (with -O0, which works) while making sure the
* ICE remained reproducible.  -- burley 1999-08-26

* Date: Mon, 26 May 1997 13:00:19 +0200 (GMT+0200)
* From: "D. O'Donoghue" <dod@da.saao.ac.za>
* To: Craig Burley <burley@gnu.ai.mit.edu>
* Cc: fortran@gnu.ai.mit.edu
* Subject: Re: g77 problems

	program dophot
	parameter (napple = 4)
        common /window/nwindo,ixwin(50),iywin(50),iboxwin(50),itype(50)
        common/io/luout,ludebg
	common/search/nstot,thresh
	common /fitparms / acc(npmax),alim(npmax),mit,mpar,mfit1,
     +                     mfit2,ind(npmax)
	common /starlist/ starpar(npmax,nsmax), imtype(nsmax),
	1shadow(npmax,nsmax),shaderr(npmax,nsmax),idstr(nsmax)
	common /aperlist/ apple(napple ,nsmax)
	common /parpred / ava(npmax)
	common /unitize / ufactor
	common /undergnd/ nfast, nslow
	common/bzero/ scale,zero
	common /ctimes / chiimp, apertime, filltime, addtime
	common / drfake / needit 
	common /mfit/ psfpar(npmax),starx(nfmax),stary(nfmax),xlim,ylim
	common /vers/ version
 	logical needit,screen,isub,loop,comd,burn,wrtres,fixedxy
	logical fixed,piped,debug,ex,clinfo
	character header*5760,rhead*2880
	character yn*1,version*40,ccd*4,infile*20
	character*30 numf,odir,record*80
	integer*2 instr(8)
	character*800 line
	external pseud0d, pseud2d, pseud4d, pseudmd, shape
C
C	Initialization
	data burn,   fixedxy,fixed,  piped 
     +     /.false.,.false.,.false.,.false./
	data needit,screen,comd,isub
     + /.true.,.false.,.true.,.false. /
	data acc / .01, -.03, -.03, .01, .03, .1, .03 /
	data alim / -1.0e8, 2*-1.0e3, -1.0e8, 3*-1.0e3 /
C
	version = 'DoPHOT Version 1.0 LINUX May 97 '
        debug=.false.
        clinfo=.false.
	line(1:800) = ' '
	odir = ' '
C
C
C	Read default tuneable parameters 
	call tuneup ( nccd, ccd, piped, debug )
	version(33:36) = ccd(1:4)
C
      
        ludebg=6
        if(piped)then
          yn='n'
        else
	  write(*,'(''****************************************'')')
	  write(*,1000) version
	  write(*,'(''****************************************''//)')           
C                                                             
          write(*,'(''Screen output (y/[n])? '',$)')             
	  read(*,1000) yn
        end if
	if(yn.eq.'y'.or.yn.eq.'Y') then
          screen=.true.                
          luout=6
        else
          luout=2
        end if
C
        if(piped)then
          yn='y'
        else
          write(*,'(''Batch mode ([y]/n)? '',$)')
          read(*,1000) yn
        end if
	if(yn.eq.'n'.or.yn.eq.'N') comd = .false.
C                                          
	if(.not.comd) then                                     
          write(*,
     *         '(''Do you want windowing ([y]/n)? '',$)')
          read(*,1000)yn
          iwindo=1
          if(yn.eq.'n'.or.yn.eq.'N')then
            nwindo=0
            iwindo=0
          end if
C
          write(*,
     *       '(''Star classification info (y/[n]) ?'',$)')
          read(*,1000)yn
          clinfo=.false.
          if(yn.eq.'y'.or.yn.eq.'Y')clinfo=.true.
C
	  write(*,
     *        '(''Create a star-subtracted frame (y/[n])? '',$)')
	  read(*,1000) yn                                     
 	  if(yn.eq.'y'.or.yn.eq.'Y') isub = .true.
C               
	  write(*,'(''Apply after-burner (y/[n])? '',$)')
	  read(*,1000) yn
	  if ( yn.eq.'y'.or.yn.eq.'Y' ) burn = .true.
	  wrtres = burn
C
	  write(*,'(''Read from fixed (X,Y) list (y/[n])? '',$)')
	  read(*,1000) yn
	  if ( yn.eq.'y'.or.yn.eq.'Y' ) then
	    fixedxy = .true.
	    fixed = .true.
	    burn = .true.
	    wrtres = .true.
	  endif
	endif         
        iopen=0
C
C       This is the start of the loop over the input files
c
        iframe=0
        open(10,file='timing',status='unknown',access='append')

1	ifit = 0
	iapr = 0
	itmn = 0
	model = 1
	xc = 0.0
	yc = 0.0
	rc = 0.0
	ibr = 0
	ixy = 0
C	
        iframe=iframe+1
        tgetpar=0.0
        tsearch=0.0
        tshape=0.0
        timprove=0.0
C
C	Batch mode ...

	if ( comd ) then
          if(iopen.eq.0)then
            iopen=1
            open(11,file='dophot.bat',status='old',err=995)
          end if
          read(11,1000,end=999)infile
c         now read in the parameter instructions. these are:
c         instr(1) : if 1, specifies uncrowded field, otherwise crowded 
c         instr(2) : if 1, specifies sequential frames of same field
c                          with a window around the stars of interest -
c                          all other objects are ignored
c         instr(3) : if 0, takes cmin from dophot.inp (via tuneup)
c                    if>0, sets cmin=instr(3)
c         instr(4) : if 0, does nothing
c                    if 1, then opens a file called classifications
c                    sets clinfo to .true. and writes out the star
c                    typing info to this file
c         instr(5) : Delete the shd.nnnnnnn file
c         instr(6) : Delete the out.nnnnnnn file
c         instr(7) : Delete the input frame
c         instr(8) : Create a star-subtracted frame 
          read(11,*)instr
          read(11,*)ifit,iapr,tmn,model,xc,yc,rc,ibr,ixy
          nocrwd = instr(1)
          iwindo=instr(2)
          if(iwindo.eq.0)nwindo=0
          itmn=tmn
          if ( instr(3).gt.0 ) cmin=instr(3)
          clinfo=.false.
          if ( instr(4).gt.0 )then
            clinfo=.true.
            open(12,file='classifications',status='unknown')
            ludebg=12
          end if
 	  if ( instr(8).ne.0 ) then
	    isub = .true.            
	  else
	    isub = .false.
	  endif
C
	  if(ibr.ne.0) burn = .true.
	  if(ixy.ne.0) then
	    fixedxy = .true.
	    fixed = .true.
	    burn = .true.
	    goto 20
          endif
          if(iwindo.eq.0)then
            write(6,10)iframe,infile(1:15)
   10       format('  ***** DoPHOT-ing frame ',i4,': ',a)
            if(ludebg.eq.12)write(ludebg,11)iframe,infile(1:15)
   11       format(////'  ',62('*')/
     *                 '  *     DoPHOT-ing frame ',i4,': ',a,
     *                 '                 *'/'  ',62('*'))
          end if
          if(iwindo.eq.1)then
            write(6,12)iframe,infile(1:15)
   12       format('  ***** DoPHOT-ing frame ',i4,': ',a,
     *             '   - Windowed *****')
            if(ludebg.eq.12)write(ludebg,13)iframe,infile(1:15)
   13       format(////'  ',62('*')/
     *                 '  *     DoPHOT-ing frame ',i4,': ',a,
     *                 '   - Windowed    *'/2x,62('*'))
          end if
C
C	Interactive...
	else
	  write(*,'(''Image name: '',$)')
	  read(*,1000) infile
	  if(infile(1:1).eq.' ') goto 999                     
1000	  format(a)                          
          write(*,'(''Crowded field mode ([y]/n) ? '',$)')
          read(*,1000)yn
          nocrwd=0
          if(yn.eq.'n'.or.yn.eq.'N')nocrwd=1
	  if(.not.fixed) then
	    write(*,1001)
1001        format('Sky model ([1]=Plane, 2=Power, 3=Hubble)? ',$)
            read(*,1000)record
            if(record.ne.' ')then
	      read(record,*) model
            else
              model=1
            end if
	  else         
	    burn=.true.
	    goto 20           
	  endif
	endif
C
C       if windowing, open the file and read the window
        if(iwindo.eq.1)then
          inquire(file='windows',exist=ex)
          if(.not.ex)go to 997
          if(iframe.eq.1)open(9,file='windows',status='old')
          nwindo=0
    2     read(9,*,end=3)intype,inx,iny,inbox
          nwindo=nwindo+1
          if(nwindo.gt.50)then
            print *,'too many windows - max = 50'
            stop
          end if
          ixwin(nwindo)=inx
          iywin(nwindo)=iny
          iboxwin(nwindo)=inbox
          itype(nwindo)=intype
          go to 2

    3     rewind 9
          if(screen)print 4,(itype(j),ixwin(j),iywin(j),iboxwin(j),
     *                       j=1,nwindo)
    4     format(' Windows: Type   X    Y   Size'/
     *           (I13,i6,i5,i5))
        end if

	t1 = cputime(0.0)
C
C	Read FITS frame.
	call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line,ccd)
C
C	Ignore frame if not the correct chip
	if(nc.lt.0) goto 900
C                     
C	Estimate starting PSF parameters.
   15   call getparams(nfast,nslow,gxwid,gywid,skyval,tmin,tmax,
     *                 iframe)
        tgetpar = cputime(t1) + tgetpar
        if(debug)write(ludebg,16)iframe,skyval,gxwid,gywid,tmin,tmax
   16   format(' Getparams on frame ',i4,'  sky ',f6.1,'  gxwid ',f5.1,
     *         '  gywid ',f5.1,'  tmin ',f5.1,'  tmax ',f5.1)
C
C	Initialize
	do j=1,nsmax
	  imtype(j) = 0
	  do i=1,npmax  
	    shadow(i,j)=0.                               
	    shaderr(i,j)=0.
	  enddo
	enddo
C              
	skyguess=skyval
	tfac = 1.0            
C	Use 4.5 X SD as fitting width       
	fitr=fitfac*(gxwid*asprat*gywid)**0.25 + 0.5                      
	i=fitr
	irect(1)=i
	irect(2)=fitr/asprat 
C	Use 4/3 X FitFac X SD as aperture width
	gmax = asprat*gywid
 	if(gxwid.gt.gmax) gmax=gxwid
	aprw = 1.33*fitfac*sqrt(gmax) + 0.5
	i = aprw
	arect(1) = i
	i = aprw/asprat + 0.1
	arect(2) = i          
C                                     
	if(irect(1).gt.50) irect(1)=50
	if(irect(2).gt.50) irect(2)=50  
	if(arect(1).gt.45.) arect(1)=45.
	if(arect(2).gt.45.) arect(2)=45.
C
	if (screen) call htype(line,skyval,.false.,fitr,ngr,ncon)
C
C       Prompt for further information         
	if ( .not.comd ) then
          write(*,1002)
 1002     format(/'The above are the inital parameters DoPHOT'/
     *            'has found. You can change them now or accept'/
     *            'the values in [ ] by pressing enter'/)

          write(*,1004)tmin
 1004     format('Enter Tmin: threshold for star detection',
     *           ' [',f5.1,']  ',$)
          read(*,1000)record
          if(record.ne.' ')read(record,*)tmin

          write(*,1005)cmin
 1005     format('Enter Cmin: threshold for PSF stars',
     *           '      [',f5.1,']  ',$)
          read(*,1000)record
          if(record.ne.' ')read(record,*)cmin

          write(*,1006)
 1006     format('Do you want to fix the aperture mag size ?',
     *           ' (y/[n]) ')
          read(*,1000)record
          if(record.eq.'y'.or.record.eq.'Y')then
            write(*,1007)
 1007       format('Enter the size in pixels: ',$)
            read(*,*)iapr
 	    if(iapr.gt.0) then          
              arect(1)=iapr
              i = iapr/asprat + 0.1
              arect(2)=i
            end if
	  endif                     
C
	  write(*,1008)
 1008     format('Satisfied with other input parameters ? ([y]/n)?',$)
	  read(*,1000) yn        
          if(yn.eq.'n'.or.yn.eq.'N')then
            yn='n'
          else
            yn='y'
          end if
	  if(.not.(yn.eq.'y'.or.yn.eq.'Y') ) call input
	else
	  if ( ifit.ne.0 ) then
	    irect(1)=ifit
	    irect(2)=(ifit/asprat + 0.1)
	  endif              
	  if ( iapr.ne.0 ) then
	    arect(1)=iapr
	    i = iapr/asprat + 0.1
	    arect(2)=i
	  endif                                       
	  if ( itmn.ne.0 ) tmin = itmn
 	  if ( .not.(xc.eq.0.0.and.yc.eq.0.0) ) then
	    xcen = xc
 	    ycen = yc
          endif
	endif          
C
C--------------------------------
C
C
	call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
     +nfast, nslow )
C
C       if the uncrowded field option has been chosen, jump
C       straight to the minimum threshold
C
        if(nocrwd.eq.1)tmax=tmin
C             
C	Adjust tfac so that thresh ends precisely on Tmin.
	if(tmin/tmax .gt. 0.999) then
	  thresh = tmin
	  tfac = 1.          
	else                                                     
	  thresh = tmax
	  xnum = alog10(tmax/tmin)/alog10(2.**tfac)
	  if(xnum.gt.1.5) then
	    xnum = float(nint(xnum))
	  else if(xnum.ge.1) then               
	    xnum = 2.0
	  else             
	    xnum = 1.0             
	  endif                                         
	  tfac = alog10(tmax/tmin)/alog10(2.)/xnum                   
	endif
C                       
C------------------------------------------------------------------------
C                            
C         This is the BIG LOOP which searches the frame for stars
C               with intensities > thresh.                
C        
C-----------------------------------------------------------------------
C                           
	loop = .true.
	nstot = 0
	do while ( loop )   
	  loop = thresh/tmin .ge. 1.01
	  write(luout,1050) thresh
1050	  format(/20('-')/'THRESHOLD: ', f10.3)
	  if(ludebg.eq.12)write(ludebg,1050) thresh
C
C         Fit given model to sky values.
C
          call varipar(nstot, nfast, nslow )            
	  t1 = cputime(0.0)
C               
C         Identifies potential objects in cleaned array IMG
 	  nstar = isearch( pseud2d, nfast, nslow , clinfo)
	  tsearch = cputime(t1) + tsearch
C                                                                   
    	  if ( (nstar .ne. 0).or.(xnum.lt.1.5) ) then
C                                           
C           Performs 7-parameter PSF fit and determines nature of object.
	    t1 = cputime(0.0)
	    call shape(pseud2d,pseud4d,nfast,nslow,clinfo)
	    tshape = cputime(t1) + tshape
C                           
C           Computes average sky values etc from star list
 	    call paravg
  	    t1 = cputime(0.0)
C                                                          
C           Computes 4-parameter fits for all stellar objects using 
C           new average shape parameters.  
  	    call improve(pseud2d,nfast,nslow,clinfo)
	    timprove = cputime(t1) + timprove
	  end if                         
C
C         Calculate aperture photometry on last pass.
	  if(.not.loop) call aper ( pseud2d, nstot, nfast, nslow )
C             
 	  totaltime = (tgetpar+tsearch+tshape+timprove)
	  write(3,1060) totaltime
 	  write(4,1060) totaltime
	  write(luout,1060) totaltime
1060	  format('Total CPU time consumed:',F10.2,' seconds.')
          write(10,1070)infile,tgetpar,tsearch,tshape,timprove,
     *                  totaltime
1070      format(a20,'   T(getp/f)',f5.1,'  T(search)',f5.1,
     *               '  T(shape)',f5.1,'  T(improve)',f5.1,
     *               '  Total',f6.1)
	  call title (line,skyval,.false.,fitr,ngr,ncon,strint,ztot,nums)
	  rewind(2)          
 	  rewind(3)                              
  	  rewind(4)
C
	  call output ( line )
C
C         Now reduce the threshold and loop back
C
	  thresh = thresh/2.**tfac
  	end do                   
C                              
C--------- END OF BIG LOOP ---------------------------------------
C                      
C	If after-burner required, residuals from analytic PSF are computed
C	and stored in RES.
C	
20	if ( burn ) then
C 	
C	If using a fixed (X,Y) coordinate list, read it.
	 if (fixed) then
C	 Read the image frame
 	  call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line)
C
C	 Initialize arrays, open files etc.
	  call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon,
     +nfast, nslow )
C
C	 Read the XY list                                 
	  write(luout,'(''Reading XY list ...'')')
	  call xylist(numf, nc, ios )
	  if(ios.ne.0) then
	   fixed = .false.
	   write(luout,'(''SXY file absent or incorrect...'')')
	   goto 15
	  endif
C
	  call htype(line,skyval,.false.,fitr,ngr,ncon)
C
C	 Remove good stars
	  write(luout,'(''Cleaning frame of stars: '',i8)') nstot
	  call clean ( pseud2d, nstot, nfast, nslow, -1)
C             
C	Calculate aperture photometry
C	  call aper ( pseud2d, nstot, nfast, nslow )
	 else            
	   rewind(3)
	   rewind(4)         
	 endif 
C             
C-----------------------
C	Flag all stars close together in groups.  Keep making the distance
C	criterion FITR smaller until the maximum number in a group is less
C	than NFMAX         
C                      
	 fitr = amax1(arect(1),arect(2))
	 fitr = fitr + 2.0      
	 nmax = 10000 
	 write(*,'(''Regrouping ...'')')
C              
	 do while ( nmax.gt.nfmax )
  	  fitr = fitr - 1.0       
	  write(luout,'(''Min distance ='',f8.1)') fitr
	  call regroup( fitr, ngr, nmax )     
	 enddo
C
	 xlim = irect(1)/2 
	 ylim = irect(2)/2
C
C	Calculate normalized PSF residual from PSEUD2D
	 call getres (pseud0d,pseud2d,strint,rmn,rmx,nfast,nslow,irect,
     +arect,ztot,nums)
	 if(nums.eq.0) then
	  write(luout,'(''No suitable PSF stars!'')')
	  goto 30
	 endif
C
	 write(luout,'(/''AFTERBURNER tuned ON!'')')
C
C	Fit multiple stars in a group with enhanced PSF using box size IRECT.
	 call mulfit( pseud2d,pseudmd,ngr,ncon,nfast,nslow,irect )        
C
C	Re-calculate aperture photometry
	 call aperm ( pseudmd, nstot, nfast, nslow )
C
	 call skyadj ( nstot )
C
 	 call title (line,skyval,.true.,fitr,ngr,ncon,strint,ztot,nums)
	 call output ( line )
	endif 
C---------------------                                                
C
C-----  This section skipped if PSF residual not written out ------
C                     
30	if( isub ) then
C
C	Write final Cleaned array.
 	 infile = 'x'//numf(1:nc)//'.fits'
	 call putfits(2,infile,header,nhead,nfast,nslow)
	 close(2)                 
C             
C	If afterburner used, then residual array also written out.
C	Find suitable scale for writing residual PSF to FITS "R" file.
C
 	 if ( wrtres ) then
	  scale=20000.0/(rmx-rmn)
	  zero=-scale*rmn   
	  do j=-nres,nres
  	   jj=nres+j+1   
	   do i=-nres,nres
  	    ii=nres+i+1                                  
	    big(ii,jj)=scale*res(i,j)+zero
	   enddo
	  enddo
	  nx=2*nres+1  
C
	  infile = 'r'//numf(1:nc)//'.fits'
  	  zer=-zero/scale
 	  scl=1.0/scale            
C                                           
C	Create a FITS header for the normalized PSF residual image    
	  call sethead(rhead,numf,nx,nx,zer,scl)
 	  scale=1.0             
	  zero=0.0    
C	Write the normalized PSF residual image
	  call putfits(2,infile,rhead,1,nx,nx)
	  close(2)
	 endif
C                        
	end if
C                     
C                     
900	close(1)
	close(3)                                 
 	close(4)
	if ( .not.screen ) close(luout)
	if(comd) then
          if(instr(5).eq.1)call system('rm shd.'//numf(1:nc))
          if(instr(6).eq.1)call system('rm out.'//numf(1:nc))
          n=1
          do while(infile(n:n).ne.' ')
            n=n+1
          end do
          if(instr(7).eq.1)call system('rm '//infile(1:n-1))
        end if
	fixed = fixedxy
	goto 1
C
995     print 996
996     format(/'*** Fatal error ***'/
     *          'You asked for batch processing but'/
     *          'I cant open the "dophot.bat" file.'/
     *          'Please make one (using batchdophot)'/
     *          'and restart DoPHOT'/)
        go to 999

C
997     print 998
998     format(/'*** Fatal error ***'/
     *          'You asked for "windowed" processing'/
     *          'but I cant open the "windows" file.'/
     *          'Please make one and restart DoPHOT'/)

999	call exit(0)
	end

* (gdb) r
* Starting program: /home3/craig/gnu/f77-e/gcc/f771 -quiet < ../../play/19990826-4.f -O
* [...]
* Breakpoint 2, fancy_abort (
*     file=0x8285220 "../../g77-e/gcc/config/i386/i386.c", line=4399,
*     function=0x82860df "output_fp_cc0_set") at ../../g77-e/gcc/rtl.c:1010
* (gdb) up
* #1  0x8222fab in output_fp_cc0_set (insn=0x8382324)
*     at ../../g77-e/gcc/config/i386/i386.c:4399
* (gdb) p insn
* $1 = 0x3a
* (gdb) up
* #2  0x8222b81 in output_float_compare (insn=0x8382324, operands=0x82acc60)
*     at ../../g77-e/gcc/config/i386/i386.c:4205
* (gdb) p insn
* $2 = 0x8382324
* (gdb) whatis insn
* type = rtx
* (gdb) pr
* (insn 2181 2180 2191 (parallel[
*             (set (cc0)
*                 (compare (reg:SF 8 %st(0))
*                     (mem:SF (plus:SI (reg:SI 6 %ebp)
*                             (const_int -9948 [0xffffd924])) 0)))
*             (clobber (reg:HI 0 %ax))
*         ] ) 29 {*cmpsf_cc_1} (insn_list 2173 (insn_list 2173 (nil)))
*     (expr_list:REG_DEAD (reg:DF 8 %st(0))
*         (expr_list:REG_UNUSED (reg:HI 0 %ax)
*             (nil))))
* (gdb)