1* Culled from 970528-1.f in Burley's g77 test suite. Copyright 2* status not clear. Feel free to chop down if the bug is still 3* reproducible (see end of test case for how bug shows up in gdb 4* run of f771). No particular reason it should be a noncompile 5* case, other than that I didn't want to spend time "fixing" it 6* to compile cleanly (with -O0, which works) while making sure the 7* ICE remained reproducible. -- burley 1999-08-26 8 9* Date: Mon, 26 May 1997 13:00:19 +0200 (GMT+0200) 10* From: "D. O'Donoghue" <dod@da.saao.ac.za> 11* To: Craig Burley <burley@gnu.ai.mit.edu> 12* Cc: fortran@gnu.ai.mit.edu 13* Subject: Re: g77 problems 14 15 program dophot 16 parameter (napple = 4) 17 common /window/nwindo,ixwin(50),iywin(50),iboxwin(50),itype(50) 18 common/io/luout,ludebg 19 common/search/nstot,thresh 20 common /fitparms / acc(npmax),alim(npmax),mit,mpar,mfit1, 21 + mfit2,ind(npmax) 22 common /starlist/ starpar(npmax,nsmax), imtype(nsmax), 23 1shadow(npmax,nsmax),shaderr(npmax,nsmax),idstr(nsmax) 24 common /aperlist/ apple(napple ,nsmax) 25 common /parpred / ava(npmax) 26 common /unitize / ufactor 27 common /undergnd/ nfast, nslow 28 common/bzero/ scale,zero 29 common /ctimes / chiimp, apertime, filltime, addtime 30 common / drfake / needit 31 common /mfit/ psfpar(npmax),starx(nfmax),stary(nfmax),xlim,ylim 32 common /vers/ version 33 logical needit,screen,isub,loop,comd,burn,wrtres,fixedxy 34 logical fixed,piped,debug,ex,clinfo 35 character header*5760,rhead*2880 36 character yn*1,version*40,ccd*4,infile*20 37 character*30 numf,odir,record*80 38 integer*2 instr(8) 39 character*800 line 40 external pseud0d, pseud2d, pseud4d, pseudmd, shape 41C 42C Initialization 43 data burn, fixedxy,fixed, piped 44 + /.false.,.false.,.false.,.false./ 45 data needit,screen,comd,isub 46 + /.true.,.false.,.true.,.false. / 47 data acc / .01, -.03, -.03, .01, .03, .1, .03 / 48 data alim / -1.0e8, 2*-1.0e3, -1.0e8, 3*-1.0e3 / 49C 50 version = 'DoPHOT Version 1.0 LINUX May 97 ' 51 debug=.false. 52 clinfo=.false. 53 line(1:800) = ' ' 54 odir = ' ' 55C 56C 57C Read default tuneable parameters 58 call tuneup ( nccd, ccd, piped, debug ) 59 version(33:36) = ccd(1:4) 60C 61 62 ludebg=6 63 if(piped)then 64 yn='n' 65 else 66 write(*,'(''****************************************'')') 67 write(*,1000) version 68 write(*,'(''****************************************''//)') 69C 70 write(*,'(''Screen output (y/[n])? '',$)') 71 read(*,1000) yn 72 end if 73 if(yn.eq.'y'.or.yn.eq.'Y') then 74 screen=.true. 75 luout=6 76 else 77 luout=2 78 end if 79C 80 if(piped)then 81 yn='y' 82 else 83 write(*,'(''Batch mode ([y]/n)? '',$)') 84 read(*,1000) yn 85 end if 86 if(yn.eq.'n'.or.yn.eq.'N') comd = .false. 87C 88 if(.not.comd) then 89 write(*, 90 * '(''Do you want windowing ([y]/n)? '',$)') 91 read(*,1000)yn 92 iwindo=1 93 if(yn.eq.'n'.or.yn.eq.'N')then 94 nwindo=0 95 iwindo=0 96 end if 97C 98 write(*, 99 * '(''Star classification info (y/[n]) ?'',$)') 100 read(*,1000)yn 101 clinfo=.false. 102 if(yn.eq.'y'.or.yn.eq.'Y')clinfo=.true. 103C 104 write(*, 105 * '(''Create a star-subtracted frame (y/[n])? '',$)') 106 read(*,1000) yn 107 if(yn.eq.'y'.or.yn.eq.'Y') isub = .true. 108C 109 write(*,'(''Apply after-burner (y/[n])? '',$)') 110 read(*,1000) yn 111 if ( yn.eq.'y'.or.yn.eq.'Y' ) burn = .true. 112 wrtres = burn 113C 114 write(*,'(''Read from fixed (X,Y) list (y/[n])? '',$)') 115 read(*,1000) yn 116 if ( yn.eq.'y'.or.yn.eq.'Y' ) then 117 fixedxy = .true. 118 fixed = .true. 119 burn = .true. 120 wrtres = .true. 121 endif 122 endif 123 iopen=0 124C 125C This is the start of the loop over the input files 126c 127 iframe=0 128 open(10,file='timing',status='unknown',access='append') 129 1301 ifit = 0 131 iapr = 0 132 itmn = 0 133 model = 1 134 xc = 0.0 135 yc = 0.0 136 rc = 0.0 137 ibr = 0 138 ixy = 0 139C 140 iframe=iframe+1 141 tgetpar=0.0 142 tsearch=0.0 143 tshape=0.0 144 timprove=0.0 145C 146C Batch mode ... 147 148 if ( comd ) then 149 if(iopen.eq.0)then 150 iopen=1 151 open(11,file='dophot.bat',status='old',err=995) 152 end if 153 read(11,1000,end=999)infile 154c now read in the parameter instructions. these are: 155c instr(1) : if 1, specifies uncrowded field, otherwise crowded 156c instr(2) : if 1, specifies sequential frames of same field 157c with a window around the stars of interest - 158c all other objects are ignored 159c instr(3) : if 0, takes cmin from dophot.inp (via tuneup) 160c if>0, sets cmin=instr(3) 161c instr(4) : if 0, does nothing 162c if 1, then opens a file called classifications 163c sets clinfo to .true. and writes out the star 164c typing info to this file 165c instr(5) : Delete the shd.nnnnnnn file 166c instr(6) : Delete the out.nnnnnnn file 167c instr(7) : Delete the input frame 168c instr(8) : Create a star-subtracted frame 169 read(11,*)instr 170 read(11,*)ifit,iapr,tmn,model,xc,yc,rc,ibr,ixy 171 nocrwd = instr(1) 172 iwindo=instr(2) 173 if(iwindo.eq.0)nwindo=0 174 itmn=tmn 175 if ( instr(3).gt.0 ) cmin=instr(3) 176 clinfo=.false. 177 if ( instr(4).gt.0 )then 178 clinfo=.true. 179 open(12,file='classifications',status='unknown') 180 ludebg=12 181 end if 182 if ( instr(8).ne.0 ) then 183 isub = .true. 184 else 185 isub = .false. 186 endif 187C 188 if(ibr.ne.0) burn = .true. 189 if(ixy.ne.0) then 190 fixedxy = .true. 191 fixed = .true. 192 burn = .true. 193 goto 20 194 endif 195 if(iwindo.eq.0)then 196 write(6,10)iframe,infile(1:15) 197 10 format(' ***** DoPHOT-ing frame ',i4,': ',a) 198 if(ludebg.eq.12)write(ludebg,11)iframe,infile(1:15) 199 11 format(////' ',62('*')/ 200 * ' * DoPHOT-ing frame ',i4,': ',a, 201 * ' *'/' ',62('*')) 202 end if 203 if(iwindo.eq.1)then 204 write(6,12)iframe,infile(1:15) 205 12 format(' ***** DoPHOT-ing frame ',i4,': ',a, 206 * ' - Windowed *****') 207 if(ludebg.eq.12)write(ludebg,13)iframe,infile(1:15) 208 13 format(////' ',62('*')/ 209 * ' * DoPHOT-ing frame ',i4,': ',a, 210 * ' - Windowed *'/2x,62('*')) 211 end if 212C 213C Interactive... 214 else 215 write(*,'(''Image name: '',$)') 216 read(*,1000) infile 217 if(infile(1:1).eq.' ') goto 999 2181000 format(a) 219 write(*,'(''Crowded field mode ([y]/n) ? '',$)') 220 read(*,1000)yn 221 nocrwd=0 222 if(yn.eq.'n'.or.yn.eq.'N')nocrwd=1 223 if(.not.fixed) then 224 write(*,1001) 2251001 format('Sky model ([1]=Plane, 2=Power, 3=Hubble)? ',$) 226 read(*,1000)record 227 if(record.ne.' ')then 228 read(record,*) model 229 else 230 model=1 231 end if 232 else 233 burn=.true. 234 goto 20 235 endif 236 endif 237C 238C if windowing, open the file and read the window 239 if(iwindo.eq.1)then 240 inquire(file='windows',exist=ex) 241 if(.not.ex)go to 997 242 if(iframe.eq.1)open(9,file='windows',status='old') 243 nwindo=0 244 2 read(9,*,end=3)intype,inx,iny,inbox 245 nwindo=nwindo+1 246 if(nwindo.gt.50)then 247 print *,'too many windows - max = 50' 248 stop 249 end if 250 ixwin(nwindo)=inx 251 iywin(nwindo)=iny 252 iboxwin(nwindo)=inbox 253 itype(nwindo)=intype 254 go to 2 255 256 3 rewind 9 257 if(screen)print 4,(itype(j),ixwin(j),iywin(j),iboxwin(j), 258 * j=1,nwindo) 259 4 format(' Windows: Type X Y Size'/ 260 * (I13,i6,i5,i5)) 261 end if 262 263 t1 = cputime(0.0) 264C 265C Read FITS frame. 266 call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line,ccd) 267C 268C Ignore frame if not the correct chip 269 if(nc.lt.0) goto 900 270C 271C Estimate starting PSF parameters. 272 15 call getparams(nfast,nslow,gxwid,gywid,skyval,tmin,tmax, 273 * iframe) 274 tgetpar = cputime(t1) + tgetpar 275 if(debug)write(ludebg,16)iframe,skyval,gxwid,gywid,tmin,tmax 276 16 format(' Getparams on frame ',i4,' sky ',f6.1,' gxwid ',f5.1, 277 * ' gywid ',f5.1,' tmin ',f5.1,' tmax ',f5.1) 278C 279C Initialize 280 do j=1,nsmax 281 imtype(j) = 0 282 do i=1,npmax 283 shadow(i,j)=0. 284 shaderr(i,j)=0. 285 enddo 286 enddo 287C 288 skyguess=skyval 289 tfac = 1.0 290C Use 4.5 X SD as fitting width 291 fitr=fitfac*(gxwid*asprat*gywid)**0.25 + 0.5 292 i=fitr 293 irect(1)=i 294 irect(2)=fitr/asprat 295C Use 4/3 X FitFac X SD as aperture width 296 gmax = asprat*gywid 297 if(gxwid.gt.gmax) gmax=gxwid 298 aprw = 1.33*fitfac*sqrt(gmax) + 0.5 299 i = aprw 300 arect(1) = i 301 i = aprw/asprat + 0.1 302 arect(2) = i 303C 304 if(irect(1).gt.50) irect(1)=50 305 if(irect(2).gt.50) irect(2)=50 306 if(arect(1).gt.45.) arect(1)=45. 307 if(arect(2).gt.45.) arect(2)=45. 308C 309 if (screen) call htype(line,skyval,.false.,fitr,ngr,ncon) 310C 311C Prompt for further information 312 if ( .not.comd ) then 313 write(*,1002) 314 1002 format(/'The above are the inital parameters DoPHOT'/ 315 * 'has found. You can change them now or accept'/ 316 * 'the values in [ ] by pressing enter'/) 317 318 write(*,1004)tmin 319 1004 format('Enter Tmin: threshold for star detection', 320 * ' [',f5.1,'] ',$) 321 read(*,1000)record 322 if(record.ne.' ')read(record,*)tmin 323 324 write(*,1005)cmin 325 1005 format('Enter Cmin: threshold for PSF stars', 326 * ' [',f5.1,'] ',$) 327 read(*,1000)record 328 if(record.ne.' ')read(record,*)cmin 329 330 write(*,1006) 331 1006 format('Do you want to fix the aperture mag size ?', 332 * ' (y/[n]) ') 333 read(*,1000)record 334 if(record.eq.'y'.or.record.eq.'Y')then 335 write(*,1007) 336 1007 format('Enter the size in pixels: ',$) 337 read(*,*)iapr 338 if(iapr.gt.0) then 339 arect(1)=iapr 340 i = iapr/asprat + 0.1 341 arect(2)=i 342 end if 343 endif 344C 345 write(*,1008) 346 1008 format('Satisfied with other input parameters ? ([y]/n)?',$) 347 read(*,1000) yn 348 if(yn.eq.'n'.or.yn.eq.'N')then 349 yn='n' 350 else 351 yn='y' 352 end if 353 if(.not.(yn.eq.'y'.or.yn.eq.'Y') ) call input 354 else 355 if ( ifit.ne.0 ) then 356 irect(1)=ifit 357 irect(2)=(ifit/asprat + 0.1) 358 endif 359 if ( iapr.ne.0 ) then 360 arect(1)=iapr 361 i = iapr/asprat + 0.1 362 arect(2)=i 363 endif 364 if ( itmn.ne.0 ) tmin = itmn 365 if ( .not.(xc.eq.0.0.and.yc.eq.0.0) ) then 366 xcen = xc 367 ycen = yc 368 endif 369 endif 370C 371C-------------------------------- 372C 373C 374 call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon, 375 +nfast, nslow ) 376C 377C if the uncrowded field option has been chosen, jump 378C straight to the minimum threshold 379C 380 if(nocrwd.eq.1)tmax=tmin 381C 382C Adjust tfac so that thresh ends precisely on Tmin. 383 if(tmin/tmax .gt. 0.999) then 384 thresh = tmin 385 tfac = 1. 386 else 387 thresh = tmax 388 xnum = alog10(tmax/tmin)/alog10(2.**tfac) 389 if(xnum.gt.1.5) then 390 xnum = float(nint(xnum)) 391 else if(xnum.ge.1) then 392 xnum = 2.0 393 else 394 xnum = 1.0 395 endif 396 tfac = alog10(tmax/tmin)/alog10(2.)/xnum 397 endif 398C 399C------------------------------------------------------------------------ 400C 401C This is the BIG LOOP which searches the frame for stars 402C with intensities > thresh. 403C 404C----------------------------------------------------------------------- 405C 406 loop = .true. 407 nstot = 0 408 do while ( loop ) 409 loop = thresh/tmin .ge. 1.01 410 write(luout,1050) thresh 4111050 format(/20('-')/'THRESHOLD: ', f10.3) 412 if(ludebg.eq.12)write(ludebg,1050) thresh 413C 414C Fit given model to sky values. 415C 416 call varipar(nstot, nfast, nslow ) 417 t1 = cputime(0.0) 418C 419C Identifies potential objects in cleaned array IMG 420 nstar = isearch( pseud2d, nfast, nslow , clinfo) 421 tsearch = cputime(t1) + tsearch 422C 423 if ( (nstar .ne. 0).or.(xnum.lt.1.5) ) then 424C 425C Performs 7-parameter PSF fit and determines nature of object. 426 t1 = cputime(0.0) 427 call shape(pseud2d,pseud4d,nfast,nslow,clinfo) 428 tshape = cputime(t1) + tshape 429C 430C Computes average sky values etc from star list 431 call paravg 432 t1 = cputime(0.0) 433C 434C Computes 4-parameter fits for all stellar objects using 435C new average shape parameters. 436 call improve(pseud2d,nfast,nslow,clinfo) 437 timprove = cputime(t1) + timprove 438 end if 439C 440C Calculate aperture photometry on last pass. 441 if(.not.loop) call aper ( pseud2d, nstot, nfast, nslow ) 442C 443 totaltime = (tgetpar+tsearch+tshape+timprove) 444 write(3,1060) totaltime 445 write(4,1060) totaltime 446 write(luout,1060) totaltime 4471060 format('Total CPU time consumed:',F10.2,' seconds.') 448 write(10,1070)infile,tgetpar,tsearch,tshape,timprove, 449 * totaltime 4501070 format(a20,' T(getp/f)',f5.1,' T(search)',f5.1, 451 * ' T(shape)',f5.1,' T(improve)',f5.1, 452 * ' Total',f6.1) 453 call title (line,skyval,.false.,fitr,ngr,ncon,strint,ztot,nums) 454 rewind(2) 455 rewind(3) 456 rewind(4) 457C 458 call output ( line ) 459C 460C Now reduce the threshold and loop back 461C 462 thresh = thresh/2.**tfac 463 end do 464C 465C--------- END OF BIG LOOP --------------------------------------- 466C 467C If after-burner required, residuals from analytic PSF are computed 468C and stored in RES. 469C 47020 if ( burn ) then 471C 472C If using a fixed (X,Y) coordinate list, read it. 473 if (fixed) then 474C Read the image frame 475 call getfits(1,infile,header,nhead,nfast,nslow,numf,nc,line) 476C 477C Initialize arrays, open files etc. 478 call setup ( numf,nc,screen,line,skyval,fitr,ngr,ncon, 479 +nfast, nslow ) 480C 481C Read the XY list 482 write(luout,'(''Reading XY list ...'')') 483 call xylist(numf, nc, ios ) 484 if(ios.ne.0) then 485 fixed = .false. 486 write(luout,'(''SXY file absent or incorrect...'')') 487 goto 15 488 endif 489C 490 call htype(line,skyval,.false.,fitr,ngr,ncon) 491C 492C Remove good stars 493 write(luout,'(''Cleaning frame of stars: '',i8)') nstot 494 call clean ( pseud2d, nstot, nfast, nslow, -1) 495C 496C Calculate aperture photometry 497C call aper ( pseud2d, nstot, nfast, nslow ) 498 else 499 rewind(3) 500 rewind(4) 501 endif 502C 503C----------------------- 504C Flag all stars close together in groups. Keep making the distance 505C criterion FITR smaller until the maximum number in a group is less 506C than NFMAX 507C 508 fitr = amax1(arect(1),arect(2)) 509 fitr = fitr + 2.0 510 nmax = 10000 511 write(*,'(''Regrouping ...'')') 512C 513 do while ( nmax.gt.nfmax ) 514 fitr = fitr - 1.0 515 write(luout,'(''Min distance ='',f8.1)') fitr 516 call regroup( fitr, ngr, nmax ) 517 enddo 518C 519 xlim = irect(1)/2 520 ylim = irect(2)/2 521C 522C Calculate normalized PSF residual from PSEUD2D 523 call getres (pseud0d,pseud2d,strint,rmn,rmx,nfast,nslow,irect, 524 +arect,ztot,nums) 525 if(nums.eq.0) then 526 write(luout,'(''No suitable PSF stars!'')') 527 goto 30 528 endif 529C 530 write(luout,'(/''AFTERBURNER tuned ON!'')') 531C 532C Fit multiple stars in a group with enhanced PSF using box size IRECT. 533 call mulfit( pseud2d,pseudmd,ngr,ncon,nfast,nslow,irect ) 534C 535C Re-calculate aperture photometry 536 call aperm ( pseudmd, nstot, nfast, nslow ) 537C 538 call skyadj ( nstot ) 539C 540 call title (line,skyval,.true.,fitr,ngr,ncon,strint,ztot,nums) 541 call output ( line ) 542 endif 543C--------------------- 544C 545C----- This section skipped if PSF residual not written out ------ 546C 54730 if( isub ) then 548C 549C Write final Cleaned array. 550 infile = 'x'//numf(1:nc)//'.fits' 551 call putfits(2,infile,header,nhead,nfast,nslow) 552 close(2) 553C 554C If afterburner used, then residual array also written out. 555C Find suitable scale for writing residual PSF to FITS "R" file. 556C 557 if ( wrtres ) then 558 scale=20000.0/(rmx-rmn) 559 zero=-scale*rmn 560 do j=-nres,nres 561 jj=nres+j+1 562 do i=-nres,nres 563 ii=nres+i+1 564 big(ii,jj)=scale*res(i,j)+zero 565 enddo 566 enddo 567 nx=2*nres+1 568C 569 infile = 'r'//numf(1:nc)//'.fits' 570 zer=-zero/scale 571 scl=1.0/scale 572C 573C Create a FITS header for the normalized PSF residual image 574 call sethead(rhead,numf,nx,nx,zer,scl) 575 scale=1.0 576 zero=0.0 577C Write the normalized PSF residual image 578 call putfits(2,infile,rhead,1,nx,nx) 579 close(2) 580 endif 581C 582 end if 583C 584C 585900 close(1) 586 close(3) 587 close(4) 588 if ( .not.screen ) close(luout) 589 if(comd) then 590 if(instr(5).eq.1)call system('rm shd.'//numf(1:nc)) 591 if(instr(6).eq.1)call system('rm out.'//numf(1:nc)) 592 n=1 593 do while(infile(n:n).ne.' ') 594 n=n+1 595 end do 596 if(instr(7).eq.1)call system('rm '//infile(1:n-1)) 597 end if 598 fixed = fixedxy 599 goto 1 600C 601995 print 996 602996 format(/'*** Fatal error ***'/ 603 * 'You asked for batch processing but'/ 604 * 'I cant open the "dophot.bat" file.'/ 605 * 'Please make one (using batchdophot)'/ 606 * 'and restart DoPHOT'/) 607 go to 999 608 609C 610997 print 998 611998 format(/'*** Fatal error ***'/ 612 * 'You asked for "windowed" processing'/ 613 * 'but I cant open the "windows" file.'/ 614 * 'Please make one and restart DoPHOT'/) 615 616999 call exit(0) 617 end 618 619* (gdb) r 620* Starting program: /home3/craig/gnu/f77-e/gcc/f771 -quiet < ../../play/19990826-4.f -O 621* [...] 622* Breakpoint 2, fancy_abort ( 623* file=0x8285220 "../../g77-e/gcc/config/i386/i386.c", line=4399, 624* function=0x82860df "output_fp_cc0_set") at ../../g77-e/gcc/rtl.c:1010 625* (gdb) up 626* #1 0x8222fab in output_fp_cc0_set (insn=0x8382324) 627* at ../../g77-e/gcc/config/i386/i386.c:4399 628* (gdb) p insn 629* $1 = 0x3a 630* (gdb) up 631* #2 0x8222b81 in output_float_compare (insn=0x8382324, operands=0x82acc60) 632* at ../../g77-e/gcc/config/i386/i386.c:4205 633* (gdb) p insn 634* $2 = 0x8382324 635* (gdb) whatis insn 636* type = rtx 637* (gdb) pr 638* (insn 2181 2180 2191 (parallel[ 639* (set (cc0) 640* (compare (reg:SF 8 %st(0)) 641* (mem:SF (plus:SI (reg:SI 6 %ebp) 642* (const_int -9948 [0xffffd924])) 0))) 643* (clobber (reg:HI 0 %ax)) 644* ] ) 29 {*cmpsf_cc_1} (insn_list 2173 (insn_list 2173 (nil))) 645* (expr_list:REG_DEAD (reg:DF 8 %st(0)) 646* (expr_list:REG_UNUSED (reg:HI 0 %ax) 647* (nil)))) 648* (gdb) 649