//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 SIMPLEX.CNTL CURRENT VERSION OF SIMPLEX C SUBROUTINE STINT (FUNK) C C COMMON INTERFACE ROUTINE TO FACILITATE SWAPPING OF STEPIT, C SIMPLEX, MARQ, ETC. C EXTERNAL FUNK CALL SMPLX (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 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="SIMPTEST.DAT") C OPEN (UNIT=6,FILE="SIMPOUT") 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 NTRACE=1 NFMAX=550 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 IF(KR.EQ.5) GO TO 90 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 90 CONTINUE C CLOSE (UNIT=6) 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 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.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 SMPLX (FUNK) C C SIMPLEX 2.12 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER C (PRESENT ADDRESS ... C COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY, C STILLWATER, OKLAHOMA 74078 C (405)-744-5676 ) C C SIMPLEX FINDS LOCAL MINIMA OF A SMOOTH FUNCTION OF SEVERAL C PARAMETERS. IT WILL ALSO HANDLE SOME NON-SMOOTH FUNCTIONS. C THE METHOD IS OFTEN VERY SLOW IF THE NUMBER OF PARAMETERS IS C AT ALL LARGE (GREATER THAN ABOUT SIX TO EIGHT). C C "A SIMPLEX METHOD FOR FUNCTION MINIMIZATION", C J. A. NELDER AND R. MEAD, C THE COMPUTER JOURNAL 7 (1965) 308-313 C C FOR APPLICATIONS, SEE C C "'DIRECT SEARCH' SOLUTION OF C NUMERICAL AND STATISTICAL PROBLEMS", C ROBERT HOOKE AND T. A. JEEVES, JOURNAL OF THE C ASSOCIATION FOR COMPUTING MACHINERY 8 (1961) 212-229 C C "THE NELDER-MEAD SIMPLEX PROCEDURE FOR FUNCTION C MINIMIZATION", C D. M. OLSSON AND L. S. NELSON, C TECHNOMETRICS 17 (1975) 45-51 C C C SIMPLEX 2.9 IS 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,MASK(*),X(*), C XMAX(*),XMIN(*),DELTX(*), C DELMIN(*),NFMAX,NFLAT,NXTRA,KW C C OUTPUT QUANTITIES.... X(*),FOBJ, KFLAG,NOREP C C FUNK -- THE NAME OF THE SUBROUTINE THAT COMPUTES C FOBJ GIVEN X(1),...,X(NV) C (EACH SUCH SUBROUTINE MUST BE NAMED IN C AN EXTERNAL STATEMENT IN THE CALLING C PROGRAM) C C NV -- THE NUMBER OF PARAMETERS, X(J) C C NTRACE -- =0 FOR NORMAL OUTPUT, C =+1 FOR TRACE OUTPUT, C =+2 FOR DEBUG OUTPUT, C =-1 FOR NO OUTPUT EXCEPT ERROR MESSAGES C C MATRX -- USED BY STEPIT PROPER BUT NOT BY SIMPLEX 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) ON C THE STEP SIZE FOR X(J) C C ERR(J,K) -- SCRATCH STORAGE C C NFMAX -- THE MAXIMUM NUMBER OF FUNCTION C COMPUTATIONS TO BE ALLOWED C C NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN C ALL TRIAL STEPS GIVE IDENTICAL FUNCTION C VALUES. C THE RECOMMENDED VALUE OF NFLAT IS C USUALLY NFLAT=1 . C C JVARY -- USED BY STEPIT AND STP BUT NOT BY SIMPLEX C (SIMPLEX SETS JVARY TO ZERO C PERMANENTLY) C C NXTRA -- NUMBER OF EXTRA POINTS TO BE ADDED TO C THE SIMPLEX C (NXTRA.GT.0 CAUSES A MORE THOROUGH C SEARCH) 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 -- NOT USED BY THIS ROUTINE 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 FORBIDDEN BY OTHERS C (THE MODCOMP II, FOR EXAMPLE). C EXTERNAL FUNK C C SIMPLEX 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 DOUBLE PRECISION FZ,Z,ZBAR,ZSTAR DOUBLE PRECISION HUGE,ALPHA,BETA,GAMMA,RZERO,UNITR,RTWO, * SCALK,XS,FSAVE,DX,DZ,ZMAX,ZMIN,FSTAR,ZBARJ,ZJK,ZJJH, * XJ,ZKJ,ZKJL,XK C INTEGER J,JDIFF,JH,JHSAV,JHTIE,JL,JLOW,JS,JUMP,JVARY, * K,KERFL,KFLAG,KW,LATER,MASK,MATRX,MAXPT,METHD,NF, * NFLAT,NFMAX,NOREP,NOW,NSSW,NTRACE,NV,NVPLUS,NXTRA, * KONTR C DIMENSION Z(20,21),ZBAR(20),ZSTAR(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 INTERNAL SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT,NVPLUS, * NSSW,NF C EQUIVALENCE (Z(1,1),ERR(1,1)) C C SIMPLEX CALLS NO FUNCTIONS, EITHER EXTERNAL OR INTRINSIC. C THE SUBROUTINES CALLED ARE FUNK, SIBEG, SIFUN, AND DATSW. C SIMPLEX 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, THE USER SHOULD SUPPLY C A DUMMY SUBROUTINE FOR DATSW. C THIS SUBROUTINE CONTAINS NO REAL OR DOUBLE PRECISION C CONSTANTS. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C CALL SIBEG (FUNK) C IF(KFLAG.LT.0) GO TO 720 FSAVE=FOBJ KERFL=0 C C SET FIXED QUANTITIES.... C C METHD=1 TO USE THE METHOD OF NELDER AND MEAD, C METHD=2 TO USE A MODIFIED METHOD. C METHD=1 CAN CAUSE THE BEST KNOWN POINT TO BE DISCARDED, C AND HENCE IS NOT RECOMMENDED. C METHD=2 IS RECOMMENDED. C METHD=2 C RZERO=0.0D0 UNITR=1.0D0 RTWO=2.0D0 C JL=NVPLUS LATER=0 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C BEGIN THE NEXT ITERATION. C FIND JH AND JL, THE INDICES OF THE POINTS WITH THE HIGHEST C AND LOWEST FUNCTION VALUES, RESPECTIVELY. C 10 JH=JL DO 20 J=1,NVPLUS IF(FZ(J).GT.FZ(JH)) JH=J IF(FZ(J).LT.FZ(JL)) JL=J 20 CONTINUE C C CHECK FOR POSSIBLE TERMINATION. C IF(KFLAG.NE.0) GO TO 690 IF(JH.EQ.JL .AND. NFLAT.GT.0) GO TO 660 C C CHECK FOR TIES -- MORE THAN ONE POINT IN THE SIMPLEX HAVING C THE HIGHEST FUNCTION VALUE. C JHTIE=0 DO 30 J=1,NVPLUS IF(J.NE.JH .AND. FZ(J).GE.FZ(JH)) JHTIE=7 30 CONTINUE C IF(JH.NE.JL .AND. JHTIE.EQ.0) GO TO 70 C C THERE IS A TIE FOR HIGHEST FUNCTION VALUE. C CHOOSE H FAR FROM L. C DX=RZERO JHSAV=JH C DO 50 J=1,NVPLUS IF(J.EQ.JL .OR. FZ(J).LT.FZ(JH)) GO TO 50 C DO 40 K=1,NV IF(MASK(K).NE.0) GO TO 40 SCALK=X(K) IF(SCALK.EQ.RZERO) SCALK=DELTX(K) DZ=(Z(K,J)-Z(K,JL))/SCALK IF(DZ.LT.RZERO) DZ=-DZ IF(DZ.LE.DX) GO TO 40 JH=J DX=DZ 40 CONTINUE C 50 CONTINUE C IF(NTRACE.GE.2) WRITE(KW,60)JHSAV,JH 60 FORMAT(' TIE BREAKING....',5X,'JHSAV =',I3,5X,'JH =',I3) IF(DX.LE.RZERO) GO TO 640 70 IF(NTRACE.LT.2) GO TO 100 WRITE(KW,80)JL,FZ(JL),(Z(J,JL),J=1,NV) 80 FORMAT(' FZ(JL=',I2,') =',1PG24.16,10X,'Z(J,JL)....'/ * (1X,3G24.16)) WRITE(KW,90)JH,FZ(JH),(Z(J,JH),J=1,NV) 90 FORMAT(' FZ(JH=',I2,') =',1PG24.16,10X,'Z(J,JH)....'/ * (1X,3G24.16)) C C ADD EXTRA POINTS TO THE SIMPLEX, IF DESIRED. C ANY EXTRA POINTS WILL BE SUPERIMPOSED ON THE HIGHEST POINT, C P(JH), SO THAT THEY WILL SPLIT AS SOON AS POSSIBLE. C 100 IF(LATER.NE.0) GO TO 140 LATER=7 IF(NVPLUS+NXTRA.GT.MAXPT) NXTRA=MAXPT-NVPLUS IF(NXTRA.LE.0) GO TO 130 JLOW=NVPLUS+1 NVPLUS=NVPLUS+NXTRA C DO 120 J=JLOW,NVPLUS FZ(J)=FZ(JH) C DO 110 K=1,NV Z(K,J)=Z(K,JH) 110 CONTINUE C 120 CONTINUE C 130 FSAVE=FZ(JL) C C CALCULATE PBAR, THE CENTROID OF ALL POINTS EXCEPT P(JH). C THE MEAN VALUE IS COMPUTED USING A STABLE UPDATE FORMULA BY C D. H. D. WEST, COMMUNICATIONS OF THE A.C.M. C 22 (1979) 532-535. C 140 DO 200 J=1,NV ZBARJ=X(J) IF(MASK(J).NE.0) GO TO 190 ZMAX=Z(J,JL) ZMIN=Z(J,JL) XS=RZERO C DO 150 K=1,NVPLUS IF(K.EQ.JH) GO TO 150 ZJK=Z(J,K) IF(ZJK.GT.ZMAX) ZMAX=ZJK IF(ZJK.LT.ZMIN) ZMIN=ZJK XS=XS+UNITR ZBARJ=ZBARJ+(ZJK-ZBARJ)/XS 150 CONTINUE C C CHECK THE ROUNDING. C ZBARJ MUST LIE IN THE CLOSED INTERVAL (ZMIN,ZMAX). C IF(ZBARJ.LE.ZMAX) GO TO 160 ZBARJ=ZMAX GO TO 170 C 160 IF(ZBARJ.GE.ZMIN) GO TO 190 ZBARJ=ZMIN C 170 NOW=1 IF(NTRACE.GE.1) WRITE(KW,180)J,NOW,NF 180 FORMAT(' ROUNDING OF X(',I3,') CORRECTED AT', * ' CHECKPOINT ',I1,' WITH NF =',I6) C 190 ZBAR(J)=ZBARJ 200 CONTINUE C IF(NTRACE.GE.2) WRITE(KW,210)(ZBAR(J),J=1,NV) 210 FORMAT(' ZBAR(J)....'/(1X,1PG24.16,2G24.16)) C C ATTEMPT A REFLECTION. FORM P* . C THE FORMS USED BELOW FOR REFLECTION, EXPANSION, AND C CONTRACTION ALL HAVE OPTIMAL ROUNDOFF PROPERTIES. C JDIFF=0 C DO 230 J=1,NV XJ=X(J) IF(MASK(J).NE.0) GO TO 220 XJ=ZBAR(J)+ALPHA*(ZBAR(J)-Z(J,JH)) IF(XJ.NE.Z(J,JH)) JDIFF=7 X(J)=XJ 220 ZSTAR(J)=XJ 230 CONTINUE C IF(JDIFF.NE.0) GO TO 250 IF(NTRACE.GE.1) WRITE(KW,240) 240 FORMAT(' ZSTAR = Z(JH). TREAT AS A CONTRACTION', * ' FAILURE.') GO TO 440 C 250 CALL SIFUN (FUNK) FSTAR=FOBJ KONTR=0 IF(FSTAR.LT.FZ(JL)) GO TO 280 C IF(NTRACE.GE.1) WRITE(KW,260)FSTAR,NF 260 FORMAT(' FOBJ =',1PG24.16,5X,'REFLECTION FAILED',8X, * 'NF =',I7) C DO 270 J=1,NVPLUS C C HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER C AND MEAD. THEY USE FSTAR.LE.FZ(J) . C THAT IS NOT CONSISTENT WITH THE TEXT OF THE ARTICLE, AND CAN C CREATE AN INFINITE LOOP OF REFLECTIONS. C IF(J.NE.JH .AND. FSTAR.LT.FZ(J)) GO TO 340 270 CONTINUE C C HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER C AND MEAD. THEY USE FSTAR.LE.FZ(JH) (SEE COMMENT ABOVE). C KONTR=1 IF(FSTAR.LT.FZ(JH)) GO TO 340 GO TO 380 C 280 IF(NTRACE.LT.1) GO TO 310 WRITE(KW,290)FSTAR,NF 290 FORMAT(' FOBJ =',1PG24.16,5X,'REFLECTION SUCCEEDED',5X, * 'NF =',I7) WRITE(KW,300)(X(J),J=1,NV) 300 FORMAT(' X =',1PG15.7,4G15.7/(4X,5G15.7)) C C THE REFLECTED VALUE FSTAR IS LESS THAN THE LOWEST VALUE, C FZ(JL), IN THE SIMPLEX. C C ATTEMPT AN EXPANSION. FORM P** . C 310 DO 320 J=1,NV IF(MASK(J).EQ.0) * X(J)=X(J)+(GAMMA-UNITR)*(X(J)-ZBAR(J)) 320 CONTINUE C CALL SIFUN (FUNK) C C CHOOSE ONE OF TWO STRATEGIES FOR ACCEPTING THE EXPANSION C POINT P** . C IF((METHD.EQ.1 .AND. FOBJ.LT.FZ(JL)) .OR. * (METHD.EQ.2 .AND. FOBJ.LT.FSTAR)) GO TO 490 C IF(NTRACE.GE.1) WRITE(KW,330)FOBJ 330 FORMAT(' FOBJ =',1PG24.16,5X,'EXPANSION FAILED') C C THE REFLECTION FAILED BUT THE REFLECTED POINT WAS BETTER C THAN SOME OTHER POINT, OR ELSE THE REFLECTION SUCCEEDED C BUT THE EXPANSION FAILED. C C ACCEPT THE REFLECTED POINT. REPLACE P(JH) BY P* . C 340 IF(NTRACE.LT.1) GO TO 360 WRITE(KW,350) 350 FORMAT(36X,'ACCEPT REFLECTED POINT') WRITE(KW,300)(ZSTAR(J),J=1,NV) C 360 DO 370 J=1,NV Z(J,JH)=ZSTAR(J) 370 CONTINUE C FZ(JH)=FSTAR C C IF THE REFLECTION FAILED AND THE REFLECTED POINT WAS BETTER C THAN ONLY P(JH), CONTRACT THE SIMPLEX. C IF(KONTR.EQ.0) GO TO 530 C C ATTEMPT A CONTRACTION. FORM P** . C 380 JDIFF=0 C DO 400 J=1,NV IF(MASK(J).NE.0) GO TO 400 ZBARJ=ZBAR(J) ZJJH=Z(J,JH) XJ=ZBARJ+BETA*(ZJJH-ZBARJ) C C CHECK THE ROUNDING. C XJ MIGHT HAVE BEEN ROUNDED UP TO ZJJH, WHICH COULD CREATE C AN INFINITE LOOP. C XJ MUST LIE IN THE CLOSED INTERVAL (ZJJH,ZBARJ), AND C XJ MUST NOT BE EQUAL TO ZJJH UNLESS ZJJH.EQ.ZBARJ . C EQUIVALENTLY, AN OPEN INTERVAL (ZJJH,ZBARJ) MUST EXIST AND C XJ MUST LIE IN IT, OR XJ MUST BE EQUAL TO ZBARJ . C IF((XJ.GT.ZJJH .AND. XJ.LT.ZBARJ) .OR. * (XJ.GT.ZBARJ .AND. XJ.LT.ZJJH) .OR. * XJ.EQ.ZBARJ) GO TO 390 NOW=2 IF(NTRACE.GE.1) WRITE(KW,180)J,NOW,NF XJ=ZBARJ C 390 IF(XJ.NE.ZJJH) JDIFF=7 X(J)=XJ 400 CONTINUE C IF(JDIFF.EQ.0) GO TO 420 CALL SIFUN (FUNK) IF(FOBJ.GT.FZ(JH)) GO TO 420 IF(NTRACE.LT.1) GO TO 420 WRITE(KW,410)FOBJ 410 FORMAT(' FOBJ =',1PG24.16,5X,'CONTRACTION SUCCEEDED') WRITE(KW,300)(X(J),J=1,NV) GO TO 510 C 420 IF(NTRACE.GE.1) WRITE(KW,430)FOBJ 430 FORMAT(' FOBJ =',1PG24.16,5X,'CONTRACTION FAILED') C 440 JS=JL C C REPLACE ALL P(J) BY (P(J)+P(JL))/2.0 . C DO 470 J=1,NVPLUS IF(J.EQ.JL) GO TO 470 C DO 460 K=1,NV IF(MASK(K).NE.0) GO TO 460 ZKJ=Z(K,J) ZKJL=Z(K,JL) XK=ZKJL+(ZKJ-ZKJL)/RTWO C C CHECK THE ROUNDING. C XK MIGHT HAVE BEEN ROUNDED UP TO ZKJ, WHICH COULD CREATE C AN INFINITE LOOP. C XK MUST LIE IN THE CLOSED INTERVAL (ZKJ,ZKJL), AND C XK MUST NOT BE EQUAL TO ZKJ UNLESS ZKJ.EQ.ZKJL . C EQUIVALENTLY, AN OPEN INTERVAL (ZKJ,ZKJL) MUST EXIST AND C XK MUST LIE IN IT, OR XK MUST BE EQUAL TO ZKJL . C IF((XK.GT.ZKJ .AND. XK.LT.ZKJL) .OR. * (XK.GT.ZKJL .AND. XK.LT.ZKJ) .OR. * XK.EQ.ZKJL) GO TO 450 NOW=3 IF(NTRACE.GE.1) WRITE(KW,180)J,NOW,NF XK=ZKJL C 450 Z(K,J)=XK X(K)=XK 460 CONTINUE C CALL SIFUN (FUNK) IF(FOBJ.LT.FZ(JS)) JS=J FZ(J)=FOBJ 470 CONTINUE C IF(JS.EQ.JL) GO TO 530 JL=JS IF(NTRACE.LT.1) GO TO 530 WRITE(KW,480)FZ(JS) 480 FORMAT(' FZ(JS) =',1PG24.16) WRITE(KW,300)(Z(K,JS),K=1,NV) GO TO 530 C 490 IF(NTRACE.LT.1) GO TO 510 WRITE(KW,500)FOBJ 500 FORMAT(' FOBJ =',1PG24.16,5X,'EXPANSION SUCCEEDED') WRITE(KW,300)(X(J),J=1,NV) C C REPLACE P(JH) BY P* OR P** . C 510 DO 520 J=1,NV Z(J,JH)=X(J) 520 CONTINUE C FZ(JH)=FOBJ C C TEST FOR CONVERGENCE. C 530 DO 550 J=1,NV IF(MASK(J).NE.0) GO TO 550 ZMAX=Z(J,1) ZMIN=Z(J,1) C DO 540 K=2,NVPLUS ZJK=Z(J,K) IF(ZJK.GT.ZMAX) ZMAX=ZJK IF(ZJK.LT.ZMIN) ZMIN=ZJK 540 CONTINUE C DZ=ZMAX-ZMIN IF(DZ.GT.DELMIN(J)) GO TO 560 550 CONTINUE C GO TO 640 C C RETURN IF THE SENSE SWITCH IS ON. C 560 JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.LE.1) GO TO 570 C IF(NF.GT.NFMAX) GO TO 590 GO TO 10 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C 570 KFLAG=-3 IF(NTRACE.GE.-1) WRITE(KW,580)NSSW 580 FORMAT(///' ABNORMAL TERMINATION...'//6X,'TERMINATED BY', * ' OPERATOR VIA SENSE SWITCH ',I2) GO TO 610 C 590 KFLAG=-2 IF(NTRACE.GE.-1) WRITE(KW,600)NFMAX 600 FORMAT(////' ABNORMAL TERMINATION...'//6X,'MORE THAN', * ' NFMAX =',I11,' CALLS TO THE FOBJ SUBROUTINE.') C 610 IF(NTRACE.LT.-1) GO TO 680 C DO 630 J=1,NVPLUS WRITE(KW,620)J,FZ(J),(Z(K,J),K=1,NV) 620 FORMAT(/' SIMPLEX POINT',I3,' .... FOBJ =',1PG24.16, * 9X,'X(J)....'/(1X,5G15.7)) 630 CONTINUE C GO TO 680 C 640 KFLAG=1 IF(NTRACE.GE.0) WRITE(KW,650) 650 FORMAT(///' TERMINATED WHEN THE DIMENSIONS OF THE'/6X, * 'SIMPLEX BECAME AS SMALL AS THE DELMIN(J).') GO TO 680 C 660 KFLAG=2 IF(NTRACE.GE.0) WRITE(KW,670) 670 FORMAT(///' TERMINATED WHEN THE FUNCTION VALUES AT'/6X, * 'ALL POINTS OF THE SIMPLEX WERE EXACTLY EQUAL.') C 680 CONTINUE C C GO BACK AND COMPUTE JL AND JH ONE LAST TIME. C GO TO 10 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C FINISH UP AND EXIT. C 690 DO 700 J=1,NV X(J)=Z(J,JL) 700 CONTINUE C CALL SIFUN (FUNK) IF(FOBJ.LE.FSAVE .AND. FOBJ.EQ.FZ(JL)) GO TO 720 NOREP=NOREP+2 IF(NTRACE.GE.-1) WRITE(KW,710)NF,FSAVE,FZ(JL),FOBJ 710 FORMAT(////' WARNING.... FOBJ IS NOT A REPRODUCIBLE', * ' FUNCTION OF X(J).',5X,'NF =',I7//5X,1PG24.16, * 2G24.16) C 720 IF(NTRACE.GE.0) WRITE(KW,730)NF,FOBJ,(X(J),J=1,NV) 730 FORMAT(///1X,I7,' FUNCTION COMPUTATIONS'/// * ' FINAL VALUE OF FOBJ =',1PG24.16/// * 9X,'FINAL VALUES OF X(J)....'//(1X,5G15.7)) C C THE FOLLOWING STATEMENT IS THE ONLY RETURN IN THIS ROUTINE. C RETURN C C END SIMPLEX C END SUBROUTINE SIBEG (FUNK) C C SIBEG 1.3 OCTOBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1965, 1975, 1990 J. P. CHANDLER, C COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C SIBEG SETS DEFAULT VALUES AND PRINTS INITIAL OUTPUT C FOR SIMPLEX. C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C INPUT QUANTITIES..... FUNK,X(*),XMAX(*),XMIN(*),DELTX(*), C DELMIN(*),NV,NTRACE,MASK(*), C NFMAX,NFLAT,NXTRA,KW C C OUTPUT QUANTITIES.... FZ,HUGE,ALPHA,BETA,GAMMA,MAXPT, C NVPLUS,NSSW,NF,KFLAG,NOREP, C DELTX(*), AND DELMIN(*) C C THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME C COMPILERS (WATFIV, FOR EXAMPLE) AND FORBIDDEN BY OTHERS C (THE MODCOMP II, FOR EXAMPLE). C EXTERNAL FUNK C DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * Z,FZ,HUGE,ALPHA,BETA,GAMMA,DELDF,XS,XPLUS,FSAVE,RZERO C INTEGER J,JUMP,JVARY,K,KERFL,KFLAG,KTYPE,KW,MASK, * MATRX,MAXPT,NACTIV,NF,NFLAT,NFMAX,NOREP,NSSW, * NTRACE,NV,NVMAX,NVPLUS,NXTRA C DIMENSION Z(20,21) 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 SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT, * NVPLUS,NSSW,NF C EQUIVALENCE (Z(1,1),ERR(1,1)) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C SET FIXED QUANTITIES .... C C KTYPE = CONSOLE TYPEWRITER UNIT NUMBER C (IRRELEVANT IF A DUMMY DATSW ROUTINE IS USED) C KTYPE=1 C C NSSW = SENSE SWITCH NUMBER FOR SEARCH TERMINATION C (IRRELEVANT IF A DUMMY DATSW ROUTINE IS USED) C NSSW=6 C C HUGE = A VERY LARGE REAL NUMBER C (DEFAULT VALUE FOR XMAX AND -XMIN, AND C OUT-OF-BOUNDS VALUE FOR FOBJ) C HUGE=1.0D35 C C NVMAX IS THE MAXIMUM VALUE OF NV. C NVMAX IS ALSO THE DIMENSION OF X, XMAX, XMIN, DELTX, DELMIN, C MASK, ZBAR, ZSTAR, AND THE FIRST DIMENSION OF ERR AND Z. C NVMAX=20 C C MAXPT IS THE MAXIMUM NUMBER OF POINTS IN THE SIMPLEX. C MAXPT IS ALSO THE DIMENSION OF FZ, AND THE SECOND DIMENSION C OF ERR AND Z. C MAXPT MUST BE .GE. (NVMAX+1). C MAXPT=21 C C DELDF = DEFAULT VALUE FOR DELTX(J) C DELDF=0.01D0 C C ALPHA = REFLECTION PARAMETER C ALPHA=1.0D0 C C BETA = CONTRACTION PARAMETER C BETA=0.5D0 C C GAMMA = EXPANSION PARAMETER C GAMMA=2.0D0 C RZERO=0.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 C MAKE SURE THAT THE SENSE SWITCH IS OFF. C JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.EQ.2) GO TO 30 C C THIS IS THE ONLY USAGE OF THE CONSOLE TYPEWRITER. C WRITE(KTYPE,10)NSSW 10 FORMAT(/' TURN OFF SENSE SWITCH ',I2//' ') 20 CALL DATSW (NSSW,JUMP) IF(JUMP.EQ.1) GO TO 20 C 30 JVARY=0 C IF(NV.GE.1 .AND. NV.LE.NVMAX .AND. NVMAX.LE.MAXPT) * GO TO 50 KFLAG=-1 WRITE(KW,40)NV,NVMAX,MAXPT 40 FORMAT(//' ERROR IN INPUT TO SUBROUTINE SIMPLEX.'/6X, * 'NV =',I14,6X,'NVMAX =',I14,6X,'MAXPT =',I14) IF(NTRACE.GE.-1) WRITE(KW,40)NV,NVMAX,MAXPT,NACTIV STOP C 50 DO 100 J=1,NV IF(MASK(J).NE.0) GO TO 100 IF(DELMIN(J).LT.RZERO) DELMIN(J)=-DELMIN(J) C C CHECK THAT DELTX(J) IS NOT NEGLIGIBLE. C XPLUS=X(J)+DELTX(J) IF(XPLUS.EQ.X(J)) GO TO 60 XPLUS=X(J)-DELTX(J) IF(XPLUS.NE.X(J)) GO TO 80 C 60 IF(X(J).EQ.RZERO) GO TO 70 DELTX(J)=DELDF*X(J) GO TO 80 C 70 DELTX(J)=DELDF C 80 IF(XMAX(J).GT.XMIN(J)) GO TO 90 XMAX(J)=HUGE XMIN(J)=-HUGE C 90 IF(X(J).GT.XMAX(J)) X(J)=XMAX(J) IF(X(J).LT.XMIN(J)) X(J)=XMIN(J) 100 CONTINUE C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C IF(NTRACE.LT.0) GO TO 180 WRITE(KW,110) 110 FORMAT('1SIMPLEX MINIMIZATION'// * ' INITIAL VALUES....'/' ') WRITE(KW,120)(MASK(J),J=1,NV) 120 FORMAT(/' MASK =',I7,4I13/(2X,5I13)) WRITE(KW,130)(X(J),J=1,NV) 130 FORMAT(/' X =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,140)(XMAX(J),J=1,NV) 140 FORMAT(/' XMAX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,150)(XMIN(J),J=1,NV) 150 FORMAT(/' XMIN =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,160)(DELTX(J),J=1,NV) 160 FORMAT(/' DELTX =',1PG13.5,4G13.5/(8X,5G13.5)) WRITE(KW,170)(DELMIN(J),J=1,NV) 170 FORMAT(/' DELMIN=',1PG13.5,4G13.5/(8X,5G13.5)) C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CALCULATE THE INITIAL P(I) AND Y(I) (Z AND FZ, C RESPECTIVELY). C C NF = NUMBER OF CALLS TO FUNK SO FAR C 180 NF=0 C C NACT = NUMBER OF ACTIVE X(J) C NACTIV=0 C DO 200 J=1,NV IF(MASK(J).NE.0) GO TO 200 NACTIV=NACTIV+1 C DO 190 K=1,NV Z(K,NACTIV)=X(K) 190 CONTINUE C Z(J,NACTIV)=Z(J,NACTIV)+DELTX(J) XS=X(J) X(J)=Z(J,NACTIV) CALL SIFUN (FUNK) X(J)=XS FZ(NACTIV)=FOBJ 200 CONTINUE C IF(NTRACE.GE.0) WRITE(KW,210)NV,NACTIV,NFMAX, * NFLAT,NXTRA,ALPHA,BETA,GAMMA 210 FORMAT(//1X,I3,' VARIABLES,',I3,' ACTIVE.'//9X,'NFMAX =', * I8,9X,'NFLAT =',I3,9X,'NXTRA =',I3//9X, * 'ALPHA =',F5.2,9X,'BETA =',F5.2,9X,'GAMMA =',F5.2) C IF(NACTIV.GT.0) GO TO 230 C KFLAG=-1 IF(NTRACE.GE.-1) WRITE(KW,220) 220 FORMAT(//' WARNING ... MASK(J).EQ.0 FOR ALL J IN', * ' SUBROUTINE SMPLX.'// * ' FOBJ WILL BE EVALUATED BUT NOT MINIMIZED.') C C SET THE BASE POINT. C 230 NVPLUS=NACTIV+1 C DO 240 K=1,NV Z(K,NVPLUS)=X(K) 240 CONTINUE C CALL SIFUN (FUNK) FZ(NVPLUS)=FOBJ C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C CHECK THE REPRODUCIBILITY OF FOBJ, AND PRINT THE REMAINING C LINES. C FSAVE=FOBJ CALL SIFUN (FUNK) IF(FOBJ.EQ.FSAVE) GO TO 260 NOREP=1 IF(NTRACE.GE.-1) WRITE(KW,250)NF,FSAVE,FOBJ 250 FORMAT(///' WARNING.... FOBJ IS NOT A REPRODUCIBLE', * ' FUNCTION OF X(J).',8X,'NF =',I5//5X,1PG24.16, * 2G24.16) C 260 IF(NTRACE.GE.0) WRITE(KW,270)FOBJ 270 FORMAT(///' FOBJ =',1PG24.16///' BEGIN MINIMIZATION....') C RETURN C C END SIBEG C END SUBROUTINE SIFUN (FUNK) C C SIFUN 1.2 DECEMBER 1991 C C A.N.S.I. STANDARD FORTRAN 77 C C COPYRIGHT (C) 1975, 1990 BY J. P. CHANDLER, C COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY C C SIFUN CALLS FUNK TO COMPUTE FOBJ, IF X IS IN BOUNDS. C SIFUN IS CALLED BY SIMPLEX AND SIBEG. C DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, * FZ,HUGE,ALPHA,BETA,GAMMA C INTEGER JF,JVARY,KERFL,KFLAG,KW,MASK,MATRX,MAXPT, * NF,NFLAT,NFMAX,NOREP,NSSW,NTRACE,NV,NVPLUS,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 SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT, * NVPLUS,NSSW,NF C DO 10 JF=1,NV C IF(MASK(JF).EQ.0 .AND. * (X(JF).GT.XMAX(JF) .OR. X(JF).LT.XMIN(JF))) * GO TO 20 C 10 CONTINUE C CALL FUNK NF=NF+1 GO TO 30 C 20 FOBJ=HUGE C 30 IF(NTRACE.GE.2) WRITE(KW,40)FOBJ,(X(JF),JF=1,NV) 40 FORMAT(' TRIAL FOBJ =',1PG24.16,9X,'X(J)....'/ * (1X,3G24.16)) C RETURN C C END SIFUN C END SUBROUTINE DATSW (NSSW,JUMP) C C DUMMY VERSION OF SUBROUTINE DATSW -- ALL SWITCHES OFF. C C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, C OKLAHOMA STATE UNIVERSITY 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 3 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 //