//U11701AA JOB (11701,123-45-6789),TIME=(,5),CLASS=1,MSGCLASS=Y, // NOTIFY=* /*JOBPARM ROOM=H,FORMS=9021 // EXEC WATFIV,REGION.WATFIV=500K $JOB ,XREF,LIBLIST C C STEPIT.CNTL CURRENT VERSIONS OF STEPIT AND FIDO C SUBROUTINE STINT (FUNK) C C COMMON INTERFACE ROUTINE TO FACILITATE SWAPPING OF STEPIT, C SIMPLEX, MARQ, ETC. C EXTERNAL FUNK CALL STEPIT (FUNK) RETURN END 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 DOUBLE PRECISION ERGSAV,FIDINT, * X,XMAX,XMIN,DELTX,DELMIN,ERR, * FIT,Y,YSIG,T, * ERFRAC,ERGUES,FLAMBD,FNU,FOBJ,RELDIF,RELMIN, * STDEVS,TOLFID, * DABS,DSQRT C INTEGER MASK, * JROUTE,JVARY,JXFID,K,KALCP,KERFL,KFLAG,KORDIF,KR,KW, * LEQU,MATRX,MAXIND,MAXIT,MAXSUB,METHD,NFLAT,NFMAX, * NOREP,NPTS,NTRACE,NTRACF,NTRSET,NV,NXTRA,MAXUPD 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 KR=5 KW=6 C C OPEN (UNIT=5,FILE="STEPTEST.DAT") 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 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 DOUBLE PRECISION FITM,FIT,Y,YSIG,T, * X,XMAX,XMIN,DELTX,DELMIN,ERR, * FOBJ C INTEGER MASK, * J,JVARY,KERFL,KFLAG,KW,MATRX,NFLAT,NFMAX,NOREP,NPTS, * NTRACE,NV,NXTRA 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.5 DECEMBER 1991 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)-YDATA(JPT))/YSIGMA(JPT))**2 C JPT=1 C C THE STANDARD ERRORS YSIGMA( ) MUST BE CORRECTLY SCALED. C IF THE YSIGMA( ) 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 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 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 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 STEPIT (FUNK) C C STEPIT 7.7 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER C (PRESENT ADDRESS .... COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY, C STILLWATER, OKLAHOMA 74078 C (405)-744-5676 ) C C STEPIT FINDS LOCAL MINIMA OF A SMOOTH FUNCTION OF SEVERAL C PARAMETERS. C C "STEPIT IS A PHLEGMATIC METHOD OF SOLVING A PROBLEM." C -- J. H. BURRILL, JR., 360 STEPIT - A USER'S MANUAL C C STEPIT 7.4 AND A WRITE-UP ARE AVAILABLE FROM THE C QUANTUM CHEMISTRY PROGRAM EXCHANGE C DEPT. OF CHEMISTRY, INDIANA UNIVERSITY C BLOOMINGTON, INDIANA 47401 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,NV,NTRACE,MATRX,MASK,X,XMAX, C XMIN,DELTX,DELMIN,NFMAX,NFLAT,KW C C OUTPUT QUANTITIES.... X,FOBJ,ERR, KFLAG,NOREP,KERFL C C FUNK -- THE NAME OF THE SUBROUTINE THAT COMPUTES C FOBJ GIVEN X(1),X(2),...,X(NV) (EACH C SUCH SUBROUTINE MUST BE NAMED IN AN C EXTERNAL STATEMENT IN THE CALLING C PROGRAM) C C NV -- THE NUMBER OF PARAMETERS, X C C NTRACE -- = 0 FOR NORMAL OUTPUT, C =+1 FOR TRACE OUTPUT, C =-1 FOR NO OUTPUT C C MATRX -- = 0 FOR NO ERROR CALCULATION, C = 100+M TO APPROXIMATE THE ERRORS IN THE C X(J) USING STEPS 10**(-M) TIMES AS C LARGE AS X(J), IF NONZERO C C FOBJ -- THE VALUE OF THE FUNCTION TO BE MINIMIZED C C MASK(J) -- NONZERO IF X(J) IS TO BE HELD FIXED C C X(J) -- THE J-TH PARAMETER C C XMAX(J) -- THE UPPER LIMIT ON X(J) C C XMIN(J) -- THE LOWER LIMIT ON X(J) C C DELTX(J) -- THE INITIAL STEP SIZE FOR X(J) C C DELMIN(J) -- THE LOWER LIMIT (CONVERGENCE TOLERANCE) C ON THE STEP SIZE FOR X(J) C C ERR(J,K) -- RETURNS THE ERROR MATRIX IF -MATRX- C IS NONZERO C (ERR IS ALSO USED FOR SCRATCH STORAGE) C C NFMAX -- THE MAXIMUM NUMBER OF FUNCTION C EVALUATIONS C C NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN C ALL TRIAL STEPS GIVE IDENTICAL FUNCTION C VALUES. THE RECOMMENDED VALUE OF NFLAT C IS USUALLY NFLAT=1 . C C JVARY -- STEPIT SETS JVARY NONZERO IF X(JVARY) IS C THE ONLY X(J) THAT HAS CHANGED SINCE C THE LAST CALL TO FUNK C (THIS CAN BE USED TO SPEED UP FUNK) C C NXTRA -- USED BY SUBROUTINE SIMPLEX BUT NOT BY C STEPIT C C KFLAG -- RETURNED .GT. ZERO FOR A NORMAL EXIT, C RETURNED .LT. ZERO FOR AN ABNORMAL EXIT C C NOREP -- RETURNED .GT. ZERO IF THE FUNCTION WAS NOT C REPRODUCIBLE C C KERFL -- RETURNED .LT. ZERO IF SUBROUTINE STERR C TERMINATED ABNORMALLY C C KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER 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 C STEPIT SHOULD USUALLY BE RUN USING A FLOATING POINT C PRECISION OF AT LEAST TEN SIGNIFICANT DIGITS. ON MOST C COMPUTERS, THIS REQUIRES THE USE OF DOUBLE PRECISION. C DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * VEC,DLX,XS,FSTORE,DX,SALVO,XOSC,FOSC,ARG,STCUT,ACK, * FACUP DOUBLE PRECISION RZERO,XPLUS, * FSAVE,FBEST,XSAVE,ABSDX,FPREV,DENOM,DEL,DXZ,DXU,DFZ, * DFU,ABSVEC,SUMV,CINDER,COXCOM,COSIN,STEPS,ZSQRT, * DSQRT C INTEGER J,JFLAT,JFLMIN,JOCK,JUMP,JVARY,JX,K,KERFL,KFLAG, * KL,KW,MASK,MATRX,MINOSC,MAXOSC,MAXSTP,NACK,NACTIV, * NAH,NCIRC,NEQUAL,NF,NFLAT,NFMAX,NFSAV,NGATE,NGIANT, * NONZER,NOREP,NOSC,NOUT,NRETRY,NSSW,NSTEPS,NT,NTRACE, * NV,NXTRA,NZIP C C THE DIMENSIONS OF ALL VECTORS AND MATRICES C (AS OPPOSED TO ARRAYS) ARE NV, EXCEPT FOR .... C ERR(NV,MAXOSC),XOSC(NV,MAXOSC),FOSC(MAXOSC). C IF ERRORS ARE TO BE CALCULATED BY SUBROUTINE STERR, HOWEVER, C THEN ERR MUST BE DIMENSIONED AT LEAST AS LARGE AS C ERR(NV,MAX(NV,MAXOSC)) . C DIMENSION VEC(20),FSTORE(20),SALVO(20),JFLAT(20) DIMENSION XOSC(20,5),FOSC(5) 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 INTERNAL STEPIT COMMON..... COMMON /STORK/ DX(20),XS(20),DLX(20),NACTIV,NSSW,NF C C THE ONLY SUBROUTINES CALLED ARE FUNK, STBEG, STERR, AND C DATSW. C C STEPIT TERMINATES IF SENSE SWITCH NUMBER NSSW IS ON. C THE STATEMENT CALL DATSW(NSSW,JUMP) RETURNS JUMP=1 IF C SENSE SWITCH NUMBER NSSW IS ON, AND JUMP=2 IF IT IS OFF. C IF NO SENSE SWITCH IS TO BE USED, SUPPLY A DUMMY SUBROUTINE C FOR DATSW. C C SET THE LIBRARY FUNCTION FOR SINGLE PRECISION (SQRT) OR FOR C DOUBLE PRECISION (DSQRT). NO OTHER FUNCTIONS ARE USED, C EITHER EXTERNAL OR INTRINSIC, EXCEPT THE ROUTINE INVOKED BY C REAL**INTEGER . C C ZSQRT(ARG)=SQRT(ARG) ZSQRT(ARG)=DSQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CALL STBEG TO SET DEFAULT VALUES AND PRINT INITIAL OUTPUT. C CALL STBEG (FUNK) C C FSAVE IS USED TO CHECK THE REPRODUCIBILITY OF FOBJ. C FSAVE=FOBJ C C SET FIXED QUANTITIES. C C MAXSTP = LOG2(MAXIMUM NUMBER OF STEPS) C MAXSTP=3 C C FACUP ... IF MORE THAN FACUP STEPS ARE TAKEN, THE STEP SIZE C IS INCREASED C FACUP=4.0D0 C C ACK = RATIO OF STEP SIZE INCREASE C ACK=2.0D0 C C STCUT = RATIO OF STEP SIZE DECREASE C STCUT=10.0D0 C C MAXOSC = MAXIMUM DEPTH OF SEARCH FOR ZIGZAGGING C MAXOSC=5 C C MINOSC = MINIMUM PERIOD OF ZIGZAGGING SEARCH C MINOSC=2 C RZERO=0.0D0 C C NO REAL OR DOUBLE PRECISION CONSTANTS ARE USED BEYOND THIS C POINT IN THIS SUBROUTINE. C KERFL=0 C C JOCK IS A FLAG USED IN SETTING JVARY. C JOCK=1 C C JUMP IS A FLAG SET BY SUBROUTINE DATSW. C JUMP=2 C C NOSC = CURRENT DEPTH OF ZIGZAGGING INFORMATION C NOSC=0 C C FBEST = BEST VALUE OF FOBJ FOUND SO FAR C FBEST=FOBJ C C DX(J) = CURRENT STEP SIZE FOR X(J) C DO 10 J=1,NV DX(J)=DELTX(J) 10 CONTINUE C IF(KFLAG.LT.0) GO TO 760 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C VARY THE PARAMETERS, ONE AT A TIME. C THIS IS THE STARTING POINT USED EACH TIME THE STEP SIZE IS C REDUCED OR A SUCCESSFUL GIANT STEP IS COMPLETED. C C NCIRC = NUMBER OF CONSECUTIVE X(JX) WITHOUT SIZABLE C CHANGES C 20 NCIRC=0 C C NZIP = NUMBER OF CONSECUTIVE CYCLES WITHOUT A GIANT STEP C NZIP=0 C C MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES.... C THE FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE. C C NACK = NUMBER OF ACTIVE X(JX) CYCLED THROUGH C 30 NACK=0 C DO 540 JX=1,NV C C NOUT = NUMBER OF TRIAL POINTS OUT OF BOUNDS C (USED IN SETTING JFLAT(JX)) C NOUT=0 C C NEQUAL = NUMBER OF TRIAL POINTS WITH FOBJ .EQ. FBEST C (USED IN SETTING JFLAG(JX)) C NEQUAL=0 C C JFLAT(JX) WILL BE NONZERO IF CHANGING X(JX) DID NOT CHANGE C FOBJ. C JFLAT(JX)=0 C C VEC(J) = CURRENT VECTOR OF NUMBER OF STEPS IN X(J) C VEC(JX)=RZERO C C DLX(JX) = CHANGE IN X(JX) C DLX(JX)=RZERO C IF(MASK(JX).EQ.0) GO TO 40 JFLAT(JX)=1 GO TO 530 C 40 NACK=NACK+1 ABSDX=DX(JX) IF(ABSDX.LT.RZERO) ABSDX=-ABSDX C C CHECK THAT DX(JX) IS NOT NEGLIGIBLE. C XSAVE=X(JX) XPLUS=XSAVE+DX(JX) IF(XPLUS.EQ.XSAVE) GO TO 50 XPLUS=XSAVE-DX(JX) IF(XPLUS.NE.XSAVE) GO TO 60 C C DX(JX) IS NEGLIGIBLE COMPARED TO X(JX) , SO THERE IS NO C REASON TO STEP X(JX) . C 50 JFLAT(JX)=2 GO TO 140 C C STEP X(JX). C 60 X(JX)=XSAVE+DX(JX) JVARY=0 IF(JOCK.LE.0) GO TO 70 JOCK=0 JVARY=JX C 70 IF(X(JX).GE.XMIN(JX) .AND. X(JX).LE.XMAX(JX)) GO TO 80 NOUT=1 GO TO 90 C 80 CALL FUNK NF=NF+1 JVARY=JX FPREV=FOBJ IF(FOBJ.LT.FBEST) GO TO 170 IF(FOBJ.EQ.FBEST) NEQUAL=1 C C STEP X(JX) THE OTHER WAY. C 90 XPLUS=X(JX) X(JX)=XSAVE-DX(JX) IF(X(JX).GE.XMIN(JX) .AND. X(JX).LE.XMAX(JX)) * GO TO 100 NOUT=NOUT+1 GO TO 110 C 100 CALL FUNK NF=NF+1 JVARY=JX IF(FOBJ.LT.FBEST) GO TO 160 IF(FOBJ.EQ.FBEST) NEQUAL=NEQUAL+1 C 110 IF(NEQUAL.EQ.2 .OR. (NOUT.EQ.1 .AND. NEQUAL.EQ.1 .AND. * (XSAVE.EQ.XMIN(JX) .OR. XSAVE.EQ.XMAX(JX)))) * GO TO 130 IF(NOUT.GT.0) GO TO 140 C C PERFORM QUADRATIC INTERPOLATION. C DENOM=(FOBJ-FBEST)-(FBEST-FPREV) IF(DENOM.LE.RZERO) GO TO 140 DLX(JX)=-DX(JX)*(FOBJ-FPREV)/(DENOM+DENOM) VEC(JX)=DLX(JX)/ABSDX X(JX)=XSAVE+DLX(JX) IF(X(JX).EQ.XSAVE) GO TO 120 CALL FUNK NF=NF+1 IF(FOBJ.GE.FBEST) GO TO 120 FBEST=FOBJ JOCK=1 GO TO 150 C 120 DLX(JX)=RZERO VEC(JX)=RZERO GO TO 140 C C BOTH TRIAL POINTS HAD FOBJ .EQ. FBEST , OR ELSE C ONE TRIAL POINT DID AND THE BASE POINT WAS ON A CONSTRAINT. C 130 JFLAT(JX)=1 C C WE WERE UNABLE TO IMPROVE FOBJ BY VARYING THIS X(JX) . C RETREAT TO THE BASE POINT. C 140 X(JX)=XSAVE C 150 NCIRC=NCIRC+1 IF(NCIRC.GE.NACTIV) GO TO 570 GO TO 210 C C FLIP DX(JX) FOR MORE EFFICIENT OPERATION NEXT TIME. C 160 DX(JX)=-DX(JX) C C A LOWER VALUE OF FOBJ HAS BEEN FOUND. C TAKE A STEP, INCREASE THE STEP SIZE, AND REPEAT AS LONG AS C FOBJ DECREASES, UP TO MAXSTP TIMES. C 170 NCIRC=0 NSTEPS=0 DEL=DX(JX) C 180 FPREV=FBEST FBEST=FOBJ VEC(JX)=VEC(JX)+DEL/ABSDX DLX(JX)=DLX(JX)+DEL NSTEPS=NSTEPS+1 IF(NSTEPS.GE.MAXSTP) GO TO 190 DEL=ACK*DEL XPLUS=XSAVE XSAVE=X(JX) X(JX)=XSAVE+DEL IF(X(JX).LT.XMIN(JX) .OR. X(JX).GT.XMAX(JX)) GO TO 200 CALL FUNK NF=NF+1 IF(FOBJ.LT.FBEST) GO TO 180 C C PERFORM QUADRATIC INTERPOLATION. C DXZ=XSAVE-XPLUS DXU=X(JX)-XSAVE DFZ=FBEST-FPREV DFU=FOBJ-FBEST DENOM=DFZ*DXU-DFU*DXZ IF(DENOM.EQ.RZERO) GO TO 200 DEL=(DFZ*DXU**2+DFU*DXZ**2)/(DENOM+DENOM) X(JX)=XSAVE+DEL IF(X(JX).EQ.XSAVE) GO TO 210 CALL FUNK NF=NF+1 IF(FOBJ.GE.FBEST) GO TO 200 FBEST=FOBJ DLX(JX)=DLX(JX)+DEL VEC(JX)=VEC(JX)+DEL/ABSDX C 190 JOCK=1 GO TO 210 C C RETREAT TO THE BEST KNOWN POINT. C 200 X(JX)=XSAVE C C CHECK WHETHER THE STEP SIZE SHOULD BE INCREASED. C 210 IF((NZIP.LE.0 .AND. NACK.LE.1) .OR. * VEC(JX).EQ.RZERO) GO TO 530 ABSVEC=VEC(JX) IF(ABSVEC.LT.RZERO) ABSVEC=-ABSVEC IF(ABSVEC.LT.FACUP) GO TO 250 C C INCREASE THE STEP SIZE. C DX(JX)=DX(JX)*ACK C C RESCALE THE NUMBERS OF STEPS STORED IN VEC(JX) AND C ERR(JX,*) . C VEC(JX)=VEC(JX)/ACK IF(NOSC.LE.0) GO TO 230 C DO 220 J=1,NOSC ERR(JX,J)=ERR(JX,J)/ACK 220 CONTINUE C 230 IF(NTRACE.GE.1) WRITE(KW,240)JX,DX(JX) 240 FORMAT(/' STEP SIZE',I3,' INCREASED TO ',1PG12.4) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C STEP ALONG A RESULTANT DIRECTION, IF POSSIBLE. C 250 IF(NZIP.LE.0) GO TO 530 NONZER=0 SUMV=RZERO C DO 260 J=1,NV IF(VEC(J).NE.RZERO) NONZER=NONZER+1 SUMV=SUMV+VEC(J)**2 260 CONTINUE C IF(NONZER.LT.2) GO TO 530 IF(SUMV.LE.RZERO) GO TO 380 C C GIANT STEPS WILL BE ATTEMPTED. C FIRST, CHECK FOR POSSIBLE GIGANTIC STEPS. C IF(MAXOSC.LT.1) GO TO 380 C C ZIGZAGGING SEARCH SECTION...... C C KL = POINTER FOR ZIGZAGGING CHECK. C KL=1 C C STORE ZIGZAGGING INFORMATION. C NOSC=NOSC+1 IF(NOSC.LE.MAXOSC) GO TO 290 NOSC=MAXOSC IF(NOSC.EQ.1) GO TO 290 C C THE QUEUE OF ZIGZAGGING INFORMATION IS FULL. C PUSH IT DOWN, THROWING AWAY THE OLDEST ITEM. C DO 280 K=2,NOSC FOSC(K-1)=FOSC(K) C DO 270 J=1,NV XOSC(J,K-1)=XOSC(J,K) ERR(J,K-1)=ERR(J,K) 270 CONTINUE C 280 CONTINUE C C ADD THE NEW ITEM TO THE QUEUE. C 290 SUMV=ZSQRT(SUMV) C DO 300 J=1,NV XOSC(J,NOSC)=X(J) ERR(J,NOSC)=VEC(J)/SUMV 300 CONTINUE C FOSC(NOSC)=FBEST IF(NOSC.LT.2) GO TO 380 C C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION C MORE NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP C THAN WAS THE IMMEDIATELY PREVIOUS STEP. C THIS MAY MEAN THAT THE DIRECTIONS OF THE GIANT STEPS ZIGZAG. C IF SO, TRY GIGANTIC (ZIGZAG) STEPS OF DECREASING PERIOD, C THEN TRY ORDINARY GIANT STEPS. C SINCE THE DIRECTIONS ARE GIVEN AS NUMBERS OF STEPS, THIS C PROCEDURE IS SCALE-INDEPENDENT. C COXCOM=RZERO C DO 310 J=1,NV COXCOM=COXCOM+ERR(J,NOSC)*ERR(J,NOSC-1) 310 CONTINUE C NAH=NOSC-MINOSC C 320 IF(KL.GT.NAH) GO TO 380 C DO 340 K=KL,NAH C C NRETRY = NUMBER OF ZIGZAGGING PERIODS YET TO BE TESTED. C NRETRY=NAH-K COSIN=RZERO C DO 330 J=1,NV COSIN=COSIN+ERR(J,NOSC)*ERR(J,K) 330 CONTINUE C IF(K.GE.NOSC-1 .OR. * (COSIN.GT.RZERO .AND. COSIN.GT.COXCOM)) * GO TO 350 340 CONTINUE C GO TO 380 C C ZIGZAGGING HAS BEEN DETECTED. C ATTEMPT TO TAKE GIGANTIC STEPS. C 350 KL=K+1 NT=NOSC-K IF(NTRACE.GE.1) WRITE(KW,360)NT,COXCOM,COSIN 360 FORMAT(/' ******** GIGANTIC STEP WITH PERIOD ',I2, * ' BEING ATTEMPTED.'/6X,'COXCOM, COSIN = ',1PG12.4, * G12.4) C DO 370 J=1,NV C C SALVO SAVES DLX DURING GIGANTIC STEPS. C SALVO(J)=DLX(J) DLX(J)=X(J)-XOSC(J,K) 370 CONTINUE C FPREV=FOSC(K) GO TO 390 C C SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE. C 380 FPREV=FSTORE(JX) C C NRETRY = -1 IF A GIANT STEP IS BEING ATTEMPTED C NRETRY=-1 C C NGIANT = NUMBER OF GIANT OR GIGANTIC STEPS COMPLETED C 390 NGIANT=0 NFSAV=NF C 400 DO 410 J=1,NV XS(J)=X(J) IF(MASK(J).NE.0) GO TO 410 X(J)=X(J)+DLX(J) IF(X(J).GT.XMAX(J)) X(J)=XMAX(J) IF(X(J).LT.XMIN(J)) X(J)=XMIN(J) 410 CONTINUE C JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(FOBJ.GE.FBEST) GO TO 470 FPREV=FBEST FBEST=FOBJ C DO 420 J=1,NV DLX(J)=DLX(J)*ACK 420 CONTINUE C NGIANT=NGIANT+1 IF(NTRACE.LT.1) GO TO 460 IF(NGIANT.GT.1) GO TO 450 WRITE(KW,430)(VEC(J),J=1,JX) 430 FORMAT(/' NO. OF STEPS = ',1PG10.2,4G10.2/ * (16X,5G10.2)) WRITE(KW,440)FPREV,NFSAV,(XS(J),J=1,NV) 440 FORMAT(/' FOBJ =',1PG15.7,7X,' NF =',I7/ * ' X =',5G15.7/(4X,5G15.7)) 450 WRITE(KW,440)FOBJ,NF,(X(J),J=1,NV) C 460 CONTINUE GO TO 400 C 470 IF(NGIANT.LE.0) GO TO 500 C C PERFORM QUADRATIC INTERPOLATION. C DENOM=ACK*(FPREV-FBEST)-(FBEST-FOBJ) IF(DENOM.LE.RZERO) GO TO 500 CINDER=((FPREV-FBEST)*ACK+(FBEST-FOBJ)/ACK)/ * (DENOM+DENOM) C DO 480 J=1,NV IF(MASK(J).NE.0) GO TO 480 X(J)=XS(J)+CINDER*DLX(J) IF(X(J).GT.XMAX(J)) X(J)=XMAX(J) IF(X(J).LT.XMIN(J)) X(J)=XMIN(J) 480 CONTINUE C JOCK=0 JVARY=0 CALL FUNK NF=NF+1 IF(FOBJ.GE.FBEST) GO TO 500 C FBEST=FOBJ JOCK=1 STEPS=NGIANT STEPS=STEPS+CINDER IF(NTRACE.GE.1) WRITE(KW,490)FBEST,STEPS,(X(J),J=1,NV) 490 FORMAT(/' FOBJ =',1PG15.7,' AFTER',G10.2, * ' GIANT STEPS.'//' X =',5G15.7/(4X,5G15.7)) GO TO 640 C 500 DO 510 J=1,NV IF(NRETRY.GE.0) DLX(J)=SALVO(J) X(J)=XS(J) 510 CONTINUE C IF(NTRACE.GE.1) WRITE(KW,520)FBEST,NGIANT, * (X(J),J=1,NV) 520 FORMAT(/' FOBJ =',1PG15.7,' AFTER',I3, * ' GIANT STEP(S).'/' X =',5G15.7/(4X,5G15.7)) C IF(NGIANT.GT.0) GO TO 640 C IF(NRETRY.GT.0) GO TO 320 C C IF ALL GIGANTIC STEPS WERE UNSUCCESSFUL, TRY A GIANT STEP. C IF(NRETRY.EQ.0) GO TO 380 C C AN UNSUCCESSFUL GIANT STEP HAS OCCURRED. C DELETE ITS ZIGZAGGING INFORMATION. C NOSC=NOSC-1 IF(NOSC.LT.0) NOSC=0 C C COMPLETE THE MAIN DO LOOP. C C FSTORE(JX) SAVES FBEST FOR INTERPOLATION IN GIANT STEPS. C 530 FSTORE(JX)=FBEST C C RETURN IF THE SENSE SWITCH IS ON. C CALL DATSW (NSSW,JUMP) IF(JUMP.LE.1) GO TO 710 C IF(NF.GT.NFMAX) GO TO 690 C 540 CONTINUE C C THIS IS THE END OF THE MAIN DO LOOP. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED. C PRINT ANOTHER LINE OF TRACES. C IF(NTRACE.GE.1) WRITE(KW,430)(VEC(J),J=1,NV) IF(NZIP.NE.0 .OR. NTRACE.LT.1) GO TO 560 WRITE(KW,440)FBEST,NF,(X(J),J=1,NV) WRITE(KW,550) 550 FORMAT(' ') C 560 NZIP=NZIP+1 GO TO 30 C C ESTABLISH THE CURRENT BEST KNOWN POINT AS THE BASE POINT. C 570 FSTORE(JX)=FBEST C C PRINT THE REMAINING TRACES. C IF(NTRACE.LT.1) GO TO 580 WRITE(KW,430)(VEC(J),J=1,JX) WRITE(KW,440)FBEST,NF,(X(J),J=1,NV) C C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES, AND C CHECK WHETHER OR NOT ALL ABS(DX(J)) .LE. DELMIN(J) . C 580 NGATE=1 C DO 600 J=1,NV IF(MASK(J).NE.0) GO TO 590 ABSDX=DX(J) IF(ABSDX.LT.RZERO) ABSDX=-ABSDX IF(ABSDX.GT.DELMIN(J)) NGATE=0 590 DX(J)=DX(J)/STCUT 600 CONTINUE C IF(NGATE.EQ.1) GO TO 650 C C CHECK THE JFLAT(J) . C IF(NFLAT.LE.0) GO TO 620 JFLMIN=5 C DO 610 J=1,NV IF(MASK(J).EQ.0 .AND. JFLAT(J).LT.JFLMIN) * JFLMIN=JFLAT(J) 610 CONTINUE C IF(JFLMIN.GE.1) GO TO 670 C 620 IF(NTRACE.GE.1) WRITE(KW,630)(DX(J),J=1,NV) 630 FORMAT(/36(' *')//' STEP SIZES REDUCED TO....'/ * (1X,1PG12.4,4G12.4)) C C RETURN IF THE SENSE SWITCH IS ON. C CALL DATSW (NSSW,JUMP) IF(JUMP.LE.1) GO TO 710 C IF(NF.GT.NFMAX) GO TO 690 C C SEARCH SOME MORE. C 640 CONTINUE GO TO 20 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C THE MINIMIZATION IS FINISHED. C 650 KFLAG=1 IF(NTRACE.GE.0) WRITE(KW,660) 660 FORMAT(//' TERMINATED WHEN THE STEP SIZES', * ' BECAME AS SMALL AS THE DELMIN(J).') GO TO 730 C 670 KFLAG=2 IF(NTRACE.GE.0) WRITE(KW,680) 680 FORMAT(//' TERMINATED WHEN THE FUNCTION VALUES', * ' AT ALL TRIAL POINTS WERE IDENTICAL.') GO TO 730 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 690 KFLAG=-2 IF(NTRACE.GE.-1) WRITE(KW,700)NFMAX 700 FORMAT(//' ABNORMAL TERMINATION....', * ' MORE THAN NFMAX =',I8, * ' CALLS TO THE FOBJ SUBROUTINE.') GO TO 740 C 710 KFLAG=-3 IF(NTRACE.GE.-1) WRITE(KW,720)NSSW 720 FORMAT(//' ABNORMAL TERMINATION.... TERMINATED BY', * ' OPERATOR VIA SENSE SWITCH ',I2) GO TO 740 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 730 IF(NTRACE.LT.0) GO TO 760 740 IF(NTRACE.LT.-1) GO TO 760 WRITE(KW,750)(DX(J),J=1,NV) 750 FORMAT(//' CURRENT STEP SIZES....'//(1X,1PG12.4,4G12.4)) C C CALL FUNK WITH THE BEST SET OF X(J) . C 760 JVARY=0 CALL FUNK NF=NF+1 IF(FBEST.LE.FSAVE .AND. FOBJ.EQ.FBEST) GO TO 780 NOREP=NOREP+2 IF(NTRACE.GE.-1) WRITE(KW,770)NF,FSAVE,FBEST,FOBJ 770 FORMAT(//' WARNING.... FOBJ IS NOT A REPRODUCIBLE', * ' FUNCTION OF X(J).',7X,' NF = ',I5//1X,1PG23.15, * 2G23.15) C 780 IF(NTRACE.GE.0) WRITE(KW,790)NF,FOBJ,(X(J),J=1,NV) 790 FORMAT(//1X,I6,' FUNCTION COMPUTATIONS'// * ' FINAL VALUE OF FOBJ =',1PG23.15// * 9X,' FINAL VALUES OF X(J)....'//(1X,5G15.7)) IF(KFLAG.LT.0) GO TO 800 C IF(MATRX.LT.70 .OR. MATRX.GT.130) GO TO 800 C C CALL STERR TO COMPUTE AN APPROXIMATE ERROR MATRIX. C CALL STERR (FUNK) C C THIS IS THE ONLY RETURN STATEMENT IN THIS SUBROUTINE.... C 800 RETURN C C END STEPIT C END SUBROUTINE STBEG (FUNK) C C STBEG 1.5 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER C COMPUTER SCIENCE DEPT., OKLAHOMA STATE UNIVERSITY C C STBEG SETS DEFAULT VALUES AND PRINTS INITIAL OUTPUT FOR C STEPIT. THE CALL TO STBEG IS THE FIRST EXECUTABLE STATEMENT C IN STEPIT, TO FACILITATE OVERLAYING IF NECESSARY. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,X,XMAX,XMIN,DELTX,DELMIN,NV, C NTRACE,MATRX,MASK,NFMAX,NFLAT,KW C C OUTPUT QUANTITIES.... NSSW,NACTIV,NF,KFLAG,NOREP, C AND SOMETIMES X,XMAX,XMIN,DELTX, C DELMIN C DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,DX,XS, * DLX DOUBLE PRECISION HUGE,DELDF,RZERO,UNITR,RTEN,RELACC, * XPLUS,FSAVE C INTEGER J,JUMP,JVARY,KERFL,KFLAG,KTYPE,KW,MASK,MATRX, * NACTIV,NF,NFLAT,NFMAX,NOREP,NSSW,NTRACE,NV,NVMAX, * NXTRA 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 INTERNAL STEPIT COMMON..... COMMON /STORK/ DX(20),XS(20),DLX(20),NACTIV,NSSW,NF C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C SET FIXED QUANTITIES .... C C KTYPE = LOGICAL UNIT NUMBER OF THE CONSOLE TYPEWRITER, IF C ANY (IRRELEVANT IF A DUMMY DATSW IS USED) C KTYPE=1 C C NSSW = SENSE SWITCH NUMBER FOR TERMINATION BY OPERATOR C (IRRELEVANT IF A DUMMY DATSW IS USED) C NSSW=6 C C HUGE = A LARGE REAL NUMBER, THE DEFAULT VALUE FOR XMAX AND C FOR (-XMIN) C HUGE=1.0D35 C C NVMAX = MAXIMUM VALUE OF NV C NVMAX=20 C C DELDF = DEFAULT VALUE FOR DELTX(J) C DELDF=0.01D0 C RZERO=0.0D0 UNITR=1.0D0 RTEN=10.0D0 C C NO REAL OR DOUBLE PRECISION CONSTANTS ARE USED BEYOND THIS C POINT IN THIS SUBROUTINE. C KFLAG=0 NOREP=0 KERFL=0 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CHECK SOME INPUT QUANTITIES, AND SET THEM TO DEFAULT VALUES C IF DESIRED. C IF(NV.GE.1 .AND. NV.LE.NVMAX) GO TO 20 C WRITE(KW,10)NV,NVMAX 10 FORMAT(///' ****** FATAL ERROR IN SUBROUTINE STEPIT...' * //6X,'NV =',I14,6X,'NVMAX =',I14) STOP C C MAKE SURE THAT THE SENSE SWITCH IS OFF. C 20 JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.GE.2) GO TO 50 C C THIS IS THE ONLY USAGE OF THE CONSOLE TYPEWRITER. C WRITE(KTYPE,30)NSSW 30 FORMAT(/' TURN OFF SENSE SWITCH ',I2//' ') 40 CALL DATSW (NSSW,JUMP) IF(JUMP.LE.1) GO TO 40 C C COMPUTE RELACC, THE RELATIVE PRECISION OF THE MACHINE AND C ARITHMETIC BEING USED. C RELACC IS USED IN SETTING DELMIN(J) TO A DEFAULT VALUE. C 50 RELACC=UNITR 60 RELACC=RELACC/RTEN XPLUS=UNITR+RELACC IF(XPLUS.GT.UNITR) GO TO 60 C C NACTIV = NUMBER OF ACTIVE X(J) C NACTIV=0 C DO 110 J=1,NV IF(MASK(J).NE.0) GO TO 110 NACTIV=NACTIV+1 C C CHECK THAT DELTX(J) IS NOT NEGLIGIBLE. C IF(DELTX(J).EQ.RZERO) GO TO 70 XPLUS=X(J)+DELTX(J) IF(XPLUS.EQ.X(J)) GO TO 70 XPLUS=X(J)-DELTX(J) IF(XPLUS.NE.X(J)) GO TO 90 C C DELTX(J) IS NEGLIGIBLE COMPARED TO X(J) . RESET DELTX(J) . C 70 IF(X(J).EQ.RZERO) GO TO 80 DELTX(J)=DELDF*X(J) GO TO 90 C 80 DELTX(J)=DELDF C 90 IF(DELMIN(J).EQ.RZERO) DELMIN(J)=DELTX(J)*RELACC IF(DELMIN(J).LT.RZERO) DELMIN(J)=-DELMIN(J) C IF(XMAX(J).GT.XMIN(J)) GO TO 100 XMAX(J)=HUGE XMIN(J)=-HUGE C 100 IF(X(J).GT.XMAX(J)) X(J)=XMAX(J) IF(X(J).LT.XMIN(J)) X(J)=XMIN(J) 110 CONTINUE C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IF(NTRACE.LT.0) GO TO 200 WRITE(KW,120) 120 FORMAT('1SUBROUTINE STEPIT'//' INITIAL VALUES....'/' ') WRITE(KW,130)(MASK(J),J=1,NV) 130 FORMAT(/' MASK =',I7,4I13/(2X,5I13)) WRITE(KW,140)(X(J),J=1,NV) 140 FORMAT(/' X =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,150)(XMAX(J),J=1,NV) 150 FORMAT(/' XMAX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,160)(XMIN(J),J=1,NV) 160 FORMAT(/' XMIN =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,170)(DELTX(J),J=1,NV) 170 FORMAT(/' DELTX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,180)(DELMIN(J),J=1,NV) 180 FORMAT(/' DELMIN=',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,190)NV,NACTIV,MATRX,NFMAX,NFLAT,RELACC 190 FORMAT(//1X,I3,' VARIABLES,',I3,' ACTIVE.',8X,' MATRX =', * I4//' NFMAX =',I8,8X,' NFLAT =',I2,9X,' RELACC =', * G11.4) C 200 JVARY=0 CALL FUNK FSAVE=FOBJ CALL FUNK C IF(NTRACE.GE.0) WRITE(KW,210)FOBJ 210 FORMAT(//' FOBJ =',1PG18.10) C C NF = NUMBER OF CALLS TO SUBROUTINE FUNK SO FAR C NF=2 C IF(FOBJ.EQ.FSAVE) GO TO 230 NOREP=1 IF(NTRACE.GE.-1) WRITE(KW,220)NF,FSAVE,FOBJ 220 FORMAT(//' WARNING.... FOBJ IS NOT A REPRODUCIBLE', * ' FUNCTION OF X(J).',7X,' NF =',I6//5X,1PG23.16, * G23.16) C 230 IF(NACTIV.GT.0) GO TO 250 KFLAG=-1 IF(NTRACE.GE.-1) WRITE(KW,240) 240 FORMAT(///' ****** WARNING... MASK(J).NE.0 FOR ALL J ,' * ,' IN A CALL TO SUBROUTINE STEPIT.'// * 6X,'FOBJ WILL BE EVALUATED BUT NOT MINIMIZED.'//' ') GO TO 270 C 250 IF(NTRACE.GE.0) WRITE(KW,260) 260 FORMAT(///' BEGIN MINIMIZATION'//' ') C 270 RETURN C C END STBEG C END SUBROUTINE STERR (FUNK) C C STERR 1.6 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER C DEPARTMENT OF COMPUTER SCIENCE, C OKLAHOMA STATE UNIVERSITY C C STERR IS CALLED BY STEPIT TO COMPUTE AN APPROXIMATE ERROR C MATRIX FOR A NONLINEAR FITTING PROBLEM. C THE VALUES COMPUTED ARE SOMETIMES POOR APPROXIMATIONS. C FOR EACH CLASS OF PROBLEMS, THE ERRORS SHOULD BE CHECKED C USING SUBROUTINE FIDO. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,KW,NSSW,DX,NF,X,NTRACE,NV C C OUTPUT QUANTITIES.... NF,ERR,KERFL, AND SOMETIMES DX C C SCRATCH STORAGE...... XS,DLX C C THE DX(J) ARE THE STEP SIZES USED IN APPROXIMATING THE C SECOND PARTIAL DERIVATIVES OF FOBJ WITH RESPECT TO THE X(J) C BY FINITE DIFFERENCES. C ERR RETURNS THE ERROR MATRIX. C XMAX, XMIN, AND MASK ARE IGNORED IN STERR. C THE REAL FORMAT SPECIFICATIONS USED ARE G13.5, G16.8, AND C G23.16 . C DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,DX,XS, * DLX,SECND,FBEST,RZERO,UNITR,RTWO,ABSERR,DENOM,RTEN, * P,Q,ARG,ZSQRT,DSQRT,ABSX,DXFAC C INTEGER J,JACTIV,JJ,JMU,JUMP,JVARY,K,KACTIV,KERFL,KFLAG, * KK,KW,L,LL,M,MASK,MATRX,NACTIV,NEGDEF,NF,NFLAT, * NFMAX,NOREP,NSSW,NTRACE,NV,NXTRA C DIMENSION SECND(2,2) 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 INTERNAL STEPIT COMMON..... COMMON /STORK/ DX(20),XS(20),DLX(20),NACTIV,NSSW,NF C C XS AND DLX ARE IN COMMON ONLY TO CONSERVE STORAGE. C C ZSQRT(ARG)=SQRT(ARG) ZSQRT(ARG)=DSQRT(ARG) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C RZERO=0.0D0 UNITR=1.0D0 RTWO=2.0D0 RTEN=10.0D0 C C NO REAL OR DOUBLE PRECISION CONSTANTS ARE USED BEYOND THIS C POINT, IN THIS SUBROUTINE. C KERFL=0 C C SET THE STEP SIZES. C DXFAC=RTEN**(100-MATRX) C DO 10 K=1,NV XS(K)=X(K) IF(MASK(K).NE.0) GO TO 10 ABSX=X(K) IF(ABSX.LT.RZERO) ABSX=-ABSX DX(K)=DXFAC*ABSX IF((ABSX+DX(K))-ABSX.GT.RZERO) DX(K)=(ABSX+DX(K))-ABSX IF(DX(K).LE.RZERO .OR. (ABSX+DX(K))-ABSX.LE.RZERO) * DX(K)=DXFAC 10 CONTINUE C CALL FUNK NF=NF+1 FBEST=FOBJ C IF(NTRACE.LT.0) GO TO 40 WRITE(KW,20) 20 FORMAT('1SUBROUTINE STERR.'//' COMPUTE AN APPROXIMATE', * ' ERROR MATRIX USING FINITE DIFFERENCES.'/// * ' INCREMENTS IN X(J) TO BE USED....') WRITE(KW,30)(DX(K),K=1,NV) 30 FORMAT(/(1X,1PG13.5,4G13.5)) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C APPROXIMATE THE (SYMMETRIC) MATRIX OF SECOND PARTIAL C DERIVATIVES OF FOBJ WITH RESPECT TO THE X(J), USING DIVIDED C DIFFERENCES. C COMPUTE THE DIAGONAL PARTIALS FIRST. C 40 DO 60 J=1,NV ERR(J,J)=RZERO IF(MASK(J).NE.0) GO TO 60 JVARY=0 C DO 50 K=1,2 X(J)=XS(J)+DX(J) CALL FUNK NF=NF+1 JVARY=J SECND(K,1)=FOBJ DX(J)=-DX(J) 50 CONTINUE C X(J)=XS(J) ERR(J,J)=((SECND(1,1)-FBEST)-(FBEST-SECND(2,1)))/ * DX(J)**2 60 CONTINUE C C COMPUTE THE OFF-DIAGONAL PARTIALS. C IF(NV.LT.2) GO TO 110 C DO 100 J=2,NV JMU=J-1 C DO 90 K=1,JMU ERR(J,K)=RZERO IF(MASK(J).NE.0 .OR. MASK(K).NE.0) GO TO 90 C DO 80 L=1,2 X(J)=XS(J)+DX(J) JVARY=0 C DO 70 M=1,2 X(K)=XS(K)+DX(K) CALL FUNK NF=NF+1 JVARY=K SECND(L,M)=FOBJ X(K)=XS(K) DX(K)=-DX(K) 70 CONTINUE C X(J)=XS(J) DX(J)=-DX(J) C C RETURN IF THE SENSE SWITCH IS ON. C JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.EQ.2) GO TO 80 KERFL=-1 GO TO 470 80 CONTINUE C ERR(J,K)=((SECND(1,1)-SECND(1,2))- * (SECND(2,1)-SECND(2,2)))/ * (RTWO*DX(J)*RTWO*DX(K)) 90 CONTINUE C 100 CONTINUE C C THIS IS THE END OF THE DERIVATIVE COMPUTATION. C 110 IF(NTRACE.LT.0) GO TO 180 WRITE(KW,120) 120 FORMAT(///' MATRIX OF THE SECOND PARTIAL DERIVATIVES...') WRITE(KW,130)(K,K=1,NV) 130 FORMAT(//11X,'K =',I3,4I12/(6X,5I12)) WRITE(KW,140)(MASK(K),K=1,NV) 140 FORMAT(/6X,'MASK(K) =',I3,4I12/(6X,5I12)) WRITE(KW,150) 150 FORMAT(/' J MASK(J)') C DO 170 J=1,NV WRITE(KW,160)J,MASK(J),(ERR(J,K),K=1,J) 160 FORMAT(/1X,I3,I6,2X,1PG12.4,4G12.4/(12X,5G12.4)) 170 CONTINUE C C PACK THE MATRIX OF SECOND DERIVATIVES. C 180 NACTIV=0 C DO 200 J=1,NV IF(MASK(J).NE.0) GO TO 200 NACTIV=NACTIV+1 KACTIV=0 C DO 190 K=1,J IF(MASK(K).NE.0) GO TO 190 KACTIV=KACTIV+1 ERR(NACTIV,KACTIV)=ERR(J,K) IF(ERR(J,K).EQ.RZERO) KERFL=1 190 CONTINUE C 200 CONTINUE C IF(KERFL.GE.1) WRITE(KW,210) 210 FORMAT(//' THE ABOVE MATRIX CONTAINS ONE OR MORE', * ' UNEXPECTED ZEROES.'/' PERHAPS A SMALLER VALUE OF', * ' -MATRX- SHOULD BE TRIED,', * ' TO SEE IF THEY ARE LEGITIMATE.') C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INVERT THE MATRIX OF SECOND PARTIAL DERIVATIVES USING THE C GAUSS-JORDAN METHOD C (F. L. BAUER AND C. REINSCH, P. 45 IN "LINEAR ALGEBRA" C BY J. H. WILKINSON AND C. REINSCH (SPRINGER-VERLAG, 1971)). C ONLY THE LOWER TRIANGLE OF ERR IS USED OR ALTERED. C C NEGDEF = 1 IF THE MATRIX IS NEGATIVE DEFINITE C NEGDEF=0 C DO 310 LL=1,NACTIV L=NACTIV+1-LL P=ERR(1,1) IF(P.LT.RZERO) GO TO 230 IF(P.GT.RZERO) GO TO 240 KERFL=-2 WRITE(KW,220) 220 FORMAT(//' A PIVOT ELEMENT OF THE MATRIX IS ZERO.', * ' PERHAPS -MATRX- SHOULD BE DECREASED.') GO TO 470 C 230 NEGDEF=1 C 240 IF(NACTIV.LT.2) GO TO 290 C DO 280 K=2,NACTIV Q=ERR(K,1) IF(K.LE.L) GO TO 250 XS(K)=Q/P GO TO 260 C 250 XS(K)=-Q/P C 260 DO 270 M=2,K ERR(K-1,M-1)=ERR(K,M)+Q*XS(M) 270 CONTINUE C 280 CONTINUE C 290 ERR(NACTIV,NACTIV)=UNITR/P IF(NACTIV.LT.2) GO TO 310 C DO 300 K=2,NACTIV ERR(NACTIV,K-1)=XS(K) 300 CONTINUE C 310 CONTINUE C IF(NEGDEF.LE.0) GO TO 330 C KERFL=-3 WRITE(KW,320) 320 FORMAT(//' THE ERROR MATRIX IS NEGATIVE DEFINITE.'/ * ' PERHAPS -MATRX- SHOULD BE INCREASED.') C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C UNPACK, DOUBLE, AND SYMMETRIZE THE INVERSE TO FORM THE C ERROR MATRIX. C 330 JACTIV=NACTIV C DO 370 JJ=1,NV J=NV+1-JJ KACTIV=JACTIV C DO 360 KK=1,J K=J+1-KK IF(MASK(J).EQ.0 .AND. MASK(K).EQ.0) GO TO 340 ERR(J,K)=RZERO GO TO 350 C 340 ERR(J,K)=ERR(JACTIV,KACTIV)*RTWO KACTIV=KACTIV-1 C 350 ERR(K,J)=ERR(J,K) 360 CONTINUE C IF(MASK(J).EQ.0) JACTIV=JACTIV-1 370 CONTINUE C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C PRINT THE STANDARD ERRORS AND THE CORRELATIONS, AND RETURN. C IF(NTRACE.GE.0) WRITE(KW,380) 380 FORMAT(///' APPROXIMATE STANDARD ERRORS....'/// * 3X,'J',6X,'MASK(J)',9X,'X(J)',14X,'ERROR') C DO 420 J=1,NV ABSERR=ERR(J,J) IF(ABSERR.EQ.RZERO) GO TO 400 IF(ABSERR.LT.RZERO) ABSERR=-ABSERR ABSERR=ZSQRT(ABSERR) IF(MASK(J).NE.0) GO TO 400 IF(ERR(J,J).GT.RZERO) GO TO 400 IF(ERR(J,J).LT.RZERO) ABSERR=-ABSERR KERFL=-4 WRITE(KW,390)ERR(J,J) 390 FORMAT(//' NEGATIVE OR ZERO MEAN SQUARE', * ' ERROR ENCOUNTERED....',3X,G16.8/ * ' PERHAPS -MATRX- SHOULD BE INCREASED.') C 400 IF(NTRACE.GE.0) WRITE(KW,410)J,MASK(J),X(J),ABSERR 410 FORMAT(/1X,I3,I10,6X,1PG16.8,4X,G13.5) XS(J)=ABSERR 420 CONTINUE C C COMPUTE AND PRINT THE CORRELATIONS. C IF(NTRACE.LT.0 .OR. NV.LT.2) GO TO 470 WRITE(KW,430) 430 FORMAT(///' LOWER TRIANGLE OF THE CORRELATION MATRIX....' * ) WRITE(KW,130)(K,K=1,NV) WRITE(KW,140)(MASK(K),K=1,NV) WRITE(KW,150) C DO 460 J=1,NV C DO 450 K=1,J DENOM=XS(J)*XS(K) IF(DENOM.NE.RZERO) GO TO 440 DLX(K)=RZERO GO TO 450 C 440 IF(DENOM.LT.RZERO) DENOM=-DENOM DLX(K)=ERR(J,K)/DENOM 450 CONTINUE C WRITE(KW,160)J,MASK(J),(DLX(K),K=1,J) 460 CONTINUE C C CALL FUNK AGAIN, SO THAT THE LAST CALL TO FUNK WILL HAVE C BEEN MADE WITH THE OPTIMAL VECTOR X(*). C 470 JVARY=0 CALL FUNK NF=NF+1 RETURN C C END STERR C END SUBROUTINE DATSW (NSSW,JUMP) C C THIS IS A DUMMY VERSION OF SUBROUTINE DATSW, WITH ALL C SWITCHES PERMANENTLY OFF. C INTEGER NSSW,JUMP C JUMP=2 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 DOUBLE PRECISION XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,FLAMBD, * FNU,RELDIF,RELMIN,RZERO,HUGE, X C INTEGER JVARY,JX,KALCP,KERFL,KFLAG,KORDIF,KW,LEQU,MASK, * MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP, * NTRACE,NV,NVMAX,NXTRA 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 $ENTRY 2 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 //