c  patrn2.f        Pattern search (Comm. ACM Algorithm 178, fixed up)
c
      subroutine patrn2(psi,delta,delmin,phi,n,Spsi,
     *   rho,konvrg,mxeval,neval,ntracp,lp)
c
c  Subroutine patrn2 searches for a local minimum of a smooth
c  function S(phi,n) of n floating point parameters phi(k),k=1,n.
c
c  The user provides a function subprogram to compute S given phi().
c
c  The user sets initial values for
c
c     n        = the number of parameters, phi()
c
c     psi()    = best available estimates for the parameters phi()
c
c     delta()  = initial step sizes
c
c     delmin() = absolute convergence tolerances for the step sizes
c
c     rho      = step size reduction factor (rho=0.2, say, or 0.1)
c
c     mxeval   = approximate maximum number of evaluations of
c                   the function S to be allowed.
c
c  phi() is a scratch vector of n floating point elements.
c
c  The user calls patrn2, which calls function S(phi,n) repeatedly
c  and finally returns to the calling program.  At this time psi()
c  contains the best point found, Spsi contains the function value
c  at psi(), and the last call to function S was made with psi() as the
c  argument.  delta() contains the last step sizes used, konvrg.eq.1
c  if convergence was achieved in neval function evaluations,
c  and neval returns the actual number of evaluations,
c  which could be as large as (mxeval+2*n), for mxeval.gt.0.
c
c  The phi() must not be modified inside of function S.
c
c  To hold any phi(k) fixed at its initial value, set delta(k) to zero.
c
c  To enforce constant constraints of inequality,
c     phimin(k) .le. phi(k) .le. phimax(k),
c  set the initial psi() within this feasible region and arrange for
c  function S to return an extremely large positive value of the
c  function whenever a constraint is violated.  (This feature is left
c  to the user to implement in detail, perhaps using arrays
c  phimin() and phimax() in common, setting the values in the main
c  program and accessing them in function S.)
c
c
      integer n,konvrg,mxeval,neval,ntracp,lp,     k,kallno,jbpok
c
      double precision psi(n),delta(n),delmin(n),phi(n),
     *   S,Spsi,rho,     Ss,theta,     dabs
c
c
      Spsi=S(psi,n)
      neval=1
c
      if(ntracp.ge.0) then
         write(lp,10)(psi(k),k=1,n)
   10    format(/' Begin pattern search (patrn2.f).'/
     *      ' psi=',1p,5g15.7:/(5x,5g15.7:))
         write(lp,20)Spsi
   20    format(5x,'Spsi =',1p,g15.7)
      endif
c
c  repeat ... until konvrg.eq.1 or neval.ge.mxeval
c
c  The best known point is (psi,Spsi).
c
   30 do 40 k=1,n
         phi(k)=psi(k)
   40 continue
      Ss=Spsi
c
c  Call procedure E to explore.
c
      kallno=1
      call E(phi,delta,n,Ss,neval,kallno,ntracp,lp)
      jbpok=1
c
c  The best known point is now (phi,Ss).
c
c  while Ss.lt.Spsi and jbpok.eq.1 and neval.lt.mxeval do
c
   50 if(Ss.ge.Spsi .or. jbpok.eq.0 .or. neval.ge.mxeval) go to 80
c
c  Make a pattern move.
c
      do 60 k=1,n
         theta=psi(k)
         psi(k)=phi(k)
         theta=phi(k)-theta
         phi(k)=psi(k)+theta
   60 continue
      Spsi=Ss
      Ss=S(phi,n)
      neval=neval+1
      kallno=2
      call E(phi,delta,n,Ss,neval,kallno,ntracp,lp)
c
c  Make the Bell-Pike test for appreciable net progress.
c
      jbpok=0
      do 70 k=1,n
         if(dabs(phi(k)-psi(k)) .gt. 0.5d0*dabs(delta(k))) jbpok=1
   70 continue
c
c  end while
c
      go to 50
c
   80 if(Ss.lt.Spsi) then
         do 90 k=1,n
            psi(k)=phi(k)
   90    continue
         Spsi=Ss
      endif
c
      konvrg=0
      if(kallno.eq.1 .or. neval.ge.mxeval) then
c
c  Test for convergence.
c
         konvrg=1
         do 100 k=1,n
            if(dabs(delta(k)) .gt. dabs(delmin(k))) konvrg=0
  100    continue
c
         if(konvrg.eq.0 .and. neval.lt.mxeval) then
c
c  Reduce the step sizes.
c
            do 110 k=1,n
               delta(k)=rho*delta(k)
  110       continue
c
            if(ntracp.ge.1) write(lp,120)rho,(delta(k),k=1,n)
  120       format(//' Step sizes reduced.   rho =',f5.3/
     *         ' delta=',1p,5g14.4:/(7x,5g14.4:))
         endif
      endif
c
c  until konvrg.eq.1 or neval.ge.mxeval
c
      if(konvrg.eq.0 .and. neval.lt.mxeval) go to 30
c
c  Evaluate S at the final point.
c
      Spsi=S(psi,n)
      neval=neval+1
c
      if(ntracp.ge.0) write(lp,130)Spsi,neval,konvrg,(psi(k),k=1,n)
  130 format(/' Spsi =',1p,g15.7,3x,'neval =',i11,3x,
     *   'Return from patrn2 with konvrg =',i2/
     *   ' psi=',5g15.7:/(5x,5g15.7:))
c
      return
c
c  End patrn2
c
      end
      subroutine E(phi,delta,n,Ss,neval,kallno,ntracp,lp)
c
      integer n,neval,kallno,ntracp,lp,     k
c
      double precision phi(n),delta(n),     Ss,
     *   phisav,S,Sphi
c
      do 10 k=1,n
         if(delta(k).ne.0.0d0) then
            phisav=phi(k)
            phi(k)=phisav+delta(k)
            Sphi=S(phi,n)
            neval=neval+1
            if(Sphi.lt.Ss) then
               Ss=Sphi
            else
               delta(k)=-delta(k)
               phi(k)=phisav+delta(k)
               Sphi=S(phi,n)
               neval=neval+1
               if(Sphi.lt.Ss) then
                  Ss=Sphi
               else
                  delta(k)=-delta(k)
                  phi(k)=phisav
               endif
            endif
         endif
   10 continue
c
      if(ntracp.ge.1) write(lp,20)Ss,neval,kallno,(phi(k),k=1,n)
   20 format(/' Ss =',1p,g15.7,6x,'neval =',i6,10x,
     *   'Return from proc E call no.',i2/
     *   ' phi=',5g15.7:/(5x,5g15.7:))
c
      return
c
c  End procedure E (patrn2)
c
      end