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 c 3/14/2020 Adam Liwo: added the condition to prevent from infinite
 162 c       iteration loop
 163       if (ncalls.gt.200) then
 164          do i = 1, nvar
 165             x(i) = x_0(i)
 166          end do
 167          f_b=f_a
 168          deallocate (x_0)
 169          deallocate (s)
 170          return
 171       endif
 172       f_b = fgvalue (x,g)
 173 c      write (2,*) "ncalls",ncalls," f_a",f_a," f_b",f_b," sg_a",sg_a
 174       sg_b = 0.0d0
 175       do i = 1, nvar
 176          sg_b = sg_b + s(i)*g(i)
 177       end do
 178 c
 179 c     scale stepsize if initial gradient change is too large
 180 c
 181       if (abs(sg_b/sg_a).ge.slpmax .and. restart) then
 182          do i = 1, nvar
 183             x(i) = x_0(i)
 184          end do
 185          step = step / 10.0d0
 186          status = 'ScaleStep'
 187          goto 10
 188       end if
 189       restart = .false.
 190 c
 191 c     return if the gradient is small and function decreases
 192 c
 193       if (abs(sg_b/sg_0).le.cappa .and. f_b.lt.f_a) then
 194          f = f_b
 195          if (status .eq. blank)  status = ' Success '
 196          deallocate (x_0)
 197          deallocate (s)
 198          return
 199       end if
 200 c
 201 c     interpolate if gradient changes sign or function increases
 202 c
 203       if (sg_b*sg_a.lt.0.0d0 .or. f_b.gt.f_a)  goto 30
 204 c
 205 c     if the finite difference curvature is negative double the step;
 206 c     or if (step < parabolic estimate < 4*step) use this estimate,
 207 c     otherwise truncate to step or 4*step, respectively
 208 c
 209       step = 2.0d0 * step
 210       if (sg_b .gt. sg_a) then
 211          parab = (f_a-f_b) / (sg_b-sg_a)
 212          if (parab .gt. 2.0d0*step)  parab = 2.0d0 * step
 213          if (parab .lt. 0.5d0*step)  parab = 0.5d0 * step
 214          step = parab
 215       end if
 216       if (step .gt. stpmax)  step = stpmax
 217       goto 20
 218 c
 219 c     beginning of the cubic interpolation procedure
 220 c
 221    30 continue
 222       intpln = intpln + 1
 223       sss = 3.0d0*(f_b-f_a)/step - sg_a - sg_b
 224       ttt = sss*sss - sg_a*sg_b
 225       if (ttt .lt. 0.0d0) then
 226          f = f_b
 227          status = 'IntplnErr'
 228          deallocate (x_0)
 229          deallocate (s)
 230          return
 231       end if
 232       ttt = sqrt(ttt)
 233       cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0d0*ttt)
 234       if (cube.lt.0.0d0 .or. cube.gt.step) then
 235          f = f_b
 236          status = 'IntplnErr'
 237          deallocate (x_0)
 238          deallocate (s)
 239          return
 240       end if
 241       do i = 1, nvar
 242          x(i) = x(i) - cube*s(i)
 243       end do
 244 c
 245 c     get new function and gradient, then test for termination
 246 c
 247       ncalls = ncalls + 1
 248       f_c = fgvalue (x,g)
 249       sg_c = 0.0d0
 250       do i = 1, nvar
 251          sg_c = sg_c + s(i)*g(i)
 252       end do
 253       if (abs(sg_c/sg_0) .le. cappa) then
 254          f = f_c
 255          if (status .eq. blank)  status = ' Success '
 256          deallocate (x_0)
 257          deallocate (s)
 258          return
 259       end if
 260 c
 261 c     get the next pair of bracketing points by replacing one
 262 c     of the current brackets with the interpolated point
 263 c
 264       if (f_c.le.f_a .or. f_c.le.f_b) then
 265          cubstp = min(abs(cube),abs(step-cube))
 266          if (cubstp.ge.stpmin .and. intpln.lt.intmax) then
 267 c
 268 c     if the current brackets have slopes of opposite sign,
 269 c     then substitute the interpolated point for the bracket
 270 c     point with slope of same sign as the interpolated point
 271 c
 272             if (sg_a*sg_b .lt. 0.0d0) then
 273                if (sg_a*sg_c .lt. 0.0d0) then
 274                   f_b = f_c
 275                   sg_b = sg_c
 276                   step = step - cube
 277                else
 278                   f_a = f_c
 279                   sg_a = sg_c
 280                   step = cube
 281                   do i = 1, nvar
 282                      x(i) = x(i) + cube*s(i)
 283                   end do
 284                end if
 285 c
 286 c     if current brackets have slope of same sign, then replace
 287 c     the far bracket if the interpolated point has a slope of
 288 c     the opposite sign or a lower function value than the near
 289 c     bracket, otherwise replace the near bracket point
 290 c
 291             else
 292                if (sg_a*sg_c.lt.0.0d0 .or. f_a.le.f_c) then
 293                   f_b = f_c
 294                   sg_b = sg_c
 295                   step = step - cube
 296                else
 297                   f_a = f_c
 298                   sg_a = sg_c
 299                   step = cube
 300                   do i = 1, nvar
 301                      x(i) = x(i) + cube*s(i)
 302                   end do
 303                end if
 304             end if
 305             goto 30
 306          end if
 307       end if
 308 c
 309 c     interpolation has failed, reset to best current point
 310 c
 311       f_1 = min(f_a,f_b,f_c)
 312       if (f_1 .eq. f_a) then
 313          sg_1 = sg_a
 314          do i = 1, nvar
 315             x(i) = x(i) + (cube-step)*s(i)
 316          end do
 317       else if (f_1 .eq. f_b) then
 318          sg_1 = sg_b
 319          do i = 1, nvar
 320             x(i) = x(i) + cube*s(i)
 321          end do
 322       else if (f_1 .eq. f_c) then
 323          sg_1 = sg_c
 324       end if
 325 c
 326 c     try to restart from best point with smaller stepsize
 327 c
 328       if (f_1 .gt. f_0) then
 329          ncalls = ncalls + 1
 330          f = fgvalue (x,g)
 331          status = 'IntplnErr'
 332          deallocate (x_0)
 333          deallocate (s)
 334          return
 335       end if
 336       f_0 = f_1
 337       sg_0 = sg_1
 338       if (sg_1 .gt. 0.0d0) then
 339          do i = 1, nvar
 340             s(i) = -s(i)
 341          end do
 342          sg_0 = -sg_1
 343       end if
 344       step = max(cube,step-cube) / 10.0d0
 345       if (step .lt. stpmin)  step = stpmin
 346 c
 347 c     if already restarted once, then return with best point
 348 c
 349       if (status .eq. ' ReSearch') then
 350          ncalls = ncalls + 1
 351          f = fgvalue (x,g)
 352          status = 'BadIntpln'
 353          deallocate (x_0)
 354          deallocate (s)
 355          return
 356       else
 357          status = ' ReSearch'
 358          goto 10
 359       end if
 360       end