source/unres/src-HCD-5D/lbfgs.F

   1 c
   2 c
   3 c     ###################################################
   4 c     ##  COPYRIGHT (C)  1999  by  Jay William Ponder  ##
   5 c     ##              All Rights Reserved              ##
   6 c     ###################################################
   7 c
   8 c     ##############################################################
   9 c     ##                                                          ##
  10 c     ##  subroutine lbfgs  --  limited memory BFGS optimization  ##
  11 c     ##                                                          ##
  12 c     ##############################################################
  13 c
  14 c
  15 c     "lbfgs" is a limited memory BFGS quasi-newton nonlinear
  16 c     optimization routine
  17 c
  18 c     literature references:
  19 c
  20 c     J. Nocedal, "Updating Quasi-Newton Matrices with Limited
  21 c     Storage", Mathematics of Computation, 35, 773-782 (1980)
  22 c
  23 c     D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
  24 c     for Large Scale Optimization", Mathematical Programming,
  25 c     45, 503-528 (1989)
  26 c
  27 c     J. Nocedal and S. J. Wright, "Numerical Optimization",
  28 c     Springer-Verlag, New York, 1999, Section 9.1
  29 c
  30 c     variables and parameters:
  31 c
  32 c     nvar      number of parameters in the objective function
  33 c     x0        contains starting point upon input, upon return
  34 c                 contains the best point found
  35 c     minimum   during optimization contains best current function
  36 c                 value; returns final best function value
  37 c     grdmin    normal exit if rms gradient gets below this value
  38 c     ncalls    total number of function/gradient evaluations
  39 c
  40 c     required external routines:
  41 c
  42 c     fgvalue    function to evaluate function and gradient values
  43 c     optsave    subroutine to write out info about current status
  44 c
  45 c
  46       subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
  47       use inform
  48       use iounit
  49       use keys
  50       use linmin
  51       use math
  52       use minima
  53       use output
  54       use scales
  55       implicit none
  56       integer i,j,k,m
  57       integer nvar,next
  58       integer msav,muse
  59       integer niter,ncalls
  60       integer nerr,maxerr
  61       real*8 f,f_old,fgvalue
  62       real*8 f_move,x_move
  63       real*8 g_norm,g_rms
  64       real*8 minimum,grdmin
  65       real*8 angle,rms,beta
  66       real*8 ys,yy,gamma
  67       real*8 x0(*)
  68       real*8, allocatable :: rho(:)
  69       real*8, allocatable :: alpha(:)
  70       real*8, allocatable :: x_old(:)
  71       real*8, allocatable :: g(:)
  72       real*8, allocatable :: g_old(:)
  73       real*8, allocatable :: p(:)
  74       real*8, allocatable :: q(:)
  75       real*8, allocatable :: r(:)
  76       real*8, allocatable :: h0(:)
  77       real*8, allocatable :: s(:,:)
  78       real*8, allocatable :: y(:,:)
  79       logical done
  80       character*9 blank,status
  81       character*20 keyword
  82       character*240 record
  83       character*240 string
  84       external fgvalue,optsave
  85 c
  86 c
  87 c     initialize some values to be used below
  88 c
  89       ncalls = 0
  90       rms = sqrt(dble(nvar))
  91       if (coordtype .eq. 'CARTESIAN') then
  92          rms = rms / sqrt(3.0d0)
  93       else if (coordtype .eq. 'RIGIDBODY') then
  94          rms = rms / sqrt(6.0d0)
  95       end if
  96       blank = '         '
  97       done = .false.
  98       nerr = 0
  99       maxerr = 2
 100 c
 101 c     perform dynamic allocation of some global arrays
 102 c
 103       if (.not. allocated(scale))  allocate (scale(nvar))
 104 c
 105 c     set default values for variable scale factors
 106 c
 107       if (.not. set_scale) then
 108          do i = 1, nvar
 109             if (scale(i) .eq. 0.0d0)  scale(i) = 1.0d0
 110          end do
 111       end if
 112 c
 113 c     set default parameters for the optimization
 114 c
 115       msav = min(nvar,20)
 116       if (fctmin .eq. 0.0d0)  fctmin = -100000000.0d0
 117       if (maxiter .eq. 0)  maxiter = 1000000
 118       if (nextiter .eq. 0)  nextiter = 1
 119       if (jprint .lt. 0)  jprint = 1
 120       if (iwrite .lt. 0)  iwrite = 1
 121 c
 122 c     set default parameters for the line search
 123 c
 124       if (stpmax .eq. 0.0d0)  stpmax = 5.0d0
 125       stpmin = 1.0d-16
 126       cappa = 0.9d0
 127       slpmax = 10000.0d0
 128       angmax = 180.0d0
 129       intmax = 5
 130 c
 131 c     search the keywords for optimization parameters
 132 c
 133 #ifdef LBFGSREAD
 134       do i = 1, nkey
 135          next = 1
 136          record = keyline(i)
 137          call gettext (record,keyword,next)
 138          call upcase (keyword)
 139          string = record(next:240)
 140          if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
 141             read (string,*,err=10,end=10)  msav
 142             msav = max(0,min(msav,nvar))
 143          else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
 144             msav = 0
 145          else if (keyword(1:7) .eq. 'FCTMIN ') then
 146             read (string,*,err=10,end=10)  fctmin
 147          else if (keyword(1:8) .eq. 'MAXITER ') then
 148             read (string,*,err=10,end=10)  maxiter
 149          else if (keyword(1:8) .eq. 'STEPMAX ') then
 150             read (string,*,err=10,end=10)  stpmax
 151          else if (keyword(1:8) .eq. 'STEPMIN ') then
 152             read (string,*,err=10,end=10)  stpmin
 153          else if (keyword(1:6) .eq. 'CAPPA ') then
 154             read (string,*,err=10,end=10)  cappa
 155          else if (keyword(1:9) .eq. 'SLOPEMAX ') then
 156             read (string,*,err=10,end=10)  slpmax
 157          else if (keyword(1:7) .eq. 'ANGMAX ') then
 158             read (string,*,err=10,end=10)  angmax
 159          else if (keyword(1:7) .eq. 'INTMAX ') then
 160             read (string,*,err=10,end=10)  intmax
 161          end if
 162    10    continue
 163       end do
 164 #endif
 165 c
 166 c     print header information about the optimization method
 167 c
 168       if (jprint .gt. 0) then
 169          if (msav .eq. 0) then
 170             write (jout,20)
 171    20       format (/,' Steepest Descent Gradient Optimization :')
 172             write (jout,30)
 173    30       format (/,' SD Iter     F Value      G RMS      F Move',
 174      &                 '   X Move   Angle  FG Call  Comment',/)
 175          else
 176             write (jout,40)
 177    40       format (/,' Limited Memory BFGS Quasi-Newton',
 178      &                 ' Optimization :')
 179             write (jout,50)
 180    50       format (/,' QN Iter     F Value      G RMS      F Move',
 181      &                 '   X Move   Angle  FG Call  Comment',/)
 182          end if
 183          flush (jout)
 184       end if
 185 c
 186 c     perform dynamic allocation of some local arrays
 187 c
 188       allocate (x_old(nvar))
 189       allocate (g(nvar))
 190       allocate (g_old(nvar))
 191       allocate (p(nvar))
 192       allocate (q(nvar))
 193       allocate (r(nvar))
 194       allocate (h0(nvar))
 195       if (msav .ne. 0) then
 196          allocate (rho(msav))
 197          allocate (alpha(msav))
 198          allocate (s(nvar,msav))
 199          allocate (y(nvar,msav))
 200       end if
 201 c
 202 c     evaluate the function and get the initial gradient
 203 c
 204       niter = nextiter - 1
 205       maxiter = niter + maxiter
 206       ncalls = ncalls + 1
 207       f = fgvalue (x0,g)
 208       f_old = f
 209       m = 0
 210       gamma = 1.0d0
 211       g_norm = 0.0d0
 212       g_rms = 0.0d0
 213       do i = 1, nvar
 214          g_norm = g_norm + g(i)*g(i)
 215          g_rms = g_rms + (g(i)*scale(i))**2
 216       end do
 217       g_norm = sqrt(g_norm)
 218       g_rms = sqrt(g_rms) / rms
 219       f_move = 0.5d0 * stpmax * g_norm
 220 c
 221 c     print initial information prior to first iteration
 222 c
 223       if (jprint .gt. 0) then
 224          if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
 225             write (jout,60)  niter,f,g_rms,ncalls
 226    60       format (i6,f14.4,f11.4,29x,i7)
 227          else
 228             write (jout,70)  niter,f,g_rms,ncalls
 229    70       format (i6,d14.4,d11.4,29x,i7)
 230          end if
 231          flush (jout)
 232       end if
 233 c
 234 c     write initial intermediate prior to first iteration
 235 c
 236       if (iwrite .gt. 0)  call optsave (niter,f,x0)
 237 c
 238 c     tests of the various termination criteria
 239 c
 240       if (niter .ge. maxiter) then
 241          status = 'IterLimit'
 242          done = .true.
 243       end if
 244       if (f .le. fctmin) then
 245          status = 'SmallFct '
 246          done = .true.
 247       end if
 248       if (g_rms .le. grdmin) then
 249          status = 'SmallGrad'
 250          done = .true.
 251       end if
 252 c
 253 c     start of a new limited memory BFGS iteration
 254 c
 255       do while (.not. done)
 256          niter = niter + 1
 257          muse = min(niter-1,msav)
 258          m = m + 1
 259          if (m .gt. msav)  m = 1
 260 c
 261 c     estimate Hessian diagonal and compute the Hg product
 262 c
 263          do i = 1, nvar
 264             h0(i) = gamma
 265             q(i) = g(i)
 266          end do
 267          k = m
 268          do j = 1, muse
 269             k = k - 1
 270             if (k .eq. 0)  k = msav
 271             alpha(k) = 0.0d0
 272             do i = 1, nvar
 273                alpha(k) = alpha(k) + s(i,k)*q(i)
 274             end do
 275             alpha(k) = alpha(k) * rho(k)
 276             do i = 1, nvar
 277                q(i) = q(i) - alpha(k)*y(i,k)
 278             end do
 279          end do
 280          do i = 1, nvar
 281             r(i) = h0(i) * q(i)
 282          end do
 283          do j = 1, muse
 284             beta = 0.0d0
 285             do i = 1, nvar
 286                beta = beta + y(i,k)*r(i)
 287             end do
 288             beta = beta * rho(k)
 289             do i = 1, nvar
 290                r(i) = r(i) + s(i,k)*(alpha(k)-beta)
 291             end do
 292             k = k + 1
 293             if (k .gt. msav)  k = 1
 294          end do
 295 c
 296 c     set search direction and store current point and gradient
 297 c
 298          do i = 1, nvar
 299             p(i) = -r(i)
 300             x_old(i) = x0(i)
 301             g_old(i) = g(i)
 302          end do
 303 c
 304 c     perform line search along the new conjugate direction
 305 c
 306          status = blank
 307          call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,status)
 308 c
 309 c     update variables based on results of this iteration
 310 c
 311          if (msav .ne. 0) then
 312             ys = 0.0d0
 313             yy = 0.0d0
 314             do i = 1, nvar
 315                s(i,m) = x0(i) - x_old(i)
 316                y(i,m) = g(i) - g_old(i)
 317                ys = ys + y(i,m)*s(i,m)
 318                yy = yy + y(i,m)*y(i,m)
 319             end do
 320             gamma = abs(ys/yy)
 321             rho(m) = 1.0d0 / ys
 322          end if
 323 c
 324 c     get the sizes of the moves made during this iteration
 325 c
 326          f_move = f_old - f
 327          f_old = f
 328          x_move = 0.0d0
 329          do i = 1, nvar
 330             x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
 331          end do
 332          x_move = sqrt(x_move) / rms
 333          if (coordtype .eq. 'INTERNAL') then
 334             x_move = radian * x_move
 335          end if
 336 c
 337 c     compute the rms gradient per optimization parameter
 338 c
 339          g_rms = 0.0d0
 340          do i = 1, nvar
 341             g_rms = g_rms + (g(i)*scale(i))**2
 342          end do
 343          g_rms = sqrt(g_rms) / rms
 344 c
 345 c     test for error due to line search problems
 346 c
 347          if (status.eq.'BadIntpln' .or. status.eq.'IntplnErr') then
 348             nerr = nerr + 1
 349             if (nerr .ge. maxerr)  done = .true.
 350          else
 351             nerr = 0
 352          end if
 353 c
 354 c     test for too many total iterations
 355 c
 356          if (niter .ge. maxiter) then
 357             status = 'IterLimit'
 358             done = .true.
 359          end if
 360 c
 361 c     test the normal termination criteria
 362 c
 363          if (f .le. fctmin) then
 364             status = 'SmallFct '
 365             done = .true.
 366          end if
 367          if (g_rms .le. grdmin) then
 368             status = 'SmallGrad'
 369             done = .true.
 370          end if
 371 c
 372 c     print intermediate results for the current iteration
 373 c
 374          if (jprint .gt. 0) then
 375             if (done .or. mod(niter,jprint).eq.0) then
 376                if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.
 377      &             g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.
 378      &             f_move.gt.-1.0d5) then
 379                   write (jout,80)  niter,f,g_rms,f_move,x_move,
 380      &                             angle,ncalls,status
 381    80             format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
 382                else
 383                   write (jout,90)  niter,f,g_rms,f_move,x_move,
 384      &                             angle,ncalls,status
 385    90             format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
 386                end if
 387             end if
 388             flush (jout)
 389          end if
 390 c
 391 c     write intermediate results for the current iteration
 392 c
 393          if (iwrite .gt. 0) then
 394             if (done .or. mod(niter,iwrite).eq.0) then
 395                call optsave (niter,f,x0)
 396             end if
 397          end if
 398       end do
 399 c
 400 c     perform deallocation of some local arrays
 401 c
 402       deallocate (x_old)
 403       deallocate (g)
 404       deallocate (g_old)
 405       deallocate (p)
 406       deallocate (q)
 407       deallocate (r)
 408       deallocate (h0)
 409       if (msav .ne. 0) then
 410          deallocate (rho)
 411          deallocate (alpha)
 412          deallocate (s)
 413          deallocate (y)
 414       end if
 415 c
 416 c     set final value of the objective function
 417 c
 418       minimum = f
 419       if (jprint .gt. 0) then
 420          if (status.eq.'SmallGrad' .or. status.eq.'SmallFct ') then
 421             write (jout,100)  status
 422   100       format (/,' LBFGS  --  Normal Termination due to ',a9)
 423          else
 424             write (jout,110)  status
 425   110       format (/,' LBFGS  --  Incomplete Convergence due to ',a9)
 426          end if
 427          flush (jout)
 428       end if
 429       return
 430       end