1 PROGRAM COLFIT 2 IMPLICIT NONE 3C 4C Program to do a nonlinear least squares fit of a two-column file in 5C the format 6C 7C x1 y1 8C x2 y2 9C .. .. 10C xN yN 11C 12C to a formula of the form y(x) = expression(x,a,b,c,...), where 13C a, b, c, ... are parameters to be found. 14C----------------------------------------------------------------------- 15C Common block XYDATA holds the data read from disk 16C 17 INTEGER NPMAX 18 PARAMETER ( NPMAX = 32768 ) 19 INTEGER NUMPTS 20 REAL*8 XX(NPMAX) , YY(NPMAX) 21 COMMON /XYDATA/ XX , YY , NUMPTS 22C....................................................................... 23C Common block EXDATA holds the parsed expression for evaluation 24C 25 INTEGER NCMAX 26 PARAMETER ( NCMAX = 512 ) 27 CHARACTER*8 CCODE(NCMAX) 28 INTEGER NUMCOD , IVEVAR,NUMPAR,IFXVAR(26) 29 COMMON /EXDATA/ CCODE , NUMCOD , IVEVAR,NUMPAR,IFXVAR 30C....................................................................... 31C local variables 32C 33 CHARACTER*64 CFILE 34 CHARACTER*128 CEXPR 35 CHARACTER*1 CHANS 36C 37 INTEGER I , IERR , ICH , NUMCH , IUPPER,ILOWER , 38 X IOPT,M,N,NPRINT,INFO , ISTART,IEND 39 REAL*8 XIN , YIN , TOL 40 LOGICAL LCHUSE(26) 41C 42 INTEGER NWORK 43 PARAMETER ( NWORK = 30*NPMAX ) 44 INTEGER IWORK(66) 45 REAL*8 WORK(NPMAX,30) , XPAR(66) , FVEC(NPMAX) , 46 X WAR1(30),WAR2(30),WAR3(30),WAR4(NPMAX) 47C 48 EXTERNAL FF , DNLS1E , DCOV 49C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 50C 51C Open data file 52C 53100 CONTINUE 54 WRITE(*,101) 55101 FORMAT(' Enter data file name : '$) 56 READ(*,111) CFILE 57111 FORMAT(A) 58 IF( CFILE(1:1) .EQ. ' ' )GOTO 9900 59C 60 OPEN( UNIT = 1 , 61 X FILE = CFILE , 62 X FORM = 'FORMATTED' , 63 X STATUS = 'OLD' , 64 X IOSTAT = IERR ) 65C 66 IF( IERR .NE. 0 )THEN 67 CLOSE( UNIT = 1 ) 68 GOTO 100 69 ENDIF 70C....................................................................... 71C Read data in 72C 73 WRITE(*,102) 74102 FORMAT(' Start, stop indices in data file [0,0=all] ? '$) 75 READ(*,112) ISTART,IEND 76112 FORMAT(2I10) 77C 78 DO 110 I=1,ISTART-1 79 READ(1,*,IOSTAT=IERR) XIN,YIN 80 IF( IERR .NE. 0 )THEN 81 WRITE(*,119) I 82119 FORMAT(' *** failure when skipping record ',I5) 83 GOTO 9900 84 ENDIF 85110 CONTINUE 86C 87 IF( IEND .LE. ISTART )THEN 88 IEND = NPMAX 89 ELSE 90 IEND = MIN( NPMAX , IEND-ISTART+1 ) 91 ENDIF 92C 93 I = 0 94200 CONTINUE 95 READ(1,*,IOSTAT=IERR) XIN,YIN 96 IF( IERR .EQ. 0 )THEN 97 I = I + 1 98 XX(I) = XIN 99 YY(I) = YIN 100 IF( I .LT. IEND )GOTO 200 101 ENDIF 102C 103 CLOSE( UNIT = 1 ) 104 NUMPTS = I 105 WRITE(*,201) NUMPTS 106201 FORMAT(' --- read',I6,' data lines from file') 107 IF( NUMPTS .LE. 1 )GOTO 9900 108C....................................................................... 109C Read in expression that will be used to fit 110C 111300 CONTINUE 112 WRITE(*,301) 113301 FORMAT(/' Expression to fit 2nd column to 1st ? '$) 114 READ(*,111,IOSTAT=IERR) CEXPR 115 IF( IERR .NE. 0 )GOTO 9900 116 CALL PARSER( CEXPR , .TRUE. , NUMCOD , CCODE ) 117 IF( NUMCOD .LE. 0 )GOTO 300 118 IF( NUMCOD .GT. NCMAX )THEN 119 WRITE(*,309) 120309 FORMAT(' *** expression too complex!' ) 121 GOTO 9900 122 ENDIF 123C....................................................................... 124C Get the names of all the variables referred to in the expression 125C 126 DO 310 I=1,26 127 LCHUSE(I) = .FALSE. 128310 CONTINUE 129C 130 IUPPER = ICHAR('A') - 1 131 ILOWER = ICHAR('a') - 1 132 DO 320 I=1,NUMCOD 133 IF( CCODE(I) .EQ. 'PUSHSYM' )THEN 134 ICH = ICHAR(CCODE(I+1)(1:1)) - IUPPER 135 LCHUSE(ICH) = .TRUE. 136 ENDIF 137320 CONTINUE 138C 139 CFILE = ' ' 140 NUMCH = 0 141 DO 330 I=1,26 142 IF( LCHUSE(I) )THEN 143 NUMCH = NUMCH + 1 144 CFILE(2*NUMCH:2*NUMCH+1) = CHAR(I+IUPPER) // ' ' 145 ENDIF 146330 CONTINUE 147C 148 IF( NUMCH .LE. 1 )THEN 149 WRITE(*,338) 150338 FORMAT(' *** not enough variable names in expression!') 151 GOTO 300 152 ENDIF 153 IF( NUMCH .GT. NUMPTS )THEN 154 WRITE(*,339) 155339 FORMAT(' *** too many variable names in expression!') 156 GOTO 300 157 ENDIF 158C....................................................................... 159C Find out which variable is the 1st data column 160C 161340 CONTINUE 162 WRITE(*,341) CFILE(1:2*NUMCH) 163341 FORMAT( 164 X ' Which variable represents 1st data column ? ' / 1X,A,' : '$) 165 READ(*,111) CHANS 166 IF( CHANS.GE.'A' .AND. CHANS.LE.'Z' )THEN 167 ICH = ICHAR(CHANS) - IUPPER 168 ELSEIF( CHANS.GE.'a' .AND. CHANS.LE.'z' )THEN 169 ICH = ICHAR(CHANS) - ILOWER 170 ELSE 171 GOTO 340 172 ENDIF 173 IF( .NOT. LCHUSE(ICH) )GOTO 340 174C....................................................................... 175C Assign 1st column variable name to VEctor VARiable, and 176C all others to be FiXed Variables; get initial values for the latter 177C 178 IVEVAR = ICH 179 NUMPAR = NUMCH - 1 180 ICH = 0 181 DO 360 I=1,26 182 IF( LCHUSE(I) .AND. I.NE.IVEVAR )THEN 183 ICH = ICH + 1 184 IFXVAR(ICH) = I 185C 186350 CONTINUE 187 WRITE(*,351) CHAR(I+IUPPER) 188351 FORMAT(' Initial value for parameter ',A,' ? '$) 189 READ(*,*,IOSTAT=IERR) XPAR(ICH) 190 IF( IERR .NE. 0 )GOTO 350 191 ENDIF 192360 CONTINUE 193C....................................................................... 194C Call the fitting routine 195C 196 IOPT = 1 197 M = NUMPTS 198 N = NUMPAR 199 TOL = 1.D-06 200 NPRINT = 0 201C 202C ************ 203 CALL DNLS1E( FF,IOPT , M,N , XPAR,FVEC , TOL,NPRINT , INFO , 204 X IWORK,WORK,NWORK ) 205C ************ 206C 207C Write results out: 208C 209 WRITE(*,501) 'DNLS1E' , INFO 210501 FORMAT(' --- ',A,' exits with INFO =',I2) 211 IF( INFO .EQ. 0 )GOTO 9900 212C 213 CALL DCOV( FF,IOPT , M,N , XPAR,FVEC , WORK,NPMAX , INFO , 214 X WAR1,WAR2,WAR3,WAR4 ) 215C 216 WRITE(*,501) 'DCOV' , INFO 217 IF( INFO .EQ. 0 )GOTO 9900 218C 219 DO 510 I=1,NUMPAR 220 IF( INFO .EQ. 1 )THEN 221 XIN = SQRT( WORK(I,I) ) 222 WRITE(*,511) CHAR( IFXVAR(I)+IUPPER ) , XPAR(I) , XIN 223511 FORMAT(3X,A,' -> ',1PG14.7,' +/- ',1PG14.7) 224 ELSE 225 WRITE(*,511) CHAR( IFXVAR(I)+IUPPER ) , XPAR(I) 226 ENDIF 227510 CONTINUE 228C 229 XIN = 0.D+00 230 DO 520 I=1,NUMPTS 231 XIN = XIN + ABS( FVEC(I) ) 232520 CONTINUE 233 XIN = XIN / NUMPTS 234 WRITE(*,521) XIN 235521 FORMAT(' --- mean absolute deviation =',1PG10.3) 236C....................................................................... 237C Write fit error curve to a file (if desired) 238C 239600 CONTINUE 240 WRITE(*,601) 241601 FORMAT(' Filename to save fit error curve in (blank=none) ? '$) 242 READ(*,111) CFILE 243 IF( CFILE(1:1) .NE. ' ' )THEN 244 OPEN( UNIT = 1 , 245 X FILE = CFILE , 246 X FORM = 'FORMATTED' , 247 X STATUS = 'UNKNOWN' , 248 X IOSTAT = IERR ) 249 IF( IERR .NE. 0 )GOTO 600 250C 251 DO 610 I=1,NUMPTS 252 WRITE(1,611) XX(I) , FVEC(I) 253611 FORMAT(1PG20.13,1X,1PG20.13) 254610 CONTINUE 255C 256 CLOSE( UNIT = 1 ) 257 ENDIF 258C----------------------------------------------------------------------- 2599900 CONTINUE 260 END 261C 262C 263C 264 SUBROUTINE FF( IFLAG , M , N , X , FVEC , FJAC , LDFJAC ) 265 IMPLICIT NONE 266C 267C Routine supplied to DNLS1E to compute the functions we are 268C are trying to fit. 269C 270 INTEGER IFLAG , LDFJAC , M , N 271 REAL*8 X(N) , FVEC(M) , FJAC 272C....................................................................... 273C Common block XYDATA holds the data read from disk 274C 275 INTEGER NPMAX 276 PARAMETER ( NPMAX = 32768 ) 277 INTEGER NUMPTS 278 REAL*8 XX(NPMAX) , YY(NPMAX) 279 COMMON /XYDATA/ XX , YY , NUMPTS 280C....................................................................... 281C Common block EXDATA holds the parsed expression for evaluation 282C 283 INTEGER NCMAX 284 PARAMETER ( NCMAX = 512 ) 285 CHARACTER*8 CCODE(NCMAX) 286 INTEGER NUMCOD , IVEVAR,NUMPAR,IFXVAR(26) 287 COMMON /EXDATA/ CCODE , NUMCOD , IVEVAR,NUMPAR,IFXVAR 288C....................................................................... 289C Local variables 290C 291 INTEGER I , J , K 292 REAL*8 VAZ(NPMAX,26) 293C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 294C 295 IF( IFLAG .EQ. 0 )THEN 296 WRITE(*,101) X 297101 FORMAT(' --- FF X=',5(1X,1PG12.5) ) 298C....................................................................... 299 ELSEIF( IFLAG .EQ. 1 )THEN 300 DO 110 I=1,N 301 K = IFXVAR(I) 302 DO 105 J=1,M 303 VAZ(J,K) = X(I) 304105 CONTINUE 305110 CONTINUE 306 DO 115 J=1,M 307 VAZ(J,IVEVAR) = XX(J) 308115 CONTINUE 309 CALL PAREVEC( NUMCOD , CCODE , 310 X VAZ(1, 1) , VAZ(1, 2) , VAZ(1, 3) , VAZ(1, 4) , 311 X VAZ(1, 5) , VAZ(1, 6) , VAZ(1, 7) , VAZ(1, 8) , 312 X VAZ(1, 9) , VAZ(1,10) , 313 X VAZ(1,11) , VAZ(1,12) , VAZ(1,13) , VAZ(1,14) , 314 X VAZ(1,15) , VAZ(1,16) , VAZ(1,17) , VAZ(1,18) , 315 X VAZ(1,19) , VAZ(1,20) , 316 X VAZ(1,21) , VAZ(1,22) , VAZ(1,23) , VAZ(1,24) , 317 X VAZ(1,25) , VAZ(1,26) , 318 X M , FVEC ) 319 DO 120 I=1,M 320 FVEC(I) = FVEC(I) - YY(I) 321120 CONTINUE 322C....................................................................... 323 ELSE 324 WRITE(*,999) IFLAG 325999 FORMAT(' *** FF has illegal IFLAG=',I5) 326 STOP 327 ENDIF 328C 329 RETURN 330 END 331