Metodo del simplesso in f77

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

//