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
85 common /lbfgstat/ status,niter,ncalls
88 c initialize some values to be used below
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)
102 c perform dynamic allocation of some global arrays
104 if (.not. allocated(scale)) allocate (scale(nvar))
106 c set default values for variable scale factors
108 if (.not. set_scale) then
110 if (scale(i) .eq. 0.0d0) scale(i) = 1.0d0
114 c set default parameters for the optimization
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
123 c set default parameters for the line search
125 if (stpmax .eq. 0.0d0) stpmax = 5.0d0
132 c search the keywords for optimization parameters
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
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
167 c print header information about the optimization method
169 if (jprint .gt. 0) then
170 if (msav .eq. 0) then
172 20 format (/,' Steepest Descent Gradient Optimization :')
174 30 format (/,' SD Iter F Value G RMS F Move',
175 & ' X Move Angle FG Call Comment',/)
178 40 format (/,' Limited Memory BFGS Quasi-Newton',
181 50 format (/,' QN Iter F Value G RMS F Move',
182 & ' X Move Angle FG Call Comment',/)
187 c perform dynamic allocation of some local arrays
189 allocate (x_old(nvar))
191 allocate (g_old(nvar))
196 if (msav .ne. 0) then
198 allocate (alpha(msav))
199 allocate (s(nvar,msav))
200 allocate (y(nvar,msav))
203 c evaluate the function and get the initial gradient
206 maxiter = niter + maxiter
215 g_norm = g_norm + g(i)*g(i)
216 g_rms = g_rms + (g(i)*scale(i))**2
218 g_norm = sqrt(g_norm)
219 g_rms = sqrt(g_rms) / rms
220 f_move = 0.5d0 * stpmax * g_norm
222 c print initial information prior to first iteration
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)
229 write (jout,70) niter,f,g_rms,ncalls
230 70 format (i6,d14.4,d11.4,29x,i7)
235 c write initial intermediate prior to first iteration
237 if (iwrite .gt. 0) call optsave (niter,f,x0)
239 c tests of the various termination criteria
241 if (niter .ge. maxiter) then
245 if (f .le. fctmin) then
249 if (g_rms .le. grdmin) then
254 c start of a new limited memory BFGS iteration
256 do while (.not. done)
258 c write (jout,*) "LBFGS niter",niter
259 muse = min(niter-1,msav)
261 if (m .gt. msav) m = 1
263 c estimate Hessian diagonal and compute the Hg product
272 if (k .eq. 0) k = msav
275 alpha(k) = alpha(k) + s(i,k)*q(i)
277 alpha(k) = alpha(k) * rho(k)
279 q(i) = q(i) - alpha(k)*y(i,k)
288 beta = beta + y(i,k)*r(i)
292 r(i) = r(i) + s(i,k)*(alpha(k)-beta)
295 if (k .gt. msav) k = 1
298 c set search direction and store current point and gradient
306 c perform line search along the new conjugate direction
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"
313 c update variables based on results of this iteration
315 if (msav .ne. 0) then
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)
328 c get the sizes of the moves made during this iteration
334 x_move = x_move + ((x0(i)-x_old(i))/scale(i))**2
336 x_move = sqrt(x_move) / rms
337 if (coordtype .eq. 'INTERNAL') then
338 x_move = radian * x_move
341 c compute the rms gradient per optimization parameter
345 g_rms = g_rms + (g(i)*scale(i))**2
347 g_rms = sqrt(g_rms) / rms
349 c test for error due to line search problems
351 if (status.eq.'BadIntpln' .or. status.eq.'IntplnErr') then
353 if (nerr .ge. maxerr) done = .true.
358 c test for too many total iterations
360 if (niter .ge. maxiter) then
365 c test the normal termination criteria
367 if (f .le. fctmin) then
371 if (g_rms .le. grdmin) then
376 c print intermediate results for the current iteration
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)
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)
395 c write intermediate results for the current iteration
397 if (iwrite .gt. 0) then
398 if (done .or. mod(niter,iwrite).eq.0) then
399 call optsave (niter,f,x0)
404 c perform deallocation of some local arrays
413 if (msav .ne. 0) then
420 c set final value of the objective function
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)
428 write (jout,110) status
429 110 format (/,' LBFGS -- Incomplete Convergence due to ',a9)