-
Notifications
You must be signed in to change notification settings - Fork 1
/
nllssrr.f
610 lines (607 loc) · 25.9 KB
/
nllssrr.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
c***********************************************************************
SUBROUTINE NLLSSRR(NDATA,NPTOT,NPMAX,CYCMAX,IROUND,ROBUST,LPRINT,
1 IFXP,YO,YU,YD,PV,PU,PS,CM,TSTPS,TSTPU,DSE)
c** Program for performing linear or non-linear least-squares fits and
c (if desired) automatically using sequential rounding and refitting
c to minimize the numbers of parameter digits which must be quoted [see
c R.J. Le Roy, J.Mol.Spectrosc. 191, 223-231 (1998)].
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c** Program uses orthogonal decomposition of the "design" (partial
c derivative) matrix for the core locally linear (steepest descent)
c step, following a method introduced (to me) by Dr. Michael Dulick.
c** If no parameters are free, simply return RMS(residuals) as
c calculated from the input parameter values {PV(j)}.
c** A user MUST SUPPLY subroutine DYIDPJ to generate the predicted
c value of each datum and the partial derivatives of each datum w.r.t.
c each parameter (see below) from the current trial parameters.
c
c** On entry:
c NDATA is the number of data to be fitted
c NPTOT the total number of parameters in the model (.le.NPMAX).
c If NPTOT.le.0 , assume YD(i)=YO(i) and calculate the (RMS
c dimensionless deviation)=DSE from them & YU(i)
c NPMAX is the maximum number of model parameters allowed by current
c external array sizes. Should set internal NPINTMX = NPMAX
c (may be freely changed by the user).
c CYCMAX is the upper bound on the allowed number of iterative cycles
c IROUND .ne. 0 causes Sequential Rounding & Refitting to be
c performed, with each parameter being rounded at the
c |IROUND|'th sig. digit of its local incertainty.
c > 0 rounding selects in turn remaining parameter with largest
c relative uncertainy
c < 0 round parameters sequentially from last to first
c = 0 simply stops after full convergence (without rounding).
c ROBUST > 0 causes fits to use Watson's ``robust'' weighting
c 1/[u^2 +{(c-o)^2}/3]. ROBUST > 1 uses normal 1/u^2 on first
c fit cycle and 'robust' on later cycles.
c LPRINT specifies the level of printing inside NLLSSRR
c if: = 0, no print except for failed convergence.
c .NE.0 also print DRMSD and convergence tests on each cycle
c and indicate nature of convergence
c >= 1 also parameters changes & uncertainties, each cycle
c >= 2 also print parameter change each rounding step
c IFXP(j) specifies whether parameter j is to be held fixed
c [IFXP > 0] or to be freely varied in the fit [IFXP= 0]
c YO(i) are the NDATA 'observed' data to be fitted
c YU(i) are the uncertainties in these YO(i) values
c PV(j) are initial trial parameter values (for non-linear fits);
c should be set at zero for initially undefined parameters.
c
c** On Exit:
c YD(i) is the array of differences [Ycalc(i) - YO(i)]
c PV(j) are the final converged parameter values
c PU(j) are 95% confidence limit uncertainties in the PV(j)'s
c PS(j) are 'parameter sensitivities' for the PV(j)'s, defined such
c that the RMS displacement of predicted data due to rounding
c off parameter-j by PS(j) is .le. DSE/10*NPTOT
c CM(j,k) is the correlation matrix obtained by normalizing variance
c /covariance matrix: CM(j,k) = CM(j,k)/SQRT[CM(j,j)*CM(k,k)]
c TSTPS = max{|delta[PV(j)]/PS(j)|} is the parameter sensitivity
c convergence test: delta[PV(j)] is last change in parameter-j
c TSTPU = max{|delta[PV(j)]/PU(j)|} is the parameter uncertainty
c convergence test: delta[PV(j)] is last change in parameter-j
c DSE is the predicted (dimensionless) standard error of the fit
c
c NOTE that the squared 95% confidence limit uncertainty in a property
c F({PV(j)}) defined in terms of the fitted parameters {PV(j)} (where
c the L.H.S. involves [row]*[matrix]*[column] multiplication) is:
c [D(F)]^2 = [PU(1)*dF/dPV(1), PU(2)*dF/dPV(2), ...]*[CM(j,k)]*
c [PU(2)*dF/dPV(1), PU(2)*dF/dPV(2), ...]
c
c** Externally dimension: YO, YU and YD .ge. NDATA
c PV, PU and PS .ge. NPTOT (say as NPMAX),
c CM as a square matrix with column & row length NPMAX
c***********************************************************************
INTEGER NPINTMX
PARAMETER (NPINTMX=8000)
INTEGER I,J,K,L,IDF,ITER,NITER,CYCMAX,IROUND,ISCAL,JROUND,LPRINT,
1 NDATA,NPTOT,NPMAX,NPARM,NPFIT,JFIX,QUIT,ROBUST,
2 IFXP(NPMAX),JFXP(NPINTMX)
REAL*8 YO(NDATA), YU(NDATA), YD(NDATA), PV(NPTOT), PU(NPTOT),
1 PS(NPTOT),PSS(NPINTMX),PC(NPINTMX),PX(NPINTMX),
2 PY(NPINTMX),CM(NPMAX,NPMAX), F95(10),
3 RMSR, RMSRB, DSE, TSTPS, TSTPSB, TSTPU, TFACT, S, YC, Zthrd
DATA F95/12.7062D0,4.3027D0,3.1824D0,2.7764D0,2.5706D0,2.4469D0,
1 2.3646D0,2.3060D0,2.2622D0,2.2281D0/
IF((NPTOT.GT.NPMAX).OR.(NPTOT.GT.NPINTMX).OR.(NPINTMX.NE.NPMAX)
1 .OR.(NPTOT.GT.NDATA)) THEN
c** If array dimensioning inadequate, print warning & then STOP
WRITE(6,602) NPTOT,NPINTMX,NPMAX,NDATA
STOP
ENDIF
Zthrd= 0.d0
IF(ROBUST.GE.2) Zthrd= 1.d0/3.d0
TSTPS= 0.d0
RMSR= 0.d0
NITER= 0
QUIT= 0
NPARM= NPTOT
DO J= 1, NPTOT
PS(J)= 0.d0
JFXP(J)= IFXP(J)
IF(IFXP(J).GT.0) NPARM= NPARM- 1
ENDDO
NPFIT= NPARM
JROUND= IABS(IROUND)
c=======================================================================
c** Beginning of loop to perform rounding (if desired). NOTE that in
c sequential rounding, NPARM is the current (iteratively shrinking)
c number of free parameters.
6 IF(NPARM.GT.0) TSTPS= 9.d99
c** TFACT is 95% student t-value for (NDATA-NPARM) degrees of freedom.
c [Approximate expression for (NDATA-NPARM).GT.10 accurate to ca. 0.002]
TFACT= 0.D0
IF(NDATA.GT.NPARM) THEN
IDF= NDATA-NPARM
IF(IDF.GT.10) TFACT= 1.960D0*DEXP(1.265D0/DBLE(IDF))
IF(IDF.LE.10) TFACT= F95(IDF)
ELSE
TFACT= 0.D0
ENDIF
c======================================================================
c** Begin iterative convergence loop: try for up to 30 cycles
DO 50 ITER= 1, CYCMAX
ISCAL= 0
NITER= NITER+ 1
DSE= 0.d0
TSTPSB= TSTPS
RMSRB= RMSR
c** Zero out various arrays
IF(NPARM.GT.0) THEN
DO I = 1,NPARM
c** PSS is the array of Saved Parameter Sensitivities from previous
c iteration to be carried into dyidpj subroutine - used in predicting
c increment for derivatives by differences.
PSS(I)= PS(I)
PS(I) = 0.D0
PU(I) = 0.D0
PX(I) = 0.D0
PY(I) = 0.D0
DO J = 1,NPARM
CM(I,J) = 0.D0
ENDDO
ENDDO
ENDIF
c
c========Beginning of core linear least-squares step====================
c
c** Begin by forming the Jacobian Matrix from partial derivative matrix
DO I = 1,NDATA
c** User-supplied subroutine DYIDPJ uses current (trial) parameter
c values {PV} to generate predicted datum # I [y(calc;I)=YC] and its
c partial derivatives w.r.t. each of the parameters, returning the
c latter in 1-D array PC. See dummy sample version at end of listing.
c* NOTE 1: if more convenient, DYIDPJ could prepare the y(calc) values
c and derivatives for all data at the same time (when I=1), but only
c returned the values here one datum at a time (for I > 1).]
c* NOTE 2: the partial derivative array PC returned by DYIDPJ must have
c an entry for every parameter in the model, though for parameters
c which are held fixed [JFXP(j)=1], those PC(j) values are ignored.
CALL DYIDPJ(I,NPTOT,YO(I),YC,PV,PC)
IF(NPARM.LT.NPTOT) THEN
c** For constrained parameter or sequential rounding, collapse partial
c derivative array here
DO J= NPTOT,1,-1
IF(JFXP(J).GT.0) THEN
c!! First ... move derivative for constrained-parameter case
cc666 FORMAT(' For IDAT=',I5,' add PC(',I3,') =',1pD15.8,
cc 1 ' to PC(',0pI3,') =',1pD15.8)
IF(JFXP(J).GT.1) THEN
cc write(6,666) I,J,PC(J),JFXP(J),PC(JFXP(J))
PC(JFXP(J))= PC(JFXP(J))+ PC(J)
ENDIF
c ... now continue collapsing partial derivative array
IF(J.LT.NPTOT) THEN
DO K= J,NPTOT-1
PC(K)= PC(K+1)
ENDDO
ENDIF
PC(NPTOT)= 0.d0
ENDIF
ENDDO
ENDIF
YD(I)= YC - YO(I)
S = 1.D0/YU(I)
cc *** For 'Robust' fitting, adjust uncertainties here
IF(Zthrd.GT.0.d0) S= 1.d0/DSQRT(YU(I)**2 + Zthrd*YD(I)**2)
YC= -YD(I)*S
DSE= DSE+ YC*YC
IF(NPARM.GT.0) THEN
DO J = 1,NPARM
PC(J) = PC(J)*S
PS(J) = PS(J)+ PC(J)**2
ENDDO
CALL QROD(NPARM,NPMAX,NPMAX,CM,PC,PU,YC,PX,PY)
ENDIF
ENDDO
RMSR= DSQRT(DSE/NDATA)
IF(NPARM.LE.0) GO TO 60
c
c** Compute the inverse of CM
CM(1,1) = 1.D0 / CM(1,1)
DO I = 2,NPARM
L = I - 1
DO J = 1,L
S = 0.D0
DO K = J,L
S = S + CM(K,I) * CM(J,K)
ENDDO
CM(J,I) = -S / CM(I,I)
ENDDO
CM(I,I) = 1.D0 / CM(I,I)
ENDDO
c
c** Solve for parameter changes PC(j)
DO I = 1,NPARM
J = NPARM - I + 1
PC(J) = 0.D0
DO K = J,NPARM
PC(J) = PC(J) + CM(J,K) * PU(K)
ENDDO
ENDDO
c
c** Get (upper triangular) "dispersion Matrix" [variance-covarience
c matrix without the sigma^2 factor].
DO I = 1,NPARM
DO J = I,NPARM
YC = 0.D0
DO K = J,NPARM
YC = YC + CM(I,K) * CM(J,K)
ENDDO
CM(I,J) = YC
ENDDO
ENDDO
c** Generate core of Parameter Uncertainties PU(j) and (symmetric)
c correlation matrix CM
DO J = 1,NPARM
PU(J) = DSQRT(CM(J,J))
DO K= J,NPARM
CM(J,K)= CM(J,K)/PU(J)
ENDDO
DO K= 1,J
CM(K,J)= CM(K,J)/PU(J)
CM(J,K)= CM(K,J)
ENDDO
ENDDO
c
c** Generate standard error DSE = sigma^2, and prepare to calculate
c Parameter Sensitivities PS
IF(NDATA.GT.NPARM) THEN
DSE= DSQRT(DSE/(NDATA-NPARM))
ELSE
DSE= 0.d0
ENDIF
c** Use DSE to get final (95% confid. limit) parameter uncertainties PU
c** Calculate 'parameter sensitivities', changes in PV(j) which would
c change predictions of input data by an RMS average of DSE*0.1/NPARM
YC= DSE*0.1d0/DBLE(NPARM)
S= DSE*TFACT
DO J = 1,NPARM
PU(J)= S* PU(J)
PS(J)= YC*DSQRT(NDATA/PS(J))
ENDDO
c========End of core linear least-squares step==========================
c ... early exit if Rounding cycle finished ...
IF(QUIT.GT.0) GO TO 54
c
c** Next test for convergence
TSTPS= 0.D0
TSTPU= 0.D0
DO J= 1, NPARM
TSTPS= MAX(TSTPS,DABS(PC(J)/PS(J)))
TSTPU= MAX(TSTPU,DABS(PC(J)/PU(J)))
ENDDO
IF(LPRINT.NE.0) WRITE(6,604) ITER,RMSR,TSTPS,TSTPU
c** Now ... update parameters (careful about rounding)
DO J= 1,NPTOT
IF(JFXP(J).GT.0) THEN
IF(JFXP(J).GT.1) THEN
c** If this parameter constrained to equal some earlier parameter ....
PV(J)= PV(JFXP(J))
WRITE(6,668) J,JFXP(J),PV(J),ITER
ENDIF
668 FORMAT(' Constrain PV(',i3,') = PV(',I3,') =',1pd15.8,
1 ' on cycle',i3)
c** If parameter held fixed (by input or rounding process), shift values
c of change, sensitivity & uncertainty to correct label.
IF(J.LT.NPTOT) THEN
DO I= NPTOT,J+1,-1
PC(I)= PC(I-1)
PS(I)= PS(I-1)
PU(I)= PU(I-1)
ENDDO
ENDIF
PC(J)= 0.d0
PS(J)= 0.d0
PU(J)= 0.d0
ELSE
PV(J)= PV(J)+ PC(J)
ENDIF
ENDDO
IF(LPRINT.GE.1) WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),
1 J=1,NPTOT)
IF(ITER.GT.1) THEN
c** New Convergence test is to require RMSD to be constant to 1 part in
c 10^7 in adjacent cycles (unlikely to occur by accident)
IF(ABS((RMSR/RMSRB)-1.d0).LT.1.d-07) THEN
IF(LPRINT.NE.0) WRITE(6,607) ITER,
1 ABS(RMSR/RMSRB-1.d0),TSTPS
GO TO 54
ENDIF
ENDIF
IF(ROBUST.GT.0) Zthrd= 1.d0/3.d0
50 CONTINUE
WRITE(6,610) NPARM,NDATA,ITER,RMSR,TSTPS,TSTPU
c** End of iterative convergence loop for (in general) non-linear case.
c======================================================================
c
54 IF(NPARM.LT.NPTOT) THEN
c** If necessary, redistribute correlation matrix elements to full
c NPTOT-element correlation matrix
DO J= 1,NPTOT
IF(JFXP(J).GT.0) THEN
c* If parameter J was held fixed
IF(J.LT.NPTOT) THEN
c ... then move every lower CM element down one row:
DO I= NPTOT,J+1,-1
c ... For K < J, just shift down or over to the right
IF(J.GT.1) THEN
DO K= 1,J-1
CM(I,K)= CM(I-1,K)
CM(K,I)= CM(I,K)
ENDDO
ENDIF
c ... while for K > J also shift elements one column to the right
DO K= NPTOT,J+1,-1
CM(I,K)= CM(I-1,K-1)
ENDDO
ENDDO
ENDIF
c ... and finally, insert appropriate row/column of zeros ....
DO I= 1,NPTOT
CM(I,J)= 0.d0
CM(J,I)= 0.d0
ENDDO
CM(J,J)= 1.d0
ENDIF
ENDDO
ENDIF
IF(QUIT.GT.0) GOTO 60
IF(NPARM.EQ.NPFIT) THEN
c** If desired, print unrounded parameters and fit properties
IF(LPRINT.NE.0) THEN
WRITE(6,616) NDATA,NPARM,RMSR,TSTPS
WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),J=1,NPTOT)
ENDIF
ENDIF
IF(IROUND.EQ.0) RETURN
c** Automated 'Sequential Rounding and Refitting' section: round
c selected parameter, fix it, and return (above) to repeat fit.
IF(IROUND.LT.0) THEN
c ... if IROUND < 0, sequentially round off 'last' remaining parameter
DO J= 1, NPTOT
IF(JFXP(J).LE.0) THEN
JFIX= J
ENDIF
ENDDO
ELSE
c ... if IROUND > 0, sequentially round off remaining parameter with
c largest relative uncertainty.
c ... First, select parameter JFIX with the largest relative uncertainty
JFIX= NPTOT
K= 0
TSTPS= 0.d0
DO J= 1,NPTOT
IF(JFXP(J).LE.0) THEN
K= K+1
TSTPSB= DABS(PU(J)/PV(J))
IF(TSTPSB.GT.TSTPS) THEN
JFIX= J
TSTPS= TSTPSB
ENDIF
ENDIF
ENDDO
ENDIF
YC= PV(JFIX)
CALL ROUND(JROUND,NPMAX,NPTOT,NPTOT,JFIX,PV,PU,PS,CM)
JFXP(JFIX)= 1
IF(LPRINT.GE.2)
1 WRITE(6,614) JFIX,YC,PU(JFIX),PS(JFIX),JFIX,PV(JFIX),RMSR
NPARM= NPARM-1
IF(NPARM.EQ.0) THEN
c** After rounding complete, make one more pass with all non-fixed
c parameters set free to get full correct final correlation matrix,
c uncertainties & sensitivities. Don't update parameters on this pass!
NPARM= NPFIT
QUIT= 1
DO J= 1,NPTOT
JFXP(J)= IFXP(J)
ENDDO
c ... reinitialize for derivative-by-differences calculation
RMSR= 0.d0
ENDIF
GO TO 6
c
c** If no parameters varied or sequential rounding completed - simply
c calculate DSE from RMS residuals and return.
60 DSE= 0.d0
IF(NDATA.GT.NPFIT) THEN
DSE= RMSR*DSQRT(DBLE(NDATA)/DBLE(NDATA-NPFIT))
ELSE
DSE= 0.d0
ENDIF
IF(NPFIT.GT.0) THEN
IF(LPRINT.GT.0) THEN
c** Print final rounded parameters with original Uncert. & Sensitivities
IF(QUIT.LT.1) WRITE(6,616) NDATA, NPFIT, RMSR, TSTPS
IF(QUIT.EQ.1) WRITE(6,616) NDATA, NPFIT, RMSR
DO J= 1, NPTOT
IF(JFXP(J).GT.0) THEN
c** If parameter held fixed (by rounding process), shift values of
c change, sensitivity & uncertainty to correct absolute number label.
DO I= NPTOT,J+1,-1
PC(I)= PC(I-1)
PS(I)= PS(I-1)
PU(I)= PU(I-1)
ENDDO
PC(J)= 0.d0
PS(J)= 0.d0
PU(J)= 0.d0
ENDIF
ENDDO
WRITE(6,612) (J,PV(J),PU(J),PS(J),PC(J),J=1,NPTOT)
ENDIF
ENDIF
RETURN
c
602 FORMAT(/' *** NLLSSRR problem: [NPTOT=',i4,'] > min{NPINTMX=',
1 i4,' NPMAX=',i4,', NDATA=',i6,'}')
604 FORMAT(' After Cycle #',i2,': DRMSD=',1PD14.7,' test(PS)=',
1 1PD8.1,' test(PU)=',D8.1)
!!606 FORMAT(/' Effective',i3,'-cycle Cgce: MAX{|change/unc.|}=',
!! 1 1PD8.1,' < 0.01 DRMSD=',D10.3)
607 FORMAT(/' Full',i3,'-cycle convergence: {ABS(RMSR/RMSRB)-1}=',
1 1PD9.2,' TSTPS=',D8.1)
610 FORMAT(/ ' !! CAUTION !! fit of',i4,' parameters to',I6,' data not
1 converged after',i3,' Cycles'/5x,'DRMS(deviations)=',1PD10.3,
2 ' test(PS) =',D9.2,' test(PU) =',D9.2/1x,31('**'))
612 FORMAT((3x,'PV(',i4,') =',1PD22.14,' (+/-',D8.1,') PS=',d8.1,
1 ' PC=',d9.1))
614 FORMAT(' =',39('==')/' Round Off PV(',i4,')=',1PD21.13,' (+/-',
1 D9.2,') PS=',d9.2/4x,'fix PV(',I4,') as ',D19.11,
2 ' & refit: DRMS(deviations)=',D12.5)
616 FORMAT(/i6,' data fit to',i5,' param. yields DRMS(devn)=',
1 1PD14.7:' tst(PS)=',D8.1)
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE QROD(N,NR,NC,A,R,F,B,GC,GS)
C** Performs ORTHOGONAL DECOMPOSITION OF THE LINEAR LEAST-SQUARES
C EQUATION J * X = F TO A * X = B(TRANSPOSE) * F WHERE
C J IS THE JACOBIAN IN WHICH THE FIRST N ROWS AND COLUMNS
C ARE TRANSFORMED TO THE UPPER TRIANGULAR MATRIX A
C (J = B * A), X IS THE INDEPENDENT VARIABLE VECTOR, AND
C F IS THE DEPENDENT VARIABLE VECTOR. THE TRANSFORMATION
C IS APPLIED TO ONE ROW OF THE JACOBIAN MATRIX AT A TIME.
C PARAMETERS :
C N - (INTEGER) DIMENSION OF A TO BE TRANSFORMED.
C NR - (INTEGER) ROW DIMENSION OF A DECLARED IN CALLING PROGRAM.
C NC - (INTEGER) Column DIMENSION OF F DECLARED IN CALLING PROGRAM.
C A - (REAL*8 ARRAY OF DIMENSIONS .GE. N*N) UPPER TRIANGULAR
C TRANSFORMATION MATRIX.
C R - (REAL*8 LINEAR ARRAY OF DIMENSION .GE. N) ROW OF
C JACOBIAN TO BE ADDED.
C F - (REAL*8 LINEAR ARRAY .GE. TO THE ROW DIMENSION OF THE
C JACOBIAN) TRANSFORMED DEPENDENT VARIABLE MATRIX.
C B - (REAL*8) VALUE OF F THAT CORRESPONDS TO THE ADDED
C JACOBIAN ROW.
C GC - (REAL*8 LINEAR ARRAY .GE. N) GIVENS COSINE TRANSFORMATIONS.
C GS - (REAL*8 LINEAR ARRAY .GE. N) GIVENS SINE TRANSFORMATIONS.
C--------------------------------------------------------------------
C AUTHOR : MICHAEL DULICK, Department of Chemistry,
C UNIVERSITY OF WATERLOO, WATERLOO, ONTARIO N2L 3G1
C--------------------------------------------------------------------
INTEGER I,J,K,N,NC,NR
REAL*8 A(NR,NC), R(N), F(NR), GC(N), GS(N), B, Z(2)
DO 10 I = 1,N
Z(1) = R(I)
J = I - 1
DO K = 1,J
Z(2) = GC(K) * A(K,I) + GS(K) * Z(1)
Z(1) = GC(K) * Z(1) - GS(K) * A(K,I)
A(K,I) = Z(2)
ENDDO
GC(I) = 1.D0
GS(I) = 0.D0
IF(DABS(Z(1)).LE.0.D0) GOTO 10
IF(DABS(A(I,I)) .LT. DABS(Z(1))) THEN
Z(2) = A(I,I) / Z(1)
GS(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2))
GC(I) = Z(2) * GS(I)
ELSE
Z(2) = Z(1) / A(I,I)
GC(I) = 1.D0 / DSQRT(1.D0 + Z(2) * Z(2))
GS(I) = Z(2) * GC(I)
ENDIF
A(I,I) = GC(I) * A(I,I) + GS(I) * Z(1)
Z(2) = GC(I) * F(I) + GS(I) * B
B = GC(I) * B - GS(I) * F(I)
F(I) = Z(2)
10 CONTINUE
RETURN
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
SUBROUTINE ROUND(IROUND,NPMAX,NPARM,NPTOT,IPAR,PV,PU,PS,CM)
c** Subroutine to round off parameter # IPAR with value PV(IPAR) at the
c |IROUND|'th significant digit of: [its uncertainty PU(IPAR)] .
c** On return, the rounded value replaced the initial value PV(IPAR).
c** Then ... use the correlation matrix CM and the uncertainties PU(I)
c in the other (NPTOT-1) [or (NPARM-1) free] parameters to calculate
c the optimum compensating changes PV(I) in their values.
INTEGER IROUND,NPMAX,NPARM,NPTOT,IPAR,I,IRND,KRND
REAL*8 PU(NPMAX),PS(NPMAX),PV(NPMAX),CM(NPMAX,NPMAX),CNST,
1 CRND,XRND,FCT,Z0
DATA Z0/0.d0/
CNST= PV(IPAR)
XRND= DLOG10(PU(IPAR))
c** If appropriate, base last rounding step on sensitivity (not uncert.)
IF((NPARM.EQ.1).AND.(PS(IPAR).LT.PU(IPAR))) XRND= DLOG10(PS(IPAR))
c** First ... fiddle with log's to perform the rounding
IRND= INT(XRND)
IF(XRND.GT.0) IRND=IRND+1
IRND= IRND- IROUND
FCT= 10.D0**IRND
CRND= PV(IPAR)/FCT
XRND= Z0
c ... if rounding goes past REAL*8 precision, retain unrounded constant
IF(DABS(CRND).GE.1.D+16) THEN
WRITE(6,601) IROUND,IPAR
RETURN
ENDIF
IF(DABS(CRND).GE.1.D+8) THEN
c ... to avoid problems from overflow of I*4 integers ...
KRND= NINT(CRND/1.D+8)
XRND= KRND*1.D+8
CRND= CRND-XRND
XRND= XRND*FCT
END IF
IRND= NINT(CRND)
CNST= IRND*FCT+ XRND !! rounded constant !!
c** Zero parameters more aggressively ... if unc. > 2* value
if(dabs(PU(IPAR)/PV(IPAR)).GT.2.d0) then
cnst= 0.d0
endif
c** Now ... combine rounding change in parameter # IPAR, together with
c correlation matrix CM and parameter uncertainties PU to predict
c changes in other parameters to optimally compensate for rounding off
c of parameter-IPAR. Method pointed out by Mary Thompson (Dept. of
c Statistics, UW),
IF(IPAR.GT.1) THEN
XRND= (CNST-PV(IPAR))/PU(IPAR)
DO I= 1,NPTOT
IF(I.NE.IPAR) THEN
PV(I)= PV(I)+ CM(IPAR,I)*PU(I)*XRND
ENDIF
ENDDO
ENDIF
PV(IPAR)= CNST
RETURN
601 FORMAT(' =',39('==')/' Caution:',i3,'-digit rounding of parameter-
1',i2,' would exceed (assumed) REAL*8'/' ******** precision overf
2low at 1.D+16, so keep unrounded constant')
END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
c***********************************************************************
c SUBROUTINE DYIDPJ(I,NPTOT,YC,PV,PD)
c** Illustrative dummy version of DYIDPJ for the case of a fit to a
c power series of order (NPTOT-1) in X(i). *** For datum number-i,
c calculate and return PD(j)=[partial derivatives of datum-i] w.r.t.
c each of the free polynomial coefficients varied in the fit
c (for j=1 to NPTOT).
c=====================================================================
c** Use COMMON block(s) to bring in values of the independent variable
c [here XX(i)] and any other parameters or variables needeed to
c calculate YC and the partial derivatives.
c=====================================================================
c INTEGER I,J,NDATA,NPTOT,MXDATA,IFXP(NPTOT)
c PARAMETER (MXDATA= 501)
c REAL*8 RMSR,YC,PV(NPTOT),PD(NPTOT),PS(NPTOT),POWER,XX(MXDATA)
c COMMON /DATABLK/XX
c=====================================================================
c** NOTE BENE(!!) for non-linear fits, need to be sure that the
c calculations of YC and PD(j) are based on the current UPDATED PV(j)
c values. If other (than PV) parameter labels are used internally
c in the calculations, UPDATE them whenever (say) I = 1 .
c=====================================================================
c POWER= 1.D0
c YC= PV(1)
c PD(1)= POWER
c DO 10 J= 2,NPTOT
c POWER= POWER*XX(I)
c YC= YC+ PV(J)*POWER
c PD(J)= POWER
c 10 CONTINUE
c RETURN
c END
c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12