Ricerca dei massimi - minimi di una funzione a più variabili assegnate le condizioni al contorno col metodo del simplesso scritta in fortran 77
C J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT,
C OKLAHOMA STATE UNIVERSITY
//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
//