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)]. 23/03/16
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c COPYRIGHT 1998-2016 by Robert J. Le Roy +
c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada +
c This software may not be sold or any other commercial use made +
c of it without the express written permission of the author. +
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++ Version of 14 August 2017 {after vble cleanup}
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/DFLOAT(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/DFLOAT(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(DFLOAT(NDATA)/DFLOAT(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.
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c COPYRIGHT 1998 by Robert J. Le Roy +
c Dept. of Chemistry, Univ. of Waterloo, Waterloo, Ontario, Canada +
c This software may not be sold or any other commercial use made +
c of it without the express written permission of the author. +
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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