L-BFGS-B  3.0
Large-scale Bound-constrained Optimization
driver1.f

This simple driver demonstrates how to call the L-BFGS-B code to solve a sample problem (the extended Rosenbrock function subject to bounds on the variables). The dimension n of this problem is variable. (Fortran-77 version)

1 c> \file driver1.f
2 
3 c
4 c L-BFGS-B is released under the "New BSD License" (aka "Modified BSD License"
5 c or "3-clause license")
6 c Please read attached file License.txt
7 c
8 c DRIVER 1 in Fortran 77
9 c --------------------------------------------------------------
10 c SIMPLE DRIVER FOR L-BFGS-B (version 3.0)
11 c --------------------------------------------------------------
12 c
13 c L-BFGS-B is a code for solving large nonlinear optimization
14 c problems with simple bounds on the variables.
15 c
16 c The code can also be used for unconstrained problems and is
17 c as efficient for these problems as the earlier limited memory
18 c code L-BFGS.
19 c
20 c This is the simplest driver in the package. It uses all the
21 c default settings of the code.
22 c
23 c
24 c References:
25 c
26 c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
27 c memory algorithm for bound constrained optimization'',
28 c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
29 c
30 c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
31 c Subroutines for Large Scale Bound Constrained Optimization''
32 c Tech. Report, NAM-11, EECS Department, Northwestern University,
33 c 1994.
34 c
35 c
36 c (Postscript files of these papers are available via anonymous
37 c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
38 c
39 c * * *
40 c
41 c March 2011 (latest revision)
42 c Optimization Center at Northwestern University
43 c Instituto Tecnologico Autonomo de Mexico
44 c
45 c Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778:
46 c L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained
47 c Optimization" (2011). To appear in ACM Transactions on
48 c Mathematical Software,
49 
50 c --------------------------------------------------------------
51 c DESCRIPTION OF THE VARIABLES IN L-BFGS-B
52 c --------------------------------------------------------------
53 c
54 c n is an INTEGER variable that must be set by the user to the
55 c number of variables. It is not altered by the routine.
56 c
57 c m is an INTEGER variable that must be set by the user to the
58 c number of corrections used in the limited memory matrix.
59 c It is not altered by the routine. Values of m < 3 are
60 c not recommended, and large values of m can result in excessive
61 c computing time. The range 3 <= m <= 20 is recommended.
62 c
63 c x is a DOUBLE PRECISION array of length n. On initial entry
64 c it must be set by the user to the values of the initial
65 c estimate of the solution vector. Upon successful exit, it
66 c contains the values of the variables at the best point
67 c found (usually an approximate solution).
68 c
69 c l is a DOUBLE PRECISION array of length n that must be set by
70 c the user to the values of the lower bounds on the variables. If
71 c the i-th variable has no lower bound, l(i) need not be defined.
72 c
73 c u is a DOUBLE PRECISION array of length n that must be set by
74 c the user to the values of the upper bounds on the variables. If
75 c the i-th variable has no upper bound, u(i) need not be defined.
76 c
77 c nbd is an INTEGER array of dimension n that must be set by the
78 c user to the type of bounds imposed on the variables:
79 c nbd(i)=0 if x(i) is unbounded,
80 c 1 if x(i) has only a lower bound,
81 c 2 if x(i) has both lower and upper bounds,
82 c 3 if x(i) has only an upper bound.
83 c
84 c f is a DOUBLE PRECISION variable. If the routine setulb returns
85 c with task(1:2)= 'FG', then f must be set by the user to
86 c contain the value of the function at the point x.
87 c
88 c g is a DOUBLE PRECISION array of length n. If the routine setulb
89 c returns with taskb(1:2)= 'FG', then g must be set by the user to
90 c contain the components of the gradient at the point x.
91 c
92 c factr is a DOUBLE PRECISION variable that must be set by the user.
93 c It is a tolerance in the termination test for the algorithm.
94 c The iteration will stop when
95 c
96 c (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
97 c
98 c where epsmch is the machine precision which is automatically
99 c generated by the code. Typical values for factr on a computer
100 c with 15 digits of accuracy in double precision are:
101 c factr=1.d+12 for low accuracy;
102 c 1.d+7 for moderate accuracy;
103 c 1.d+1 for extremely high accuracy.
104 c The user can suppress this termination test by setting factr=0.
105 c
106 c pgtol is a double precision variable.
107 c On entry pgtol >= 0 is specified by the user. The iteration
108 c will stop when
109 c
110 c max{|proj g_i | i = 1, ..., n} <= pgtol
111 c
112 c where pg_i is the ith component of the projected gradient.
113 c The user can suppress this termination test by setting pgtol=0.
114 c
115 c wa is a DOUBLE PRECISION array of length
116 c (2mmax + 5)nmax + 11mmax^2 + 8mmax used as workspace.
117 c This array must not be altered by the user.
118 c
119 c iwa is an INTEGER array of length 3nmax used as
120 c workspace. This array must not be altered by the user.
121 c
122 c task is a CHARACTER string of length 60.
123 c On first entry, it must be set to 'START'.
124 c On a return with task(1:2)='FG', the user must evaluate the
125 c function f and gradient g at the returned value of x.
126 c On a return with task(1:5)='NEW_X', an iteration of the
127 c algorithm has concluded, and f and g contain f(x) and g(x)
128 c respectively. The user can decide whether to continue or stop
129 c the iteration.
130 c When
131 c task(1:4)='CONV', the termination test in L-BFGS-B has been
132 c satisfied;
133 c task(1:4)='ABNO', the routine has terminated abnormally
134 c without being able to satisfy the termination conditions,
135 c x contains the best approximation found,
136 c f and g contain f(x) and g(x) respectively;
137 c task(1:5)='ERROR', the routine has detected an error in the
138 c input parameters;
139 c On exit with task = 'CONV', 'ABNO' or 'ERROR', the variable task
140 c contains additional information that the user can print.
141 c This array should not be altered unless the user wants to
142 c stop the run for some reason. See driver2 or driver3
143 c for a detailed explanation on how to stop the run
144 c by assigning task(1:4)='STOP' in the driver.
145 c
146 c iprint is an INTEGER variable that must be set by the user.
147 c It controls the frequency and type of output generated:
148 c iprint<0 no output is generated;
149 c iprint=0 print only one line at the last iteration;
150 c 0<iprint<99 print also f and |proj g| every iprint iterations;
151 c iprint=99 print details of every iteration except n-vectors;
152 c iprint=100 print also the changes of active set and final x;
153 c iprint>100 print details of every iteration including x and g;
154 c When iprint > 0, the file iterate.dat will be created to
155 c summarize the iteration.
156 c
157 c csave is a CHARACTER working array of length 60.
158 c
159 c lsave is a LOGICAL working array of dimension 4.
160 c On exit with task = 'NEW_X', the following information is
161 c available:
162 c lsave(1) = .true. the initial x did not satisfy the bounds;
163 c lsave(2) = .true. the problem contains bounds;
164 c lsave(3) = .true. each variable has upper and lower bounds.
165 c
166 c isave is an INTEGER working array of dimension 44.
167 c On exit with task = 'NEW_X', it contains information that
168 c the user may want to access:
169 c isave(30) = the current iteration number;
170 c isave(34) = the total number of function and gradient
171 c evaluations;
172 c isave(36) = the number of function value or gradient
173 c evaluations in the current iteration;
174 c isave(38) = the number of free variables in the current
175 c iteration;
176 c isave(39) = the number of active constraints at the current
177 c iteration;
178 c
179 c see the subroutine setulb.f for a description of other
180 c information contained in isave
181 c
182 c dsave is a DOUBLE PRECISION working array of dimension 29.
183 c On exit with task = 'NEW_X', it contains information that
184 c the user may want to access:
185 c dsave(2) = the value of f at the previous iteration;
186 c dsave(5) = the machine precision epsmch generated by the code;
187 c dsave(13) = the infinity norm of the projected gradient;
188 c
189 c see the subroutine setulb.f for a description of other
190 c information contained in dsave
191 c
192 c --------------------------------------------------------------
193 c END OF THE DESCRIPTION OF THE VARIABLES IN L-BFGS-B
194 c --------------------------------------------------------------
195 
196  program driver
197  use json_writer, only: json_write_aggregate
198 
199 c This simple driver demonstrates how to call the L-BFGS-B code to
200 c solve a sample problem (the extended Rosenbrock function
201 c subject to bounds on the variables). The dimension n of this
202 c problem is variable.
203 
204  integer nmax, mmax
205  parameter(nmax=1024, mmax=17)
206 c nmax is the dimension of the largest problem to be solved.
207 c mmax is the maximum number of limited memory corrections.
208 
209 c Declare the variables needed by the code.
210 c A description of all these variables is given at the end of
211 c the driver.
212 
213  character*60 task, csave
214  logical lsave(4)
215  integer n, m, iprint,
216  + nbd(nmax), iwa(3*nmax), isave(44)
217  double precision f, factr, pgtol,
218  + x(nmax), l(nmax), u(nmax), g(nmax), dsave(29),
219  + wa(2*mmax*nmax + 5*nmax + 11*mmax*mmax + 8*mmax)
220 
221 c Declare a few additional variables for this sample problem.
222 
223  double precision t1, t2
224  integer i
225 
226 c Test instrumentation: dump per-iteration state to JSON when the
227 c environment variable LBFGSB_JSON_OUTPUT is set. No-op otherwise.
228  character*512 lbfgsb_json
229  logical json_active
230 
231 c We wish to have output at every iteration.
232 
233  iprint = 1
234 
235 c We specify the tolerances in the stopping criteria.
236 
237  factr=1.0d+7
238  pgtol=1.0d-5
239 
240 c We specify the dimension n of the sample problem and the number
241 c m of limited memory corrections stored. (n and m should not
242 c exceed the limits nmax and mmax respectively.)
243 
244  n=25
245  m=5
246 
247 c We now provide nbd which defines the bounds on the variables:
248 c l specifies the lower bounds,
249 c u specifies the upper bounds.
250 
251 c First set bounds on the odd-numbered variables.
252 
253  do 10 i=1,n,2
254  nbd(i)=2
255  l(i)=1.0d0
256  u(i)=1.0d2
257  10 continue
258 
259 c Next set bounds on the even-numbered variables.
260 
261  do 12 i=2,n,2
262  nbd(i)=2
263  l(i)=-1.0d2
264  u(i)=1.0d2
265  12 continue
266 
267 c We now define the starting point.
268 
269  do 14 i=1,n
270  x(i)=3.0d0
271  14 continue
272 
273  write (6,16)
274  16 format(/,5x, 'Solving sample problem.',
275  + /,5x, ' (f = 0.0 at the optimal solution.)',/)
276 
277 c We start the iteration by initializing task.
278 c
279  task = 'START'
280 
281 c Test instrumentation: enabled when LBFGSB_JSON_OUTPUT is set.
282  call get_environment_variable('LBFGSB_JSON_OUTPUT', lbfgsb_json)
283  json_active = (len_trim(lbfgsb_json) .gt. 0)
284 
285 c ------- the beginning of the loop ----------
286 
287  111 continue
288 
289 c This is the call to the L-BFGS-B code.
290 
291  call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa,task,iprint,
292  + csave,lsave,isave,dsave)
293 
294  if (task(1:2) .eq. 'FG') then
295 c the minimization routine has returned to request the
296 c function f and gradient g values at the current x.
297 
298 c Compute function value f for the sample problem.
299 
300  f=.25d0*(x(1)-1.d0)**2
301  do 20 i=2,n
302  f=f+(x(i)-x(i-1)**2)**2
303  20 continue
304  f=4.d0*f
305 
306 c Compute gradient g for the sample problem.
307 
308  t1=x(2)-x(1)**2
309  g(1)=2.d0*(x(1)-1.d0)-1.6d1*x(1)*t1
310  do 22 i=2,n-1
311  t2=t1
312  t1=x(i+1)-x(i)**2
313  g(i)=8.d0*t2-1.6d1*x(i)*t1
314  22 continue
315  g(n)=8.d0*t1
316 
317 c go back to the minimization routine.
318  goto 111
319  endif
320 c
321  if (task(1:5) .eq. 'NEW_X') goto 111
322 c the minimization routine has returned with a new iterate,
323 c and we have opted to continue the iteration.
324 
325 c ---------- the end of the loop -------------
326 
327 c If task is neither FG nor NEW_X we terminate execution.
328 
329  if (json_active) then
330  call json_write_aggregate(trim(lbfgsb_json), task, f,
331  + dsave(13), n, x)
332  endif
333 
334  stop
335 
336  end
337 
338 c======================= The end of driver1 ============================
subroutine setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)
This subroutine partitions the working arrays wa and iwa, and then uses the limited memory BFGS metho...
Definition: setulb.f:190