//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