source/maxlik/src_MD_T_maxlik-NEWCORR.safe/amebsa.f

   1       SUBROUTINE amebsa(p,y,mp,np,ndim,pb,yb,ftol,funk,iter,temptr)
   2       INTEGER iter,mp,ndim,np,NMAX
   3       REAL ftol,temptr,yb,p(mp,np),pb(np),y(mp),funk
   4       PARAMETER (NMAX=200)
   5       EXTERNAL funk
   6 CU    USES amotsa,funk,ran1
   7       INTEGER i,idum,ihi,ilo,inhi,j,m,n
   8       REAL rtol,sum,swap,tt,yhi,ylo,ynhi,ysave,yt,ytry,psum(NMAX),
   9      *amotsa
  10       real ran1
  11       COMMON /ambsa/ tt,idum
  12 c      write (2,*) "amebsa: p and y"
  13 c      do i=1,n+1
  14 c        write (iout,'(20e15.5)') (p(i,j),j=1,n),y(i)
  15 c      enddo
  16 c      write (2,*) "amebsa: ftol",ftol," iter",iter," temptr",temptr
  17       tt=-temptr
  18 1     do 12 n=1,ndim
  19         sum=0.
  20         do 11 m=1,ndim+1
  21           sum=sum+p(m,n)
  22 11      continue
  23         psum(n)=sum
  24 12    continue
  25 2     ilo=1
  26       inhi=1
  27       ihi=2
  28       ylo=y(1)+tt*log(ran1(idum))
  29       ynhi=ylo
  30       yhi=y(2)+tt*log(ran1(idum))
  31       if (ylo.gt.yhi) then
  32         ihi=1
  33         inhi=2
  34         ilo=2
  35         ynhi=yhi
  36         yhi=ylo
  37         ylo=ynhi
  38       endif
  39       do 13 i=3,ndim+1
  40         yt=y(i)+tt*log(ran1(idum))
  41         if(yt.le.ylo) then
  42           ilo=i
  43           ylo=yt
  44         endif
  45         if(yt.gt.yhi) then
  46           inhi=ihi
  47           ynhi=yhi
  48           ihi=i
  49           yhi=yt
  50         else if(yt.gt.ynhi) then
  51           inhi=i
  52           ynhi=yt
  53         endif
  54 13    continue
  55       rtol=2.*abs(yhi-ylo)/(abs(yhi)+abs(ylo))
  56       if (rtol.lt.ftol.or.iter.lt.0) then
  57         swap=y(1)
  58         y(1)=y(ilo)
  59         y(ilo)=swap
  60         do 14 n=1,ndim
  61           swap=p(1,n)
  62           p(1,n)=p(ilo,n)
  63           p(ilo,n)=swap
  64 14      continue
  65         return
  66       endif
  67       iter=iter-2
  68       ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,-1.0)
  69       if (ytry.le.ylo) then
  70         ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,2.0)
  71       else if (ytry.ge.ynhi) then
  72         ysave=yhi
  73         ytry=amotsa(p,y,psum,mp,np,ndim,pb,yb,funk,ihi,yhi,0.5)
  74         if (ytry.ge.ysave) then
  75           do 16 i=1,ndim+1
  76             if(i.ne.ilo)then
  77               do 15 j=1,ndim
  78                 psum(j)=0.5*(p(i,j)+p(ilo,j))
  79                 p(i,j)=psum(j)
  80 15            continue
  81               y(i)=funk(psum)
  82             endif
  83 16        continue
  84           iter=iter-ndim
  85           goto 1
  86         endif
  87       else
  88         iter=iter+1
  89       endif
  90       goto 2
  91       END
  92 C  (C) Copr. 1986-92 Numerical Recipes Software 0!5,.