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