//U11701AA JOB (*),TIME=(,5),CLASS=1,MSGCLASS=Y, // NOTIFY=* /*JOBPARM ROOM=H,FORMS=9021 // EXEC WATFIV,REGION.WATFIV=500K $JOB ,XREF,LIBLIST C C MARQ.CNTL MARQ SOURCE CODE AND SIMPLE TEST C C C STMQDRIV 1.0 DECEMBER 1991 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C C THIS PROGRAM FITS THE MODEL C C FIT(J) = X(1)/(1.0 + X(2)*T(J,1) + X(3)*T(J,2)) C C TO DATA POINTS (T(J,1),T(J,2),Y(J)) . C T(J,1) AND T(J,2) ARE THE TWO INDEPENDENT VARIABLES IN DATA C SPACE. Y(J) IS THE MEASURED DEPENDENT VARIABLE AND FIT(J) C IS THE VALUE TO BE FITTED TO Y(J) BY MINIMIZING THE C WEIGHTED SUM OF SQUARES, C C J=NPTS C PHI = SUM ((FIT(J)-Y(J))/YSIG(J))**2 C J=1 C C THE VALUE OF PHI IS STORED IN THE VARIABLE FOBJ. C IMPLICIT REAL*8 (A-H,O-Z) C EXTERNAL FUNK C DIMENSION ERGSAV(20),FIDINT(2) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CSTOR/ T(300,2) C C OPEN (UNIT=5,FILE="MARQTEST.DAT") C KR=5 KW=6 C C READ IN THE VALUE OF JROUTE... C JROUTE=1 TO USE MARQ, C JROUTE=2 TO USE STEPIT, C JROUTE=3 TO USE SIMPLEX C READ (KR,10)JROUTE 10 FORMAT(I1) C C READ IN THE VALUE OF NPTS. THEN READ IN THE DATA POINTS. C READ(KR,20)NPTS 20 FORMAT(I5) WRITE(KW,30)NPTS 30 FORMAT(///' NPTS =',I5) C DO 60 K=1,NPTS READ(KR,40)T(K,1),T(K,2),Y(K),YSIG(K) 40 FORMAT(4F10.5) WRITE(KW,50)K,T(K,1),T(K,2),Y(K),YSIG(K) 50 FORMAT(1X,I5,4F10.3) 60 CONTINUE C C INITIALIZE FOR THE FIT. C CALL STSET NV=3 X(1)=10.0 X(2)=1.0 X(3)=1.0 C C FIT THE MODEL TO THE DATA. C CALL STINT (FUNK) C C USE SUBROUTINE FIDO TO CHECK THE ERRORS FROM MARQ C OR STEPIT, OR TO COMPUTE ERRORS USING SIMPLEX. C ERFRAC=0.1D0 C DO 70 JXFID=1,NV C IF(JROUTE.EQ.1 .OR. JROUTE.EQ.2) THEN ERGSAV(JXFID)=DSQRT(DABS(ERR(JXFID,JXFID))) ELSE ERGSAV(JXFID)=ERFRAC*DABS(X(JXFID)) ENDIF C IF(ERGSAV(JXFID).EQ.0.0D0) ERGSAV(JXFID)=ERFRAC 70 CONTINUE C STDEVS=1.0D0 TOLFID=0.05D0 MAXIND=2 NTRACF=1 C IF(JROUTE.EQ.1) THEN NTRSET=-2 ELSE NTRSET=-1 ENDIF C DO 80 JXFID=1,NV ERGUES=ERGSAV(JXFID) CALL FIDO (FUNK,JXFID,STDEVS,TOLFID,ERGUES,MAXIND, * NTRSET,NTRACF,JROUTE,FIT,FIDINT) 80 CONTINUE C STOP C C END STMQDRIV C END SUBROUTINE FUNK (FITM) C C J. P. CHANDLER, COMPUTER SCIENCE DEPT., C OKLAHOMA STATE UNIVERSITY C C THIS VERSION OF SUBROUTINE FUNK COMPUTES THE FITTED VALUES C AND THE WEIGHTED SUM OF SQUARES FOR THE PROBLEM OF FITTING C THE MODEL C C FIT(J) = X(1)/(1.0 + X(2)*T(J,1) + X(3)*T(J,2)) C C TO DATA POINTS (T(J,1),T(J,2),Y(J)) . C C TO RUN THIS ROUTINE WITH MARQ, USE THE STATEMENT C SUBROUTINE FUNK (FITM) . C C TO RUN THIS ROUTINE WITH STEPIT OR SIMPLEX, USE C SUBROUTINE FUNK . C IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION FITM(300) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CSTOR/ T(300,2) C C FOBJ WILL CONTAIN THE WEIGHTED SUM OF SQUARES, PHI. C MARQ REQUIRES THAT FUNK COMPUTE FITM. C WHETHER OR NOT FUNK COMPUTES FOBJ IS IRRELEVANT TO MARQ. C STEPIT, SIMPLEX, ETC. REQUIRE THAT FUNK COMPUTE FOBJ. C THOSE ROUTINES NEVER SEE THE FITTED VALUES AT ALL. C TO ALLOW MARQ TO BE INTERCHANGED WITH STEPIT, SIMPLEX, ETC. C WITH A MINIMUM NUMBER OF CHANGES, WE WILL COMPUTE FOBJ C AS WELL AS FITM. C FOBJ=0.0D0 C DO 10 J=1,NPTS FITM(J)=X(1)/(1.0D0+X(2)*T(J,1)+X(3)*T(J,2)) FOBJ=FOBJ+((FITM(J)-Y(J))/YSIG(J))**2 10 CONTINUE C RETURN C C END FUNK C END SUBROUTINE FIDO (FUNK,JXFID,STDEVS,TOLFID,ERGUES,MAXIND, * NTRSET,NTRACF,JROUTE,FITM, FIDINT) C C FIDO 4.6 FEBRUARY 1993 C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA UNIVERSITY, STILLWATER, OKLAHOMA 74078 C C SUBROUTINE FIDO IS USED IN CONJUNCTION WITH MARQ, STEPIT, OR C SIMPLEX TO COMPUTE CONFIDENCE HALF-INTERVALS ("ERRORS") FOR C THE PARAMETERS IN A FITTING PROBLEM, USING THE METHOD OF C SUPPORT PLANES. C THE QUANTITY BEING MINIMIZED (FOBJ) MUST EITHER BE C CHI-SQUARE OR IT MUST BE TWICE THE NEGATIVE OF THE NATURAL C LOGARITHM OF THE LIKELIHOOD FUNCTION. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C USAGE...... C FIRST CALL STEPIT TO FIND THE OPTIMUM VALUES OF THE C PARAMETERS. THEN CALL FIDO ONCE FOR EACH PARAMETER FOR C WHICH THE CONFIDENCE INTERVALS ARE DESIRED. C C INPUT QUANTITIES..... FUNK,JXFID,STDEVS,TOLFID,ERGUES, C MAXIND,NTRSET,NTRACF,JROUTE,X(*), C KW C C OUTPUT QUANTITY...... FIDINT( ) C C SCRATCH ARRAY (FOR USE WITH MARQ ONLY)... FITM(*) C C FUNK -- THE NAME OF THE SUBROUTINE WHICH COMPUTES C FOBJ GIVEN THE VALUES OF THE X(J) C C JXFID -- THE INDEX OF THE PARAMETER, X(JXFID), FOR C WHICH CONFIDENCE HALF-INTERVALS ARE TO C BE COMPUTED C C STDEVS -- THE NUMBER OF STANDARD DEVIATIONS TO WHICH C THE HALF-INTERVALS ARE TO CORRESPOND C C TOLFID -- CONVERGENCE TOLERANCE (USE TOLFID=0.05) C C ERGUES -- ROUGH ESTIMATE OF THE LENGTH OF THE C HALF-INTERVALS C C MAXIND -- =1 TO COMPUTE THE HALF-INTERVAL WITH THE C SAME SIGN AS ERGUES, C =2 TO COMPUTE BOTH HALF-INTERVALS C C NTRSET -- VALUE OF NTRACE TO BE USED IN THE C FITTING SUBROUTINE (STEPIT OR MARQ, C ETC.) WHEN CALLED BY FIDO C C NTRACF -- A SWITCH THAT CONTROLS PRINTING... C C =-1 TO PRINT NOTHING IN FIDO EXCEPT ANY C WARNING OR ERROR MESSAGES, C C = 0 FOR INITIAL AND SUMMARY PRINTOUT, C C =+1 TO PRINT EVERY FIDO ITERATION C C JROUTE -- =1 TO USE FIDO WITH SUBROUTINE MARQ, C =2 TO USE FIDO WITH STEPIT OR SIMPLEX, ETC. C C FITM(*) -- A SCRATCH ARRAY, USED WHEN FIDO IS USED C WITH MARQ C C X( ) -- THE POSITION OF THE MINIMUM OF FOBJ C C KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER C C FIDINT(J) -- RETURN THE COMPUTED LENGTHS OF THE C CONFIDENCE HALF-INTERVALS C C IF STEPIT COMPUTES THE ERROR MATRIX ERR, ERGUES CAN BE SET C TO PLUS OR MINUS STDEVS*SQRT(ERR(JX,JX)) (SEE BELOW). C C NOTE: C FIDO CALLS STEPIT, WHICH MAY USE THE ARRAY ERR FOR SCRATCH C STORAGE, DESTROYING ITS CONTENTS. C C FOR A LEAST SQUARES PROBLEM, FOBJ IS EQUAL TO CHI-SQUARE, C THE WEIGHTED SUM OF SQUARES... C C NPTS C FOBJ = SUM ((FIT(JPT)-Y(JPT))/YSIG(JPT))**2 C JPT=1 C C THE STANDARD ERRORS YSIG( ) MUST BE CORRECTLY SCALED. C IF THE YSIG( ) ARE NOT KNOWN ACCURATELY, NORMALIZE THEM C SO THAT THE VALUE OF FOBJ AT THE MINIMUM IS EQUAL TO THE C NUMBER OF DEGREES OF FREEDOM, C N.D.F. = (NO. DATA POINTS) - (NO. ADJUSTABLE PARAMETERS) C C FIDO IS ESSENTIALLY A ROOT FINDING ROUTINE. IT USES INVERSE C QUADRATIC INTERPOLATION OR EXTRAPOLATION, INVERSE LINEAR C INTERPOLATION, AND BISECTION (INTERVAL HALVING), AS C SUCCESSIVELY MORE DIFFICULTIES ARE ENCOUNTERED IN FINDING C THE ROOT. C C PRIMARY REFERENCES.... C C W. T. EADIE ET AL., "STATISTICAL METHODS IN EXPERIMENTAL C EXPERIMENTAL PHYSICS" (AMERICAN ELSEVIER, 1971), C CHAPTER 9 C C P. R. BEVINGTON, "DATA REDUCTION AND ERROR ANALYSIS IN C THE PHYSICAL SCIENCES" (MCGRAW-HILL, 1969), C PAGES 242-245 C C OTHER REFERENCES.... C C H. SCHEFFE, "THE ANALYSIS OF VARIANCE" (WILEY, 1959) C C H. STONE, DISCUSSION ON PAPER BY E. M. L. BEALE, C J. ROY. STATIS. SOC. B., V. 22 (1960), P. 41, PP. 84-5 C C G. E. P. BOX AND G. A. COUTIE, PROC. INST. ELEC. ENGRS. C V. 103, PART B, SUPPL. NO. 1 (1956). C C G. W. BOOTH, G. E. P. BOX, M. E. MULLER, AND C T. I. PETERSON, "FORECASTING BY GENERALIZED REGRESSION C METHODS; NON-LINEAR ESTIMATION" (PRINCETON/IBM), C FEBRUARY 1959, IBM MANUAL C C D. W. MARQUARDT, "LEAST-SQUARES ESTIMATION OF NONLINEAR C PARAMETERS", SHARE DISTRIBUTION 3094 C C D. W. MARQUARDT, R. G. BENNETT, AND E. J. BURRELL, "LEAST C SQUARES ANALYSIS OF ELECTRON PARAMAGNETIC RESONANCE C SPECTRA", J. OF MOLECULAR SPECTROSCOPY 7 (1961) 269 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME C COMPILERS (WATFIV, FOR EXAMPLE) AND IT IS FORBIDDEN BY SOME C OTHERS (MODCOMP II, FOR EXAMPLE). C EXTERNAL FUNK C INTEGER NV,NTRACE,MATRX,MASK,NFMAX,NFLAT,JVARY,NXTRA, * KFLAG,NOREP,KERFL,KW, * JXFID,MAXIND,NTRACF,JROUTE, * ITER,J,JJ,K,KLOSB,MATSAV,MAXTRY,MLIN,MSKSAV,NACTIV, * NTRSAV,NTRSET C DOUBLE PRECISION STDEVS,TOLFID,ERGUES,FITM,FIDINT,XSAVE, * XKA,XKB DOUBLE PRECISION FOBJ,X,XMAX,XMIN,DELTX,DELMIN,ERR DOUBLE PRECISION FSAVE,T,DABS,DSQRT,ARG,ZSQRT DOUBLE PRECISION RSMALL,CONVRG,BRACK,FACMAX,FRMIN,RZERO, * RHALF,UNITR,FPSGSQ,SGN,FB,FC,DX,DC,DIFF,FKA,FKB,FAC, * XX,FRAC,XBEST,FBEST C DIMENSION FIDINT(2),FITM(1) DIMENSION XSAVE(20),XKA(20),XKB(20),XBEST(20) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C C SET THE LIBRARY FUNCTIONS FOR REAL PRECISION (SQRT, ETC.) C OR FOR DOUBLE PRECISION (DSQRT, ETC.). C ZSQRT(ARG)= SQRT(ARG) ZSQRT(ARG)=DSQRT(ARG) C T(ARG)= ABS( SQRT(ARG-FSAVE)-STDEVS) T(ARG)=DABS(DSQRT(ARG-FSAVE)-STDEVS) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C SET SOME PARAMETERS. C C RSMALL IS USED IN SETTING A PHONY VALUE IF A VALUE OF FOBJ C IS FOUND THAT IS LESS THAN THE "GLOBAL MINIMUM" VALUE. C RSMALL=0.0001D0 C C CONVRG = SATISFACTORY RATE OF CONVERGENCE C CONVRG=0.5D0 C C BRACK = FACTOR FOR BRACKETTING SHOTS C BRACK=2.0D0 C C MAXTRY = MAXIMUM NUMBER OF ITERATIONS C MAXTRY=20 C C FACMAX = MAXIMUM FACTOR FOR QUADRATIC INTERPOLATION C FACMAX=4.0D0 C C FRMIN = MINIMUM FRACTION FOR LINEAR INTERPOLATION C FRMIN=0.1D0 C RZERO=0.0D0 RHALF=0.5D0 UNITR=1.0D0 C C SAVE SOME INPUT QUANTITIES, AND INITIALIZE EVERYTHING. C NTRSAV=NTRACE NTRACE=NTRSET C MATSAV=MATRX MATRX=0 C MSKSAV=MASK(JXFID) MASK(JXFID)=1 C NACTIV=0 C DO 10 J=1,NV XBEST(J)=X(J) XSAVE(J)=X(J) IF(MASK(J).EQ.0) NACTIV=NACTIV+1 10 CONTINUE C IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF C FSAVE=FOBJ FBEST=FOBJ FPSGSQ=FOBJ+STDEVS**2 X(JXFID)=XSAVE(JXFID)+STDEVS*ERGUES C C LOOP OVER THE FIDINT(JJ). C JJ=1 C 20 IF(ERGUES.NE.RZERO) GO TO 40 WRITE(KW,30)ERGUES 30 FORMAT(//' ERROR IN INPUT TO FIDO...', * ' ERGUES IS EQUAL TO ZERO.') STOP C 40 SGN=UNITR IF(ERGUES.LT.RZERO) SGN=-UNITR CONTINUE IF(NTRACF.LT.0) GO TO 70 WRITE(KW,50)JXFID,STDEVS,TOLFID,FSAVE, * XSAVE(JXFID),ERGUES 50 FORMAT(////' SUBROUTINE FIDO.'//' JXFID =',I3,8X, * 'STDEVS =',1PG12.5,4X,' TOLFID = ',G12.5// * ' GLOBAL MINIMUM OF FOBJ =',G15.7, * ' IS AT X(JXFID) = ',G15.7// * ' FIRST GUESS FOR LENGTH OF HALF-INTERVAL = ',E12.5) WRITE(KW,60)NV,NACTIV 60 FORMAT(/10X,' NV =',I11,8X,'NACTIV =',I11//' ') C 70 KLOSB=0 MLIN=0 C DO 80 K=1,NV XKA(K)=XSAVE(K) 80 CONTINUE C FB=FSAVE FKA=FSAVE C C BEGIN THE ITERATION LOOP FOR LOCATING THE PLANE C PERPENDICULAR TO THE X(JXIFD) AXIS IN WHICH THE MINIMUM C VALUE OF FOBJ IS EQUAL TO FPSGSQ = FSAVE + STDEVS**2 , C WHERE FSAVE IS THE VALUE OF FOBJ AT THE GLOBAL MINIMUM. C DO 310 ITER=1,MAXTRY FC=FB C IF(NACTIV.EQ.0) THEN C IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF C ELSE CALL STINT (FUNK) ENDIF C IF(FOBJ.GE.FBEST) GO TO 100 FBEST=FOBJ C DO 90 K=1,NV XBEST(K)=X(K) 90 CONTINUE C 100 DX=X(JXFID)-XSAVE(JXFID) DC=FOBJ-FPSGSQ IF(NTRACF.GE.1) WRITE(KW,110)ITER,JXFID,X(JXFID),DX, * FOBJ,FPSGSQ,DC 110 FORMAT(//' ITERATION ',I2//5X,'X(',I3,') = ',1PG18.10/ * 5X,'DISTANCE FROM GLOBAL MINIMUM = ',G12.5// * 5X,'MINIMUM FOBJ IN THIS PLANE = ',G15.7//5X, * 'VALUE SOUGHT = ',G15.7,6X,'DIFFERENCE = ',G12.5) C IF(FOBJ.EQ.FSAVE) GO TO 140 IF(FOBJ.LT.FSAVE) GO TO 120 C C TEST FOR CONVERGENCE. C C IF(ABS(FOBJ-FPSGSQ).LE.TOL) ............................ C DIFF=FOBJ-FPSGSQ IF(DIFF.LT.RZERO) DIFF=-DIFF IF(DIFF.LE.TOLFID) GO TO 330 IF(FOBJ.GE.FPSGSQ) GO TO 170 GO TO 150 C 120 DIFF=FSAVE-FOBJ WRITE(KW,130)DIFF,(X(K),K=1,NV) 130 FORMAT(//' ***** FOBJ LESS THAN VALUE AT INPUT POINT', * ' BY',1PG13.5//9X,' X(J) AT THIS POINT....'// * (1X,3G23.15)) 140 FOBJ=FSAVE+RSMALL*(FPSGSQ-FSAVE) C C CHECK FOR TIGHTER BRACKETTING OF FPSGSQ FROM BELOW. C 150 IF((X(JXFID)-XKA(JXFID))*SGN.LE.RZERO) GO TO 190 C DO 160 K=1,NV XKA(K)=X(K) 160 CONTINUE C FKA=FOBJ GO TO 190 C C CHECK FOR BRACKETTING OF FPSGSQ FROM ABOVE. C 170 IF(KLOSB.EQ.1) THEN IF((X(JXFID)-XKB(JXFID))*SGN.GE.RZERO) GO TO 190 ENDIF KLOSB=1 C DO 180 K=1,NV XKB(K)=X(K) 180 CONTINUE C FKB=FOBJ C 190 IF(T(FOBJ).LT.T(FC)) FB=FOBJ C C CHECK THE RATE OF CONVERGENCE. IF IT IS SATISFACTORY, AND C IF LINEAR INTERPOLATION HAS BEEN USED, USE IT AGAIN. C IF(ITER.GE.2 .AND. T(FOBJ).GT.CONVRG*T(FC)) GO TO 210 IF(MLIN.EQ.1) GO TO 250 C C USE INVERSE QUADRATIC INTERPOLATION. C FAC=STDEVS/ZSQRT(FOBJ-FSAVE) C C FAC=AMIN1(FACMAX,AMAX1(FAC,1./FACMAX)) ................... C IF(FAC.LT.UNITR/FACMAX) FAC=UNITR/FACMAX IF(FAC.GT.FACMAX) FAC=FACMAX XX=XSAVE(JXFID)+FAC*(X(JXFID)-XSAVE(JXFID)) C C CHECK THAT THE PROPOSED POINT IS INSIDE THE BRACKETTED C INTERVAL. C IF((XX-XKA(JXFID))*SGN.LE.RZERO) GO TO 210 IF(KLOSB.EQ.1) THEN IF((XX-XKB(JXFID))*SGN.GE.RZERO) GO TO 240 ENDIF C DO 200 K=1,NV IF(K.EQ.JXFID .OR. MASK(K).EQ.0) * X(K)=XSAVE(K)+FAC*(X(K)-XSAVE(K)) 200 CONTINUE C GO TO 310 C 210 IF(KLOSB.EQ.1) GO TO 240 C C CONVERGENCE IS POOR, AND FPSGSQ HAS NOT YET BEEN BRACKETTED. C TRY TO BRACKET IT. C FRAC=BRACK IF(FOBJ.GE.FPSGSQ) FRAC=UNITR/BRACK C DO 220 K=1,NV IF(K.EQ.JXFID .OR. MASK(K).EQ.0) * X(K)=XSAVE(K)+FRAC*(X(K)-XSAVE(K)) 220 CONTINUE C IF(NTRACF.GE.1) WRITE(KW,230) 230 FORMAT(//' BRACKETTING SHOT....') GO TO 310 C 240 IF(MLIN.EQ.1) GO TO 270 C C TRY LINEAR INTERPOLATION BETWEEN THE TWO BRACKETTING POINTS. C 250 MLIN=1 FRAC=(FPSGSQ-FKA)/(FKB-FKA) C C FRAC=AMAX1(FRMIN,AMIN1(1.0-FRMIN,FRAC)) ............. C IF(FRAC.GT.UNITR-FRMIN) FRAC=UNITR-FRMIN IF(FRAC.LT.FRMIN) FRAC=FRMIN C CONTINUE IF(NTRACF.GE.1) WRITE(KW,260) 260 FORMAT(//' LINEAR INTERPOLATION....') GO TO 290 C C CONVERGENCE IS POOR, AND LINEAR INTERPOLATION HAS BEEN USED. C BISECT THE BRACKETTED INTERVAL. C 270 FRAC=RHALF IF(NTRACF.GE.1) WRITE(KW,280) 280 FORMAT(//' INTERVAL BISECTED....') C 290 DO 300 K=1,NV IF(K.EQ.JXFID .OR. MASK(K).EQ.0) * X(K)=XKA(K)+FRAC*(XKB(K)-XKA(K)) 300 CONTINUE C 310 CONTINUE C C END OF ITERATION LOOP. THE ITERATION FAILED TO CONVERGE. C WRITE(KW,320) 320 FORMAT(' CONVERGENCE FAILURE IN SUBROUTINE FIDO.') FIDINT(JJ)=RZERO GO TO 350 C C CONVERGENCE ACHIEVED ... A SATISFACTORY PLANE HAS BEEN C LOCATED. C 330 FIDINT(JJ)=X(JXFID)-XSAVE(JXFID) IF(NTRACF.GE.0) WRITE(KW,340)JXFID,FIDINT(JJ),STDEVS 340 FORMAT(///' THE LENGTH OF THE CONFIDENCE HALF-INTERVAL', * ' FOR X(',I3,') IS ',1PG13.5/9X,'(STDEVS = ',G12.5, * ')'//' ') C 350 JJ=JJ+1 IF(JJ.GT.MAXIND) GO TO 370 C C REFLECT X THROUGH XSAVE, AND SEARCH FOR THE OTHER C HALF-INTERVAL. C DO 360 K=1,NV IF(K.EQ.JXFID .OR. MASK(K).EQ.0) * X(K)=XSAVE(K)+(XSAVE(K)-X(K)) 360 CONTINUE C ERGUES=X(JXFID)-XSAVE(JXFID) GO TO 20 C C END OF THE JJ LOOP. PRINT SUMMARY RESULTS. C 370 IF(MAXIND.LT.2 .OR. FIDINT(1).EQ.RZERO) GO TO 400 IF(FIDINT(1).LT.RZERO .AND. FIDINT(2).LE.RZERO) GO TO 400 IF(FIDINT(1).GT.RZERO .AND. FIDINT(2).GE.RZERO) GO TO 400 IF(FIDINT(1).GT.RZERO .AND. FIDINT(2).LT.RZERO) GO TO 380 XX=FIDINT(1) FIDINT(1)=FIDINT(2) FIDINT(2)=XX 380 XX=-FIDINT(2) IF(NTRACF.GE.0) WRITE(KW,390)FIDINT(1),JXFID, * XSAVE(JXFID),STDEVS,XX 390 FORMAT(///' SUMMARY OF RESULTS FROM FIDO....'///24X, * '+',1PG12.5/' X(',I3,') = ',G15.8,25X, * ' STDEVS = ',G12.5/24X,'-',G12.5///' ') C C RESTORE THE SAVED ENTRY VALUES. C 400 DO 410 J=1,NV X(J)=XBEST(J) 410 CONTINUE C IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF C IF(FOBJ.GE.FSAVE) GO TO 430 WRITE(KW,420)FSAVE,FOBJ,(X(J),J=1,NV) 420 FORMAT(///' SUBROUTINE FIDO FOUND A BETTER MINIMUM.'// * ' OLD FOBJ = ',1PG18.10,6X,' NEW FOBJ = ',G18.10,8X, * 'X(J)....'/(1X,4G18.10)) 430 MASK(JXFID)=MSKSAV NTRACE=NTRSAV MATRX=MATSAV C RETURN C C END FIDO C END SUBROUTINE STEPT (FUNK) C C INTERFACE TO MAKE MARQ LOOK LIKE STEPT. C C J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, C OKLAHOMA STATE UNIVERSITY C C TO USE THIS ROUTINE, SET THE VALUES OF LPDMA AND LPDMB, AND SET C THE DIMENSIONS OF THE ARRAYS P, FITSV, FIT, Y, AND YSIG. C NORMALLY THE DIMENSIONS ARE.... C P(LPDMA,LPDMB), FITSV(LPDMA), FIT(LPDMA), Y(LPDMA), C YSIG(LPDMA), WHERE LPDMA IS .GE. NPTS AND LPDMB IS .GE. NV. C C COMMON/CDAT/ DOES NOT APPEAR IN ANY ROUTINE OF THE MARQ PACKAGE C OTHER THAN THIS ONE, SO THIS IS THE ONLY ROUTINE THAT MUST BE C RECOMPILED WHEN THE DIMENSIONS OF THE ARRAYS ARE CHANGED. C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME COMPILERS C (WATFIV, FOR EXAMPLE) AND FORBIDDEN BY OTHERS (MODCOMP II). EXTERNAL FUNK C INTEGER LPDMA,LPDMB,NPTS C DOUBLE PRECISION FIT,FITSV ,Y,YSIG,P C DIMENSION P(300,20),FITSV(300) C COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS C LPDMA=300 LPDMB=20 C CALL MARQ (FUNK,Y,YSIG,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C RETURN C C END STEPT C END SUBROUTINE STINT (FUNK) C C INTERFACE ROUTINE TO FACILITATE SWAPPING OF STEPIT, C SIMPLEX, MARQ, ETC. C THIS VERSION OF STINT CALLS MARQ. C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C TO USE THIS ROUTINE, SET THE VALUES OF LPDMA AND LPDMB, AND C SET THE DIMENSIONS OF THE ARRAYS P, FITSV, FIT, Y, AND YSIG. C NORMALLY THE DIMENSIONS ARE.... C P(LPDMA,LPDMB), FITSV(LPDMA), FIT(LPDMA), Y(LPDMA), C YSIG(LPDMA), WHERE LPDMA IS .GE. NPTS AND LPDMB IS .GE. NV. C C COMMON/CDAT/ DOES NOT APPEAR IN ANY ROUTINE OF THE MARQ C PACKAGE OTHER THAN THIS ROUTINE, SO THIS IS THE ONLY ROUTINE C THAT MUST BE RECOMPILED WHEN THE DIMENSIONS OF THE DATA C ARRAYS ARE CHANGED. C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME C COMPILERS (WATFIV, FOR EXAMPLE) AND IT IS FORBIDDEN BY SOME C OTHER COMPILERS (FOR THE MODCOMP II, FOR EXAMPLE). EXTERNAL FUNK C INTEGER LPDMA,LPDMB,NPTS C DOUBLE PRECISION FIT,FITSV ,Y,YSIG,P C DIMENSION P(300,20),FITSV(300) C COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS C LPDMA=300 LPDMB=20 C CALL MARQ (FUNK,Y,YSIG,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C RETURN C C END STINT C END SUBROUTINE MARQ (FUNK,Y,YSIG,NPTS,FIT,FITSV,P, * LPDMA,LPDMB) C C MARQ 3.3 FEBRUARY 1993 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1981, 1991 J. P. CHANDLER C C J. P. CHANDLER AND LEON W. JACKSON, C COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY, C STILLWATER, OKLAHOMA 74078 C C MARQ PERFORMS A NONLINEAR LEAST SQUARES FIT OF A C USER-SUPPLIED FUNCTION TO A GIVEN SET OF DATA, USING C MARQUARDT'S METHOD, THE GAUSS-NEWTON METHOD, OR THE C MODIFIED GAUSS-NEWTON METHOD. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,X(*),XMAX(*),XMIN(*),DELMIN(*), C NV,NTRACE,MASK(*),MATRX,NFMAX, C NFLAT,KW, C Y(*),YSIG(*),NPTS,LPCOL, C FLAMBD,FNU,RELDIF,METHD,KALCP, C KORDIF,MAXIT,LEQU,MAXSUB C C OUTPUT QUANTITIES.... X(*),FOBJ,ERR(*,*),KFLAG,FIT(*) C C SCRATCH VECTOR....... DELTX(*) C C UNUSED QUANTITIES (INCLUDED FOR COMPATIBILITY WITH C STEPIT COMMON).... JVARY,NXTRA,NOREP,KERFL C C C FUNK -- THE NAME OF THE SUBROUTINE CALLED TO C OBTAIN THE FITTED VALUES IF KALCP=0 C C X(JX) -- THE JX-TH PARAMETER C C XMAX(JX) -- THE UPPER LIMIT FOR X(JX) C C XMIN(JX) -- THE LOWER LIMIT FOR X(JX) C C DELMIN(JX) -- THE CONVERGENCE TOLERANCE FOR X(JX) C C FOBJ -- RETURNS THE FINAL VALUE OF PHI, THE C WEIGHTED SUM OF SQUARES (PHI IS BEING C MINIMIZED AS A FUNCTION OF THE X(JX)) C C NV -- THE NUMBER OF PARAMETERS C C NTRACE -- USER PRINT CONTROL ... C C =-4 FOR NO OUTPUT EXCEPT FATAL ERROR C MESSAGE, C C =-3 FOR NO OUTPUT EXCEPT ERROR C MESSAGES, C C =-2 FOR NO OUTPUT EXCEPT DIAGNOSTIC C MESSAGES, C C =-1 FOR STANDARD OUTPUT EXCEPT NO FINAL C FIT VALUES, C C = 0 FOR STANDARD OUTPUT, C C = 1 TO PRINT ALSO THE RESULTS OF EACH C ITERATION, C C = 2 TO PRINT ALSO THE COEFFICIENT C MATRIX, QSAV(*,*), C C = 3 TO PRINT ALSO THE JACOBIAN C MATRIX, P(*,*) C C MATRX -- NONZERO IF STATISTICAL ERRORS ARE TO BE C COMPUTED C C MASK(JX) -- NONZERO IF X(JX) IS TO BE HELD FIXED C C NFMAX -- THE MAXIMUM NUMBER OF FUNCTION C COMPUTATIONS C C NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN C TWO ITERATIONS GIVE IDENTICAL VALUES C OF PHI. C THE RECOMMENDED VALUE OF NFLAT IS C USUALLY NFLAT=1 . C C KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER C C ERR(JX,KX) -- RETURNS THE ERROR MATRIX C C KFLAG -- RETURNED .GT. ZERO FOR A NORMAL EXIT, C RETURNED .LT. ZERO FOR AN ABNORMAL EXIT C C Y(JPT) -- THE JPT-TH DATA ORDINATE C C YSIG(JPT) -- THE EXPECTED ERROR IN Y(JPT) C C NPTS -- THE NUMBER OF DATA OBSERVATIONS C C FIT(JPT) -- THE JPT-TH FITTED VALUE C C FITSV -- A SCRATCH VECTOR OF NPTS VALUES C C P(JPT,JX) -- THE FIRST PARTIAL DERIVATIVE C OF FIT(JPT) WITH RESPECT TO X(JX) C C LPDMA -- THE FIRST DIMENSION OF THE ARRAY C CONTAINING P C (LPCOL MUST BE .GE. NPTS IF KALCP IS C .GE. ZERO) C C LPDMB -- THE SECOND DIMENSION OF THE ARRAY C CONTAINING P (LPDMB MUST BE .GE. NV) C C FLAMBD -- MARQUARDT'S LAMBDA, THE RELATIVE AMOUNT BY C WHICH THE DIAGONAL COEFFICIENTS OF THE C NORMAL EQUATIONS ARE AUGMENTED. C MARQUARDT RECOMMENDED FLAMBD=0.01 BUT C SOME USERS PREFER FLAMBD=1.0, SLOWER C BUT SAFER. C C FNU -- MARQUARDT'S NU, THE FACTOR BY WHICH FLAMBD C IS CHANGED. C MARQUARDT RECOMMENDED FNU=10.0 . C C RELDIF -- DETERMINES THE MAGNITUDE OF THE C DIFFERENCING STEP. C IF KORDIF.EQ.1 , RELDIF SHOULD USUALLY C BE OF THE ORDER OF C SQRT(MACHINE EPSILON). C IF KORDIF.EQ.2 , RELDIF SHOULD USUALLY C BE OF THE ORDER OF C (MACHINE EPSILON)**(1.0/3.0) . C FOR DOUBLE PRECISION, YOU CAN USE C MACHINE EPSILON = 1.0D-16 . C C RELMIN -- RELATIVE CONVERGENCE CRITERION USED FOR C EACH X(J) FOR WHICH DELMIN(J) IS ZERO C C METHD -- DETERMINES THE METHOD USED... C C =-1, MODIFIED GAUSS-NEWTON METHOD, C C = 0, GAUSS-NEWTON METHOD, C C = 1, MODIFIED FORM OF MARQUARDT'S C METHOD, C C = 2, MARQUARDT'S METHOD C C THE RECOMMENDED VALUE OF METHD IS C USUALLY METHD=1 . C C KALCP -- DETERMINES WHICH ROUTINE IS CALLED C TO COMPUTE THE JACOBIAN MATRIX C OF FIRST PARTIAL DERIVATIVES... C C =-1, ONE ROW OF P RETURNED BY DERIV C OR CALCD. DERIV CALLS FOFX. C C = 0, ALL OF P RETURNED BY DERIV OR C CALCD. DERIV CALLS FUNK. C C = 1, ALL OF P RETURNED BY DERIV OR C CALCD. DERIV CALLS FOFX. C C THE RECOMMENDED VALUE OF KALCP IS C USUALLY KALCP=0 . C C KORDIF -- DETERMINES HOW P IS TO BE CALCULATED... C C = 1, FIRST ORDER APPROXIMATION USED BY C DERIV, C C = 2, SECOND ORDER APPROXIMATION USED BY C DERIV, C C = 3, ANALYTICAL DERIVATIVE SUPPLIED BY C CALCD C C THE RECOMMENDED VALUE OF KORDIF IS C KORDIF=1 FOR SPEED OR KORDIF=2 FOR C ACCURACY. KORDIF=3 IS ACCURATE AND C FAIRLY FAST, BUT A LOT MORE WORK. C C MAXIT -- THE MAXIMUM NUMBER OF ITERATIONS TO BE C PERFORMED C C LEQU -- SAVES STORAGE IF ALL YSIG(JPT) ARE C THE SAME... C C =1 IF ALL YSIG(JPT) ARE EQUAL C (IN THIS CASE ONLY YSIG(1) IS C REFERENCED), C C =0 OTHERWISE C C MAXSUB -- THE MAXIMUM NUMBER OF SUBITERATIONS C PERMITTED C C MAXUPD -- NOT USED AT PRESENT C C C MODIFICATIONS TO MARQUARDT'S METHOD... C C 1. THE QUANTITY (1+LAMBDA) IS NOT ALLOWED TO INCREASE BY C MORE THAN A FACTOR OF FNULM, IN ANY SINGLE CHANGE. C C 2. WHEN CUTSTEPS ARE USED IN A LINE SEARCH (WHEN THE C ANGLE GAMMA .LT. GAMMA-SUB-ZERO), THE VALUE OF LAMBDA C IS INCREASED. C C IN BOTH THE MODIFIED MARQUARDT'S METHOD AND THE MODIFIED C GAUSS-NEWTON METHOD, A HALF STEP IS ATTEMPTED FOLLOWING EACH C SUCCESSFUL STEP. THIS AVOIDS ONE FORM OF VERY SLOW C CONVERGENCE. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME C COMPILERS (WATFIV, FOR EXAMPLE) AND IT IS FORBIDDEN BY C SOME OTHERS (MODCOMP II, FOR EXAMPLE). C EXTERNAL FUNK C INTEGER ITER,J,JBND,JINV,JPT,JPU,JQ,JSMAX,JSUB, * JT,JVARY,JX,JXLIM,K,KALCP,KBND,KERFL,KFLAG, * KORDIF,KPT,KQ,KRANK,KT,KW,KX,L,LEQU,LPDMA, * LPDMB,MASK,MASKT,MATRX,MAXIT,METHD,MRANK, * MAXSUB,MAXUPD,NACSV,NACT,NACTIV,NF,NFLAT,NFMAX, * NLOOP,NMU,NOREP,NPTS,NSMAL,NTRACE,NV,NVMAX, * NXTRA C DOUBLE PRECISION P,Y,YSIG,H,XMAX,XMIN,DELTX,DELMIN, * ERR,FOBJ,FLAMBD,FNU,RELDIF,ZSQRT,ARG,CRIT,FLDEF, * RELMIN,RELTOL,FACCL,HUGE,RZERO,UNITR,RTWO,DELN,FMGN, * FCUT,STFAC,PTERM,SCALJ,SA,PIVOT,EM,SUM,COSIN,SB,SC,HH, * FRMIN,XMX,XMN,DENOM,FRAC,UPFAC,RLFAC,RSFAC,RMSDV, * SDVMX,ZMAX1,ZMIN1,ZABS,ARGB,DOWNU,FNULM,UPNU,YY, * RLAMBD DOUBLE PRECISION X,XSAVE,XTEMP,GRAD,FIT,FITSV, * XSAV,SIG,RTERM,PHI,PHNEW,PHALF,XLIM, * DABS,DMAX1,DMIN1,DSQRT C C THE DIMENSIONS OF THE VECTORS AND MATRICES ARE.... C P(NPTS,NACTIV) (OR P(1,NACTIV) IF KALCP.EQ.-1), C FITSV(NPTS) (OR FITSV(1) IF KALCP.EQ.-1), C X(NV),XMAX(NV),XMIN(NV),DELMIN(NV),ERR(NV,NV+1), C XSAVE(NV),H(NV),MASKT(NV), C GRAD(NACTIV),SCALE(NACTIV), C Y(NPTS),FIT(NPTS),YSIG(NPTS) (OR YSIG(1) IF LEQU.NE.0), C WHERE NACTIV IS THE NUMBER OF ACTIVE (UNMASKED) X(J). C DIMENSION P(LPDMA,LPDMB) DIMENSION Y(1),YSIG(1),FIT(1),FITSV(1) DIMENSION XSAVE(20),H(20),GRAD(20),MASKT(20),XTEMP(20) C C USER COMMON..... COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C C MARQ COMMON..... COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C C SET THE LIBRARY FUNCTIONS FOR DOUBLE PRECISION C (DSQRT, DABS, DMAX1, DMIN1) C OR FOR SINGLE PRECISION C (SQRT, ABS, AMAX1, AMIN1). C IT IS RECOMMENDED THAT AN ARITHMETIC WITH AT LEAST C TWELVE SIGNIFICANT DECIMAL DIGITS BE USED. C ZSQRT(ARG)=DSQRT(ARG) ZABS(ARG)=DABS(ARG) ZMAX1(ARG,ARGB)=DMAX1(ARG,ARGB) ZMIN1(ARG,ARGB)=DMIN1(ARG,ARGB) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C SET FIXED QUANTITIES .... C C NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV. IT IS ALSO C THE DIMENSION OF THE ARRAYS X, XMAX, XMIN, MASK, DELMIN, C XSAVE, H, MASKT, GRAD, AND SCALE, AND THE FIRST DIMENSION C OF ERR. THE SECOND DIMENSION OF ERR IS NVMAX+1. C NVMAX=20 C C CRIT = COSINE OF MARQUARDT'S CRITICAL ANGLE, GAMMA-SUB-ZERO C CRIT=0.70711D0 C C FLDEF = DEFAULT VALUE FOR LAMBDA C FLDEF=1.0D0 C C RELTOL = TOLERANCE FOR A MESSAGE WARNING OF ILL-CONDITIONING C RELTOL=1.0D-4 C C FACCL = ACCELERATION FACTOR FOR FCUT C FACCL=2.0D0 C C FNULM = LIMIT ON THE FACTOR BY WHICH (1+FLAMBD) MAY CHANGE C FNULM=2.0D0 C C HUGE = A LARGE REAL NUMBER, THE DEFAULT VALUE FOR XMAX C AND (-XMIN) C HUGE=1.0D30 C RZERO=0.0D0 UNITR=1.0D0 RTWO=2.0D0 C C NO FLOATING POINT CONSTANTS ARE USED BEYOND THIS POINT. C KFLAG=0 NOREP=0 KERFL=0 C RLAMBD=FLAMBD ITER=0 PHI=HUGE SC=HUGE C IF(NTRACE.GE.-1) WRITE(KW,10) 10 FORMAT('1MARQ.... BEGIN NONLINEAR LEAST SQUARES', * ' SOLUTION') C IF(NV.GT.0 .AND. NV.LE.NVMAX .AND. NV.LE.LPDMB .AND. * NPTS.GE.1 .AND. (KALCP.LT.0 .OR. NPTS.LE.LPDMA)) * GO TO 30 C WRITE(KW,20)NV,NVMAX,NPTS,LPDMA,LPDMB,KALCP 20 FORMAT(//' ****** ILLEGAL INPUT VALUE IN SUBROUTINE', * ' MARQ.'/6X,'NV =',I14,6X,'NVMAX =',I14,6X,'NPTS =', * I14/6X,'LPDMA =',I14,6X,'LPDMA =',I14,6X,'KALCP =', * I14) STOP C C CHECK SOME INPUT QUANTITIES, AND SET THEM TO DEFAULT VALUES C IF DESIRED. C C NACTIV = NUMBER OF ACTIVE X(J) C 30 NACTIV=0 C DO 60 JX=1,NV DELN=ZABS(DELMIN(JX)) IF(MASK(JX).NE.0) GO TO 50 NACTIV=NACTIV+1 IF(DELN.EQ.RZERO) DELN=ZABS(RELMIN*X(JX)) IF(DELN.EQ.RZERO) DELN=RELMIN IF(XMAX(JX).GT.XMIN(JX)) GO TO 40 XMAX(JX)=HUGE XMIN(JX)=-HUGE C 40 X(JX)=ZMAX1(XMIN(JX),ZMIN1(XMAX(JX),X(JX))) C 50 DELMIN(JX)=DELN 60 CONTINUE C IF(NTRACE.LT.-1) GO TO 130 WRITE(KW,70)(MASK(J),J=1,NV) 70 FORMAT(/' MASK =',I7,4I13/(2X,5I13)) WRITE(KW,80)(X(J),J=1,NV) 80 FORMAT(/' X =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,90)(XMAX(J),J=1,NV) 90 FORMAT(/' XMAX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,100)(XMIN(J),J=1,NV) 100 FORMAT(/' XMIN =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,110)(DELMIN(J),J=1,NV) 110 FORMAT(/' DELMIN=',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,120)NV,NPTS,LPDMA,LPDMB,NTRACE,METHD,KALCP, * KORDIF,NFLAT,NFMAX,MAXIT,MAXSUB,CRIT,RELDIF,RELMIN 120 FORMAT(//' NV =',I11,8X,'NPTS =',I11,8X,'LPDMA =',I11// * ' LPDMB =',I11,8X,'NTRACE =',I4,8X,'METHD =',I3// * ' KALCP =',I3,8X,'KORDIF =',I2,8X,'NFLAT =',I2,8X, * 'NFMAX =',I11//' MAXIT =',I8,8X,'MAXSUB =',I6// * ' CRIT =',1PG12.5,8X,'RELDIF =',G12.5,8X,'RELMIN =', * G12.5) C 130 JVARY=0 C C SET FMGN. IF NECESSARY, RESET RLAMBD AND/OR FACCL. C FMGN=UNITR UPNU=FNU DOWNU=FNU IF(RLAMBD.LE.RZERO) RLAMBD=FLDEF IF(METHD.LE.0) RLAMBD=RZERO IF(METHD.NE.1) FACCL=UNITR C C COMPUTE THE INITIAL GOODNESS OF FIT OF THE MODEL TO THE C DATA. CALL FUNC TO CALCULATE THE VECTOR OF FITTED VALUES. C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) C C NF = EQUIVALENT NUMBER OF CALLS TO FUNC C NF=1 IF(NTRACE.GE.-1) WRITE(KW,140)PHI,RLAMBD 140 FORMAT(//' PHI (THE SUM OF SQUARES) =',1PG15.8,10X, * 'LAMBDA =',G12.5//' ') C IF(NACTIV.GT.0) GO TO 160 KFLAG=-1 IF(NTRACE.GE.-3) WRITE(KW,150) 150 FORMAT(//' ****** WARNING... MASK(J).NE.0 FOR ALL J,', * ' IN SUBROUTINE MARQ.'//6X, * 'PHI WILL BE EVALUATED BUT NOT MINIMIZED.') GO TO 1210 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C BEGIN THE NEXT ITERATION. C C THIS IS THE ENTRY POINT AFTER A SUCCESSFUL STEP, IF THE C CONVERGENCE CRITERION IS NOT MET. C 160 JSUB=0 FCUT=RTWO ITER=ITER+1 IF(NTRACE.GE.1) WRITE(KW,170)ITER,FMGN,RLAMBD 170 FORMAT(//' BEGIN ITERATION',I5,8X,'FMGN =',1PG9.2, * 8X,'LAMBDA =',G9.2) IF(NTRACE.GE.3) WRITE(KW,180) 180 FORMAT(/' P(*,*) (THE JACOBIAN MATRIX) ='/' ') C C INITIALIZE FOR THIS ITERATION. C STFAC=UNITR C DO 200 JX=1,NACTIV GRAD(JX)=RZERO C DO 190 KX=1,JX ERR(JX,KX)=RZERO 190 CONTINUE C 200 CONTINUE C C CALL DERIV (OR CALCD) TO COMPUTE THE JACOBIAN MATRIX, C P(*,*). DERIV IS CALLED NPTS TIMES IF KALCP=-1 . C SIG=YSIG(1) C DO 270 JPT=1,NPTS KPT=JPT IF(KALCP.LT.0) GO TO 210 IF(JPT.NE.1) GO TO 230 C 210 KPT=1 IF(KORDIF.LE.2) GO TO 220 CALL CALCD (JPT,NPTS,P,LPDMA,LPDMB) GO TO 230 C 220 CALL DERIV (JPT,FUNK,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C 230 IF(NTRACE.GE.3) WRITE(KW,240)JPT, * (P(KPT,JX),JX=1,NACTIV) 240 FORMAT(1X,I3,2X,1PE13.5,4E13.5/(6X,5E13.5)) C C COMPUTE QSAV(*,*) AND GRAD(*). C QSAV(*,*), WHICH IS STORED IN THE UPPER TRIANGLE OF C THE ARRAY ERR(*,*), IS PT*P . C GRAD(*) IS HALF THE GRADIENT OF PHI. C IF(LEQU.EQ.0) SIG=YSIG(JPT) C IF(SIG.LE.RZERO) THEN WRITE(KW,245)LEQU,JPT,SIG 245 FORMAT(/' ILLEGAL YSIG IN MARQ: LEQU =',I11/ * 5X,'JPT =',I11,5X,'SIG =',1PG15.7) STOP ENDIF C RTERM=(FIT(JPT)-Y(JPT))/SIG**2 C DO 260 JX=1,NACTIV GRAD(JX)=GRAD(JX)+P(KPT,JX)*RTERM PTERM=P(KPT,JX)/SIG**2 C DO 250 KX=1,JX ERR(JX,KX)=ERR(JX,KX)+P(KPT,KX)*PTERM 250 CONTINUE C 260 CONTINUE C 270 CONTINUE C C RESTORE FIT(*) IF IT WAS DESTROYED IN DERIV. C IF(KORDIF.NE.2) GO TO 290 IF(KALCP.NE.0) GO TO 280 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) NF=NF+1 C 280 NF=NF+NACTIV C 290 NF=NF+NACTIV C C COMPUTE THE SCALE FACTORS AND STORE THEM IN DELTX(*). C SCALE GRAD(*). C DO 300 JX=1,NACTIV SCALJ=ZSQRT(ERR(JX,JX)) IF(SCALJ.EQ.RZERO) SCALJ=UNITR DELTX(JX)=SCALJ GRAD(JX)=GRAD(JX)/SCALJ 300 CONTINUE C IF(NTRACE.GE.1) WRITE(KW,310)(GRAD(JX),JX=1,NACTIV) 310 FORMAT(/' SCALED GRADIENT VECTOR GRAD(*) ='/ * (1X,1PG15.5,4G15.5)) C C SCALE QSAV(*,*). C THE DIAGONAL ELEMENTS OF QSAV(*,*) ARE SCALED TO UNITY. C DO 370 JX=1,NACTIV C DO 360 KX=1,JX SA=ERR(JX,KX)/(DELTX(JX)*DELTX(KX)) IF(KX.NE.JX) GO TO 320 IF(SA.EQ.RZERO) GO TO 330 SA=UNITR GO TO 350 C 320 IF(ZABS(SA).LT.UNITR-RELTOL) GO TO 350 C 330 IF(NTRACE.GE.-2) WRITE(KW,340)JX,KX,SA,ITER 340 FORMAT(' ******** POSSIBLY DANGEROUS VALUE OF', * ' COEFFICIENT....'/6X,'QSAV(',I3,',',I3,') =', * 1PG15.7,10X,'ITERATION',I5) C 350 ERR(JX,KX)=SA 360 CONTINUE C 370 CONTINUE C IF(NTRACE.LT.2) GO TO 400 WRITE(KW,380) 380 FORMAT(/' QSAV(*,*) (PT*P, SCALED, WHERE P IS THE', * ' JACOBIAN MATRIX) ='/' ') C DO 390 JX=1,NACTIV WRITE(KW,240)JX,(ERR(JX,KX),KX=1,JX) 390 CONTINUE C 400 DO 410 JX=1,NV XSAVE(JX)=X(JX) 410 CONTINUE C C INITIALIZE MASKT AND NACT. C 420 NACT=NACTIV C DO 430 JX=1,NV MASKT(JX)=MASK(JX) 430 CONTINUE C C COPY QSAV(*,*) INTO Q(*,*) AND COPY GRAD(*) INTO H(*), C AND SET THE DIAGONAL ELEMENTS OF Q(*,*). C Q(*,*) IS STORED IN THE UPPER TRIANGLE OF ERR(*,*) . C C THIS IS THE ENTRY POINT FOR SUBITERATIONS IN WHICH RLAMBD IS C INCREASED OR CONSTRAINTS ARE IMPOSED. C 440 KRANK=0 JQ=0 JT=0 C DO 470 JX=1,NV IF(MASK(JX).NE.0) GO TO 470 JQ=JQ+1 IF(MASKT(JX).NE.0) GO TO 470 JT=JT+1 H(JT)=-GRAD(JQ) KQ=0 KT=0 C DO 460 KX=1,JX IF(MASK(KX).NE.0) GO TO 460 KQ=KQ+1 IF(MASKT(KX).NE.0) GO TO 460 KT=KT+1 SA=ERR(JQ,KQ) IF(KX.NE.JX .OR. SA.EQ.RZERO) GO TO 450 SA=UNITR+RLAMBD KRANK=KRANK+1 C 450 ERR(KT,JT+1)=SA 460 CONTINUE C 470 CONTINUE C C SOLVE THE NORMAL EQUATIONS FOR H(*), THE CORRECTION VECTOR. C NSMAL=0 NMU=NACT-1 IF(NMU.EQ.0) GO TO 520 C C REDUCE THE SYSTEM TO TRIANGULAR FORM, UTILIZING THE SYMMETRY C OF THE MATRIX. C DO 510 J=1,NMU PIVOT=ERR(J,J+1) IF(PIVOT.EQ.RZERO) GO TO 510 JPU=J+1 C DO 500 K=JPU,NACT EM=ERR(J,K+1)/PIVOT IF(EM.EQ.RZERO) GO TO 490 C DO 480 L=K,NACT ERR(K,L+1)=ERR(K,L+1)-ERR(J,L+1)*EM 480 CONTINUE C 490 H(K)=H(K)-H(J)*EM 500 CONTINUE C 510 CONTINUE C C DO THE BACK SOLUTION. C 520 DO 560 JINV=1,NACT J=(NACT+1)-JINV PIVOT=ERR(J,J+1) IF(PIVOT.LE.RZERO) NSMAL=NSMAL+1 IF(PIVOT.NE.RZERO) GO TO 530 H(J)=RZERO GO TO 560 C 530 SUM=RZERO IF(J.EQ.NACT) GO TO 550 JPU=J+1 C DO 540 K=JPU,NACT SUM=SUM+ERR(J,K+1)*H(K) 540 CONTINUE C 550 H(J)=(H(J)-SUM)/PIVOT 560 CONTINUE C C IF THE COEFFICIENT MATRIX WAS RANK DEFICIENT, PRINT A C MESSAGE. C MRANK=NACT-NSMAL IF(MRANK.EQ.NACT) GO TO 590 COSIN=HUGE IF(NTRACE.GE.-2) WRITE(KW,570)MRANK,NACT,ITER 570 FORMAT(/' RANK-DEFICIENT NORMAL EQUATIONS IN MARQ.'/8X, * 'RANK =',I4,8X,'ORDER OF MATRIX =',I4,8X, * 'ITERATION',I5) IF(MRANK.GT.0) GO TO 580 KFLAG=-4 GO TO 1210 C 580 IF(METHD.GT.0 .AND. MRANK.LT.KRANK) GO TO 870 C C UNPACK AND DE-SCALE THE CORRECTION VECTOR H(*). C COMPUTE THE INNER PRODUCTS SA, SB, AND SC. C 590 SA=RZERO SB=RZERO SC=RZERO KX=NV KQ=NACTIV KT=NACT C DO 620 JX=1,NV HH=RZERO IF(MASK(KX).NE.0) GO TO 610 IF(MASKT(KX).NE.0) GO TO 600 HH=H(KT) SA=SA-HH*GRAD(KQ) SB=SB+HH*HH SC=SC+GRAD(KQ)**2 HH=HH*FMGN/DELTX(KQ) KT=KT-1 C 600 KQ=KQ-1 C 610 H(KX)=HH KX=KX-1 620 CONTINUE C C ADD THE CORRECTION VECTOR TO THE PARAMETER VECTOR AND C CHECK FOR CONSTRAINT VIOLATIONS. C C THIS IS THE ENTRY POINT FOLLOWING A CUTSTEP. C 630 IF(NTRACE.GE.1) WRITE(KW,640)(H(JX),JX=1,NV) 640 FORMAT(/' CORRECTION VECTOR H(*) ='/(1X,1PG15.5,4G15.5)) NACSV=NACT FRMIN=UNITR NLOOP=0 JXLIM=0 C 650 DO 740 JX=1,NV IF(MASKT(JX).NE.0) GO TO 740 XSAV=XSAVE(JX) XMX=XMAX(JX) XMN=XMIN(JX) HH=H(JX) IF((XSAV.LT.XMX .OR. HH.LE.RZERO) .AND. * (XSAV.GT.XMN .OR. HH.GE.RZERO)) GO TO 670 MASKT(JX)=1 NACT=NACT-1 X(JX)=XSAV IF(NTRACE.GE.-1) WRITE(KW,660)JX,XSAV,ITER 660 FORMAT(/' FIX X(',I3,') = ',1PG15.8, * ' TEMPORARILY'/5X,'TO AVOID VIOLATING A', * ' CONSTRAINT.',8X,'ITERATION',I5) GO TO 740 C 670 XLIM=XSAV+HH*FRMIN IF(JX.NE.JXLIM) GO TO 680 IF(KBND.LT.0) GO TO 700 C 680 X(JX)=XLIM IF(XLIM.LE.XMX) GO TO 690 C X(JX)=XMX JBND=1 GO TO 710 C 690 IF(XLIM.GE.XMN) GO TO 740 C 700 X(JX)=XMN JBND=-1 C 710 IF(JX.NE.JXLIM) GO TO 730 IF(NTRACE.LT.-1) GO TO 740 WRITE(KW,720)JX,X(JX),FRMIN,ITER 720 FORMAT(/' CONSTRAINT VIOLATED BY X(',I3, * ').'/' VALUE RESET TO',1PG15.7/6X, * ' USING CUTSTEP FACTOR =',G12.5,8X,'ITERATION',I5) GO TO 740 C 730 IF(NLOOP.NE.0) GO TO 740 DENOM=XLIM-XSAV IF(DENOM.EQ.RZERO) GO TO 740 FRAC=(X(JX)-XSAV)/DENOM IF(FRAC.GE.FRMIN) GO TO 740 FRMIN=FRAC JXLIM=JX KBND=JBND 740 CONTINUE C C IF THE PROPOSED STEP WOULD VIOLATE ANY ALREADY ACTIVE C CONTRAINTS, FIX THOSE COMPONENTS OF H(*) EQUAL TO ZERO AND C RECOMPUTE THE OTHER COMPONENTS. C IF(NACT.GT.0) GO TO 760 KFLAG=3 IF(NTRACE.LT.-2) GO TO 1210 WRITE(KW,750) 750 FORMAT(//' APPARENT CONSTRAINED OPTIMUM LIES IN A', * ' CORNER.') GO TO 1210 C 760 IF(NACT.LT.NACSV) GO TO 440 IF(NLOOP.NE.0) GO TO 770 NLOOP=1 IF(JXLIM.NE.0) GO TO 650 C 770 IF(NTRACE.GE.1) WRITE(KW,780)(X(JX),JX=1,NV) 780 FORMAT(/' X(*) ='/(1X,1PG15.7,4G15.7)) C C CALCULATE THE NEW FITTED VALUES. C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHNEW) NF=NF+1 IF(PHNEW.LT.PHI) GO TO 940 IF(PHNEW.GT.PHI) GO TO 800 C C THE NEW VALUE OF PHI IS EXACTLY EQUAL TO THE OLD VALUE. C CHECK FOR CONVERGENCE UNDER THE NFLAT OPTION. C IF(NFLAT.EQ.0) GO TO 940 KFLAG=2 IF(NTRACE.LT.-1) GO TO 840 WRITE(KW,790) 790 FORMAT(//' CONVERGENCE ACHIEVED UNDER THE NFLAT OPTION', * ' IN MARQ.') GO TO 840 C C THE NEW VALUE OF PHI IS GREATER THAN THE OLD VALUE. C 800 IF(NTRACE.GE.1) WRITE(KW,810)PHI,PHNEW 810 FORMAT(/' OLD PHI =',1PG15.8,8X,'NEW PHI =',G15.8) C C CHECK WHETHER JSUB HAS EXCEEDED MAXSUB. C JSUB=JSUB+1 IF(JSUB.GT.MAXSUB) GO TO 820 IF(METHD.EQ.0) GO TO 1100 IF(METHD.LT.0) GO TO 890 GO TO 860 C 820 KFLAG=-1 IF(NTRACE.GE.-1) WRITE(KW,830)MAXSUB 830 FORMAT(//' EXCEEDED MAXIMUM NUMBER OF SUBITERATIONS =', * I3,' IN MARQ.') C C RESTORE X(*) TO THE BASE POINT. C 840 DO 850 JX=1,NV X(JX)=XSAVE(JX) 850 CONTINUE C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) GO TO 1210 C C THE NEW FIT IS WORSE THAN THE OLD FIT. C COMPUTE COSIN, THE COSINE OF THE ANGLE BETWEEN THE SCALED C GRADIENT AND THE SCALED CORRECTION VECTOR. C COMPUTE DENOM IN SUCH A WAY AS TO AVOID OVERFLOW. C 860 DENOM=ZSQRT(SB)*ZSQRT(SC) IF(DENOM.LE.RZERO) GO TO 870 COSIN=SA/DENOM IF(COSIN.GT.CRIT) GO TO 890 C C COSIN IS NOT GREATER THAN CRIT. C INCREASE THE VALUE OF LAMBDA. C 870 UPFAC=UPNU UPNU=ZMIN1(UPNU*RTWO,FNU) IF(METHD.EQ.1) UPFAC=ZMIN1(UPFAC,FNULM+UNITR/RLAMBD) RLAMBD=RLAMBD*UPFAC IF(NTRACE.GE.1) WRITE(KW,880)JSUB,COSIN,RLAMBD 880 FORMAT(/' **** SUBITERATION',I3,8X,'INCREASE LAMBDA.'/8X, * 'COS(GAMMA) =',1PG12.5,8X,'LAMBDA =',G12.5) C C GO BACK AND FORM THE NORMAL EQUATIONS USING A LARGER VALUE C OF LAMBDA. THIS WILL CAUSE A SHORTER STEP VECTOR H(*) TO BE C COMPUTED. C GO TO 420 C C COSIN IS GREATER THAN CRIT. C CUT THE MAGNITUDE OF THE STEP VECTOR H(*). C 890 STFAC=STFAC/FCUT IF(METHD.GE.0) GO TO 900 FMGN=FMGN/FCUT GO TO 910 C 900 IF(METHD.EQ.1) RLAMBD=RLAMBD*RTWO C 910 DO 920 JX=1,NV H(JX)=(X(JX)-XSAVE(JX))/FCUT 920 CONTINUE C FCUT=FCUT*FACCL IF(NTRACE.GE.1) WRITE(KW,930)JSUB,COSIN,STFAC 930 FORMAT(/' **** SUBITERATION',I3,8X,'TAKE CUT STEPS.'/5X, * 'COS(GAMMA) =',1PG12.5,8X,'CUTSTEP FACTOR =',G12.5) C C GO BACK AND TRY A SMALLER CUTSTEP. C GO TO 630 C C THE VALUE OF PHI HAS DECREASED. TRY A HALF STEP. C 940 IF(METHD.EQ.0 .OR. METHD.EQ.2) GO TO 1100 C DO 950 JX=1,NV XTEMP(JX)=X(JX) IF(MASK(JX).NE.0) GO TO 950 X(JX)=XSAVE(JX)+(X(JX)-XSAVE(JX))/RTWO X(JX)=ZMAX1(XMIN(JX),ZMIN1(XMAX(JX),X(JX))) 950 CONTINUE C DO 960 JPT=1,NPTS FITSV(JPT)=FIT(JPT) 960 CONTINUE C CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHALF) NF=NF+1 C C USE QUADRATIC INTERPOLATION, IN ORDER TO TRY TO REFINE THE C POSITION OF THE MINIMUM OF PHI. C RLFAC=UNITR DENOM=RTWO*((PHNEW-PHALF)-(PHALF-PHI)) STFAC=RZERO IF(DENOM.LE.RZERO) GO TO 970 STFAC=(PHI-PHNEW)/DENOM RSFAC=(UNITR+STFAC)/RTWO C C DO NOT EXTRAPOLATE. C IF(STFAC.GE.UNITR) STFAC=RZERO C 970 DO 980 JX=1,NV H(JX)=X(JX) X(JX)=X(JX)+(XTEMP(JX)-X(JX))*STFAC 980 CONTINUE C IF(PHALF.GE.PHNEW) GO TO 1020 RLFAC=UNITR/RTWO JSUB=JSUB+1 C DO 990 JX=1,NV XTEMP(JX)=H(JX) 990 CONTINUE C DO 1000 JPT=1,NPTS FITSV(JPT)=FIT(JPT) 1000 CONTINUE C IF(NTRACE.GE.1) WRITE(KW,1010)PHNEW,PHALF 1010 FORMAT(/' HALF STEP SUCCEEDED.'/8X,'PHNEW =',1PG15.8,8X, * 'PHALF =',G15.8) PHNEW=PHALF C 1020 IF(STFAC.EQ.RZERO) GO TO 1030 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) NF=NF+1 IF(PHI.LT.PHNEW) GO TO 1070 C 1030 DO 1040 JX=1,NV X(JX)=XTEMP(JX) 1040 CONTINUE C DO 1050 JPT=1,NPTS FIT(JPT)=FITSV(JPT) 1050 CONTINUE C IF(STFAC.NE.RZERO .AND. NTRACE.GE.1) WRITE(KW,1060)RSFAC, * PHI 1060 FORMAT(/' QUADRATIC INTERPOLATION.'/8X,'RSFAC =',1PG12.5, * 8X,'PHI =',G15.8) GO TO 1090 C 1070 RLFAC=RSFAC PHNEW=PHI IF(NTRACE.GE.1) WRITE(KW,1080)RLFAC,PHI 1080 FORMAT(/' QUADRATIC INTERPOLATION SUCCEEDED.'/8X, * 'RLFAC =',1PG12.5,8X,'PHI =',G15.8) C 1090 IF(RLFAC.LE.RZERO) GO TO 1100 RLAMBD=RLAMBD/RLFAC IF(METHD.LT.0) FMGN=FMGN*RLFAC C C THE STEP IS NOW ACCEPTED. C C TEST FOR CONVERGENCE, PROVIDED THAT NO CONSTRAINT BECAME C ACTIVE DURING THIS ITERATION. C 1100 IF(NTRACE.GE.1) WRITE(KW,1110)ITER,PHNEW 1110 FORMAT(/' END OF ITERATION',I5,8X,'PHI =',1PG15.8) PHI=PHNEW IF(JXLIM.GT.0) GO TO 1140 C DO 1120 JX=1,NV IF(MASK(JX).NE.0) GO TO 1120 IF(ZABS(X(JX)-XSAVE(JX)).GT.DELMIN(JX)) GO TO 1140 1120 CONTINUE C KFLAG=1 IF(NTRACE.LT.-1) GO TO 1210 WRITE(KW,1130) 1130 FORMAT(//' MARQ CONVERGED WHEN THE STEP BECAME SMALL.') GO TO 1210 C C THE ITERATION HAS NOT YET CONVERGED. C 1140 IF(ITER.LT.MAXIT) GO TO 1160 KFLAG=-6 IF(NTRACE.GE.-3) WRITE(KW,1150)MAXIT 1150 FORMAT(//' MAXIMUM NUMBER OF ITERATIONS REACHED IN', * ' MARQ.',8X,'MAXIT =',I5) GO TO 1210 C C IF SUBITERATIONS WERE NOT PERFORMED THIS ITERATION, DECREASE C THE VALUE OF LAMBDA. C 1160 IF(NF.GE.NFMAX) GO TO 1190 IF(JSUB.GT.0) GO TO 1170 FMGN=ZMIN1(FMGN*RTWO,UNITR) SCALJ=UNITR+RLAMBD IF(SCALJ.GT.UNITR) RLAMBD=RLAMBD/DOWNU UPNU=DOWNU DOWNU=ZMIN1(DOWNU*RTWO,FNU) GO TO 1180 C 1170 IF(METHD.NE.1) GO TO 1180 UPNU=FNULM DOWNU=FNULM C 1180 CONTINUE C C GO BACK AND DO ANOTHER ITERATION. C GO TO 160 C 1190 KFLAG=-7 IF(NTRACE.GE.-3) WRITE(KW,1200)NFMAX 1200 FORMAT(//' NF HAS REACHED NFMAX =',I8,'IN MARQ.') C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE ITERATION HAS TERMINATED. C C PRINT OUT THE DATA, FITTED VALUES, AND RESIDUALS. C COMPUTE AND PRINT THE STANDARD DEVIATION OF THE DATA FROM C THE FIT. C 1210 SC=ZSQRT(SC) IF(NTRACE.LT.-1) GO TO 1300 WRITE(KW,1220)ITER,NF,PHI,FMGN,RLAMBD,SC 1220 FORMAT(/1X,I4,' ITERATIONS',8X,'NF =',I5,8X,'PHI =', * 1PG15.8//8X,'FMGN =',G12.5,8X,'LAMBDA =',G12.5// * ' NORM OF SCALED GRADIENT VECTOR =',G12.5) WRITE(KW,780)(X(JX),JX=1,NV) IF(NTRACE.GE.0) WRITE(KW,1230) 1230 FORMAT(///4X,'J',5X,'Y(J)',10X,'FIT(J)',6X, * 'Y(J)-FIT(J)',2X,'YSIG(J)',3X,'(Y-FIT)/YSIG'/' ') SIG=YSIG(1) RMSDV=RZERO SDVMX=RZERO JSMAX=1 C DO 1250 JPT=1,NPTS IF(LEQU.EQ.0) SIG=YSIG(JPT) YY=Y(JPT) PTERM=YY-FIT(JPT) RTERM=PTERM/SIG IF(NTRACE.GE.0) WRITE(KW,1240)JPT,YY,FIT(JPT),PTERM, * SIG,RTERM 1240 FORMAT(1X,I4,1PE15.7,E15.7,3E12.4) RMSDV=RMSDV+RTERM**2 IF(ZABS(RTERM).LE.SDVMX) GO TO 1250 SDVMX=ZABS(RTERM) JSMAX=JPT 1250 CONTINUE C DENOM=NPTS-NACTIV WRITE(KW,1260)DENOM 1260 FORMAT(//' NUMBER OF DEGREES OF FREEDOM = ',1PG12.5) IF(DENOM.LE.RZERO) GO TO 1280 RMSDV=ZSQRT(RMSDV/DENOM) WRITE(KW,1270)RMSDV 1270 FORMAT(/' R.M.S. SCALED DEVIATION OF DATA FROM FIT =', * 1PG12.5) C 1280 WRITE(KW,1290)SDVMX,JSMAX 1290 FORMAT(/' MAXIMUM SCALED DEVIATION 0F',1PG12.5/ * 5X,'OCCURRED AT DATA POINT NUMBER',I6) C IF(NTRACE.GE.0) CALL RPLOT (KW,NPTS,Y,YSIG,LEQU,FIT,0) C C CALL FUNC TO SET THE FINAL VALUES OF FIT(*). C 1300 CALL FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) FOBJ=PHI C C CALL MQERR TO PRINT THE PARAMETER ERRORS AND CORRELATIONS. C A DUMMY ROUTINE MAY BE SUBSTITUTED FOR MQERR IF THESE ARE C NOT NEEDED. C IF(MATRX.GT.0) CALL MQERR (NACTIV,NPTS) C RETURN C C END MARQ C END SUBROUTINE MQERR (NACTIV,NPTS) C C MQERR 3.2 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C MQERR IS CALLED BY MARQ TO COMPUTE AND PRINT C APPROXIMATE VALUES OF THE PARAMETER ERRORS AND CORRELATIONS. C C FOR THE MEANING OF THE "MAXIMUM VARIANCE INFLATION FACTOR", C SEE... D. W. MARQUARDT AND R. D. SNEE, C "RIDGE REGRESSION IN PRACTICE", C THE AMERICAN STATISTICIAN 29 (1975) 3-20 C C INPUT QUANTITIES..... KW,ERR(*,*),NACTIV,DELTX(*),NPTS,NV, C NTRACE,MASK(*),FOBJ C C OUTPUT QUANTITIES.... ERR(*,*) C INTEGER JV,JVARY,JX,K,KERFL,KFLAG,KV,KW,KX,L,LINV, * M,MASK,MATRX,NACTIV,NDF,NFLAT,NFMAX,NOREP,NPTS, * NRANK,NSMAL,NTRACE,NV,NVPLU,NXTRA C DOUBLE PRECISION SCALJ,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * ZSQRT,ARG,RZERO,UNITR,HUGE,PIVOT,Q,VIFMX,ER,TEMP, * DENOM,SCFAC,RESCL ,X, DSQRT C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C ZSQRT(ARG)=DSQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C RZERO=0.0D0 UNITR=1.0D0 HUGE=1.0D30 C C PRINT QSAV(*,*). C IF(NTRACE.LT.-1) GO TO 40 WRITE(KW,10) 10 FORMAT(///' SUBROUTINE MQERR.'//' QSAV(*,*) (PT*P,', * ' SCALED, WHERE P IS THE JACOBIAN) =') DO 30 JX=1,NACTIV WRITE(KW,20)JX,(ERR(JX,KX),KX=1,JX) 20 FORMAT(/1X,I3,2X,1PG13.5,4G13.5/(6X,5G13.5)) 30 CONTINUE C C COMPUTE THE SCALED ERROR MATRIX, WHICH IS THE INVERSE OF C QSAV. C C INVERT QSAV USING THE GAUSS-JORDAN METHOD WITHOUT PIVOTING. C F. L. BAUER AND C. REINSCH, P. 45 IN "LINEAR ALGEBRA" C BY J. H. WILKINSON AND C. REINSCH (SPRINGER-VERLAG, 1971). C ERR(*,NV+1) IS USED AS A SCRATCH VECTOR. C 40 NVPLU=NV+1 NSMAL=0 C DO 140 LINV=1,NACTIV L=(NACTIV+1)-LINV PIVOT=ERR(1,1) IF(PIVOT.LE.RZERO) NSMAL=NSMAL+1 IF(NACTIV.LT.2) GO TO 100 C DO 90 K=2,NACTIV Q=ERR(K,1) IF(PIVOT.NE.RZERO) GO TO 50 ERR(K,NVPLU)=RZERO GO TO 70 C 50 IF(K.GT.L) GO TO 60 ERR(K,NVPLU)=-Q/PIVOT GO TO 70 C 60 ERR(K,NVPLU)=Q/PIVOT C 70 DO 80 M=2,K C ERR(K-1,M-1)=ERR(K,M)+Q*ERR(M,NVPLU) 80 CONTINUE C 90 CONTINUE C 100 IF(PIVOT.NE.RZERO) GO TO 110 ERR(NACTIV,NACTIV)=RZERO GO TO 120 C 110 ERR(NACTIV,NACTIV)=UNITR/PIVOT C 120 IF(NACTIV.LT.2) GO TO 140 C DO 130 K=2,NACTIV ERR(NACTIV,K-1)=ERR(K,NVPLU) 130 CONTINUE C 140 CONTINUE C NRANK=NACTIV-NSMAL IF(NRANK.LT.NACTIV .AND. NTRACE.GE.-3) * WRITE(KW,150)NRANK,NACTIV 150 FORMAT(//' THE SECOND DERIVATIVE MATRIX IS SINGULAR', * ' IN MQERR.'/8X,'RANK =',I3,8X,'ORDER =',I3/ * ' THEREFORE ALL PARAMETER ERRORS ARE INFINITE.') C C UNPACK THE ERROR MATRIX INTO THE UPPER C TRIANGLE OF ERR(*,*), DE-SCALING IT. C JV=0 VIFMX=RZERO C DO 210 JX=1,NV IF(MASK(JX).EQ.0) JV=JV+1 KV=0 C DO 200 KX=1,JX ER=RZERO IF(MASK(JX).NE.0 .OR. MASK(KX).NE.0) GO TO 190 KV=KV+1 TEMP=RZERO DENOM=DELTX(JV)*DELTX(KV) IF(DENOM.EQ.RZERO) GO TO 160 TEMP=ERR(JV,KV) ER=TEMP/DENOM C 160 IF(JX.NE.KX) GO TO 190 IF(TEMP.GT.RZERO) GO TO 180 IF(NTRACE.GE.-2) WRITE(KW,170)JX,KX,TEMP 170 FORMAT(/' THE (',I3,',',I3, * ') ELEMENT OF QSAV**-1 =',1PG13.5/ * ' . THEREFORE ALL PARAMETER ERRORS ARE', * ' INFINITE.') TEMP=-TEMP C 180 IF(TEMP.GT.VIFMX) VIFMX=TEMP C 190 ERR(KX,JX+1)=ER C 200 CONTINUE C 210 CONTINUE C C COMPUTE AND PRINT THE STANDARD ERRORS. C NDF=NPTS-NACTIV SCFAC=HUGE IF(NDF.LE.0) GO TO 220 SCFAC=NDF SCFAC=ZSQRT(FOBJ/SCFAC) C 220 RESCL=NDF+NDF IF(RESCL.GT.RZERO) RESCL=ZSQRT(RESCL) IF(NTRACE.LT.-1) GO TO 350 WRITE(KW,230)NDF,NDF,RESCL,FOBJ,SCFAC 230 FORMAT(///' NUMBER OF DEGREES OF FREEDOM (N.D.F.) =', * ' (NPTS-NACTIV) =',I5//' EXPECTED VALUE OF PHI =', * ' N.D.F. PLUS OR MINUS SQRT(2*N.D.F.)'/23X,'=',I5, * ' PLUS OR MINUS',1PG12.5// * ' ACTUAL VALUE OF PHI =',G12.5// * ' RESCALING FACTOR = SQRT(PHI/N.D.F.) =',G12.5) WRITE(KW,240)VIFMX 240 FORMAT(///' MAXIMUM VARIANCE INFLATION FACTOR =',1PG12.5, * ' = CONDITION NUMBER'/// * ' APPROXIMATE STANDARD ERRORS....'/ * 57X,'RESCALED'/3X,'J',6X,'MASK(J)',7X,'X(J)', * 13X,'ERROR',12X,'ERROR') RESCL=HUGE C DO 280 JX=1,NV SCALJ=UNITR ER=ERR(JX,JX+1) IF(ER.GE.RZERO) GO TO 250 ER=-ZSQRT(-ER) SCALJ=-ER GO TO 260 C 250 ER=ZSQRT(ER) SCALJ=ER C 260 IF(NDF.GT.0 .AND. NRANK.EQ.NACTIV) RESCL=SCFAC*ER DELTX(JX)=SCALJ WRITE(KW,270)JX,MASK(JX),X(JX),ER,RESCL 270 FORMAT(/1X,I3,I10,4X,1PG16.8,4X,G13.5,4X,G13.5) 280 CONTINUE C C COMPUTE AND PRINT THE CORRELATIONS. C IF(NV.LT.2) GO TO 350 WRITE(KW,290)(K,K=1,NV) 290 FORMAT(///' LOWER TRIANGLE OF THE CORRELATION MATRIX....' * //11X,'K =',I3,4I12/(6X,5I12)) WRITE(KW,300)(MASK(K),K=1,NV) 300 FORMAT(/6X,'MASK(K) =',I3,4I12/(6X,5I12)) WRITE(KW,310) 310 FORMAT(/' J MASK(J)') C DO 340 JX=1,NV C DO 320 KX=1,JX DENOM=DELTX(JX)*DELTX(KX) ERR(KX,1)=RZERO IF(DENOM.NE.RZERO) ERR(KX,1)=ERR(KX,JX+1)/DENOM 320 CONTINUE C WRITE(KW,330)JX,MASK(JX),(ERR(KX,1),KX=1,JX) 330 FORMAT(/1X,I3,I6,2X,1PG12.4,4G12.4/(12X,5G12.4)) 340 CONTINUE C C RESCALE ERR AND SYMMETRIZE IT. C 350 SCFAC=SCFAC**2 IF(SCFAC.LT.UNITR) SCFAC=UNITR C DO 370 JX=1,NV C DO 360 KX=1,JX ERR(JX,KX)=ERR(KX,JX+1)*SCFAC ERR(KX,JX)=ERR(JX,KX) 360 CONTINUE C 370 CONTINUE C RETURN C C END MQERR C END SUBROUTINE FUNC (FUNK,Y,YSIG,NPTS,FIT,PHI) C C FUNC CALLS FUNK OR FOFX TO COMPUTE THE ARRAY OF FITTED C VALUES, FITM(*), FOR THE MARQ PACKAGE. . C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C INTEGER JPT,JVARY,KALCP,KERFL,KFLAG,KORDIF,KW, * LEQU,MASK,MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT, * NFMAX,NOREP,NPTS,NTRACE,NV,NXTRA C DOUBLE PRECISION Y,YSIG,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * FLAMBD,FNU,RELDIF,RELMIN,RZERO DOUBLE PRECISION X,FIT,F,PHI,SIG C DIMENSION Y(1),YSIG(1),FIT(1) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C RZERO=0.0D0 C IF(KALCP.NE.0) GO TO 10 C CALL FUNK (FIT) GO TO 30 C 10 DO 20 JPT=1,NPTS CALL FOFX(JPT,NV,X,F) FIT(JPT)=F 20 CONTINUE C 30 PHI=RZERO SIG=YSIG(1) C DO 60 JPT=1,NPTS IF(LEQU.EQ.0) SIG=YSIG(JPT) C C CHECK FOR AN ILLEGAL VALUE OF SIG. C IF(SIG.GT.RZERO) GO TO 50 WRITE(KW,40)LEQU,JPT,SIG 40 FORMAT(/' ERROR IN FUNC.... LEQU =',I2/ * 6X,'JPT =',I5,6X,'SIG =',1PG13.5,' IS .LE. ZERO.') STOP C 50 PHI=PHI+((FIT(JPT)-Y(JPT))/SIG)**2 60 CONTINUE C RETURN C C END FUNC C END SUBROUTINE DERIV (JPT,FUNK,NPTS,FIT,FITSV,P,LPDMA,LPDMB) C C DERIV 4.3 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C DERIV APPROXIMATES THE JACOBIAN MATRIX P(*,*) FOR THE C MARQ PACKAGE, USING DIFFERENCE QUOTIENTS. C C P(J,K) IS THE PARTIAL DERIVATIVE OF FIT(J) WITH RESPECT TO C X(K) . C IF KORDIF.EQ.1, DERIV USES A NONCENTRAL DIFFERENCE FORMULA. C IF KORDIF.EQ.2, DERIV USES A CENTRAL DIFFERENCE FORMULA. C KORDIF.EQ.1 IS ABOUT TWICE AS FAST AS KORDIF.EQ.2, BUT LESS C ACCURATE. C INTEGER J,JPT,JVARY,JX,KALCP,KERFL,KFLAG,KORDIF, * KW,KX,LEQU,LPDMA,LPDMB,MASK,MATRX,MAXIT,METHD, * MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP,NPTS,NTRACE,NV, * NXTRA C DOUBLE PRECISION P,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * FLAMBD,FNU,RELDIF,RELMIN,RZERO DOUBLE PRECISION X,FIT,FITSV,DEL,ABSX,XSAVE,FX0,FX1, * XPLUS,XMINUS C DIMENSION FIT(1),FITSV(1),P(LPDMA,LPDMB) C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C RZERO=0.0D0 C JVARY=0 C C SAVE FIT IF KALCP.GE.0 . C IF(KALCP.LT.0) GO TO 20 C DO 10 J=1,NPTS FITSV(J)=FIT(J) 10 CONTINUE C C LOOP OVER THE ACTIVE PARAMETERS X(JX). C 20 KX=0 C DO 160 JX=1,NV IF(MASK(JX).NE.0) GO TO 160 KX=KX+1 XSAVE=X(JX) ABSX=XSAVE IF(ABSX.LT.RZERO) ABSX=-ABSX DEL=RELDIF*ABSX IF(DEL.EQ.RZERO) DEL=RELDIF XPLUS=XSAVE+DEL X(JX)=XPLUS XMINUS=XSAVE-DEL IF(XPLUS.NE.XSAVE .AND. XPLUS.NE.XMINUS) GO TO 40 WRITE(KW,30)RELDIF,JX,XSAVE,DEL,XPLUS,XMINUS 30 FORMAT(/' ****** UNABLE TO COMPUTE A SUITABLE VALUE', * ' FOR DEL IN SUBROUTINE DERIV.'/6X,'RELDIF =', * 1PG15.7,6X,'JX =',I11,6X,'XSAVE =',G15.7,6X,'DEL =' * ,G15.7/6X,'XPLUS =',G15.7,6X,'XMINUS =',G15.7) STOP C 40 IF(KALCP.LT.0) GO TO 120 IF(KALCP.GT.0) GO TO 90 C C KALCP.EQ.0 . COMPUTE P(*,*), ONE COLUMN AT A TIME. C CALL FUNK (FIT) IF(KORDIF.EQ.2) GO TO 60 C DO 50 J=1,NPTS P(J,KX)=(FIT(J)-FITSV(J))/(XPLUS-XSAVE) 50 CONTINUE C GO TO 150 C C KALCP.EQ.0 AND KORDIF.EQ.2 . C IN THIS CASE, THE INPUT VALUES OF FIT(J) WILL BE DESTROYED. C 60 X(JX)=XMINUS C DO 70 J=1,NPTS FITSV(J)=FIT(J) 70 CONTINUE C JVARY=JX CALL FUNK (FIT) JVARY=0 C DO 80 J=1,NPTS P(J,KX)=(FITSV(J)-FIT(J))/(XPLUS-XMINUS) 80 CONTINUE C GO TO 150 C C KALCP.GT.0 . COMPUTE P(*,*), ONE ELEMENT AT A TIME. C 90 DO 110 J=1,NPTS CALL FOFX (J,NV,X,FX1) IF(KORDIF.EQ.2) GO TO 100 P(J,KX)=(FX1-FITSV(J))/(XPLUS-XSAVE) GO TO 110 C 100 X(JX)=XSAVE-DEL CALL FOFX (J,NV,X,FX0) P(J,KX)=(FX1-FX0)/(XPLUS-XMINUS) X(JX)=XPLUS 110 CONTINUE C GO TO 150 C C KALCP.LT.0 . C COMPUTE ONE ROW OF P(*,*), ONE ELEMENT AT A TIME. C 120 FITSV(1)=FIT(JPT) CALL FOFX (JPT,NV,X,FX1) IF(KORDIF.EQ.2) GO TO 130 P(1,KX)=(FX1-FITSV(1))/(XPLUS-XSAVE) GO TO 140 C 130 X(JX)=XMINUS CALL FOFX (JPT,NV,X,FX0) P(1,KX)=(FX1-FX0)/(XPLUS-XMINUS) C 140 FIT(JPT)=FITSV(1) C C RESTORE X(JX). C 150 X(JX)=XSAVE 160 CONTINUE C IF(KALCP.LT.0) RETURN C DO 170 J=1,NPTS FIT(J)=FITSV(J) 170 CONTINUE C RETURN C C END DERIV C END SUBROUTINE FOFX (JPT,NV,X,F) C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C THIS IS A DUMMY VERSION OF SUBROUTINE FOFX, FOR THE MARQ C PACKAGE. C A NON-DUMMY VERSION OF FOFX MAY BE USED (OPTIONALLY) TO C SUPPLY TO SUBROUTINE MARQ THE VALUES OF THE FUNCTION BEING C FITTED, INSTEAD OF USING A 'FUNK' SUBROUTINE TO DO THIS. C THE USE OF FOFX REQUIRES SUBSTANTIALLY MORE OVERHEAD TIME C DURING EXECUTION, BUT SAVES CONSIDERABLE STORAGE BY NOT C REQUIRING THAT THE JACOPIAN MATRIX P(*,*) BE STORED. C INTEGER JPT,NV C DOUBLE PRECISION X,F C DIMENSION X(20) C RETURN END SUBROUTINE CALCD (JPT,NPTS,P,LPDMA,LPDMB) C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C THIS IS A DUMMY VERSION OF SUBROUTINE CALCD, FOR THE MARQ C FOR THE MARQ PACKAGE. C A NON-DUMMY VERSION OF CALCD MAY BE USED (OPTIONALLY) TO C SUPPLY TO SUBROUTINE MARQ ANALYTIC VALUES OF THE ELEMENTS C OF THE JACOBIAN MATRIX, INSTEAD OF APPROXIMATION THE C ELEMENTS USING FINITE DIFFERENCES. C HOWEVER, MOST USERS PREFER TO USE FINITE DIFFERENCES. C INTEGER JPT,LPDMA,LPDMB C DOUBLE PRECISION P C DIMENSION P(LPDMA,LPDMB) C RETURN END SUBROUTINE STSET C C STSET 3.2 DECEMBER 1991 C C STSET SETS SOME INPUT QUANTITIES TO DEFAULT VALUES, FOR C SUBROUTINES STEPIT, SIMPLEX, MARQ, STP, MINF, OR KAUPE. C C J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, C OKLAHOMA STATE UNIVERSITY C C C USAGE..... C C CALL STSET. C THEN SET SOME INPUT QUANTITIES (NV, AT LEAST) AND RESET ANY C OF THOSE SET IN STSET (BETTER VALUES OF X(J), ETC.) BEFORE C CALLING STEPIT OR SIMPLEX OR MARQ, ETC. C C INTEGER JVARY,JX,KALCP,KERFL,KFLAG,KORDIF,KW,LEQU,MASK, * MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP, * NTRACE,NV,NVMAX,NXTRA C DOUBLE PRECISION XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,FLAMBD, * FNU,RELDIF,RELMIN,RZERO,HUGE, X C COMMON /CSTEP/ X(20),XMAX(20),XMIN(20),DELTX(20), * DELMIN(20),ERR(20,21),FOBJ,NV,NTRACE,MATRX,MASK(20), * NFMAX,NFLAT,JVARY,NXTRA,KFLAG,NOREP,KERFL,KW C COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, * KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD C C HUGE=1.0D30 C RZERO=0.0D0 C C KW = LOGICAL UNIT NUMBER OF THE PRINTER C KW=6 C C NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV, GIVEN THE C PRESENT DIMENSIONS OF ARRAYS. C NVMAX IS THE DIMENSION OF THE ARRAYS X(*), XMAX(*), XMIN(*), C DELTX(*), DELMIN(*), AND MASK(*). NVMAX IS ALSO THE FIRST C DIMENSION OF ERR(*,*). THE SECOND DIMENSION OF ERR(*,*) C IS NVMAX+1. C NVMAX=20 C C THE USER MUST SET NV AFTER CALLING STSET. C NV=-1 C NTRACE=0 NFMAX=1000000 MAXIT=50 MAXSUB=30 METHD=1 KALCP=0 LEQU=0 NFLAT=1 MATRX=105 NXTRA=0 FLAMBD=1.0D0 FNU=10.0D0 C KORDIF=1 RELDIF=1.0D-8 RELMIN=1.0D-6 C C FOR GREATER ACCURACY BUT LESS SPEED IN MARQ, SET... C C KORDIF=2 C RELDIF=O.4D-5 C RELMIN=1.0D-9 C DO 10 JX=1,NVMAX X(JX)=RZERO XMAX(JX)=HUGE XMIN(JX)=-HUGE DELTX(JX)=RZERO DELMIN(JX)=RZERO MASK(JX)=0 10 CONTINUE C RETURN C C END STSET C END SUBROUTINE RPLOT (LP,NPTS,Y,YSIG,LEQU,FIT,MLARGE) C C RPLOT 2.2 JANUARY 1992 C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C RPLOT PRINTS A PLOT OF DATA AND FITTED VALUES, SIDE BY SIDE C WITH A PLOT OF NORMALIZED RESIDUALS. C C THE X AXIS RUNS FROM TOP TO BOTTOM ON THE PAPER, AND THE Y C AXIS RUNS FROM LEFT TO RIGHT. C THE POINTS ARE PLOTTED ONE PER LINE, REGARDLESS OF THE C VALUES OF THE ABSCISSAS. C C C INPUT QUANTITIES.... C C LP -- LOGICAL UNIT NUMBER OF THE PRINTER C C NPTS -- NUMBER OF POINTS C C Y(J) -- ORDINATE OF THE J-TH POINT C C YSIG(J) -- EXPECTED ERROR IN Y(J) C C LEQU -- =1 IF ALL YSIG(J) ARE EQUAL (IN THIS CASE C ONLY YSIG(1) IS REFERENCED), C C =0 OTHERWISE C C FIT(J) -- FITTED VALUE AT THE J-TH POINT C C MLARGE -- =0 TO PRINT IN 72 COLUMNS, C =1 TO PRINT IN 120 COLUMNS C C INTEGER LP,NPTS,LEQU,MLARGE, * NR,NY,JPT,J,JYZERO,JRZERO,K,L,JY,JF,JR, * JJ,MAXMSG,NBADSG C CHARACTER*1 KBL,KY,KF,KBOTH,KR,KI,KPLUS,KMINUS,KYLINE, * KRLINE C DOUBLE PRECISION Y,YSIG,FIT, * RZERO,UNITR,R15, * YMAX,YMIN,RMAX,RMIN,YY,F,SIG,R,DY,DR C DIMENSION Y(1),YSIG(1),FIT(1) DIMENSION KYLINE(81),KRLINE(21) C KBL=' ' KY='Y' KF='F' KBOTH='2' KR='R' KI='I' KPLUS='+' KMINUS='-' C RZERO=0.0D0 UNITR=1.0D0 R15=1.5D0 C MAXMSG=3 NR=21 C IF(MLARGE.EQ.0) THEN NY=41 ELSE NY=81 ENDIF C C COMPUTE THE EXTREMES OF Y(*) AND FIT(*), AND OF ABS(R(*)) . C YMAX=Y(1) YMIN=YMAX RMAX=RZERO JPT=1 NBADSG=0 C DO 10 J=1,NPTS YY=Y(J) IF(YY.GT.YMAX) YMAX=YY IF(YY.LT.YMIN) YMIN=YY F=FIT(J) IF(F.GT.YMAX) YMAX=F IF(F.LT.YMIN) YMIN=F IF(LEQU.EQ.0) JPT=J SIG=YSIG(JPT) C IF(SIG.EQ.RZERO) THEN SIG=UNITR NBADSG=NBADSG+1 IF(NBADSG.LE.MAXMSG) WRITE(LP,110)J,JPT 110 FORMAT(/' ERROR IN SUBROUTINE RPLOT ...'/ * 5X,'SIG = ZERO FOR J =',I6,' AND JPT =',I6/ * 5X,'SIG = 1.0 WILL BE USED INSTEAD.') IF(NBADSG.EQ.MAXMSG+1) WRITE(LP,115) 115 FORMAT(/' FURTHER ERROR MESSAGES SUPPRESSED.') ENDIF C R=(YY-FIT(J))/SIG IF(R.GT.RMAX) RMAX=R IF(-R.GT.RMAX) RMAX=-R 10 CONTINUE C IF(YMAX.LT.RZERO) YMAX=RZERO IF(YMIN.GT.RZERO) YMIN=RZERO C IF(YMAX.EQ.RZERO .AND. YMIN.EQ.RZERO) THEN YMAX=UNITR YMIN=-UNITR ENDIF C IF(RMAX.EQ.RZERO) RMAX=UNITR RMIN=-RMAX C C PRINT THE HEADING. C IF(MLARGE.EQ.0) THEN WRITE(LP,20) 20 FORMAT(/13X,'Y = DATA',8X,'F = FIT',14X, * 'NORMALIZED RESIDUALS') WRITE(LP,30)YMIN,YMAX,RMIN,RMAX 30 FORMAT(1X,1PE9.2,28X,E9.2,1X,E8.1,9X,E8.1) ELSE WRITE(LP,40) 40 FORMAT(31X,'Y = DATA',12X,'F = FIT',32X, * 'NORMALIZED RESIDUALS') WRITE(LP,50)YMIN,YMAX,RMIN,RMAX 50 FORMAT(1X,1PE9.2,68X,E9.2,1X,E8.1,9X,E8.1) ENDIF C DY=NY-1 DY=(YMAX-YMIN)/DY DR=NR-1 DR=(RMAX-RMIN)/DR JYZERO=-YMIN/DY+R15 C DO 150 JJ=1,2 C DO 60 K=2,NY KYLINE(K)=KMINUS 60 CONTINUE C KYLINE(1)=KPLUS KYLINE(NY)=KPLUS KYLINE(JYZERO)=KPLUS C IF(MLARGE.EQ.0) THEN WRITE(LP,70)(KYLINE(K),K=1,NY) 70 FORMAT(4X,41A1,5X,'+---------+---------+') ELSE WRITE(LP,80)(KYLINE(K),K=1,NY) 80 FORMAT(4X,81A1,5X,'+---------+---------+') ENDIF C IF(JJ.EQ.2) GO TO 150 C C PLOT THE TWO GRAPHS SIDE BY SIDE. C JRZERO=(NR+1)/2 C DO 90 K=1,NY KYLINE(K)=KBL 90 CONTINUE C DO 100 K=1,NR KRLINE(K)=KBL 100 CONTINUE C JPT=1 L=0 C DO 140 J=1,NPTS L=L+1 IF(L.GT.999) L=0 YY=Y(J) JY=(YY-YMIN)/DY+R15 IF(JY.LT.1) JY=1 IF(JY.GT.NY) JY=NY KYLINE(JYZERO)=KI F=FIT(J) JF=(F-YMIN)/DY+R15 IF(JF.LT.1) JF=1 IF(JF.GT.NY) JF=NY IF(LEQU.EQ.0) JPT=J SIG=YSIG(JPT) IF(SIG.EQ.RZERO) SIG=UNITR R=(YY-F)/SIG JR=JRZERO IF(DR.NE.RZERO) JR=(R-RMIN)/DR+R15 KYLINE(JY)=KY KYLINE(JF)=KF IF(JY.EQ.JF) KYLINE(JY)=KBOTH KRLINE(JRZERO)=KI KRLINE(JR)=KR C IF(MLARGE.EQ.0) THEN WRITE(LP,120)(KYLINE(K),K=1,NY),L, * (KRLINE(K),K=1,NR) 120 FORMAT(4X,41A1,I4,1X,21A1) ELSE WRITE(LP,130)(KYLINE(K),K=1,NY),L, * (KRLINE(K),K=1,NR) 130 FORMAT(4X,81A1,I4,1X,21A1) ENDIF C C RESET THE PRINT LINE TO BLANKS. C KYLINE(JF)=KBL KYLINE(JY)=KBL KRLINE(JR)=KBL C 140 CONTINUE C 150 CONTINUE C IF(MLARGE.EQ.0) THEN WRITE(LP,30)YMIN,YMAX,RMIN,RMAX ELSE WRITE(LP,50)YMIN,YMAX,RMIN,RMAX ENDIF C WRITE(LP,160) 160 FORMAT(' ') C RETURN C C END RPLOT C END $ENTRY 1 12 0.0 1.0 16.6 0.2 1.0 0.0 12.4 0.2 1.0 1.0 8.2 0.2 0.0 2.0 10.1 0.2 2.0 0.0 7.3 0.2 2.0 2.0 4.7 0.2 1.0 3.0 5.1 0.2 3.0 1.0 4.3 0.2 3.0 3.0 3.0 0.2 5.0 1.0 2.5 0.2 1.0 5.0 3.4 0.2 5.0 5.0 2.1 0.2 //