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