3 c ###################################################
4 c ## COPYRIGHT (C) 1999 by Jay William Ponder ##
5 c ## All Rights Reserved ##
6 c ###################################################
8 c ##############################################################
10 c ## subroutine lbfgs -- limited memory BFGS optimization ##
12 c ##############################################################
15 c "lbfgs" is a limited memory BFGS quasi-newton nonlinear
16 c optimization routine
18 c literature references:
20 c J. Nocedal, "Updating Quasi-Newton Matrices with Limited
21 c Storage", Mathematics of Computation, 35, 773-782 (1980)
23 c D. C. Lui and J. Nocedal, "On the Limited Memory BFGS Method
24 c for Large Scale Optimization", Mathematical Programming,
27 c J. Nocedal and S. J. Wright, "Numerical Optimization",
28 c Springer-Verlag, New York, 1999, Section 9.1
30 c variables and parameters:
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
40 c required external routines:
42 c fgvalue function to evaluate function and gradient values
43 c optsave subroutine to write out info about current status
46 subroutine lbfgs (nvar,x0,minimum,grdmin,fgvalue,optsave)
61 real*8 f,f_old,fgvalue
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(:,:)
80 character*9 blank,status
84 external fgvalue,optsave
87 c initialize some values to be used below
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)
101 c perform dynamic allocation of some global arrays
103 if (.not. allocated(scale)) allocate (scale(nvar))
105 c set default values for variable scale factors
107 if (.not. set_scale) then
109 if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0
113 c set default parameters for the optimization
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
122 c set default parameters for the line search
124 if (stpmax .eq. 0.0d0) stpmax = 5.0d0
131 c search the keywords for optimization parameters
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
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
166 c print header information about the optimization method
168 if (jprint .gt. 0) then
169 if (msav .eq. 0) then
171 20 format (/,' Steepest Descent Gradient Optimization :')
173 30 format (/,' SD Iter F Value G RMS F Move',
174 & ' X Move Angle FG Call Comment',/)
177 40 format (/,' Limited Memory BFGS Quasi-Newton',
180 50 format (/,' QN Iter F Value G RMS F Move',
181 & ' X Move Angle FG Call Comment',/)
186 c perform dynamic allocation of some local arrays
188 allocate (x_old(nvar))
190 allocate (g_old(nvar))
195 if (msav .ne. 0) then
197 allocate (alpha(msav))
198 allocate (s(nvar,msav))
199 allocate (y(nvar,msav))
202 c evaluate the function and get the initial gradient
205 maxiter = niter + maxiter
214 g_norm = g_norm + g(i)*g(i)
215 g_rms = g_rms + (g(i)*scale(i))**2
217 g_norm = sqrt(g_norm)
218 g_rms = sqrt(g_rms) / rms
219 f_move = 0.5d0 * stpmax * g_norm
221 c print initial information prior to first iteration
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)
228 write (jout,70) niter,f,g_rms,ncalls
229 70 format (i6,d14.4,d11.4,29x,i7)
234 c write initial intermediate prior to first iteration
236 if (iwrite .gt. 0) call optsave (niter,f,x0)
238 c tests of the various termination criteria
240 if (niter .ge. maxiter) then
244 if (f .le. fctmin) then
248 if (g_rms .le. grdmin) then
253 c start of a new limited memory BFGS iteration
255 do while (.not. done)
257 muse = min(niter-1,msav)
259 if (m .gt. msav) m = 1
261 c estimate Hessian diagonal and compute the Hg product
270 if (k .eq. 0) k = msav
273 alpha(k) = alpha(k) + s(i,k)*q(i)
275 alpha(k) = alpha(k) * rho(k)
277 q(i) = q(i) - alpha(k)*y(i,k)
286 beta = beta + y(i,k)*r(i)
290 r(i) = r(i) + s(i,k)*(alpha(k)-beta)
293 if (k .gt. msav) k = 1
296 c set search direction and store current point and gradient
304 c perform line search along the new conjugate direction
307 call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,status)
309 c update variables based on results of this iteration
311 if (msav .ne. 0) then
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)
324 c get the sizes of the moves made during this iteration
330 x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
332 x_move = sqrt(x_move) / rms
333 if (coordtype .eq. 'INTERNAL') then
334 x_move = radian * x_move
337 c compute the rms gradient per optimization parameter
341 g_rms = g_rms + (g(i)*scale(i))**2
343 g_rms = sqrt(g_rms) / rms
345 c test for error due to line search problems
347 if (status.eq.'BadIntpln' .or. status.eq.'IntplnErr') then
349 if (nerr .ge. maxerr) done = .true.
354 c test for too many total iterations
356 if (niter .ge. maxiter) then
361 c test the normal termination criteria
363 if (f .le. fctmin) then
367 if (g_rms .le. grdmin) then
372 c print intermediate results for the current iteration
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)
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)
391 c write intermediate results for the current iteration
393 if (iwrite .gt. 0) then
394 if (done .or. mod(niter,iwrite).eq.0) then
395 call optsave (niter,f,x0)
400 c perform deallocation of some local arrays
409 if (msav .ne. 0) then
416 c set final value of the objective function
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)
424 write (jout,110) status
425 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)