!****************************************************************** !* MULTI-DIMENSIONAL CURVE FITTING BY SIMPLEX METHOD * !* -------------------------------------------------------------- * !* Reference: * !* * !* J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, * !* OKLAHOMA STATE UNIVERSITY * !* -------------------------------------------------------------- * !* Explanations: * !* * !* * !* * !* THIS PROGRAM FITS THE MODEL: * !* * !* * !* A * !* Z = ---------------- * !* 1 + BX + CY * !* * !* WRITTEN AS: * !* * !* FIT(J) = X(1)/(1.0 + X(2)*T(J,1) + X(3)*T(J,2)) * !* * !* TO DATA POINTS (T(J,1),T(J,2),Y(J)). * !* * !* * !* * !* T(J,1) AND T(J,2) ARE THE TWO INDEPENDENT VARIABLES X AND Y * !* IN DATA SPACE. Y(J) IS THE MEASURED DEPENDENT VARIABLE Z AND * !* FIT(J) IS THE VALUE TO BE FITTED TO Y(J) BY MINIMIZING THE * !* * !* * !* WEIGHTED SUM OF SQUARES, * !* * !* J=NPTS * !* * !* PHI = SUM ((FIT(J)-Y(J))/YSIG(J))**2 * !* * !* J=1 * !* * !* * !* THE VALUE OF PHI IS STORED IN THE VARIABLE FOBJ. * !* THE MODEL Z=F(X,Y) IS IMPLEMENTED IN FUNCTION FUNK. * !* THE RESULTS ARE BEST COEFFICIENTS A, B, C, STORED IN X(1), * !* X(2), X(3) TO FIT THE DATA PROVIDED IN FILE SIMPTEST.DAT. * !* -------------------------------------------------------------- * !* SAMPLE RUN: * !* File Simptest.dat contains: * !* 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 * !* * !* Output file Simpout.lst contains at the end: * !* --/-- * !* * !* TERMINATED WHEN THE DIMENSIONS OF THE SIMPLEX BECAME AS SMALL * !* AS THE DELMIN(J) *. !* * !* * !* 543 FUNCTION COMPUTATIONS * !* * !* FINAL VALUE OF FOBJ = 6.504033783066415 * !* * !* FINAL VALUES OF X(J).... * !* * !* 48.20289 2.876328 1.903036 * !* * !* * !* + 3.6498 * !* X( 1) = 48.202888 STDEVS = 1.0000 * !* - 3.0447 * !* * !* * !* + 0.26378 * !* X( 2) = 2.8763282 STDEVS = 1.0000 * !* - 0.22393 * !* * !* * !* + 0.19636 * !* X( 3) = 1.9030361 STDEVS = 1.0000 * !* - 0.16812 * !* * !* -------------------------------------------------------------- * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !****************************************************************** !Note: To change example, adjust function FUNK to your problem and check ! NV, NPTS and size of T(300,3). PROGRAM TEST_SMPLX IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL FUNK DIMENSION ERGSAV(20),FIDINT(2) 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 COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, & KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CSTOR/ T(300,2) OPEN (UNIT=5,FILE='SIMPTEST.DAT') OPEN (UNIT=6,FILE='SIMPOUT.LST') KR=5 KW=6 WRITE(KW,*) ' ' WRITE(KW,*) ' ** SIMPLEX PROGRAM **' ! READ IN THE VALUE OF JROUTE... ! JROUTE=1 TO USE MARQ, ! JROUTE=2 TO USE STEPIT, ! JROUTE=3 TO USE SIMPLEX READ (KR,10) JROUTE 10 FORMAT(I1) ! READ IN THE VALUE OF NPTS, THEN READ IN THE DATA POINTS. READ(KR,20)NPTS 20 FORMAT(I5) WRITE(KW,30)NPTS 30 FORMAT(//' NPTS =',I5) 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 ! INITIALIZE FOR THE FIT. CALL STSET NV=3 NTRACE=0 !Normal output ! NTRACE=1 For additional output NFMAX=550 X(1)=10.0 X(2)=1.0 X(3)=1.0 ! FIT THE MODEL TO THE DATA. CALL STINT (FUNK) ! GO TO 90 (TO SKIP FIDO) ! USE SUBROUTINE FIDO TO CHECK THE ERRORS FROM MARQ ! OR STEPIT, OR TO COMPUTE ERRORS USING SIMPLEX. ERFRAC=0.1D0 DO JXFID=1,NV IF(JROUTE.EQ.1 .OR. JROUTE.EQ.2) THEN ERGSAV(JXFID)=DSQRT(DABS(ERR(JXFID,JXFID))) ELSE ERGSAV(JXFID)=ERFRAC*DABS(X(JXFID)) ENDIF IF(ERGSAV(JXFID).EQ.0.0D0) ERGSAV(JXFID)=ERFRAC END DO STDEVS=1.0D0 TOLFID=0.05D0 MAXIND=2 NTRACF=-1 !0 OR 1 TO PRINT MORE IN FIDO IF(JROUTE.EQ.1) THEN NTRSET=-2 ELSE NTRSET=-1 ENDIF DO JXFID=1,NV ERGUES=ERGSAV(JXFID) CALL FIDO (FUNK,JXFID,STDEVS,TOLFID,ERGUES,MAXIND, & NTRSET,NTRACF,JROUTE,FIT,FIDINT) END DO 90 CLOSE (UNIT=KW) print *,' ' print *,' Results in file Simpout.lst.' Pause ' Program terminated.' STOP END !of main program SUBROUTINE STINT (FUNK) !-------------------------------------------------------------- ! COMMON INTERFACE ROUTINE TO FACILITATE SWAPPING OF STEPIT, ! SIMPLEX, MARQ, ETC. !-------------------------------------------------------------- EXTERNAL FUNK CALL SMPLX (FUNK) RETURN END SUBROUTINE FUNK !------------------------------------------------------------- ! J. P. CHANDLER, COMPUTER SCIENCE DEPT., ! OKLAHOMA STATE UNIVERSITY ! ! THIS VERSION OF SUBROUTINE FUNK COMPUTES THE FITTED VALUES ! AND THE WEIGHTED SUM OF SQUARES FOR THE PROBLEM OF FITTING ! THE MODEL ! ! FIT(J) = X(1)/(1.0 + X(2)*T(J,1) + X(3)*T(J,2)) ! ! TO DATA POINTS (T(J,1),T(J,2),Y(J)) . ! ! TO RUN THIS ROUTINE WITH MARQ, USE THE STATEMENT ! SUBROUTINE FUNK (FITM) . ! ! TO RUN THIS ROUTINE WITH STEPIT OR SIMPLEX, USE ! SUBROUTINE FUNK . !-------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FITM(300) 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 COMMON /CDAT/ FIT(300),Y(300),YSIG(300),NPTS COMMON /CSTOR/ T(300,2) ! FOBJ WILL CONTAIN THE WEIGHTED SUM OF SQUARES, PHI. ! MARQ REQUIRES THAT FUNK COMPUTE FITM. ! WHETHER OR NOT FUNK COMPUTES FOBJ IS IRRELEVANT TO MARQ. ! STEPIT, SIMPLEX, ETC. REQUIRE THAT FUNK COMPUTE FOBJ. ! THOSE ROUTINES NEVER SEE THE FITTED VALUES AT ALL. ! TO ALLOW MARQ TO BE INTERCHANGED WITH STEPIT, SIMPLEX, ETC. ! WITH A MINIMUM NUMBER OF CHANGES, WE WILL COMPUTE FOBJ ! AS WELL AS FITM. FOBJ=0.0D0 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 RETURN END SUBROUTINE FIDO (FUNK,JXFID,STDEVS,TOLFID,ERGUES,MAXIND, & NTRSET,NTRACF,JROUTE,FITM, FIDINT) !--------------------------------------------------------------- ! FIDO 4.5 DECEMBER 1991 ! A.N.S.I. STANDARD FORTRAN 77 ! ! J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, ! OKLAHOMA UNIVERSITY, STILLWATER, OKLAHOMA 74078 ! ! SUBROUTINE FIDO IS USED IN CONJUNCTION WITH MARQ, STEPIT, OR ! SIMPLEX TO COMPUTE CONFIDENCE HALF-INTERVALS ("ERRORS") FOR ! THE PARAMETERS IN A FITTING PROBLEM, USING THE METHOD OF ! SUPPORT PLANES. ! THE QUANTITY BEING MINIMIZED (FOBJ) MUST EITHER BE ! CHI-SQUARE OR IT MUST BE TWICE THE NEGATIVE OF THE NATURAL ! LOGARITHM OF THE LIKELIHOOD FUNCTION. ! ! ------------------------------------------------------------ ! USAGE...... ! FIRST CALL STEPIT TO FIND THE OPTIMUM VALUES OF THE ! PARAMETERS. THEN CALL FIDO ONCE FOR EACH PARAMETER FOR ! WHICH THE CONFIDENCE INTERVALS ARE DESIRED. ! ! INPUT QUANTITIES..... FUNK,JXFID,STDEVS,TOLFID,ERGUES, ! MAXIND,NTRSET,NTRACF,JROUTE,X(*), ! KW ! ! OUTPUT QUANTITY...... FIDINT( ) ! ! SCRATCH ARRAY (FOR USE WITH MARQ ONLY)... FITM(*) ! ! FUNK -- THE NAME OF THE SUBROUTINE WHICH COMPUTES ! FOBJ GIVEN THE VALUES OF THE X(J) ! ! JXFID -- THE INDEX OF THE PARAMETER, X(JXFID), FOR ! WHICH CONFIDENCE HALF-INTERVALS ARE TO ! BE COMPUTED ! ! STDEVS -- THE NUMBER OF STANDARD DEVIATIONS TO WHICH ! THE HALF-INTERVALS ARE TO CORRESPOND ! ! TOLFID -- CONVERGENCE TOLERANCE (USE TOLFID=0.05) ! ! ERGUES -- ROUGH ESTIMATE OF THE LENGTH OF THE ! HALF-INTERVALS ! ! MAXIND -- =1 TO COMPUTE THE HALF-INTERVAL WITH THE ! SAME SIGN AS ERGUES, ! =2 TO COMPUTE BOTH HALF-INTERVALS ! ! NTRSET -- VALUE OF NTRACE TO BE USED IN THE ! FITTING SUBROUTINE (STEPIT OR MARQ, ! ETC.) WHEN CALLED BY FIDO ! ! NTRACF -- A SWITCH THAT CONTROLS PRINTING... ! ! =-1 TO PRINT NOTHING IN FIDO EXCEPT ANY ! WARNING OR ERROR MESSAGES, ! ! = 0 FOR INITIAL AND SUMMARY PRINTOUT, ! ! =+1 TO PRINT EVERY FIDO ITERATION ! ! JROUTE -- =1 TO USE FIDO WITH SUBROUTINE MARQ, ! =2 TO USE FIDO WITH STEPIT OR SIMPLEX, ETC. ! ! FITM(*) -- A SCRATCH ARRAY, USED WHEN FIDO IS USED ! WITH MARQ ! ! X( ) -- THE POSITION OF THE MINIMUM OF FOBJ ! ! KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER ! ! FIDINT(J) -- RETURN THE COMPUTED LENGTHS OF THE ! CONFIDENCE HALF-INTERVALS ! ! IF STEPIT COMPUTES THE ERROR MATRIX ERR, ERGUES CAN BE SET ! TO PLUS OR MINUS STDEVS*SQRT(ERR(JX,JX)) (SEE BELOW). ! ! NOTE: ! FIDO CALLS STEPIT, WHICH MAY USE THE ARRAY ERR FOR SCRATCH ! STORAGE, DESTROYING ITS CONTENTS. ! ! FOR A LEAST SQUARES PROBLEM, FOBJ IS EQUAL TO CHI-SQUARE, ! THE WEIGHTED SUM OF SQUARES... ! ! NPTS ! FOBJ = SUM ((FIT(JPT)-YDATA(JPT))/YSIGMA(JPT))**2 ! JPT=1 ! ! THE STANDARD ERRORS YSIGMA( ) MUST BE CORRECTLY SCALED. ! IF THE YSIGMA( ) ARE NOT KNOWN ACCURATELY, NORMALIZE THEM ! SO THAT THE VALUE OF FOBJ AT THE MINIMUM IS EQUAL TO THE ! NUMBER OF DEGREES OF FREEDOM, ! N.D.F. = (NO. DATA POINTS) - (NO. ADJUSTABLE PARAMETERS) ! ! FIDO IS ESSENTIALLY A ROOT FINDING ROUTINE. IT USES INVERSE ! QUADRATIC INTERPOLATION OR EXTRAPOLATION, INVERSE LINEAR ! INTERPOLATION, AND BISECTION (INTERVAL HALVING), AS ! SUCCESSIVELY MORE DIFFICULTIES ARE ENCOUNTERED IN FINDING ! THE ROOT. ! ! PRIMARY REFERENCES.... ! ! W. T. EADIE ET AL., "STATISTICAL METHODS IN EXPERIMENTAL ! EXPERIMENTAL PHYSICS" (AMERICAN ELSEVIER, 1971), ! CHAPTER 9 ! ! P. R. BEVINGTON, "DATA REDUCTION AND ERROR ANALYSIS IN ! THE PHYSICAL SCIENCES" (MCGRAW-HILL, 1969), ! PAGES 242-245 ! ! OTHER REFERENCES.... ! ! H. SCHEFFE, "THE ANALYSIS OF VARIANCE" (WILEY, 1959) ! ! H. STONE, DISCUSSION ON PAPER BY E. M. L. BEALE, ! J. ROY. STATIS. SOC. B., V. 22 (1960), P. 41, PP. 84-5 ! ! G. E. P. BOX AND G. A. COUTIE, PROC. INST. ELEC. ENGRS. ! V. 103, PART B, SUPPL. NO. 1 (1956). ! ! G. W. BOOTH, G. E. P. BOX, M. E. MULLER, AND ! T. I. PETERSON, "FORECASTING BY GENERALIZED REGRESSION ! METHODS; NON-LINEAR ESTIMATION" (PRINCETON/IBM), ! FEBRUARY 1959, IBM MANUAL ! ! D. W. MARQUARDT, "LEAST-SQUARES ESTIMATION OF NONLINEAR ! PARAMETERS", SHARE DISTRIBUTION 3094 ! ! D. W. MARQUARDT, R. G. BENNETT, AND E. J. BURRELL, "LEAST ! SQUARES ANALYSIS OF ELECTRON PARAMAGNETI! RESONANCE ! SPECTRA", J. OF MOLECULAR SPECTROSCOPY 7 (1961) 269 ! !---------------------------------------------------------------- ! ! THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME ! COMPILERS (WATFIV, FOR EXAMPLE) AND IT IS FORBIDDEN BY SOME ! OTHERS (MODCOMP II, FOR EXAMPLE). ! EXTERNAL FUNK 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 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 DIMENSION FIDINT(2),FITM(1) DIMENSION XSAVE(20),XKA(20),XKB(20),XBEST(20) 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 ! SET THE LIBRARY FUNCTIONS FOR REAL PRECISION (SQRT, ETC.) ! OR FOR DOUBLE PRECISION (DSQRT, ETC.). ! ZSQRT(ARG)= SQRT(ARG) ZSQRT(ARG)=DSQRT(ARG) ! T(ARG)= ABS( SQRT(ARG-FSAVE)-STDEVS) T(ARG)=DABS(DSQRT(ARG-FSAVE)-STDEVS) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ! SET SOME PARAMETERS. ! ! RSMALL IS USED IN SETTING A PHONY VALUE IF A VALUE OF FOBJ ! IS FOUND THAT IS LESS THAN THE "GLOBAL MINIMUM" VALUE. ! RSMALL=0.0001D0 ! ! CONVRG = SATISFACTORY RATE OF CONVERGENCE ! CONVRG=0.5D0 ! ! BRACK = FACTOR FOR BRACKETTING SHOTS ! BRACK=2.0D0 ! ! MAXTRY = MAXIMUM NUMBER OF ITERATIONS ! MAXTRY=20 ! ! FACMAX = MAXIMUM FACTOR FOR QUADRATIC INTERPOLATION ! FACMAX=4.0D0 ! ! FRMIN = MINIMUM FRACTION FOR LINEAR INTERPOLATION ! FRMIN=0.1D0 RZERO=0.0D0 RHALF=0.5D0 UNITR=1.0D0 ! ! SAVE SOME INPUT QUANTITIES, AND INITIALIZE EVERYTHING. ! NTRSAV=NTRACE NTRACE=NTRSET MATSAV=MATRX MATRX=0 MSKSAV=MASK(JXFID) MASK(JXFID)=1 NACTIV=0 DO 10 J=1,NV XBEST(J)=X(J) XSAVE(J)=X(J) IF(MASK(J).EQ.0) NACTIV=NACTIV+1 10 CONTINUE IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF FSAVE=FOBJ FBEST=FOBJ FPSGSQ=FOBJ+STDEVS**2 X(JXFID)=XSAVE(JXFID)+STDEVS*ERGUES ! ! LOOP OVER THE FIDINT(JJ). ! JJ=1 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 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) 70 KLOSB=0 MLIN=0 DO 80 K=1,NV XKA(K)=XSAVE(K) 80 CONTINUE FB=FSAVE FKA=FSAVE ! ! BEGIN THE ITERATION LOOP FOR LOCATING THE PLANE ! PERPENDICULAR TO THE X(JXIFD) AXIS IN WHICH THE MINIMUM ! VALUE OF FOBJ IS EQUAL TO FPSGSQ = FSAVE + STDEVS**2 , ! WHERE FSAVE IS THE VALUE OF FOBJ AT THE GLOBAL MINIMUM. ! DO 310 ITER=1,MAXTRY FC=FB IF(NACTIV.EQ.0) THEN IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF ELSE CALL STINT (FUNK) ENDIF IF(FOBJ.GE.FBEST) GO TO 100 FBEST=FOBJ DO 90 K=1,NV XBEST(K)=X(K) 90 CONTINUE 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) IF(FOBJ.EQ.FSAVE) GO TO 140 IF(FOBJ.LT.FSAVE) GO TO 120 ! TEST FOR CONVERGENCE. ! ! IF(ABS(FOBJ-FPSGSQ).LE.TOL) ............................ ! 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 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) ! CHECK FOR TIGHTER BRACKETTING OF FPSGSQ FROM BELOW. 150 IF((X(JXFID)-XKA(JXFID))*SGN.LE.RZERO) GO TO 190 DO 160 K=1,NV XKA(K)=X(K) 160 CONTINUE FKA=FOBJ GO TO 190 ! CHECK FOR BRACKETTING OF FPSGSQ FROM ABOVE. 170 IF(KLOSB.EQ.1) THEN IF((X(JXFID)-XKB(JXFID))*SGN.GE.RZERO) GO TO 190 ENDIF KLOSB=1 DO 180 K=1,NV XKB(K)=X(K) 180 CONTINUE FKB=FOBJ 190 IF(T(FOBJ).LT.T(FC)) FB=FOBJ ! CHECK THE RATE OF CONVERGENCE. IF IT IS SATISFACTORY, AND ! IF LINEAR INTERPOLATION HAS BEEN USED, USE IT AGAIN. IF(ITER.GE.2 .AND. T(FOBJ).GT.CONVRG*T(FC)) GO TO 210 IF(MLIN.EQ.1) GO TO 250 ! USE INVERSE QUADRATIC INTERPOLATION. FAC=STDEVS/ZSQRT(FOBJ-FSAVE) ! FAC=AMIN1(FACMAX,AMAX1(FAC,1./FACMAX)) ................... IF(FAC.LT.UNITR/FACMAX) FAC=UNITR/FACMAX IF(FAC.GT.FACMAX) FAC=FACMAX XX=XSAVE(JXFID)+FAC*(X(JXFID)-XSAVE(JXFID)) ! CHECK THAT THE PROPOSED POINT IS INSIDE THE BRACKETTED INTERVAL. 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 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 GO TO 310 210 IF(KLOSB.EQ.1) GO TO 240 ! CONVERGENCE IS POOR, AND FPSGSQ HAS NOT YET BEEN BRACKETTED. ! TRY TO BRACKET IT. FRAC=BRACK IF(FOBJ.GE.FPSGSQ) FRAC=UNITR/BRACK 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 IF(NTRACF.GE.1) WRITE(KW,230) 230 FORMAT(//' BRACKETTING SHOT....') GO TO 310 240 IF(MLIN.EQ.1) GO TO 270 ! TRY LINEAR INTERPOLATION BETWEEN THE TWO BRACKETTING POINTS. 250 MLIN=1 FRAC=(FPSGSQ-FKA)/(FKB-FKA) ! FRAC=AMAX1(FRMIN,AMIN1(1.0-FRMIN,FRAC)) ............. IF(FRAC.GT.UNITR-FRMIN) FRAC=UNITR-FRMIN IF(FRAC.LT.FRMIN) FRAC=FRMIN CONTINUE IF(NTRACF.GE.1) WRITE(KW,260) 260 FORMAT(//' LINEAR INTERPOLATION....') GO TO 290 ! CONVERGENCE IS POOR, AND LINEAR INTERPOLATION HAS BEEN USED. ! BISECT THE BRACKETTED INTERVAL. 270 FRAC=RHALF IF(NTRACF.GE.1) WRITE(KW,280) 280 FORMAT(//' INTERVAL BISECTED....') 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 310 CONTINUE ! END OF ITERATION LOOP. THE ITERATION FAILED TO CONVERGE. WRITE(KW,320) 320 FORMAT(' CONVERGENCE FAILURE IN SUBROUTINE FIDO.') FIDINT(JJ)=RZERO GO TO 350 ! CONVERGENCE ACHIEVED ... A SATISFACTORY PLANE HAS BEEN LOCATED. 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, & ')') 350 JJ=JJ+1 IF(JJ.GT.MAXIND) GO TO 370 ! REFLECT X THROUGH XSAVE, AND SEARCH FOR THE OTHER ! HALF-INTERVAL. DO 360 K=1,NV IF(K.EQ.JXFID .OR. MASK(K).EQ.0) X(K)=XSAVE(K)+(XSAVE(K)-X(K)) 360 CONTINUE ERGUES=X(JXFID)-XSAVE(JXFID) GO TO 20 ! END OF THE JJ LOOP. PRINT SUMMARY RESULTS. 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) ! PRINT A SUMMARY OF RESULTS FROM FIDO... WRITE(KW,390)FIDINT(1),JXFID, XSAVE(JXFID),STDEVS,XX 390 FORMAT(24X,'+',1PG12.5/' X(',I3,') = ',G15.8,25X, & ' STDEVS = ',G12.5/24X,'-',G12.5) ! RESTORE THE SAVED ENTRY VALUES. 400 DO 410 J=1,NV X(J)=XBEST(J) 410 CONTINUE IF(JROUTE.EQ.1) THEN CALL FUNK (FITM) ELSE CALL FUNK ENDIF 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 RETURN END SUBROUTINE SMPLX (FUNK) !------------------------------------------------------------- ! SIMPLEX 2.12 DECEMBER 1991 ! ! A.N.S.I. STANDARD FORTRAN 77 ! ! COPYRIGHT (C) 1965, 1975, 1991 J. P. CHANDLER ! (PRESENT ADDRESS ... ! COMPUTER SCIENCE DEPARTMENT, ! OKLAHOMA STATE UNIVERSITY, ! STILLWATER, OKLAHOMA 74078 ! (405)-744-5676 ) ! ! SIMPLEX FINDS LOCAL MINIMA OF A SMOOTH FUNCTION OF SEVERAL ! PARAMETERS. IT WILL ALSO HANDLE SOME NON-SMOOTH FUNCTIONS. ! THE METHOD IS OFTEN VERY SLOW IF THE NUMBER OF PARAMETERS IS ! AT ALL LARGE (GREATER THAN ABOUT SIX TO EIGHT). ! ! "A SIMPLEX METHOD FOR FUNCTION MINIMIZATION", ! J. A. NELDER AND R. MEAD, ! THE COMPUTER JOURNAL 7 (1965) 308-313 ! ! FOR APPLICATIONS, SEE ! ! "'DIRECT SEARCH' SOLUTION OF ! NUMERICAL AND STATISTICAL PROBLEMS", ! ROBERT HOOKE AND T. A. JEEVES, JOURNAL OF THE ! ASSOCIATION FOR COMPUTING MACHINERY 8 (1961) 212-229 ! ! "THE NELDER-MEAD SIMPLEX PROCEDURE FOR FUNCTION ! MINIMIZATION", ! D. M. OLSSON AND L. S. NELSON, ! TECHNOMETRICS 17 (1975) 45-51 ! ! ! SIMPLEX 2.9 IS AVAILABLE FROM THE ! QUANTUM CHEMISTRY PROGRAM EXCHANGE ! DEPT. OF CHEMISTRY, INDIANA UNIVERSITY ! BLOOMINGTON, INDIANA 47401 ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ! INPUT QUANTITIES..... FUNK,NV,NTRACE,MASK(*),X(*), ! XMAX(*),XMIN(*),DELTX(*), ! DELMIN(*),NFMAX,NFLAT,NXTRA,KW ! ! OUTPUT QUANTITIES.... X(*),FOBJ, KFLAG,NOREP ! ! FUNK -- THE NAME OF THE SUBROUTINE THAT COMPUTES ! FOBJ GIVEN X(1),...,X(NV) ! (EACH SUCH SUBROUTINE MUST BE NAMED IN ! AN EXTERNAL STATEMENT IN THE CALLING ! PROGRAM) ! ! NV -- THE NUMBER OF PARAMETERS, X(J) ! ! NTRACE -- =0 FOR NORMAL OUTPUT, ! =+1 FOR TRACE OUTPUT, ! =+2 FOR DEBUG OUTPUT, ! =-1 FOR NO OUTPUT EXCEPT ERROR MESSAGES ! ! MATRX -- USED BY STEPIT PROPER BUT NOT BY SIMPLEX ! ! FOBJ -- THE VALUE OF THE FUNCTION TO BE MINIMIZED ! ! MASK(J) -- NONZERO IF X(J) IS TO BE HELD FIXED ! ! X(J) -- THE J-TH PARAMETER ! ! XMAX(J) -- THE UPPER LIMIT ON X(J) ! ! XMIN(J) -- THE LOWER LIMIT ON X(J) ! ! DELTX(J) -- THE INITIAL STEP SIZE FOR X(J) ! ! DELMIN(J) -- THE LOWER LIMIT (CONVERGENCE TOLERANCE) ON ! THE STEP SIZE FOR X(J) ! ! ERR(J,K) -- SCRATCH STORAGE ! ! NFMAX -- THE MAXIMUM NUMBER OF FUNCTION ! COMPUTATIONS TO BE ALLOWED ! ! NFLAT -- NONZERO IF THE SEARCH IS TO TERMINATE WHEN ! ALL TRIAL STEPS GIVE IDENTICAL FUNCTION ! VALUES. ! THE RECOMMENDED VALUE OF NFLAT IS ! USUALLY NFLAT=1 . ! ! JVARY -- USED BY STEPIT AND STP BUT NOT BY SIMPLEX ! (SIMPLEX SETS JVARY TO ZERO ! PERMANENTLY) ! ! NXTRA -- NUMBER OF EXTRA POINTS TO BE ADDED TO ! THE SIMPLEX ! (NXTRA.GT.0 CAUSES A MORE THOROUGH ! SEARCH) ! ! KFLAG -- RETURNED .GT. ZERO FOR A NORMAL EXIT, ! RETURNED .LT. ZERO FOR AN ABNORMAL EXIT ! ! NOREP -- RETURNED .GT. ZERO IF THE FUNCTION WAS NOT ! REPRODUCIBLE ! ! KERFL -- NOT USED BY THIS ROUTINE ! ! KW -- THE LOGICAL UNIT NUMBER OF THE PRINTER ! ! --------------------------------------------------------------- ! ! THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME ! COMPILERS (WATFIV, FOR EXAMPLE) AND FORBIDDEN BY OTHERS ! (THE MODCOMP II, FOR EXAMPLE). ! EXTERNAL FUNK ! SIMPLEX SHOULD USUALLY BE RUN USING A FLOATING POINT ! PRECISION OF AT LEAST TEN SIGNIFICANT DIGITS. ON MOST ! COMPUTERS, THIS REQUIRES THE USE OF DOUBLE PRECISION. 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 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 DIMENSION Z(20,21),ZBAR(20),ZSTAR(20) ! 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 ! INTERNAL SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT,NVPLUS,NSSW,NF EQUIVALENCE (Z(1,1),ERR(1,1)) ! SIMPLEX CALLS NO FUNCTIONS, EITHER EXTERNAL OR INTRINSIC. ! THE SUBROUTINES CALLED ARE FUNK, SIBEG, SIFUN, AND DATSW. ! SIMPLEX TERMINATES IF SENSE SWITCH NUMBER -NSSW- IS ON. ! THE STATEMENT CALL DATSW(NSSW,JUMP) RETURNS JUMP=1 IF ! SENSE SWITCH NUMBER -NSSW- IS ON, AND JUMP=2 IF IT IS OFF. ! IF NO SENSE SWITCH IS TO BE USED, THE USER SHOULD SUPPLY ! A DUMMY SUBROUTINE FOR DATSW. ! THIS SUBROUTINE CONTAINS NO REAL OR DOUBLE PRECISION ! CONSTANTS. CALL SIBEG (FUNK) IF(KFLAG.LT.0) GO TO 720 FSAVE=FOBJ KERFL=0 ! SET FIXED QUANTITIES.... ! ! METHD=1 TO USE THE METHOD OF NELDER AND MEAD, ! METHD=2 TO USE A MODIFIED METHOD. ! METHD=1 CAN CAUSE THE BEST KNOWN POINT TO BE DISCARDED, ! AND HENCE IS NOT RECOMMENDED. ! METHD=2 IS RECOMMENDED. METHD=2 RZERO=0.0D0 UNITR=1.0D0 RTWO=2.0D0 JL=NVPLUS LATER=0 ! BEGIN THE NEXT ITERATION. ! FIND JH AND JL, THE INDICES OF THE POINTS WITH THE HIGHEST ! AND LOWEST FUNCTION VALUES, RESPECTIVELY. 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 ! CHECK FOR POSSIBLE TERMINATION. IF(KFLAG.NE.0) GO TO 690 IF(JH.EQ.JL .AND. NFLAT.GT.0) GO TO 660 ! CHECK FOR TIES -- MORE THAN ONE POINT IN THE SIMPLEX HAVING ! THE HIGHEST FUNCTION VALUE. JHTIE=0 DO 30 J=1,NVPLUS IF(J.NE.JH .AND. FZ(J).GE.FZ(JH)) JHTIE=7 30 CONTINUE IF(JH.NE.JL .AND. JHTIE.EQ.0) GO TO 70 ! THERE IS A TIE FOR HIGHEST FUNCTION VALUE. ! CHOOSE H FAR FROM L. DX=RZERO JHSAV=JH DO 50 J=1,NVPLUS IF(J.EQ.JL .OR. FZ(J).LT.FZ(JH)) GO TO 50 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 50 CONTINUE 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)) ! ADD EXTRA POINTS TO THE SIMPLEX, IF DESIRED. ! ANY EXTRA POINTS WILL BE SUPERIMPOSED ON THE HIGHEST POINT, ! P(JH), SO THAT THEY WILL SPLIT AS SOON AS POSSIBLE. 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 DO 120 J=JLOW,NVPLUS FZ(J)=FZ(JH) DO 110 K=1,NV Z(K,J)=Z(K,JH) 110 CONTINUE 120 CONTINUE 130 FSAVE=FZ(JL) ! CALCULATE PBAR, THE CENTROID OF ALL POINTS EXCEPT P(JH). ! THE MEAN VALUE IS COMPUTED USING A STABLE UPDATE FORMULA BY ! D. H. D. WEST, COMMUNICATIONS OF THE A.C.M. ! 22 (1979) 532-535. 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 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 ! CHECK THE ROUNDING. ! ZBARJ MUST LIE IN THE CLOSED INTERVAL (ZMIN,ZMAX). IF(ZBARJ.LE.ZMAX) GO TO 160 ZBARJ=ZMAX GO TO 170 160 IF(ZBARJ.GE.ZMIN) GO TO 190 ZBARJ=ZMIN 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) 190 ZBAR(J)=ZBARJ 200 CONTINUE IF(NTRACE.GE.2) WRITE(KW,210)(ZBAR(J),J=1,NV) 210 FORMAT(' ZBAR(J)....'/(1X,1PG24.16,2G24.16)) ! ATTEMPT A REFLECTION. FORM P* . ! THE FORMS USED BELOW FOR REFLECTION, EXPANSION, AND ! CONTRACTION ALL HAVE OPTIMAL ROUNDOFF PROPERTIES. JDIFF=0 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 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 250 CALL SIFUN (FUNK) FSTAR=FOBJ KONTR=0 IF(FSTAR.LT.FZ(JL)) GO TO 280 IF(NTRACE.GE.1) WRITE(KW,260)FSTAR,NF 260 FORMAT(' FOBJ =',1PG24.16,5X,'REFLECTION FAILED',8X,'NF =',I7) DO 270 J=1,NVPLUS ! HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER ! AND MEAD. THEY USE FSTAR.LE.FZ(J) . ! THAT IS NOT CONSISTENT WITH THE TEXT OF THE ARTICLE, AND CAN ! CREATE AN INFINITE LOOP OF REFLECTIONS. IF(J.NE.JH .AND. FSTAR.LT.FZ(J)) GO TO 340 270 CONTINUE ! HERE WE DEVIATE FROM THE FLOWCHART IN THE ARTICLE BY NELDER ! AND MEAD. THEY USE FSTAR.LE.FZ(JH) (SEE COMMENT ABOVE). KONTR=1 IF(FSTAR.LT.FZ(JH)) GO TO 340 GO TO 380 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)) ! THE REFLECTED VALUE FSTAR IS LESS THAN THE LOWEST VALUE, ! FZ(JL), IN THE SIMPLEX. ! ATTEMPT AN EXPANSION. FORM P** . 310 DO 320 J=1,NV IF(MASK(J).EQ.0) X(J)=X(J)+(GAMMA-UNITR)*(X(J)-ZBAR(J)) 320 CONTINUE CALL SIFUN (FUNK) ! CHOOSE ONE OF TWO STRATEGIES FOR ACCEPTING THE EXPANSION ! POINT P** . IF((METHD.EQ.1 .AND. FOBJ.LT.FZ(JL)) .OR. & (METHD.EQ.2 .AND. FOBJ.LT.FSTAR)) GO TO 490 IF(NTRACE.GE.1) WRITE(KW,330)FOBJ 330 FORMAT(' FOBJ =',1PG24.16,5X,'EXPANSION FAILED') ! THE REFLECTION FAILED BUT THE REFLECTED POINT WAS BETTER ! THAN SOME OTHER POINT, OR ELSE THE REFLECTION SUCCEEDED ! BUT THE EXPANSION FAILED. ! ACCEPT THE REFLECTED POINT. REPLACE P(JH) BY P* . 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) 360 DO 370 J=1,NV Z(J,JH)=ZSTAR(J) 370 CONTINUE FZ(JH)=FSTAR ! IF THE REFLECTION FAILED AND THE REFLECTED POINT WAS BETTER ! THAN ONLY P(JH), CONTRACT THE SIMPLEX. IF(KONTR.EQ.0) GO TO 530 ! ATTEMPT A CONTRACTION. FORM P** . 380 JDIFF=0 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) ! CHECK THE ROUNDING. ! XJ MIGHT HAVE BEEN ROUNDED UP TO ZJJH, WHICH COULD CREATE ! AN INFINITE LOOP. ! XJ MUST LIE IN THE CLOSED INTERVAL (ZJJH,ZBARJ), AND ! XJ MUST NOT BE EQUAL TO ZJJH UNLESS ZJJH.EQ.ZBARJ . ! EQUIVALENTLY, AN OPEN INTERVAL (ZJJH,ZBARJ) MUST EXIST AND ! XJ MUST LIE IN IT, OR XJ MUST BE EQUAL TO ZBARJ . 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 390 IF(XJ.NE.ZJJH) JDIFF=7 X(J)=XJ 400 CONTINUE 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 420 IF(NTRACE.GE.1) WRITE(KW,430)FOBJ 430 FORMAT(' FOBJ =',1PG24.16,5X,'CONTRACTION FAILED') 440 JS=JL ! REPLACE ALL P(J) BY (P(J)+P(JL))/2.0 . DO 470 J=1,NVPLUS IF(J.EQ.JL) GO TO 470 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 ! CHECK THE ROUNDING. ! XK MIGHT HAVE BEEN ROUNDED UP TO ZKJ, WHICH COULD CREATE ! AN INFINITE LOOP. ! XK MUST LIE IN THE CLOSED INTERVAL (ZKJ,ZKJL), AND ! XK MUST NOT BE EQUAL TO ZKJ UNLESS ZKJ.EQ.ZKJL . ! EQUIVALENTLY, AN OPEN INTERVAL (ZKJ,ZKJL) MUST EXIST AND ! XK MUST LIE IN IT, OR XK MUST BE EQUAL TO ZKJL . 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 450 Z(K,J)=XK X(K)=XK 460 CONTINUE CALL SIFUN (FUNK) IF(FOBJ.LT.FZ(JS)) JS=J FZ(J)=FOBJ 470 CONTINUE 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 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) ! REPLACE P(JH) BY P* OR P** . 510 DO 520 J=1,NV Z(J,JH)=X(J) 520 CONTINUE FZ(JH)=FOBJ ! TEST FOR CONVERGENCE. 530 DO 550 J=1,NV IF(MASK(J).NE.0) GO TO 550 ZMAX=Z(J,1) ZMIN=Z(J,1) DO 540 K=2,NVPLUS ZJK=Z(J,K) IF(ZJK.GT.ZMAX) ZMAX=ZJK IF(ZJK.LT.ZMIN) ZMIN=ZJK 540 CONTINUE DZ=ZMAX-ZMIN IF(DZ.GT.DELMIN(J)) GO TO 560 550 CONTINUE GO TO 640 ! RETURN IF THE SENSE SWITCH IS ON. 560 JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.LE.1) GO TO 570 IF(NF.GT.NFMAX) GO TO 590 GO TO 10 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 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.') 610 IF(NTRACE.LT.-1) GO TO 680 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 GO TO 680 640 KFLAG=1 IF(NTRACE.GE.0) WRITE(KW,650) 650 FORMAT(///' TERMINATED WHEN THE DIMENSIONS OF THE SIMPLEX BECAME AS SMALL AS THE DELMIN(J).') GO TO 680 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.') 680 CONTINUE ! GO BACK AND COMPUTE JL AND JH ONE LAST TIME. GO TO 10 ! ! FINISH UP AND EXIT. ! 690 DO 700 J=1,NV X(J)=Z(J,JL) 700 CONTINUE 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) 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)) WRITE(KW,*) ' ' WRITE(KW,*) ' ' ! THE FOLLOWING STATEMENT IS THE ONLY RETURN IN THIS ROUTINE. RETURN END SUBROUTINE SIBEG (FUNK) !-------------------------------------------------------------- ! SIBEG 1.3 OCTOBER 1991 ! ! A.N.S.I. STANDARD FORTRAN 77 ! ! COPYRIGHT (C) 1965, 1975, 1990 J. P. CHANDLER, ! COMPUTER SCIENCE DEPARTMENT, ! OKLAHOMA STATE UNIVERSITY ! ! SIBEG SETS DEFAULT VALUES AND PRINTS INITIAL OUTPUT ! FOR SIMPLEX. ! ! ------------------------------------------------------------ ! ! INPUT QUANTITIES..... FUNK,X(*),XMAX(*),XMIN(*),DELTX(*), ! DELMIN(*),NV,NTRACE,MASK(*), ! NFMAX,NFLAT,NXTRA,KW ! ! OUTPUT QUANTITIES.... FZ,HUGE,ALPHA,BETA,GAMMA,MAXPT, ! NVPLUS,NSSW,NF,KFLAG,NOREP, ! DELTX(*), AND DELMIN(*) !--------------------------------------------------------------- ! THE FOLLOWING EXTERNAL STATEMENT IS REQUIRED BY SOME ! COMPILERS (WATFIV, FOR EXAMPLE) AND FORBIDDEN BY OTHERS ! (THE MODCOMP II, FOR EXAMPLE). EXTERNAL FUNK DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, & Z,FZ,HUGE,ALPHA,BETA,GAMMA,DELDF,XS,XPLUS,FSAVE,RZERO INTEGER J,JUMP,JVARY,K,KERFL,KFLAG,KTYPE,KW,MASK, & MATRX,MAXPT,NACTIV,NF,NFLAT,NFMAX,NOREP,NSSW, & NTRACE,NV,NVMAX,NVPLUS,NXTRA DIMENSION Z(20,21) ! 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 ! INTERNAL SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT, & NVPLUS,NSSW,NF EQUIVALENCE (Z(1,1),ERR(1,1)) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ! SET FIXED QUANTITIES .... ! ! KTYPE = CONSOLE TYPEWRITER UNIT NUMBER ! (IRRELEVANT IF A DUMMY DATSW ROUTINE IS USED) ! KTYPE=1 ! ! NSSW = SENSE SWITCH NUMBER FOR SEARCH TERMINATION ! (IRRELEVANT IF A DUMMY DATSW ROUTINE IS USED) ! NSSW=6 ! ! HUGE = A VERY LARGE REAL NUMBER ! (DEFAULT VALUE FOR XMAX AND -XMIN, AND ! OUT-OF-BOUNDS VALUE FOR FOBJ) ! HUGE=1.0D35 ! ! NVMAX IS THE MAXIMUM VALUE OF NV. ! NVMAX IS ALSO THE DIMENSION OF X, XMAX, XMIN, DELTX, DELMIN, ! MASK, ZBAR, ZSTAR, AND THE FIRST DIMENSION OF ERR AND Z. ! NVMAX=20 ! ! MAXPT IS THE MAXIMUM NUMBER OF POINTS IN THE SIMPLEX. ! MAXPT IS ALSO THE DIMENSION OF FZ, AND THE SECOND DIMENSION ! OF ERR AND Z. ! MAXPT MUST BE .GE. (NVMAX+1). ! MAXPT=21 ! ! DELDF = DEFAULT VALUE FOR DELTX(J) ! DELDF=0.01D0 ! ! ALPHA = REFLECTION PARAMETER ! ALPHA=1.0D0 ! ! BETA = CONTRACTION PARAMETER ! BETA=0.5D0 ! ! GAMMA = EXPANSION PARAMETER ! GAMMA=2.0D0 RZERO=0.0D0 ! NO REAL OR DOUBLE PRECISION CONSTANTS ARE USED BEYOND THIS ! POINT, IN THIS SUBROUTINE. KFLAG=0 NOREP=0 KERFL=0 ! ! CHECK SOME INPUT QUANTITIES, AND SET THEM TO DEFAULT VALUES ! IF DESIRED. ! ! MAKE SURE THAT THE SENSE SWITCH IS OFF. ! JUMP=2 CALL DATSW (NSSW,JUMP) IF(JUMP.EQ.2) GO TO 30 ! ! THIS IS THE ONLY USAGE OF THE CONSOLE TYPEWRITER. ! WRITE(KTYPE,10) NSSW 10 FORMAT(/' TURN OFF SENSE SWITCH ',I2//' ') 20 CALL DATSW (NSSW,JUMP) IF(JUMP.EQ.1) GO TO 20 30 JVARY=0 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 50 DO 100 J=1,NV IF(MASK(J).NE.0) GO TO 100 IF(DELMIN(J).LT.RZERO) DELMIN(J)=-DELMIN(J) ! CHECK THAT DELTX(J) IS NOT NEGLIGIBLE. 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 60 IF(X(J).EQ.RZERO) GO TO 70 DELTX(J)=DELDF*X(J) GO TO 80 70 DELTX(J)=DELDF 80 IF(XMAX(J).GT.XMIN(J)) GO TO 90 XMAX(J)=HUGE XMIN(J)=-HUGE 90 IF(X(J).GT.XMAX(J)) X(J)=XMAX(J) IF(X(J).LT.XMIN(J)) X(J)=XMIN(J) 100 CONTINUE 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)) ! ! CALCULATE THE INITIAL P(I) AND Y(I) (Z AND FZ, ! RESPECTIVELY). ! ! NF = NUMBER OF CALLS TO FUNK SO FAR 180 NF=0 ! ! NACT = NUMBER OF ACTIVE X(J) ! NACTIV=0 DO 200 J=1,NV IF(MASK(J).NE.0) GO TO 200 NACTIV=NACTIV+1 DO 190 K=1,NV Z(K,NACTIV)=X(K) 190 CONTINUE 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 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) IF(NACTIV.GT.0) GO TO 230 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.') ! SET THE BASE POINT. 230 NVPLUS=NACTIV+1 DO 240 K=1,NV Z(K,NVPLUS)=X(K) 240 CONTINUE CALL SIFUN (FUNK) FZ(NVPLUS)=FOBJ ! ! CHECK THE REPRODUCIBILITY OF FOBJ, AND PRINT THE REMAINING ! LINES. ! 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) 260 IF(NTRACE.GE.0) WRITE(KW,270)FOBJ 270 FORMAT(///' FOBJ =',1PG24.16///' BEGIN MINIMIZATION....') RETURN END SUBROUTINE SIFUN (FUNK) !----------------------------------------------------- ! SIFUN 1.2 DECEMBER 1991 ! ! A.N.S.I. STANDARD FORTRAN 77 ! ! COPYRIGHT (C) 1975, 1990 BY J. P. CHANDLER, ! COMPUTER SCIENCE DEPARTMENT, ! OKLAHOMA STATE UNIVERSITY ! ! SIFUN CALLS FUNK TO COMPUTE FOBJ, IF X IS IN BOUNDS. ! SIFUN IS CALLED BY SIMPLEX AND SIBEG. !------------------------------------------------------ DOUBLE PRECISION X,XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ, & FZ,HUGE,ALPHA,BETA,GAMMA INTEGER JF,JVARY,KERFL,KFLAG,KW,MASK,MATRX,MAXPT, & NF,NFLAT,NFMAX,NOREP,NSSW,NTRACE,NV,NVPLUS,NXTRA ! 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 ! INTERNAL SIMPLEX COMMON..... COMMON /STICK/ FZ(21),HUGE,ALPHA,BETA,GAMMA,MAXPT,NVPLUS,NSSW,NF DO 10 JF=1,NV IF(MASK(JF).EQ.0 .AND. & (X(JF).GT.XMAX(JF) .OR. X(JF).LT.XMIN(JF))) & GO TO 20 10 CONTINUE CALL FUNK NF=NF+1 GO TO 30 20 FOBJ=HUGE 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)) RETURN END SUBROUTINE DATSW (NSSW,JUMP) !----------------------------------------------------------- ! DUMMY VERSION OF SUBROUTINE DATSW -- ALL SWITCHES OFF. ! ! J. P. CHANDLER, COMPUTER SCIENCE DEPARTMENT, ! OKLAHOMA STATE UNIVERSITY !----------------------------------------------------------- INTEGER NSSW,JUMP JUMP=2 RETURN END SUBROUTINE STSET !------------------------------------------------------------ ! STSET 3.2 DECEMBER 1991 ! ! STSET SETS SOME INPUT QUANTITIES TO DEFAULT VALUES, FOR ! SUBROUTINES STEPIT, SIMPLEX, MARQ, STP, MINF, OR KAUPE. ! ! J. P. CHANDLER, DEPARTMENT OF COMPUTER SCIENCE, ! OKLAHOMA STATE UNIVERSITY ! ! USAGE..... ! ! CALL STSET. ! THEN SET SOME INPUT QUANTITIES (NV, AT LEAST) AND RESET ANY ! OF THOSE SET IN STSET (BETTER VALUES OF X(J), ETC.) BEFORE ! CALLING STEPIT OR SIMPLEX OR MARQ, ETC. !------------------------------------------------------------- DOUBLE PRECISION XMAX,XMIN,DELTX,DELMIN,ERR,FOBJ,FLAMBD, & FNU,RELDIF,RELMIN,RZERO,HUGE, X INTEGER JVARY,JX,KALCP,KERFL,KFLAG,KORDIF,KW,LEQU,MASK, & MATRX,MAXIT,METHD,MAXSUB,MAXUPD,NFLAT,NFMAX,NOREP, & NTRACE,NV,NVMAX,NXTRA 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 COMMON /NLLS4/ FLAMBD,FNU,RELDIF,RELMIN,METHD,KALCP, & KORDIF,MAXIT,LEQU,MAXSUB,MAXUPD HUGE=1.0D30 RZERO=0.0D0 ! KW = LOGICAL UNIT NUMBER OF THE PRINTER KW=6 ! NVMAX IS THE MAXIMUM PERMISSIBLE VALUE OF NV, GIVEN THE ! PRESENT DIMENSIONS OF ARRAYS. ! NVMAX IS THE DIMENSION OF THE ARRAYS X(*), XMAX(*), XMIN(*), ! DELTX(*), DELMIN(*), AND MASK(*). NVMAX IS ALSO THE FIRST ! DIMENSION OF ERR(*,*). THE SECOND DIMENSION OF ERR(*,*) ! IS NVMAX+1. NVMAX=20 ! THE USER MUST SET NV AFTER CALLING STSET. NV=-1 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 KORDIF=1 RELDIF=1.0D-8 RELMIN=1.0D-6 ! FOR GREATER ACCURACY BUT LESS SPEED IN MARQ, SET... ! KORDIF=2 ! RELDIF=O.4D-5 ! RELMIN=1.0D-9 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 RETURN END !SAMPLE DATA (File Simptest.dat must contain:) !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 !end of file tsmplx.f90