source/unres/Tinker-minimize/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 (iprint .lt. 0)  iprint = 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       do i = 1, nkey
 134          next = 1
 135          record = keyline(i)
 136          call gettext (record,keyword,next)
 137          call upcase (keyword)
 138          string = record(next:240)
 139          if (keyword(1:14) .eq. 'LBFGS-VECTORS ') then
 140             read (string,*,err=10,end=10)  msav
 141             msav = max(0,min(msav,nvar))
 142          else if (keyword(1:17) .eq. 'STEEPEST-DESCENT ') then
 143             msav = 0
 144          else if (keyword(1:7) .eq. 'FCTMIN ') then
 145             read (string,*,err=10,end=10)  fctmin
 146          else if (keyword(1:8) .eq. 'MAXITER ') then
 147             read (string,*,err=10,end=10)  maxiter
 148          else if (keyword(1:8) .eq. 'STEPMAX ') then
 149             read (string,*,err=10,end=10)  stpmax
 150          else if (keyword(1:8) .eq. 'STEPMIN ') then
 151             read (string,*,err=10,end=10)  stpmin
 152          else if (keyword(1:6) .eq. 'CAPPA ') then
 153             read (string,*,err=10,end=10)  cappa
 154          else if (keyword(1:9) .eq. 'SLOPEMAX ') then
 155             read (string,*,err=10,end=10)  slpmax
 156          else if (keyword(1:7) .eq. 'ANGMAX ') then
 157             read (string,*,err=10,end=10)  angmax
 158          else if (keyword(1:7) .eq. 'INTMAX ') then
 159             read (string,*,err=10,end=10)  intmax
 160          end if
 161    10    continue
 162       end do
 163 c
 164 c     print header information about the optimization method
 165 c
 166       if (iprint .gt. 0) then
 167          if (msav .eq. 0) then
 168             write (jout,20)
 169    20       format (/,' Steepest Descent Gradient Optimization :')
 170             write (jout,30)
 171    30       format (/,' SD Iter     F Value      G RMS      F Move',
 172      &                 '   X Move   Angle  FG Call  Comment',/)
 173          else
 174             write (jout,40)
 175    40       format (/,' Limited Memory BFGS Quasi-Newton',
 176      &                 ' Optimization :')
 177             write (jout,50)
 178    50       format (/,' QN Iter     F Value      G RMS      F Move',
 179      &                 '   X Move   Angle  FG Call  Comment',/)
 180          end if
 181          flush (jout)
 182       end if
 183 c
 184 c     perform dynamic allocation of some local arrays
 185 c
 186       allocate (x_old(nvar))
 187       allocate (g(nvar))
 188       allocate (g_old(nvar))
 189       allocate (p(nvar))
 190       allocate (q(nvar))
 191       allocate (r(nvar))
 192       allocate (h0(nvar))
 193       if (msav .ne. 0) then
 194          allocate (rho(msav))
 195          allocate (alpha(msav))
 196          allocate (s(nvar,msav))
 197          allocate (y(nvar,msav))
 198       end if
 199 c
 200 c     evaluate the function and get the initial gradient
 201 c
 202       niter = nextiter - 1
 203       maxiter = niter + maxiter
 204       ncalls = ncalls + 1
 205       f = fgvalue (x0,g)
 206       f_old = f
 207       m = 0
 208       gamma = 1.0d0
 209       g_norm = 0.0d0
 210       g_rms = 0.0d0
 211       do i = 1, nvar
 212          g_norm = g_norm + g(i)*g(i)
 213          g_rms = g_rms + (g(i)*scale(i))**2
 214       end do
 215       g_norm = sqrt(g_norm)
 216       g_rms = sqrt(g_rms) / rms
 217       f_move = 0.5d0 * stpmax * g_norm
 218 c
 219 c     print initial information prior to first iteration
 220 c
 221       if (iprint .gt. 0) then
 222          if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and. g_rms.lt.1.0d5) then
 223             write (jout,60)  niter,f,g_rms,ncalls
 224    60       format (i6,f14.4,f11.4,29x,i7)
 225          else
 226             write (jout,70)  niter,f,g_rms,ncalls
 227    70       format (i6,d14.4,d11.4,29x,i7)
 228          end if
 229          flush (jout)
 230       end if
 231 c
 232 c     write initial intermediate prior to first iteration
 233 c
 234       if (iwrite .gt. 0)  call optsave (niter,f,x0)
 235 c
 236 c     tests of the various termination criteria
 237 c
 238       if (niter .ge. maxiter) then
 239          status = 'IterLimit'
 240          done = .true.
 241       end if
 242       if (f .le. fctmin) then
 243          status = 'SmallFct '
 244          done = .true.
 245       end if
 246       if (g_rms .le. grdmin) then
 247          status = 'SmallGrad'
 248          done = .true.
 249       end if
 250 c
 251 c     start of a new limited memory BFGS iteration
 252 c
 253       do while (.not. done)
 254          niter = niter + 1
 255          muse = min(niter-1,msav)
 256          m = m + 1
 257          if (m .gt. msav)  m = 1
 258 c
 259 c     estimate Hessian diagonal and compute the Hg product
 260 c
 261          do i = 1, nvar
 262             h0(i) = gamma
 263             q(i) = g(i)
 264          end do
 265          k = m
 266          do j = 1, muse
 267             k = k - 1
 268             if (k .eq. 0)  k = msav
 269             alpha(k) = 0.0d0
 270             do i = 1, nvar
 271                alpha(k) = alpha(k) + s(i,k)*q(i)
 272             end do
 273             alpha(k) = alpha(k) * rho(k)
 274             do i = 1, nvar
 275                q(i) = q(i) - alpha(k)*y(i,k)
 276             end do
 277          end do
 278          do i = 1, nvar
 279             r(i) = h0(i) * q(i)
 280          end do
 281          do j = 1, muse
 282             beta = 0.0d0
 283             do i = 1, nvar
 284                beta = beta + y(i,k)*r(i)
 285             end do
 286             beta = beta * rho(k)
 287             do i = 1, nvar
 288                r(i) = r(i) + s(i,k)*(alpha(k)-beta)
 289             end do
 290             k = k + 1
 291             if (k .gt. msav)  k = 1
 292          end do
 293 c
 294 c     set search direction and store current point and gradient
 295 c
 296          do i = 1, nvar
 297             p(i) = -r(i)
 298             x_old(i) = x0(i)
 299             g_old(i) = g(i)
 300          end do
 301 c
 302 c     perform line search along the new conjugate direction
 303 c
 304          status = blank
 305          call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,status)
 306 c
 307 c     update variables based on results of this iteration
 308 c
 309          if (msav .ne. 0) then
 310             ys = 0.0d0
 311             yy = 0.0d0
 312             do i = 1, nvar
 313                s(i,m) = x0(i) - x_old(i)
 314                y(i,m) = g(i) - g_old(i)
 315                ys = ys + y(i,m)*s(i,m)
 316                yy = yy + y(i,m)*y(i,m)
 317             end do
 318             gamma = abs(ys/yy)
 319             rho(m) = 1.0d0 / ys
 320          end if
 321 c
 322 c     get the sizes of the moves made during this iteration
 323 c
 324          f_move = f_old - f
 325          f_old = f
 326          x_move = 0.0d0
 327          do i = 1, nvar
 328             x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
 329          end do
 330          x_move = sqrt(x_move) / rms
 331          if (coordtype .eq. 'INTERNAL') then
 332             x_move = radian * x_move
 333          end if
 334 c
 335 c     compute the rms gradient per optimization parameter
 336 c
 337          g_rms = 0.0d0
 338          do i = 1, nvar
 339             g_rms = g_rms + (g(i)*scale(i))**2
 340          end do
 341          g_rms = sqrt(g_rms) / rms
 342 c
 343 c     test for error due to line search problems
 344 c
 345          if (status.eq.'BadIntpln' .or. status.eq.'IntplnErr') then
 346             nerr = nerr + 1
 347             if (nerr .ge. maxerr)  done = .true.
 348          else
 349             nerr = 0
 350          end if
 351 c
 352 c     test for too many total iterations
 353 c
 354          if (niter .ge. maxiter) then
 355             status = 'IterLimit'
 356             done = .true.
 357          end if
 358 c
 359 c     test the normal termination criteria
 360 c
 361          if (f .le. fctmin) then
 362             status = 'SmallFct '
 363             done = .true.
 364          end if
 365          if (g_rms .le. grdmin) then
 366             status = 'SmallGrad'
 367             done = .true.
 368          end if
 369 c
 370 c     print intermediate results for the current iteration
 371 c
 372          if (iprint .gt. 0) then
 373             if (done .or. mod(niter,iprint).eq.0) then
 374                if (f.lt.1.0d8 .and. f.gt.-1.0d7 .and.
 375      &             g_rms.lt.1.0d5 .and. f_move.lt.1.0d6 .and.
 376      &             f_move.gt.-1.0d5) then
 377                   write (jout,80)  niter,f,g_rms,f_move,x_move,
 378      &                             angle,ncalls,status
 379    80             format (i6,f14.4,f11.4,f12.4,f9.4,f8.2,i7,3x,a9)
 380                else
 381                   write (jout,90)  niter,f,g_rms,f_move,x_move,
 382      &                             angle,ncalls,status
 383    90             format (i6,d14.4,d11.4,d12.4,f9.4,f8.2,i7,3x,a9)
 384                end if
 385             end if
 386             flush (jout)
 387          end if
 388 c
 389 c     write intermediate results for the current iteration
 390 c
 391          if (iwrite .gt. 0) then
 392             if (done .or. mod(niter,iwrite).eq.0) then
 393                call optsave (niter,f,x0)
 394             end if
 395          end if
 396       end do
 397 c
 398 c     perform deallocation of some local arrays
 399 c
 400       deallocate (x_old)
 401       deallocate (g)
 402       deallocate (g_old)
 403       deallocate (p)
 404       deallocate (q)
 405       deallocate (r)
 406       deallocate (h0)
 407       if (msav .ne. 0) then
 408          deallocate (rho)
 409          deallocate (alpha)
 410          deallocate (s)
 411          deallocate (y)
 412       end if
 413 c
 414 c     set final value of the objective function
 415 c
 416       minimum = f
 417       if (iprint .gt. 0) then
 418          if (status.eq.'SmallGrad' .or. status.eq.'SmallFct ') then
 419             write (jout,100)  status
 420   100       format (/,' LBFGS  --  Normal Termination due to ',a9)
 421          else
 422             write (jout,110)  status
 423   110       format (/,' LBFGS  --  Incomplete Convergence due to ',a9)
 424          end if
 425          flush (jout)
 426       end if
 427       return
 428       end