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