C SUBROUTINE SMPLX(FUNC,XVLO,XVEC,XVHI,NVEC) C C Subroutine SMPLX: Simplex-method minimization, using C slight modification of Press et al routines C C Update history: C V1.0 Original version 15-Oct-1992 C IMPLICIT NONE C C Arguments: C FUNC = R*4 function = quantity to minimize C XVLO = R*4(1) = array of lower limits to XVEC C XVEC = R*4(1) = array of quantities to find at minimum C XVHI = R*4(1) = array of upper limits to XVEC C NVEC = I*4 = number of parameters C REAL*4 FUNC EXTERNAL FUNC C REAL*4 XVLO(1),XVEC(1),XVHI(1) INTEGER*4 NVEC C C Subprogrammes: C GETOP = subroutine = get values of controlling variables C EXTERNAL GETOP C C Local parameter: C NVMAX = I*4 = maximum length of vectors C INTEGER*4 NVMAX PARAMETER (NVMAX=128) C C Local variables: C ALPHA = R*4 = algorithm convergence factor C BETA = R*4 = algorithm convergence factor C EPSSM = R*4 = fractional accuracy demanded C FMIN = R*4 = minimum of function C GAMMA = R*4 = algorithm convergence factor C I1 = I*4 = counting index C I2 = I*4 = counting index C ITER = I*4 = number of iterations used C NIMAX = I*4 = maximum number of iterations C PTT(i) = R*4(NVMAX) = previous parameters C XINI(i) = R*4(NVMAX) = initial search direction C XIT(i) = R*4(NVMAX) = amount of shift C XMAT(i,j) = R*4(NVMAX,NVMAX) = matrix used by POWELL C REAL*4 ALPHA,BETA,EPSSM,FMIN,GAMMA,RTOL,YPR,YPRR REAL*4 FVEC(NVMAX) REAL*4 XMAT(NVMAX,NVMAX) INTEGER*4 I1,I2,ITER,MPTS,IHI,ILO,INHI,NIMAX REAL*4 PR(NVMAX),PRR(NVMAX),PBAR(NVMAX) C C Calculate size of matrix needed for minimization C MPTS=NVEC+1 C C Check that size isn't too large C IF (MPTS .GT. NVMAX) THEN WRITE (6,699) RETURN ENDIF C C Get controlling variables C CALL GETOP('SMPLX.ALPHA',ALPHA) CALL GETOP('SMPLX.BETA' ,BETA) CALL GETOP('SMPLX.GAMMA',GAMMA) CALL GETOP('SMPLX.EPS' ,EPSSM) CALL GETOP('SMPLX.NIMAX',RTOL) NIMAX=IFIX(RTOL+0.5) C C Initialise matrix of starting vertices of simplex: the initial C vector of parameters and the unit vectors in the other variables C DO I2=1,NVEC XMAT(1,I2)=XVEC(I2) ENDDO DO I1=2,NVEC+1 DO I2=1,NVEC XMAT(I1,I2)=0.0 ENDDO XMAT(I1,I1-1)=1.0 ENDDO C C Initialise vector of starting values of function on vertices C of simplex C DO I1=1,NVEC+1 DO I2=1,NVEC PRR(I2)=AMAX0(XVLO(I2),AMIN0(XMAT(I1,I2),XVHI(I2))) ENDDO FVEC(I1)=FUNC(PRR,NVEC) ENDDO C C Print value at start C WRITE (6,999) FVEC(1) C C Recoding of the AMOEBA routine from Press et al: C main loop covered once per iteration C ITER=0 1 CONTINUE C C Determine vectors to highest and lowest values of function, C as marked by indices ILO and IHI C ILO=1 IF (FVEC(1) .GT. FVEC(2)) THEN IHI =1 INHI=2 ELSE IHI =2 INHI=1 ENDIF DO I1=1,MPTS IF (FVEC(I1) .LT. FVEC(ILO)) THEN ILO=I1 ENDIF IF (FVEC(I1) .GT. FVEC(IHI)) THEN INHI=IHI IHI =I1 ELSEIF (FVEC(I1) .GT. FVEC(INHI)) THEN IF (I1 .NE. IHI) THEN INHI=I1 ENDIF ENDIF ENDDO C C Copy best result so far to vector XVEC C DO I2=1,NVEC XVEC(I2)=AMAX0(XVLO(I2),AMIN0(XMAT(ILO,I2),XVHI(I2))) ENDDO FMIN=FVEC(ILO) C C Compute fractional range from highest to lowest point of simplex. C RTOL=2.0*ABS(FVEC(IHI)-FVEC(ILO))/ + (ABS(FVEC(IHI))+ABS(FVEC(ILO))) C C If the range is small enough, return C IF (RTOL .LT. EPSSM) THEN WRITE (6,998) FMIN,ITER RETURN ENDIF C C Check number of iterations used; if at maximum return C IF (ITER .EQ. NIMAX) THEN WRITE (6,997) FMIN,ITER RETURN ENDIF C C Increment iteration counter C ITER=ITER+1 C C Find vector average of all points except the worst, PBAR, and the C extrapolation by a factor ALPHA from the high point through C the mean (a reflection about the mean), PR C DO I2=1,NVEC PBAR(I2)=0.0 ENDDO DO I1=1,MPTS IF (I1 .NE. IHI) THEN DO I2=1,NVEC PBAR(I2)=PBAR(I2)+XMAT(I1,I2) ENDDO ENDIF ENDDO DO I2=1,NVEC PBAR(I2)=PBAR(I2)/FLOAT(NVEC) PBAR(I2)=AMAX0(XVLO(I2),AMIN0(PBAR(I2),XVHI(I2))) C PR(I2) =(1.0+ALPHA)*PBAR(I2)-ALPHA*XMAT(IHI,I2) PR(I2) =AMAX0(XVLO(I2),AMIN0(PR(I2),XVHI(I2))) ENDDO C C Find the function value at the reflected point C YPR=FUNC(PR,NVEC) C C If this has improved matters to make PR the best simplex C vertex, try a further extrapolation by an amount GAMMA C to make the vector PRR C IF (YPR .LE. FVEC(ILO)) THEN DO I2=1,NVEC PRR(I2)=GAMMA*PR(I2)+(1.0-GAMMA)*PBAR(I2) PRR(I2)=AMAX0(XVLO(I2),AMIN0(PRR(I2),XVHI(I2))) ENDDO YPRR=FUNC(PRR,NVEC) C C If the additional extrapolation worked, use it to replace the C worst point in the simplex. Otherwise, replace the high C point by the PR result. C IF (YPRR .LT. FVEC(ILO)) THEN DO I2=1,NVEC XMAT(IHI,I2)=PRR(I2) ENDDO FVEC(IHI)=YPRR WRITE (6,991) FMIN,ITER ELSE DO I2=1,NVEC XMAT(IHI,I2)=PR(I2) ENDDO FVEC(IHI)=YPR WRITE (6,992) FMIN,ITER ENDIF C C If the reflection did not improve matters, look to see if the C reflected point is better than the second highest C ELSEIF (YPR .GE. FVEC(INHI)) THEN C C Reflected point is worse than second highest. Replace the highest C point if it's better than the highest C IF (YPR .LT. FVEC(IHI)) THEN DO I2=1,NVEC XMAT(IHI,I2)=PR(I2) ENDDO FVEC(IHI)=YPR WRITE (6,993) FMIN,ITER ENDIF C C Look for an intermediate better point, using it if it C is good enough C DO I2=1,NVEC PRR(I2)=BETA*XMAT(IHI,I2)+(1.0-BETA)*PBAR(I2) PRR(I2)=AMAX0(XVLO(I2),AMIN0(PRR(I2),XVHI(I2))) ENDDO YPRR=FUNC(PRR,NVEC) C IF (YPRR .LT. FVEC(IHI)) THEN DO I2=1,NVEC XMAT(IHI,I2)=PRR(I2) ENDDO FVEC(IHI)=YPRR WRITE (6,996) FMIN,ITER ELSE C C If highest point cannot be improved on in this way, contract C simplex about lowest point by a factor 0.5, recalculating C all function values C DO I1=1,MPTS IF (I1 .NE. ILO) THEN DO I2=1,NVEC PR(I2) =0.5*(XMAT(I1,I2)+XMAT(ILO,I2)) PR(I2) =AMAX0(XVLO(I2),AMIN0(PR(I2),XVHI(I2))) XMAT(I1,I2)=PR(I2) ENDDO FVEC(I1)=FUNC(PR,NVEC) ENDIF ENDDO WRITE (6,995) FMIN,ITER C ENDIF C C Original reflection gives a middling point. Replace old high C point and continue C ELSE DO I2=1,NVEC XMAT(IHI,I2)=PR(I2) ENDDO FVEC(IHI)=YPR WRITE (6,994) FMIN,ITER ENDIF C C End of main loop. Go back and try again. C GOTO 1 C C Error messages C 699 FORMAT (' '/ + ' *** Error: parameter vector too large ***'/) C C Debug progress messages C 991 FORMAT (' SMPLX: (2-extrapolate) at value = ', + 1PE15.5,' at iteration = ',I4) 992 FORMAT (' SMPLX: (extrapolated) at value = ', + 1PE15.5,' at iteration = ',I4) 993 FORMAT (' SMPLX: (nearly worst) at value = ', + 1PE15.5,' at iteration = ',I4) 994 FORMAT (' SMPLX: (mid reflect) at value = ', + 1PE15.5,' at iteration = ',I4) 995 FORMAT (' SMPLX: (contraction) at value = ', + 1PE15.5,' at iteration = ',I4) 996 FORMAT (' SMPLX: (intermediate) at value = ', + 1PE15.5,' at iteration = ',I4) 997 FORMAT (' SMPLX: iteration limit at value = ', + 1PE15.5,' at iteration = ',I4) 998 FORMAT (' SMPLX: converged to minimum = ', + 1PE15.5,' at iteration = ',I4) 999 FORMAT (' '/ + ' SMPLX: initial function value = ', + 1PE15.5) C END