source/unres/src-HCD-5D/search.f

   1 c
   2 c
   3 c     ###################################################
   4 c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
   5 c     ##              All Rights Reserved              ##
   6 c     ###################################################
   7 c
   8 c     #################################################################
   9 c     ##                                                             ##
  10 c     ##  subroutine search  --  perform unidimensional line search  ##
  11 c     ##                                                             ##
  12 c     #################################################################
  13 c
  14 c
  15 c     "search" is a unidimensional line search based upon parabolic
  16 c     extrapolation and cubic interpolation using both function and
  17 c     gradient values
  18 c
  19 c     variables used by the routine :
  20 c
  21 c     f       function value at the best line search point
  22 c     x       current values of variables during line search
  23 c     g       gradient at the current point during line search
  24 c     p       initial search vector, unchanged by this routine
  25 c     s       scaled search vector at current line search point
  26 c     angle   angle between search and negative gradient vector
  27 c
  28 c     parameters used by the routine :
  29 c
  30 c     stpmin   minimum step length in current line search direction
  31 c     stpmax   maximum step length in current line search direction
  32 c     cappa    stringency of line search (0=tight < cappa < 1=loose)
  33 c     slpmax   projected gradient above which stepsize is reduced
  34 c     angmax   maximum angle between search direction and -gradient
  35 c     intmax   maximum number of interpolations during line search
  36 c
  37 c     status codes upon return :
  38 c
  39 c     Success     normal termination after satisfying "cappa" test
  40 c     ScaleStep   normal termination after a step size rescaling
  41 c     ReSearch    normal termination after a reinterpolation
  42 c     WideAngle   large angle between search direction and -gradient
  43 c     BadIntpln   unsatisfied "cappa" test after two searches
  44 c     IntplnErr   function value increase or serious gradient error
  45 c
  46 c
  47       subroutine search (nvar,f,g,x,p,f_move,angle,ncalls,
  48      &                          fgvalue,status)
  49       use linmin
  50       use math
  51       implicit none
  52       integer i,nvar
  53       integer ncalls
  54       integer intpln
  55       real*8 fgvalue
  56       real*8 f,f_move
  57       real*8 s_norm,g_norm
  58       real*8 cosang,angle
  59       real*8 step,parab
  60       real*8 cube,cubstp
  61       real*8 sss,ttt
  62       real*8 f_0,f_1
  63       real*8 f_a,f_b,f_c
  64       real*8 sg_0,sg_1
  65       real*8 sg_a,sg_b,sg_c
  66       real*8 x(*)
  67       real*8 g(*)
  68       real*8 p(*)
  69       real*8, allocatable :: x_0(:)
  70       real*8, allocatable :: s(:)
  71       logical restart
  72       character*9 status
  73       character*9 blank
  74       external fgvalue
  75 c
  76 c
  77 c     use default parameters for the line search if needed
  78 c
  79       blank = '         '
  80       if (stpmin .eq. 0.0d0)  stpmin = 1.0d-16
  81       if (stpmax .eq. 0.0d0)  stpmax = 2.0d0
  82       if (cappa .eq. 0.0d0)  cappa = 0.1d0
  83       if (slpmax .eq. 0.0d0)  slpmax = 10000.0d0
  84       if (angmax .eq. 0.0d0)  angmax = 180.0d0
  85       if (intmax .eq. 0)  intmax = 5
  86 c
  87 c     perform dynamic allocation of some local arrays
  88 c
  89       allocate (x_0(nvar))
  90       allocate (s(nvar))
  91 c
  92 c     copy the search direction into a new vector
  93 c
  94       do i = 1, nvar
  95          s(i) = p(i)
  96       end do
  97 c
  98 c     compute the length of gradient and search direction
  99 c
 100       g_norm = 0.0d0
 101       s_norm = 0.0d0
 102       do i = 1, nvar
 103          g_norm = g_norm + g(i)*g(i)
 104          s_norm = s_norm + s(i)*s(i)
 105       end do
 106       g_norm = sqrt(g_norm)
 107       s_norm = sqrt(s_norm)
 108 c
 109 c     store initial function, then normalize the
 110 c     search vector and find projected gradient
 111 c
 112       f_0 = f
 113       sg_0 = 0.0d0
 114       do i = 1, nvar
 115          x_0(i) = x(i)
 116          s(i) = s(i) / s_norm
 117          sg_0 = sg_0 + s(i)*g(i)
 118       end do
 119 c
 120 c     check the angle between the search direction
 121 c     and the negative gradient vector
 122 c
 123       cosang = -sg_0 / g_norm
 124       cosang = min(1.0d0,max(-1.0d0,cosang))
 125       angle = radian * acos(cosang)
 126       if (angle .gt. angmax) then
 127          status = 'WideAngle'
 128          deallocate (x_0)
 129          deallocate (s)
 130          return
 131       end if
 132 c
 133 c     set the initial stepsize to the length of the passed
 134 c     search vector, or based on previous function decrease
 135 c
 136       step = 2.0d0 * abs(f_move/sg_0)
 137       step = min(step,s_norm)
 138       if (step .gt. stpmax)  step = stpmax
 139       if (step .lt. stpmin)  step = stpmin
 140 c
 141 c     beginning of the parabolic extrapolation procedure
 142 c
 143    10 continue
 144       restart = .true.
 145       intpln = 0
 146       f_b = f_0
 147       sg_b = sg_0
 148 c
 149 c     replace last point by latest and take another step
 150 c
 151    20 continue
 152       f_a = f_b
 153       sg_a = sg_b
 154       do i = 1, nvar
 155          x(i) = x(i) + step*s(i)
 156       end do
 157 c
 158 c     get new function and projected gradient following a step
 159 c
 160       ncalls = ncalls + 1
 161       f_b = fgvalue (x,g)
 162       sg_b = 0.0d0
 163       do i = 1, nvar
 164          sg_b = sg_b + s(i)*g(i)
 165       end do
 166 c
 167 c     scale stepsize if initial gradient change is too large
 168 c
 169       if (abs(sg_b/sg_a).ge.slpmax .and. restart) then
 170          do i = 1, nvar
 171             x(i) = x_0(i)
 172          end do
 173          step = step / 10.0d0
 174          status = 'ScaleStep'
 175          goto 10
 176       end if
 177       restart = .false.
 178 c
 179 c     return if the gradient is small and function decreases
 180 c
 181       if (abs(sg_b/sg_0).le.cappa .and. f_b.lt.f_a) then
 182          f = f_b
 183          if (status .eq. blank)  status = ' Success '
 184          deallocate (x_0)
 185          deallocate (s)
 186          return
 187       end if
 188 c
 189 c     interpolate if gradient changes sign or function increases
 190 c
 191       if (sg_b*sg_a.lt.0.0d0 .or. f_b.gt.f_a)  goto 30
 192 c
 193 c     if the finite difference curvature is negative double the step;
 194 c     or if (step < parabolic estimate < 4*step) use this estimate,
 195 c     otherwise truncate to step or 4*step, respectively
 196 c
 197       step = 2.0d0 * step
 198       if (sg_b .gt. sg_a) then
 199          parab = (f_a-f_b) / (sg_b-sg_a)
 200          if (parab .gt. 2.0d0*step)  parab = 2.0d0 * step
 201          if (parab .lt. 0.5d0*step)  parab = 0.5d0 * step
 202          step = parab
 203       end if
 204       if (step .gt. stpmax)  step = stpmax
 205       goto 20
 206 c
 207 c     beginning of the cubic interpolation procedure
 208 c
 209    30 continue
 210       intpln = intpln + 1
 211       sss = 3.0d0*(f_b-f_a)/step - sg_a - sg_b
 212       ttt = sss*sss - sg_a*sg_b
 213       if (ttt .lt. 0.0d0) then
 214          f = f_b
 215          status = 'IntplnErr'
 216          deallocate (x_0)
 217          deallocate (s)
 218          return
 219       end if
 220       ttt = sqrt(ttt)
 221       cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0d0*ttt)
 222       if (cube.lt.0.0d0 .or. cube.gt.step) then
 223          f = f_b
 224          status = 'IntplnErr'
 225          deallocate (x_0)
 226          deallocate (s)
 227          return
 228       end if
 229       do i = 1, nvar
 230          x(i) = x(i) - cube*s(i)
 231       end do
 232 c
 233 c     get new function and gradient, then test for termination
 234 c
 235       ncalls = ncalls + 1
 236       f_c = fgvalue (x,g)
 237       sg_c = 0.0d0
 238       do i = 1, nvar
 239          sg_c = sg_c + s(i)*g(i)
 240       end do
 241       if (abs(sg_c/sg_0) .le. cappa) then
 242          f = f_c
 243          if (status .eq. blank)  status = ' Success '
 244          deallocate (x_0)
 245          deallocate (s)
 246          return
 247       end if
 248 c
 249 c     get the next pair of bracketing points by replacing one
 250 c     of the current brackets with the interpolated point
 251 c
 252       if (f_c.le.f_a .or. f_c.le.f_b) then
 253          cubstp = min(abs(cube),abs(step-cube))
 254          if (cubstp.ge.stpmin .and. intpln.lt.intmax) then
 255 c
 256 c     if the current brackets have slopes of opposite sign,
 257 c     then substitute the interpolated point for the bracket
 258 c     point with slope of same sign as the interpolated point
 259 c
 260             if (sg_a*sg_b .lt. 0.0d0) then
 261                if (sg_a*sg_c .lt. 0.0d0) then
 262                   f_b = f_c
 263                   sg_b = sg_c
 264                   step = step - cube
 265                else
 266                   f_a = f_c
 267                   sg_a = sg_c
 268                   step = cube
 269                   do i = 1, nvar
 270                      x(i) = x(i) + cube*s(i)
 271                   end do
 272                end if
 273 c
 274 c     if current brackets have slope of same sign, then replace
 275 c     the far bracket if the interpolated point has a slope of
 276 c     the opposite sign or a lower function value than the near
 277 c     bracket, otherwise replace the near bracket point
 278 c
 279             else
 280                if (sg_a*sg_c.lt.0.0d0 .or. f_a.le.f_c) then
 281                   f_b = f_c
 282                   sg_b = sg_c
 283                   step = step - cube
 284                else
 285                   f_a = f_c
 286                   sg_a = sg_c
 287                   step = cube
 288                   do i = 1, nvar
 289                      x(i) = x(i) + cube*s(i)
 290                   end do
 291                end if
 292             end if
 293             goto 30
 294          end if
 295       end if
 296 c
 297 c     interpolation has failed, reset to best current point
 298 c
 299       f_1 = min(f_a,f_b,f_c)
 300       if (f_1 .eq. f_a) then
 301          sg_1 = sg_a
 302          do i = 1, nvar
 303             x(i) = x(i) + (cube-step)*s(i)
 304          end do
 305       else if (f_1 .eq. f_b) then
 306          sg_1 = sg_b
 307          do i = 1, nvar
 308             x(i) = x(i) + cube*s(i)
 309          end do
 310       else if (f_1 .eq. f_c) then
 311          sg_1 = sg_c
 312       end if
 313 c
 314 c     try to restart from best point with smaller stepsize
 315 c
 316       if (f_1 .gt. f_0) then
 317          ncalls = ncalls + 1
 318          f = fgvalue (x,g)
 319          status = 'IntplnErr'
 320          deallocate (x_0)
 321          deallocate (s)
 322          return
 323       end if
 324       f_0 = f_1
 325       sg_0 = sg_1
 326       if (sg_1 .gt. 0.0d0) then
 327          do i = 1, nvar
 328             s(i) = -s(i)
 329          end do
 330          sg_0 = -sg_1
 331       end if
 332       step = max(cube,step-cube) / 10.0d0
 333       if (step .lt. stpmin)  step = stpmin
 334 c
 335 c     if already restarted once, then return with best point
 336 c
 337       if (status .eq. ' ReSearch') then
 338          ncalls = ncalls + 1
 339          f = fgvalue (x,g)
 340          status = 'BadIntpln'
 341          deallocate (x_0)
 342          deallocate (s)
 343          return
 344       else
 345          status = ' ReSearch'
 346          goto 10
 347       end if
 348       end