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 (iprint .lt. 0) iprint = 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
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
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
164 c print header information about the optimization method
166 if (iprint .gt. 0) then
167 if (msav .eq. 0) then
169 20 format (/,' Steepest Descent Gradient Optimization :')
171 30 format (/,' SD Iter F Value G RMS F Move',
172 & ' X Move Angle FG Call Comment',/)
175 40 format (/,' Limited Memory BFGS Quasi-Newton',
178 50 format (/,' QN Iter F Value G RMS F Move',
179 & ' X Move Angle FG Call Comment',/)
184 c perform dynamic allocation of some local arrays
186 allocate (x_old(nvar))
188 allocate (g_old(nvar))
193 if (msav .ne. 0) then
195 allocate (alpha(msav))
196 allocate (s(nvar,msav))
197 allocate (y(nvar,msav))
200 c evaluate the function and get the initial gradient
203 maxiter = niter + maxiter
212 g_norm = g_norm + g(i)*g(i)
213 g_rms = g_rms + (g(i)*scale(i))**2
215 g_norm = sqrt(g_norm)
216 g_rms = sqrt(g_rms) / rms
217 f_move = 0.5d0 * stpmax * g_norm
219 c print initial information prior to first iteration
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)
226 write (jout,70) niter,f,g_rms,ncalls
227 70 format (i6,d14.4,d11.4,29x,i7)
232 c write initial intermediate prior to first iteration
234 if (iwrite .gt. 0) call optsave (niter,f,x0)
236 c tests of the various termination criteria
238 if (niter .ge. maxiter) then
242 if (f .le. fctmin) then
246 if (g_rms .le. grdmin) then
251 c start of a new limited memory BFGS iteration
253 do while (.not. done)
255 muse = min(niter-1,msav)
257 if (m .gt. msav) m = 1
259 c estimate Hessian diagonal and compute the Hg product
268 if (k .eq. 0) k = msav
271 alpha(k) = alpha(k) + s(i,k)*q(i)
273 alpha(k) = alpha(k) * rho(k)
275 q(i) = q(i) - alpha(k)*y(i,k)
284 beta = beta + y(i,k)*r(i)
288 r(i) = r(i) + s(i,k)*(alpha(k)-beta)
291 if (k .gt. msav) k = 1
294 c set search direction and store current point and gradient
302 c perform line search along the new conjugate direction
305 call search (nvar,f,g,x0,p,f_move,angle,ncalls,fgvalue,status)
307 c update variables based on results of this iteration
309 if (msav .ne. 0) then
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)
322 c get the sizes of the moves made during this iteration
328 x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
330 x_move = sqrt(x_move) / rms
331 if (coordtype .eq. 'INTERNAL') then
332 x_move = radian * x_move
335 c compute the rms gradient per optimization parameter
339 g_rms = g_rms + (g(i)*scale(i))**2
341 g_rms = sqrt(g_rms) / rms
343 c test for error due to line search problems
345 if (status.eq.'BadIntpln' .or. status.eq.'IntplnErr') then
347 if (nerr .ge. maxerr) done = .true.
352 c test for too many total iterations
354 if (niter .ge. maxiter) then
359 c test the normal termination criteria
361 if (f .le. fctmin) then
365 if (g_rms .le. grdmin) then
370 c print intermediate results for the current iteration
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)
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)
389 c write intermediate results for the current iteration
391 if (iwrite .gt. 0) then
392 if (done .or. mod(niter,iwrite).eq.0) then
393 call optsave (niter,f,x0)
398 c perform deallocation of some local arrays
407 if (msav .ne. 0) then
414 c set final value of the objective function
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)
422 write (jout,110) status
423 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)