c---------------------------------------------------------------------- ccc AMOEBA c---------------------------------------------------------------------- ccc Multidimensional minimization routine of the function supplied by ccc user in FUNC(ndim,x,fofx), where x is an ndim-dimensional ccc vector and fofx is the functional value at x. That is, f(x). ccc ccc Input is a matrix simplx whose ndim+1 rows are ndim-dimensional ccc vectors which are the vertices of the starting simplex. Dimensions ccc are simplx(ndim+1,ndim). Also input is the vector fvertx of length ccc ndim+1, whose components must be pre-initialized to the values of ccc FUNC evaluated at the ndim+1 vertices (rows) of simplx. ftol is ccc the fractional convergence tolerance to be achieved in the ccc function value. Note: on output, simplx and fvertx will have been ccc reset to ndim+1 new points all within ftol of the minimum function ccc value. lambda is the simplex side length parameter, indicative of ccc the "size" of the search area. A counter icalls is used to keep ccc track of the number of function evaluations. ccc User must supply subroutine FUNC(ndim,x,fofx) and routines for ccc initializing the simplex on the first call. ccc Other definitions: ccc nmax: the maximum number of degrees of freedom. This can ccc be changed and the code recompiled to tackle bigger problems. ccc Because this is an inefficient method, nmax > 40 can be slow. ccc alpha,beta,gamma: size of the reflections, contractions and ccc expansions ccc FUNC - subroutine supplied by user to evaluate the function to ccc be minimized. It is declared as follows: ccc subroutine func(ndim,x,fofx) ccc double precision x(ndim),fofx ccc ftol - convergence tolerance. Recommended to be 1.0e-05 or so. ccc itmax - maximum number of iterations until the job is killed. ccc upon output, record x, the vector containing the coordinates of ccc the minimum and fofx, the value at the minimum. ccc Method uses the downhill simplex method of Nelder and Mead ccc (Computer Journal, 7, 308 (1965). Code contains elements from ccc Numerical Recipes by Press et al. and code written by Gary Sowers ccc (Univ. of Deleware). Modified for ChEg 698D by E. Maginn, Univ. of ccc Notre Dame, Tue Jan 25, 2000 10:33PM.. ccc Please report bugs to ed@nd.edu. c---------------------------------------------------------------------- subroutine amoeba(simplx,fvertx,ndim,ftol,itmax,lambda,x > ,fofx) implicit none c---------------------------------------------------------------------- ccc Parameters - must agree with parameters in user-supplied routines! integer nmax parameter (nmax=20) double precision alpha,beta,gamma parameter(alpha=1.0,beta=0.5,gamma=2.0) integer icalls common /savethis/ icalls integer i,j,nplus1 integer ndim,itmax,ilow,ihigh,inxthi double precision simplx(nmax+1,nmax) double precision fvertx(nmax+1),preflc(nmax),pexpan(nmax) double precision pcontr(nmax),pcentr(nmax) double precision x(nmax),fofx,rtol double precision ftol,lambda,fofx,fexpan,freflc,fcontr c---------------------------------------------------------------------- nplus1 = ndim+1 100 continue ! begin minimization ccc Determine the highest, next highest and lowest vertex points. ccc Initialize... ilow = 1 if (fvertx(1) .gt. fvertx(2)) then ihigh=1 inxthi=2 else ihigh=2 inxthi=1 endif ccc loop over all points of the simplex, and sort according to highest ccc next highest and lowest points. do i=1,nplus1 if (fvertx(i) .lt. fvertx(ilow) ) ilow=i if ( fvertx(i) .gt. fvertx(ihigh) ) then inxthi=ihigh ihigh=i elseif (fvertx(i) .gt. fvertx(inxthi)) then if (i .ne. ihigh) inxthi=i endif enddo c---------------------------------------------------------------------- ccc compute the fractional range from highest to lowest rtol = 2.0*ABS(fvertx(ihigh)-fvertx(ilow))/ > (ABS(fvertx(ihigh))+ABS(fvertx(ilow))) ccc Check to see if we continue if (rtol .lt. ftol) then write(*,*) 'Converged normally' goto 998 endif if (icalls .ge. itmax) then write(*,*) 'Exceeded maximum number of calls' goto 998 endif c---------------------------------------------------------------------- ccc Now, with the best and worst points picked out, we calculate the ccc centroid of the simplex. Through this point, the simplex will be ccc expanded and contracted and it "slithers" around searching for ccc a minimum. This is done by exploring along the ray from the high ccc point through the center. ccc zero all counters do j=1,ndim pcentr(j) = 0.0d0 enddo ccc Here, we compute the vector average of all points except the ccc highest. That is, the center of the "face" across from the highest ccc point. do i=1,nplus1 if (i .ne. ihigh) then do j=1,ndim pcentr(j) = pcentr(j) + simplx(i,j) enddo endif enddo do j=1,ndim pcentr(j) = pcentr(j) / (ndim) ccc Reflect the worst point through the centroid of the simplex by ccc an amount alpha. preflc(j) = (1.0d0 + alpha)*pcentr(j) - alpha*simplx(ihigh,j) enddo c---------------------------------------------------------------------- ccc Evaluate the function at the reflected point CALL func(ndim,preflc,freflc) ccc See how it compares ot the best previous point if (freflc .le. fvertx(ilow)) then ccc Gives a better point than the previous lowest. So try ccc additional extrapolation by gamma do j=1,ndim pexpan(j) = gamma*preflc(j) + (1.0d0-gamma)*pcentr(j) enddo ccc evaluate the function at the expanded point CALL func(ndim,pexpan,fexpan) c---------------------------------------------------------------------- ccc if the expansion was successful, replace the previous worst ccc vertex with the expanded point if (fexpan .lt. fvertx(ilow)) then do j=1,ndim simplx(ihigh,j) = pexpan(j) enddo fvertx(ihigh) = fexpan else ccc otherwise, the expansion by gamma didn't improve ccc things. Still use the reflected point, since it was better ccc than the previous best vertex do j=1,ndim simplx(ihigh,j) = preflc(j) enddo fvertx(ihigh) = freflc endif c---------------------------------------------------------------------- ccc else if the reflected point is actually worse than the 2nd ccc worst vertex... elseif (freflc .ge. fvertx(inxthi)) then ccc check to see if it is at least better than the worst point, ccc and replace the worst point with the reflected point if this ccc is the case. c---------------------------------------------------------------------- if (freflc .lt. fvertx(ihigh)) then ccc The reflected point is better than worst point, so replace ccc worst point. do j=1,ndim simplx(ihigh,j) = preflc(j) enddo fvertx(ihigh) = freflc endif c---------------------------------------------------------------------- ccc Then contract the simplex looking for an intermediate lower ccc point do j=1,ndim pcontr(j) = beta*simplx(ihigh,j) + (1.0d0-beta)*pcentr(j) enddo ccc evaluate function at contracted point call func(ndim,pcontr,fcontr) c---------------------------------------------------------------------- ccc check to see if contraction is an improvement over worst point if (fcontr .lt. fvertx(ihigh)) then ccc it is an improvement, so assign it as the new worst point do j=1,ndim simplx(ihigh,j) = pcontr(j) enddo fvertx(ihigh) = fcontr else ccc else the contraction is worse than all current vertices. In ccc this case, contract around the best vertex. That is, move ccc the entire simplex in the direction of the best point. do i=1,nplus1 if (i .ne. ilow) then do j=1,ndim pcontr(j) = 0.5d0*(simplx(i,j) + simplx(ilow,j)) simplx(i,j) = pcontr(j) enddo CALL func(ndim,pcontr,fvertx(i)) endif enddo endif c---------------------------------------------------------------------- else ccc arrive here if the reflected point gives a function value ccc intermediate between the 2nd worst and the best. Call the ccc reflected point worst and move on. do j=1,ndim simplx(ihigh,j) = preflc(j) enddo fvertx(ihigh) = freflc endif ! if (freflc .le. fvertx(ilow)) then ccc End of the while loop. goto 100 998 continue ccc assign the best parameters and the function value at the best ccc vertex so they can be returned to calling routine ccc minimum function value fofx = fvertx(ilow) ccc coordinates of minimum do j=1,ndim x(j) = simplx(ilow,j) enddo return end