1 /* lbfgsb.f -- translated by f2c (version 20100827).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17 #include <float.h>
18 // Timing features disabled #include <sys/time.h>
19
20 #undef abs
21 #undef min
22 #undef max
23 #include <math.h>
24 #include <stdio.h>
25 #define abs(x) ((x) >= 0 ? (x) : -(x))
26 #define min(a,b) ((a) <= (b) ? (a) : (b))
27 #define max(a,b) ((a) >= (b) ? (a) : (b))
28
29 /* Code below must be fixed to pass a real FILE* value to fprintf if
30 this macro is defined. */
31 #undef LBFGSB_ENABLE_ITERATE_FILE
32
lbfgsb_printf_vec(const char * name,double const * v,integer len)33 static void lbfgsb_printf_vec(const char* name,
34 double const* v,
35 integer len)
36 {
37 /*
38 1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
39 */
40 integer i;
41 printf("%s =", name);
42 for(i = 1; i <= len; ++i)
43 {
44 printf(" %11.4g", v[i]);
45 }
46 printf("\n");
47 }
48
49 /* Table of constant values */
50 static doublereal c_b9 = 0.;
51 static integer c__1 = 1;
52 static integer c__11 = 11;
53 static doublereal c_b280 = .001;
54 static doublereal c_b281 = .9;
55 static doublereal c_b282 = .1;
56
epsilon1()57 static doublereal epsilon1()
58 {
59 //Fortran code used intrinsic epsilon(1.) function
60 //which is the same as DBL_EPSILON.
61 return DBL_EPSILON;
62 } /* dpmeps_ */
63
dpmeps_()64 extern doublereal dpmeps_()
65 {
66 return epsilon1();
67 } /* dpmeps_ */
68
69
70 /* =========== L-BFGS-B (version 3.0. April 25, 2011 =================== */
71
72 /* This is a modified version of L-BFGS-B. Minor changes in the updated */
73 /* code appear preceded by a line comment as follows */
74
75 /* c-jlm-jn */
76
77 /* Major changes are described in the accompanying paper: */
78
79 /* Jorge Nocedal and Jose Luis Morales, Remark on "Algorithm 778: */
80 /* L-BFGS-B: Fortran Subroutines for Large-Scale Bound Constrained */
81 /* Optimization" (2011). To appear in ACM Transactions on */
82 /* Mathematical Software, */
83
84 /* The paper describes an improvement and a correction to Algorithm 778. */
85 /* It is shown that the performance of the algorithm can be improved */
86 /* significantly by making a relatively simple modication to the subspace */
87 /* minimization phase. The correction concerns an error caused by the use */
88 /* of routine dpmeps to estimate machine precision. */
89
90 /* The total work space **wa** required by the new version is */
91
92 /* 2*m*n + 11m*m + 5*n + 8*m */
93
94 /* the old version required */
95
96 /* 2*m*n + 12m*m + 4*n + 12*m */
97
98
99 /* J. Nocedal Department of Electrical Engineering and */
100 /* Computer Science. */
101 /* Northwestern University. Evanston, IL. USA */
102
103
104 /* J.L Morales Departamento de Matematicas, */
105 /* Instituto Tecnologico Autonomo de Mexico */
106 /* Mexico D.F. Mexico. */
107
108 /* March 2011 */
109
110 /* ============================================================================= */
setulb_(integer * n,integer * m,doublereal * x,doublereal * l,doublereal * u,integer * nbd,doublereal * f,doublereal * g,doublereal * factr,doublereal * pgtol,doublereal * wa,integer * iwa,char * task,integer * iprint,char * csave,logical * lsave,integer * isave,doublereal * dsave)111 /* Subroutine */ int setulb_(integer *n, integer *m, doublereal *x,
112 doublereal *l, doublereal *u, integer *nbd, doublereal *f, doublereal
113 *g, doublereal *factr, doublereal *pgtol, doublereal *wa, integer *
114 iwa, char *task, integer *iprint, char *csave, logical *lsave,
115 integer *isave, doublereal *dsave /*f2c autogenerated , ftnlen task_len, ftnlen csave_len */)
116 {
117 /* System generated locals */
118 integer i__1;
119
120 /* Builtin functions */
121 integer s_cmp(char *, char *, ftnlen, ftnlen);
122
123 /* Local variables */
124 static integer ld, lr, lt, lz, lwa, lwn, lss, lxp, lws, lwt, lsy, lwy,
125 lsnd;
126 extern /* Subroutine */ int mainlb_(integer *, integer *, doublereal *,
127 doublereal *, doublereal *, integer *, doublereal *, doublereal *,
128 doublereal *, doublereal *, doublereal *, doublereal *,
129 doublereal *, doublereal *, doublereal *, doublereal *,
130 doublereal *, doublereal *, doublereal *, doublereal *,
131 doublereal *, doublereal *, doublereal *, integer *, integer *,
132 integer *, char *, integer *, char *, logical *, integer *,
133 doublereal *, ftnlen, ftnlen);
134
135
136 /* -jlm-jn */
137 /* ************ */
138
139 /* Subroutine setulb */
140
141 /* This subroutine partitions the working arrays wa and iwa, and */
142 /* then uses the limited memory BFGS method to solve the bound */
143 /* constrained optimization problem by calling mainlb. */
144 /* (The direct method will be used in the subspace minimization.) */
145
146 /* n is an integer variable. */
147 /* On entry n is the dimension of the problem. */
148 /* On exit n is unchanged. */
149
150 /* m is an integer variable. */
151 /* On entry m is the maximum number of variable metric corrections */
152 /* used to define the limited memory matrix. */
153 /* On exit m is unchanged. */
154
155 /* x is a double precision array of dimension n. */
156 /* On entry x is an approximation to the solution. */
157 /* On exit x is the current approximation. */
158
159 /* l is a double precision array of dimension n. */
160 /* On entry l is the lower bound on x. */
161 /* On exit l is unchanged. */
162
163 /* u is a double precision array of dimension n. */
164 /* On entry u is the upper bound on x. */
165 /* On exit u is unchanged. */
166
167 /* nbd is an integer array of dimension n. */
168 /* On entry nbd represents the type of bounds imposed on the */
169 /* variables, and must be specified as follows: */
170 /* nbd(i)=0 if x(i) is unbounded, */
171 /* 1 if x(i) has only a lower bound, */
172 /* 2 if x(i) has both lower and upper bounds, and */
173 /* 3 if x(i) has only an upper bound. */
174 /* On exit nbd is unchanged. */
175
176 /* f is a double precision variable. */
177 /* On first entry f is unspecified. */
178 /* On final exit f is the value of the function at x. */
179
180 /* g is a double precision array of dimension n. */
181 /* On first entry g is unspecified. */
182 /* On final exit g is the value of the gradient at x. */
183
184 /* factr is a double precision variable. */
185 /* On entry factr >= 0 is specified by the user. The iteration */
186 /* will stop when */
187
188 /* (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
189
190 /* where epsmch is the machine precision, which is automatically */
191 /* generated by the code. Typical values for factr: 1.d+12 for */
192 /* low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely */
193 /* high accuracy. */
194 /* On exit factr is unchanged. */
195
196 /* pgtol is a double precision variable. */
197 /* On entry pgtol >= 0 is specified by the user. The iteration */
198 /* will stop when */
199
200 /* max{|proj g_i | i = 1, ..., n} <= pgtol */
201
202 /* where pg_i is the ith component of the projected gradient. */
203 /* On exit pgtol is unchanged. */
204
205 /* wa is a double precision working array of length */
206 /* (2mmax + 5)nmax + 12mmax^2 + 12mmax. */
207
208 /* iwa is an integer working array of length 3nmax. */
209
210 /* task is a working string of characters of length 60 indicating */
211 /* the current job when entering and quitting this subroutine. */
212
213 /* iprint is an integer variable that must be set by the user. */
214 /* It controls the frequency and type of output generated: */
215 /* iprint<0 no output is generated; */
216 /* iprint=0 print only one line at the last iteration; */
217 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
218 /* iprint=99 print details of every iteration except n-vectors; */
219 /* iprint=100 print also the changes of active set and final x; */
220 /* iprint>100 print details of every iteration including x and g; */
221 /* When iprint > 0, the file iterate.dat will be created to */
222 /* summarize the iteration. */
223
224 /* csave is a working string of characters of length 60. */
225
226 /* lsave is a logical working array of dimension 4. */
227 /* On exit with 'task' = NEW_X, the following information is */
228 /* available: */
229 /* If lsave(1) = .true. then the initial X has been replaced by */
230 /* its projection in the feasible set; */
231 /* If lsave(2) = .true. then the problem is constrained; */
232 /* If lsave(3) = .true. then each variable has upper and lower */
233 /* bounds; */
234
235 /* isave is an integer working array of dimension 44. */
236 /* On exit with 'task' = NEW_X, the following information is */
237 /* available: */
238 /* isave(22) = the total number of intervals explored in the */
239 /* search of Cauchy points; */
240 /* isave(26) = the total number of skipped BFGS updates before */
241 /* the current iteration; */
242 /* isave(30) = the number of current iteration; */
243 /* isave(31) = the total number of BFGS updates prior the current */
244 /* iteration; */
245 /* isave(33) = the number of intervals explored in the search of */
246 /* Cauchy point in the current iteration; */
247 /* isave(34) = the total number of function and gradient */
248 /* evaluations; */
249 /* isave(36) = the number of function value or gradient */
250 /* evaluations in the current iteration; */
251 /* if isave(37) = 0 then the subspace argmin is within the box; */
252 /* if isave(37) = 1 then the subspace argmin is beyond the box; */
253 /* isave(38) = the number of free variables in the current */
254 /* iteration; */
255 /* isave(39) = the number of active constraints in the current */
256 /* iteration; */
257 /* n + 1 - isave(40) = the number of variables leaving the set of */
258 /* active constraints in the current iteration; */
259 /* isave(41) = the number of variables entering the set of active */
260 /* constraints in the current iteration. */
261
262 /* dsave is a double precision working array of dimension 29. */
263 /* On exit with 'task' = NEW_X, the following information is */
264 /* available: */
265 /* dsave(1) = current 'theta' in the BFGS matrix; */
266 /* dsave(2) = f(x) in the previous iteration; */
267 /* dsave(3) = factr*epsmch; */
268 /* dsave(4) = 2-norm of the line search direction vector; */
269 /* dsave(5) = the machine precision epsmch generated by the code; */
270 /* dsave(7) = the accumulated time spent on searching for */
271 /* Cauchy points; */
272 /* dsave(8) = the accumulated time spent on */
273 /* subspace minimization; */
274 /* dsave(9) = the accumulated time spent on line search; */
275 /* dsave(11) = the slope of the line search function at */
276 /* the current point of line search; */
277 /* dsave(12) = the maximum relative step length imposed in */
278 /* line search; */
279 /* dsave(13) = the infinity norm of the projected gradient; */
280 /* dsave(14) = the relative step length in the line search; */
281 /* dsave(15) = the slope of the line search function at */
282 /* the starting point of the line search; */
283 /* dsave(16) = the square of the 2-norm of the line search */
284 /* direction vector. */
285
286 /* Subprograms called: */
287
288 /* L-BFGS-B Library ... mainlb. */
289
290
291 /* References: */
292
293 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
294 /* memory algorithm for bound constrained optimization'', */
295 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
296
297 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
298 /* limited memory FORTRAN code for solving bound constrained */
299 /* optimization problems'', Tech. Report, NAM-11, EECS Department, */
300 /* Northwestern University, 1994. */
301
302 /* (Postscript files of these papers are available via anonymous */
303 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
304
305 /* * * * */
306
307 /* NEOS, November 1994. (Latest revision June 1996.) */
308 /* Optimization Technology Center. */
309 /* Argonne National Laboratory and Northwestern University. */
310 /* Written by */
311 /* Ciyou Zhu */
312 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
313
314
315 /* ************ */
316 /* -jlm-jn */
317 /* Parameter adjustments */
318 --iwa;
319 --g;
320 --nbd;
321 --u;
322 --l;
323 --x;
324 --wa;
325 --lsave;
326 --isave;
327 --dsave;
328
329 /* Function Body */
330 if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
331 isave[1] = *m * *n;
332 /* Computing 2nd power */
333 i__1 = *m;
334 isave[2] = i__1 * i__1;
335 /* Computing 2nd power */
336 i__1 = *m;
337 isave[3] = i__1 * i__1 << 2;
338 isave[4] = 1;
339 /* ws m*n */
340 isave[5] = isave[4] + isave[1];
341 /* wy m*n */
342 isave[6] = isave[5] + isave[1];
343 /* wsy m**2 */
344 isave[7] = isave[6] + isave[2];
345 /* wss m**2 */
346 isave[8] = isave[7] + isave[2];
347 /* wt m**2 */
348 isave[9] = isave[8] + isave[2];
349 /* wn 4*m**2 */
350 isave[10] = isave[9] + isave[3];
351 /* wsnd 4*m**2 */
352 isave[11] = isave[10] + isave[3];
353 /* wz n */
354 isave[12] = isave[11] + *n;
355 /* wr n */
356 isave[13] = isave[12] + *n;
357 /* wd n */
358 isave[14] = isave[13] + *n;
359 /* wt n */
360 isave[15] = isave[14] + *n;
361 /* wxp n */
362 isave[16] = isave[15] + *n;
363 /* wa 8*m */
364 }
365 lws = isave[4];
366 lwy = isave[5];
367 lsy = isave[6];
368 lss = isave[7];
369 lwt = isave[8];
370 lwn = isave[9];
371 lsnd = isave[10];
372 lz = isave[11];
373 lr = isave[12];
374 ld = isave[13];
375 lt = isave[14];
376 lxp = isave[15];
377 lwa = isave[16];
378 mainlb_(n, m, &x[1], &l[1], &u[1], &nbd[1], f, &g[1], factr, pgtol, &wa[
379 lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lwt], &wa[lwn], &wa[lsnd],
380 &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lxp], &wa[lwa], &iwa[1],
381 &iwa[*n + 1], &iwa[(*n << 1) + 1], task, iprint, csave, &lsave[1],
382 &isave[22], &dsave[1], (ftnlen)60, (ftnlen)60);
383 return 0;
384 } /* setulb_ */
385
386 /* ======================= The end of setulb ============================= */
mainlb_(integer * n,integer * m,doublereal * x,doublereal * l,doublereal * u,integer * nbd,doublereal * f,doublereal * g,doublereal * factr,doublereal * pgtol,doublereal * ws,doublereal * wy,doublereal * sy,doublereal * ss,doublereal * wt,doublereal * wn,doublereal * snd,doublereal * z__,doublereal * r__,doublereal * d__,doublereal * t,doublereal * xp,doublereal * wa,integer * index,integer * iwhere,integer * indx2,char * task,integer * iprint,char * csave,logical * lsave,integer * isave,doublereal * dsave,ftnlen task_len,ftnlen csave_len)387 /* Subroutine */ int mainlb_(integer *n, integer *m, doublereal *x,
388 doublereal *l, doublereal *u, integer *nbd, doublereal *f, doublereal
389 *g, doublereal *factr, doublereal *pgtol, doublereal *ws, doublereal *
390 wy, doublereal *sy, doublereal *ss, doublereal *wt, doublereal *wn,
391 doublereal *snd, doublereal *z__, doublereal *r__, doublereal *d__,
392 doublereal *t, doublereal *xp, doublereal *wa, integer *index,
393 integer *iwhere, integer *indx2, char *task, integer *iprint, char *
394 csave, logical *lsave, integer *isave, doublereal *dsave, ftnlen
395 task_len, ftnlen csave_len)
396 {
397 /* System generated locals */
398 integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
399 ss_dim1, ss_offset, wt_dim1, wt_offset, wn_dim1, wn_offset,
400 snd_dim1, snd_offset, i__1;
401 doublereal d__1, d__2;
402
403 /* Builtin functions */
404 integer s_cmp(char *, char *, ftnlen, ftnlen);
405 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
406
407 /* Local variables */
408 static integer i__, k;
409 static doublereal gd, dr, rr, dtd;
410 static integer col;
411 static doublereal tol;
412 static logical wrk;
413 static doublereal stp, cpu1, cpu2;
414 static integer head;
415 static doublereal fold;
416 static integer nact;
417 static doublereal ddum;
418 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
419 integer *);
420 static integer info, nseg;
421 static doublereal time;
422 static integer nfgv, ifun, iter;
423 static char word[3];
424 static doublereal time1, time2;
425 static integer iback;
426 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
427 integer *);
428 static doublereal gdold;
429 static integer nfree;
430 static logical boxed;
431 static integer itail;
432 static doublereal theta;
433 extern /* Subroutine */ int freev_(integer *, integer *, integer *,
434 integer *, integer *, integer *, integer *, logical *, logical *,
435 logical *, integer *, integer *), dcopy_(integer *, doublereal *,
436 integer *, doublereal *, integer *);
437 static doublereal dnorm;
438 extern /* Subroutine */ int timer_(doublereal *), formk_(integer *,
439 integer *, integer *, integer *, integer *, integer *, integer *,
440 logical *, doublereal *, doublereal *, integer *, doublereal *,
441 doublereal *, doublereal *, doublereal *, integer *, integer *,
442 integer *);
443 static integer nskip, iword;
444 extern /* Subroutine */ int formt_(integer *, doublereal *, doublereal *,
445 doublereal *, integer *, doublereal *, integer *), subsm_(integer
446 *, integer *, integer *, integer *, doublereal *, doublereal *,
447 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
448 doublereal *, doublereal *, doublereal *, doublereal *, integer *
449 , integer *, integer *, doublereal *, doublereal *, integer *,
450 integer *);
451 static doublereal xstep, stpmx;
452 extern /* Subroutine */ int prn1lb_(integer *, integer *, doublereal *,
453 doublereal *, doublereal *, integer *, integer *, doublereal *),
454 prn2lb_(integer *, doublereal *, doublereal *, doublereal *,
455 integer *, integer *, integer *, integer *, integer *, doublereal
456 *, integer *, char *, integer *, integer *, doublereal *,
457 doublereal *, ftnlen), prn3lb_(integer *, doublereal *,
458 doublereal *, char *, integer *, integer *, integer *, integer *,
459 integer *, integer *, integer *, integer *, doublereal *,
460 doublereal *, integer *, char *, integer *, doublereal *,
461 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
462 ftnlen, ftnlen);
463 static integer ileave;
464 extern /* Subroutine */ int errclb_(integer *, integer *, doublereal *,
465 doublereal *, doublereal *, integer *, char *, integer *, integer
466 *, ftnlen);
467 static doublereal cachyt;
468 static integer itfile;
469 extern /* Subroutine */ int active_(integer *, doublereal *, doublereal *,
470 integer *, doublereal *, integer *, integer *, logical *,
471 logical *, logical *), cauchy_(integer *, doublereal *,
472 doublereal *, doublereal *, integer *, doublereal *, integer *,
473 integer *, doublereal *, doublereal *, doublereal *, integer *,
474 doublereal *, doublereal *, doublereal *, doublereal *,
475 doublereal *, integer *, integer *, doublereal *, doublereal *,
476 doublereal *, doublereal *, integer *, integer *, doublereal *,
477 integer *, doublereal *);
478 static doublereal epsmch;
479 extern /* Subroutine */ int cmprlb_(integer *, integer *, doublereal *,
480 doublereal *, doublereal *, doublereal *, doublereal *,
481 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
482 doublereal *, integer *, integer *, integer *, logical *,
483 integer *);
484 static logical updatd;
485 static doublereal sbtime;
486 static logical prjctd;
487 static integer iupdat;
488 static doublereal sbgnrm;
489 static logical cnstnd;
490 extern /* Subroutine */ int matupd_(integer *, integer *, doublereal *,
491 doublereal *, doublereal *, doublereal *, doublereal *,
492 doublereal *, integer *, integer *, integer *, integer *,
493 doublereal *, doublereal *, doublereal *, doublereal *,
494 doublereal *);
495 static integer nenter;
496 static doublereal lnscht;
497 extern /* Subroutine */ int lnsrlb_(integer *, doublereal *, doublereal *,
498 integer *, doublereal *, doublereal *, doublereal *, doublereal *
499 , doublereal *, doublereal *, doublereal *, doublereal *,
500 doublereal *, doublereal *, doublereal *, doublereal *,
501 doublereal *, doublereal *, doublereal *, integer *, integer *,
502 integer *, integer *, integer *, char *, logical *, logical *,
503 char *, integer *, doublereal *, ftnlen, ftnlen);
504 static integer nintol;
505 extern /* Subroutine */ int projgr_(integer *, doublereal *, doublereal *,
506 integer *, doublereal *, doublereal *, doublereal *);
507
508 /* -jlm-jn */
509 /* ************ */
510
511 /* Subroutine mainlb */
512
513 /* This subroutine solves bound constrained optimization problems by */
514 /* using the compact formula of the limited memory BFGS updates. */
515
516 /* n is an integer variable. */
517 /* On entry n is the number of variables. */
518 /* On exit n is unchanged. */
519
520 /* m is an integer variable. */
521 /* On entry m is the maximum number of variable metric */
522 /* corrections allowed in the limited memory matrix. */
523 /* On exit m is unchanged. */
524
525 /* x is a double precision array of dimension n. */
526 /* On entry x is an approximation to the solution. */
527 /* On exit x is the current approximation. */
528
529 /* l is a double precision array of dimension n. */
530 /* On entry l is the lower bound of x. */
531 /* On exit l is unchanged. */
532
533 /* u is a double precision array of dimension n. */
534 /* On entry u is the upper bound of x. */
535 /* On exit u is unchanged. */
536
537 /* nbd is an integer array of dimension n. */
538 /* On entry nbd represents the type of bounds imposed on the */
539 /* variables, and must be specified as follows: */
540 /* nbd(i)=0 if x(i) is unbounded, */
541 /* 1 if x(i) has only a lower bound, */
542 /* 2 if x(i) has both lower and upper bounds, */
543 /* 3 if x(i) has only an upper bound. */
544 /* On exit nbd is unchanged. */
545
546 /* f is a double precision variable. */
547 /* On first entry f is unspecified. */
548 /* On final exit f is the value of the function at x. */
549
550 /* g is a double precision array of dimension n. */
551 /* On first entry g is unspecified. */
552 /* On final exit g is the value of the gradient at x. */
553
554 /* factr is a double precision variable. */
555 /* On entry factr >= 0 is specified by the user. The iteration */
556 /* will stop when */
557
558 /* (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
559
560 /* where epsmch is the machine precision, which is automatically */
561 /* generated by the code. */
562 /* On exit factr is unchanged. */
563
564 /* pgtol is a double precision variable. */
565 /* On entry pgtol >= 0 is specified by the user. The iteration */
566 /* will stop when */
567
568 /* max{|proj g_i | i = 1, ..., n} <= pgtol */
569
570 /* where pg_i is the ith component of the projected gradient. */
571 /* On exit pgtol is unchanged. */
572
573 /* ws, wy, sy, and wt are double precision working arrays used to */
574 /* store the following information defining the limited memory */
575 /* BFGS matrix: */
576 /* ws, of dimension n x m, stores S, the matrix of s-vectors; */
577 /* wy, of dimension n x m, stores Y, the matrix of y-vectors; */
578 /* sy, of dimension m x m, stores S'Y; */
579 /* ss, of dimension m x m, stores S'S; */
580 /* yy, of dimension m x m, stores Y'Y; */
581 /* wt, of dimension m x m, stores the Cholesky factorization */
582 /* of (theta*S'S+LD^(-1)L'); see eq. */
583 /* (2.26) in [3]. */
584
585 /* wn is a double precision working array of dimension 2m x 2m */
586 /* used to store the LEL^T factorization of the indefinite matrix */
587 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
588 /* [L_a -R_z theta*S'AA'S ] */
589
590 /* where E = [-I 0] */
591 /* [ 0 I] */
592
593 /* snd is a double precision working array of dimension 2m x 2m */
594 /* used to store the lower triangular part of */
595 /* N = [Y' ZZ'Y L_a'+R_z'] */
596 /* [L_a +R_z S'AA'S ] */
597
598 /* z(n),r(n),d(n),t(n), xp(n),wa(8*m) are double precision working arrays. */
599 /* z is used at different times to store the Cauchy point and */
600 /* the Newton point. */
601 /* xp is used to safeguard the projected Newton direction */
602
603 /* sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. */
604
605 /* index is an integer working array of dimension n. */
606 /* In subroutine freev, index is used to store the free and fixed */
607 /* variables at the Generalized Cauchy Point (GCP). */
608
609 /* iwhere is an integer working array of dimension n used to record */
610 /* the status of the vector x for GCP computation. */
611 /* iwhere(i)=0 or -3 if x(i) is free and has bounds, */
612 /* 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) */
613 /* 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) */
614 /* 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) */
615 /* -1 if x(i) is always free, i.e., no bounds on it. */
616
617 /* indx2 is an integer working array of dimension n. */
618 /* Within subroutine cauchy, indx2 corresponds to the array iorder. */
619 /* In subroutine freev, a list of variables entering and leaving */
620 /* the free set is stored in indx2, and it is passed on to */
621 /* subroutine formk with this information. */
622
623 /* task is a working string of characters of length 60 indicating */
624 /* the current job when entering and leaving this subroutine. */
625
626 /* iprint is an INTEGER variable that must be set by the user. */
627 /* It controls the frequency and type of output generated: */
628 /* iprint<0 no output is generated; */
629 /* iprint=0 print only one line at the last iteration; */
630 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
631 /* iprint=99 print details of every iteration except n-vectors; */
632 /* iprint=100 print also the changes of active set and final x; */
633 /* iprint>100 print details of every iteration including x and g; */
634 /* When iprint > 0, the file iterate.dat will be created to */
635 /* summarize the iteration. */
636
637 /* csave is a working string of characters of length 60. */
638
639 /* lsave is a logical working array of dimension 4. */
640
641 /* isave is an integer working array of dimension 23. */
642
643 /* dsave is a double precision working array of dimension 29. */
644
645
646 /* Subprograms called */
647
648 /* L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, */
649
650 /* errclb, prn1lb, prn2lb, prn3lb, active, projgr, */
651
652 /* freev, cmprlb, matupd, formt. */
653
654 /* Minpack2 Library ... timer */
655
656 /* Linpack Library ... dcopy, ddot. */
657
658
659 /* References: */
660
661 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
662 /* memory algorithm for bound constrained optimization'', */
663 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
664
665 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
666 /* Subroutines for Large Scale Bound Constrained Optimization'' */
667 /* Tech. Report, NAM-11, EECS Department, Northwestern University, */
668 /* 1994. */
669
670 /* [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of */
671 /* Quasi-Newton Matrices and their use in Limited Memory Methods'', */
672 /* Mathematical Programming 63 (1994), no. 4, pp. 129-156. */
673
674 /* (Postscript files of these papers are available via anonymous */
675 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
676
677 /* * * * */
678
679 /* NEOS, November 1994. (Latest revision June 1996.) */
680 /* Optimization Technology Center. */
681 /* Argonne National Laboratory and Northwestern University. */
682 /* Written by */
683 /* Ciyou Zhu */
684 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
685
686
687 /* ************ */
688 /* Parameter adjustments */
689 --indx2;
690 --iwhere;
691 --index;
692 --xp;
693 --t;
694 --d__;
695 --r__;
696 --z__;
697 --g;
698 --nbd;
699 --u;
700 --l;
701 --x;
702 --wa;
703 snd_dim1 = 2 * *m;
704 snd_offset = 1 + snd_dim1;
705 snd -= snd_offset;
706 wn_dim1 = 2 * *m;
707 wn_offset = 1 + wn_dim1;
708 wn -= wn_offset;
709 wt_dim1 = *m;
710 wt_offset = 1 + wt_dim1;
711 wt -= wt_offset;
712 ss_dim1 = *m;
713 ss_offset = 1 + ss_dim1;
714 ss -= ss_offset;
715 sy_dim1 = *m;
716 sy_offset = 1 + sy_dim1;
717 sy -= sy_offset;
718 wy_dim1 = *n;
719 wy_offset = 1 + wy_dim1;
720 wy -= wy_offset;
721 ws_dim1 = *n;
722 ws_offset = 1 + ws_dim1;
723 ws -= ws_offset;
724 --lsave;
725 --isave;
726 --dsave;
727
728 /* Function Body */
729 if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
730 epsmch = epsilon1();
731 timer_(&time1);
732 /* Initialize counters and scalars when task='START'. */
733 /* for the limited memory BFGS matrices: */
734 col = 0;
735 head = 1;
736 theta = 1.;
737 iupdat = 0;
738 updatd = FALSE_;
739 iback = 0;
740 itail = 0;
741 iword = 0;
742 nact = 0;
743 ileave = 0;
744 nenter = 0;
745 fold = 0.;
746 dnorm = 0.;
747 cpu1 = 0.;
748 gd = 0.;
749 stpmx = 0.;
750 sbgnrm = 0.;
751 stp = 0.;
752 gdold = 0.;
753 dtd = 0.;
754 /* for operation counts: */
755 iter = 0;
756 nfgv = 0;
757 nseg = 0;
758 nintol = 0;
759 nskip = 0;
760 nfree = *n;
761 ifun = 0;
762 /* for stopping tolerance: */
763 tol = *factr * epsmch;
764 /* for measuring running time: */
765 cachyt = 0.;
766 sbtime = 0.;
767 lnscht = 0.;
768 /* 'word' records the status of subspace solutions. */
769 s_copy(word, "---", (ftnlen)3, (ftnlen)3);
770 /* 'info' records the termination information. */
771 info = 0;
772 itfile = 8;
773 if (*iprint >= 1) {
774 /* open a summary file 'iterate.dat' */
775 /*
776 o__1.oerr = 0;
777 o__1.ounit = 8;
778 o__1.ofnmlen = 11;
779 o__1.ofnm = "iterate.dat";
780 o__1.orl = 0;
781 o__1.osta = "unknown";
782 o__1.oacc = 0;
783 o__1.ofm = 0;
784 o__1.oblnk = 0;
785 f_open(&o__1);
786 */
787 itfile = 8;
788 }
789 /* Check the input arguments for errors. */
790 errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (ftnlen)
791 60);
792 if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
793 prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &
794 nintol, &nskip, &nact, &sbgnrm, &c_b9, &nseg, word, &
795 iback, &stp, &xstep, &k, &cachyt, &sbtime, &lnscht, (
796 ftnlen)60, (ftnlen)3);
797 return 0;
798 }
799 prn1lb_(n, m, &l[1], &u[1], &x[1], iprint, &itfile, &epsmch);
800 /* Initialize iwhere & project x onto the feasible set. */
801 active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd,
802 &cnstnd, &boxed);
803 /* The end of the initialization. */
804 } else {
805 /* restore local variables. */
806 prjctd = lsave[1];
807 cnstnd = lsave[2];
808 boxed = lsave[3];
809 updatd = lsave[4];
810 nintol = isave[1];
811 itfile = isave[3];
812 iback = isave[4];
813 nskip = isave[5];
814 head = isave[6];
815 col = isave[7];
816 itail = isave[8];
817 iter = isave[9];
818 iupdat = isave[10];
819 nseg = isave[12];
820 nfgv = isave[13];
821 info = isave[14];
822 ifun = isave[15];
823 iword = isave[16];
824 nfree = isave[17];
825 nact = isave[18];
826 ileave = isave[19];
827 nenter = isave[20];
828 theta = dsave[1];
829 fold = dsave[2];
830 tol = dsave[3];
831 dnorm = dsave[4];
832 epsmch = dsave[5];
833 cpu1 = dsave[6];
834 cachyt = dsave[7];
835 sbtime = dsave[8];
836 lnscht = dsave[9];
837 time1 = dsave[10];
838 gd = dsave[11];
839 stpmx = dsave[12];
840 sbgnrm = dsave[13];
841 stp = dsave[14];
842 gdold = dsave[15];
843 dtd = dsave[16];
844 /* After returning from the driver go to the point where execution */
845 /* is to resume. */
846 if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
847 goto L666;
848 }
849 if (s_cmp(task, "NEW_X", (ftnlen)5, (ftnlen)5) == 0) {
850 goto L777;
851 }
852 if (s_cmp(task, "FG_ST", (ftnlen)5, (ftnlen)5) == 0) {
853 goto L111;
854 }
855 if (s_cmp(task, "STOP", (ftnlen)4, (ftnlen)4) == 0) {
856 if (s_cmp(task + 6, "CPU", (ftnlen)3, (ftnlen)3) == 0) {
857 /* restore the previous iterate. */
858 dcopy_(n, &t[1], &c__1, &x[1], &c__1);
859 dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
860 *f = fold;
861 }
862 goto L999;
863 }
864 }
865 /* Compute f0 and g0. */
866 s_copy(task, "FG_START", (ftnlen)60, (ftnlen)(8+1));
867 /* return to the driver to calculate f and g; reenter at 111. */
868 goto L1000;
869 L111:
870 nfgv = 1;
871 /* Compute the infinity norm of the (-) projected gradient. */
872 projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
873 if (*iprint >= 1) {
874 /*
875 1002 format
876 + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
877 */
878 printf("At iterate %5ld f= %12.5g |proj g|= %12.5g\n",
879 iter, *f, sbgnrm);
880 /*
881 1003 format (2(1x,i4),5x,'-',5x,'-',3x,'-',5x,'-',5x,'-',8x,'-',3x,
882 + 1p,2(1x,d10.3))
883 */
884 #ifdef LBFGSB_ENABLE_ITERATE_FILE
885 fprintf(0,
886 " %4ld %4ld - - - - - - %10.3g %10.3g\n",
887 iter, nfgv, sbgnrm, *f);
888 #endif
889 }
890 if (sbgnrm <= *pgtol) {
891 /* terminate the algorithm. */
892 s_copy(task, "CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL", (
893 ftnlen)60, (ftnlen)(48 + 1));
894 goto L999;
895 }
896 /* ----------------- the beginning of the loop -------------------------- */
897 L222:
898 if (*iprint >= 99) {
899 i__1 = iter + 1;
900 printf("ITERATION %5ld\n", i__1);
901 }
902 iword = -1;
903
904 if (! cnstnd && col > 0) {
905 /* skip the search for GCP. */
906 dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
907 wrk = updatd;
908 nseg = 0;
909 goto L333;
910 }
911 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
912
913 /* Compute the Generalized Cauchy Point (GCP). */
914
915 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
916 timer_(&cpu1);
917 cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[
918 1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[
919 sy_offset], &wt[wt_offset], &theta, &col, &head, &wa[1], &wa[(*m
920 << 1) + 1], &wa[(*m << 2) + 1], &wa[*m * 6 + 1], &nseg, iprint, &
921 sbgnrm, &info, &epsmch);
922 if (info != 0) {
923 /* singular triangular system detected; refresh the lbfgs memory. */
924 if (*iprint >= 1) {
925 printf(" Singular triangular system detected;\n"
926 " refresh the lbfgs memory and restart the iteration.\n");
927 }
928 info = 0;
929 col = 0;
930 head = 1;
931 theta = 1.;
932 iupdat = 0;
933 updatd = FALSE_;
934 timer_(&cpu2);
935 cachyt = cachyt + cpu2 - cpu1;
936 goto L222;
937 }
938 timer_(&cpu2);
939 cachyt = cachyt + cpu2 - cpu1;
940 nintol += nseg;
941 /* Count the entering and leaving variables for iter > 0; */
942 /* find the index set of free and active variables at the GCP. */
943 freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], &
944 wrk, &updatd, &cnstnd, iprint, &iter);
945 nact = *n - nfree;
946 L333:
947 /* If there are no free variables or B=theta*I, then */
948 /* skip the subspace minimization. */
949 if (nfree == 0 || col == 0) {
950 goto L555;
951 }
952 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
953
954 /* Subspace minimization. */
955
956 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
957 timer_(&cpu1);
958 /* Form the LEL^T factorization of the indefinite */
959 /* matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
960 /* [L_a -R_z theta*S'AA'S ] */
961 /* where E = [-I 0] */
962 /* [ 0 I] */
963 if (wrk) {
964 formk_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iupdat, &
965 updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], &
966 wy[wy_offset], &sy[sy_offset], &theta, &col, &head, &info);
967 }
968 if (info != 0) {
969 /* nonpositive definiteness in Cholesky factorization; */
970 /* refresh the lbfgs memory and restart the iteration. */
971 if (*iprint >= 1) {
972 printf(" Nonpositive definiteness in Cholesky factorization in formk;\n"
973 " refresh the lbfgs memory and restart the iteration.\n");
974 }
975 info = 0;
976 col = 0;
977 head = 1;
978 theta = 1.;
979 iupdat = 0;
980 updatd = FALSE_;
981 timer_(&cpu2);
982 sbtime = sbtime + cpu2 - cpu1;
983 goto L222;
984 }
985 /* compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */
986 /* from 'cauchy'). */
987 cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset]
988 , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], &theta, &
989 col, &head, &nfree, &cnstnd, &info);
990 if (info != 0) {
991 goto L444;
992 }
993 /* -jlm-jn call the direct method. */
994 subsm_(n, m, &nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], &
995 xp[1], &ws[ws_offset], &wy[wy_offset], &theta, &x[1], &g[1], &col,
996 &head, &iword, &wa[1], &wn[wn_offset], iprint, &info);
997 L444:
998 if (info != 0) {
999 /* singular triangular system detected; */
1000 /* refresh the lbfgs memory and restart the iteration. */
1001 if (*iprint >= 1) {
1002 printf(" Singular triangular system detected;\n"
1003 " refresh the lbfgs memory and restart the iteration.\n");
1004 }
1005 info = 0;
1006 col = 0;
1007 head = 1;
1008 theta = 1.;
1009 iupdat = 0;
1010 updatd = FALSE_;
1011 timer_(&cpu2);
1012 sbtime = sbtime + cpu2 - cpu1;
1013 goto L222;
1014 }
1015 timer_(&cpu2);
1016 sbtime = sbtime + cpu2 - cpu1;
1017 L555:
1018 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1019
1020 /* Line search and optimality tests. */
1021
1022 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1023 /* Generate the search direction d:=z-x. */
1024 i__1 = *n;
1025 for (i__ = 1; i__ <= i__1; ++i__) {
1026 d__[i__] = z__[i__] - x[i__];
1027 /* L40: */
1028 }
1029 timer_(&cpu1);
1030 L666:
1031 lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], &
1032 d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, &
1033 stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd,
1034 csave, &isave[22], &dsave[17], (ftnlen)60, (ftnlen)60);
1035 if (info != 0 || iback >= 20) {
1036 /* restore the previous iterate. */
1037 dcopy_(n, &t[1], &c__1, &x[1], &c__1);
1038 dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
1039 *f = fold;
1040 if (col == 0) {
1041 /* abnormal termination. */
1042 if (info == 0) {
1043 info = -9;
1044 /* restore the actual number of f and g evaluations etc. */
1045 --nfgv;
1046 --ifun;
1047 --iback;
1048 }
1049 s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (ftnlen)60, (
1050 ftnlen)(30+1));
1051 ++iter;
1052 goto L999;
1053 } else {
1054 /* refresh the lbfgs memory and restart the iteration. */
1055 if (*iprint >= 1) {
1056 printf(" Bad direction in the line search;\n"
1057 " refresh the lbfgs memory and restart the iteration.\n");
1058 }
1059 if (info == 0) {
1060 --nfgv;
1061 }
1062 info = 0;
1063 col = 0;
1064 head = 1;
1065 theta = 1.;
1066 iupdat = 0;
1067 updatd = FALSE_;
1068 s_copy(task, "RESTART_FROM_LNSRCH", (ftnlen)60, (ftnlen)(19+1));
1069 timer_(&cpu2);
1070 lnscht = lnscht + cpu2 - cpu1;
1071 goto L222;
1072 }
1073 } else if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
1074 /* return to the driver for calculating f and g; reenter at 666. */
1075 goto L1000;
1076 } else {
1077 /* calculate and print out the quantities related to the new X. */
1078 timer_(&cpu2);
1079 lnscht = lnscht + cpu2 - cpu1;
1080 ++iter;
1081 /* Compute the infinity norm of the projected (-)gradient. */
1082 projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
1083 /* Print iteration information. */
1084 prn2lb_(n, &x[1], f, &g[1], iprint, &itfile, &iter, &nfgv, &nact, &
1085 sbgnrm, &nseg, word, &iword, &iback, &stp, &xstep, (ftnlen)3);
1086 goto L1000;
1087 }
1088 L777:
1089 /* Test for termination. */
1090 if (sbgnrm <= *pgtol) {
1091 /* terminate the algorithm. */
1092 s_copy(task, "CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL", (
1093 ftnlen)60, (ftnlen)(48+1));
1094 goto L999;
1095 }
1096 /* Computing MAX */
1097 d__1 = abs(fold), d__2 = abs(*f), d__1 = max(d__1,d__2);
1098 ddum = max(d__1,1.);
1099 if (fold - *f <= tol * ddum) {
1100 /* terminate the algorithm. */
1101 s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", (
1102 ftnlen)60, (ftnlen)(47+1));
1103 if (iback >= 10) {
1104 info = -5;
1105 }
1106 /* i.e., to issue a warning if iback>10 in the line search. */
1107 goto L999;
1108 }
1109 /* Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
1110 i__1 = *n;
1111 for (i__ = 1; i__ <= i__1; ++i__) {
1112 r__[i__] = g[i__] - r__[i__];
1113 /* L42: */
1114 }
1115 rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
1116 if (stp == 1.) {
1117 dr = gd - gdold;
1118 ddum = -gdold;
1119 } else {
1120 dr = (gd - gdold) * stp;
1121 dscal_(n, &stp, &d__[1], &c__1);
1122 ddum = -gdold * stp;
1123 }
1124 if (dr <= epsmch * ddum) {
1125 /* skip the L-BFGS update. */
1126 ++nskip;
1127 updatd = FALSE_;
1128 if (*iprint >= 1) {
1129 printf(" ys=%10.3g -gs=%10.3g BFGS update SKIPPED\n",
1130 dr, ddum);
1131 }
1132 goto L888;
1133 }
1134 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1135
1136 /* Update the L-BFGS matrix. */
1137
1138 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1139 updatd = TRUE_;
1140 ++iupdat;
1141 /* Update matrices WS and WY and form the middle matrix in B. */
1142 matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[
1143 ss_offset], &d__[1], &r__[1], &itail, &iupdat, &col, &head, &
1144 theta, &rr, &dr, &stp, &dtd);
1145 /* Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */
1146 /* Store T in the upper triangular of the array wt; */
1147 /* Cholesky factorize T to J*J' with */
1148 /* J' stored in the upper triangular of wt. */
1149 formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], &col, &theta, &
1150 info);
1151 if (info != 0) {
1152 /* nonpositive definiteness in Cholesky factorization; */
1153 /* refresh the lbfgs memory and restart the iteration. */
1154 if (*iprint >= 1) {
1155 printf(" Nonpositive definiteness in Cholesky factorization in formt;\n"
1156 " refresh the lbfgs memory and restart the iteration.\n");
1157 }
1158 info = 0;
1159 col = 0;
1160 head = 1;
1161 theta = 1.;
1162 iupdat = 0;
1163 updatd = FALSE_;
1164 goto L222;
1165 }
1166 /* Now the inverse of the middle matrix in B is */
1167 /* [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] */
1168 /* [ -L*D^(-1/2) J ] [ 0 J' ] */
1169 L888:
1170 /* -------------------- the end of the loop ----------------------------- */
1171 goto L222;
1172 L999:
1173 timer_(&time2);
1174 time = time2 - time1;
1175 prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &nintol,
1176 &nskip, &nact, &sbgnrm, &time, &nseg, word, &iback, &stp, &xstep,
1177 &k, &cachyt, &sbtime, &lnscht, (ftnlen)60, (ftnlen)3);
1178 L1000:
1179 /* Save local variables. */
1180 lsave[1] = prjctd;
1181 lsave[2] = cnstnd;
1182 lsave[3] = boxed;
1183 lsave[4] = updatd;
1184 isave[1] = nintol;
1185 isave[3] = itfile;
1186 isave[4] = iback;
1187 isave[5] = nskip;
1188 isave[6] = head;
1189 isave[7] = col;
1190 isave[8] = itail;
1191 isave[9] = iter;
1192 isave[10] = iupdat;
1193 isave[12] = nseg;
1194 isave[13] = nfgv;
1195 isave[14] = info;
1196 isave[15] = ifun;
1197 isave[16] = iword;
1198 isave[17] = nfree;
1199 isave[18] = nact;
1200 isave[19] = ileave;
1201 isave[20] = nenter;
1202 dsave[1] = theta;
1203 dsave[2] = fold;
1204 dsave[3] = tol;
1205 dsave[4] = dnorm;
1206 dsave[5] = epsmch;
1207 dsave[6] = cpu1;
1208 dsave[7] = cachyt;
1209 dsave[8] = sbtime;
1210 dsave[9] = lnscht;
1211 dsave[10] = time1;
1212 dsave[11] = gd;
1213 dsave[12] = stpmx;
1214 dsave[13] = sbgnrm;
1215 dsave[14] = stp;
1216 dsave[15] = gdold;
1217 dsave[16] = dtd;
1218 return 0;
1219 } /* mainlb_ */
1220
1221 /* ======================= The end of mainlb ============================= */
active_(integer * n,doublereal * l,doublereal * u,integer * nbd,doublereal * x,integer * iwhere,integer * iprint,logical * prjctd,logical * cnstnd,logical * boxed)1222 /* Subroutine */ int active_(integer *n, doublereal *l, doublereal *u,
1223 integer *nbd, doublereal *x, integer *iwhere, integer *iprint,
1224 logical *prjctd, logical *cnstnd, logical *boxed)
1225 {
1226 /* System generated locals */
1227 integer i__1;
1228
1229 /* Local variables */
1230 static integer i__, nbdd;
1231
1232 /* ************ */
1233
1234 /* Subroutine active */
1235
1236 /* This subroutine initializes iwhere and projects the initial x to */
1237 /* the feasible set if necessary. */
1238
1239 /* iwhere is an integer array of dimension n. */
1240 /* On entry iwhere is unspecified. */
1241 /* On exit iwhere(i)=-1 if x(i) has no bounds */
1242 /* 3 if l(i)=u(i) */
1243 /* 0 otherwise. */
1244 /* In cauchy, iwhere is given finer gradations. */
1245
1246
1247 /* * * * */
1248
1249 /* NEOS, November 1994. (Latest revision June 1996.) */
1250 /* Optimization Technology Center. */
1251 /* Argonne National Laboratory and Northwestern University. */
1252 /* Written by */
1253 /* Ciyou Zhu */
1254 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
1255
1256
1257 /* ************ */
1258 /* Initialize nbdd, prjctd, cnstnd and boxed. */
1259 /* Parameter adjustments */
1260 --iwhere;
1261 --x;
1262 --nbd;
1263 --u;
1264 --l;
1265
1266 /* Function Body */
1267 nbdd = 0;
1268 *prjctd = FALSE_;
1269 *cnstnd = FALSE_;
1270 *boxed = TRUE_;
1271 /* Project the initial x to the easible set if necessary. */
1272 i__1 = *n;
1273 for (i__ = 1; i__ <= i__1; ++i__) {
1274 if (nbd[i__] > 0) {
1275 if (nbd[i__] <= 2 && x[i__] <= l[i__]) {
1276 if (x[i__] < l[i__]) {
1277 *prjctd = TRUE_;
1278 x[i__] = l[i__];
1279 }
1280 ++nbdd;
1281 } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) {
1282 if (x[i__] > u[i__]) {
1283 *prjctd = TRUE_;
1284 x[i__] = u[i__];
1285 }
1286 ++nbdd;
1287 }
1288 }
1289 /* L10: */
1290 }
1291 /* Initialize iwhere and assign values to cnstnd and boxed. */
1292 i__1 = *n;
1293 for (i__ = 1; i__ <= i__1; ++i__) {
1294 if (nbd[i__] != 2) {
1295 *boxed = FALSE_;
1296 }
1297 if (nbd[i__] == 0) {
1298 /* this variable is always free */
1299 iwhere[i__] = -1;
1300 /* otherwise set x(i)=mid(x(i), u(i), l(i)). */
1301 } else {
1302 *cnstnd = TRUE_;
1303 if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) {
1304 /* this variable is always fixed */
1305 iwhere[i__] = 3;
1306 } else {
1307 iwhere[i__] = 0;
1308 }
1309 }
1310 /* L20: */
1311 }
1312 if (*iprint >= 0) {
1313 if (*prjctd) {
1314 printf("The initial X is infeasible. Restart with its projection.\n");
1315 }
1316 if (! (*cnstnd)) {
1317 printf("This problem is unconstrained.\n");
1318 }
1319 }
1320 if (*iprint > 0) {
1321 printf("At X0 %9ld variables are exactly at the bounds\n", nbdd);
1322 }
1323 return 0;
1324 } /* active_ */
1325
1326 /* ======================= The end of active ============================= */
bmv_(integer * m,doublereal * sy,doublereal * wt,integer * col,doublereal * v,doublereal * p,integer * info)1327 /* Subroutine */ int bmv_(integer *m, doublereal *sy, doublereal *wt, integer
1328 *col, doublereal *v, doublereal *p, integer *info)
1329 {
1330 /* System generated locals */
1331 integer sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2;
1332
1333 /* Builtin functions */
1334 double sqrt(doublereal);
1335
1336 /* Local variables */
1337 static integer i__, k, i2;
1338 static doublereal sum;
1339 extern /* Subroutine */ int dtrsl_(doublereal *, integer *, integer *,
1340 doublereal *, integer *, integer *);
1341
1342 /* ************ */
1343
1344 /* Subroutine bmv */
1345
1346 /* This subroutine computes the product of the 2m x 2m middle matrix */
1347 /* in the compact L-BFGS formula of B and a 2m vector v; */
1348 /* it returns the product in p. */
1349
1350 /* m is an integer variable. */
1351 /* On entry m is the maximum number of variable metric corrections */
1352 /* used to define the limited memory matrix. */
1353 /* On exit m is unchanged. */
1354
1355 /* sy is a double precision array of dimension m x m. */
1356 /* On entry sy specifies the matrix S'Y. */
1357 /* On exit sy is unchanged. */
1358
1359 /* wt is a double precision array of dimension m x m. */
1360 /* On entry wt specifies the upper triangular matrix J' which is */
1361 /* the Cholesky factor of (thetaS'S+LD^(-1)L'). */
1362 /* On exit wt is unchanged. */
1363
1364 /* col is an integer variable. */
1365 /* On entry col specifies the number of s-vectors (or y-vectors) */
1366 /* stored in the compact L-BFGS formula. */
1367 /* On exit col is unchanged. */
1368
1369 /* v is a double precision array of dimension 2col. */
1370 /* On entry v specifies vector v. */
1371 /* On exit v is unchanged. */
1372
1373 /* p is a double precision array of dimension 2col. */
1374 /* On entry p is unspecified. */
1375 /* On exit p is the product Mv. */
1376
1377 /* info is an integer variable. */
1378 /* On entry info is unspecified. */
1379 /* On exit info = 0 for normal return, */
1380 /* = nonzero for abnormal return when the system */
1381 /* to be solved by dtrsl is singular. */
1382
1383 /* Subprograms called: */
1384
1385 /* Linpack ... dtrsl. */
1386
1387
1388 /* * * * */
1389
1390 /* NEOS, November 1994. (Latest revision June 1996.) */
1391 /* Optimization Technology Center. */
1392 /* Argonne National Laboratory and Northwestern University. */
1393 /* Written by */
1394 /* Ciyou Zhu */
1395 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
1396
1397
1398 /* ************ */
1399 /* Parameter adjustments */
1400 wt_dim1 = *m;
1401 wt_offset = 1 + wt_dim1;
1402 wt -= wt_offset;
1403 sy_dim1 = *m;
1404 sy_offset = 1 + sy_dim1;
1405 sy -= sy_offset;
1406 --p;
1407 --v;
1408
1409 /* Function Body */
1410 if (*col == 0) {
1411 return 0;
1412 }
1413 /* PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] */
1414 /* [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. */
1415 /* solve Jp2=v2+LD^(-1)v1. */
1416 p[*col + 1] = v[*col + 1];
1417 i__1 = *col;
1418 for (i__ = 2; i__ <= i__1; ++i__) {
1419 i2 = *col + i__;
1420 sum = 0.;
1421 i__2 = i__ - 1;
1422 for (k = 1; k <= i__2; ++k) {
1423 sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1];
1424 /* L10: */
1425 }
1426 p[i2] = v[i2] + sum;
1427 /* L20: */
1428 }
1429 /* Solve the triangular system */
1430 dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info);
1431 if (*info != 0) {
1432 return 0;
1433 }
1434 /* solve D^(1/2)p1=v1. */
1435 i__1 = *col;
1436 for (i__ = 1; i__ <= i__1; ++i__) {
1437 p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
1438 /* L30: */
1439 }
1440 /* PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] */
1441 /* [ 0 J' ] [ p2 ] [ p2 ]. */
1442 /* solve J^Tp2=p2. */
1443 dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info);
1444 if (*info != 0) {
1445 return 0;
1446 }
1447 /* compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */
1448 /* =-D^(-1/2)p1+D^(-1)L'p2. */
1449 i__1 = *col;
1450 for (i__ = 1; i__ <= i__1; ++i__) {
1451 p[i__] = -p[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
1452 /* L40: */
1453 }
1454 i__1 = *col;
1455 for (i__ = 1; i__ <= i__1; ++i__) {
1456 sum = 0.;
1457 i__2 = *col;
1458 for (k = i__ + 1; k <= i__2; ++k) {
1459 sum += sy[k + i__ * sy_dim1] * p[*col + k] / sy[i__ + i__ *
1460 sy_dim1];
1461 /* L50: */
1462 }
1463 p[i__] += sum;
1464 /* L60: */
1465 }
1466 return 0;
1467 } /* bmv_ */
1468
1469 /* ======================== The end of bmv =============================== */
cauchy_(integer * n,doublereal * x,doublereal * l,doublereal * u,integer * nbd,doublereal * g,integer * iorder,integer * iwhere,doublereal * t,doublereal * d__,doublereal * xcp,integer * m,doublereal * wy,doublereal * ws,doublereal * sy,doublereal * wt,doublereal * theta,integer * col,integer * head,doublereal * p,doublereal * c__,doublereal * wbp,doublereal * v,integer * nseg,integer * iprint,doublereal * sbgnrm,integer * info,doublereal * epsmch)1470 /* Subroutine */ int cauchy_(integer *n, doublereal *x, doublereal *l,
1471 doublereal *u, integer *nbd, doublereal *g, integer *iorder, integer *
1472 iwhere, doublereal *t, doublereal *d__, doublereal *xcp, integer *m,
1473 doublereal *wy, doublereal *ws, doublereal *sy, doublereal *wt,
1474 doublereal *theta, integer *col, integer *head, doublereal *p,
1475 doublereal *c__, doublereal *wbp, doublereal *v, integer *nseg,
1476 integer *iprint, doublereal *sbgnrm, integer *info, doublereal *
1477 epsmch)
1478 {
1479 /* System generated locals */
1480 integer wy_dim1, wy_offset, ws_dim1, ws_offset, sy_dim1, sy_offset,
1481 wt_dim1, wt_offset, i__1, i__2;
1482 doublereal d__1;
1483
1484 /* Local variables */
1485 static integer i__, j;
1486 static doublereal f1, f2, dt, tj, tl, tu, tj0;
1487 static integer ibp;
1488 static doublereal dtm;
1489 extern /* Subroutine */ int bmv_(integer *, doublereal *, doublereal *,
1490 integer *, doublereal *, doublereal *, integer *);
1491 static doublereal wmc, wmp, wmw;
1492 static integer col2;
1493 static doublereal dibp;
1494 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
1495 integer *);
1496 static integer iter;
1497 static doublereal zibp, tsum, dibp2;
1498 static logical bnded;
1499 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
1500 integer *);
1501 static doublereal neggi;
1502 static integer nfree;
1503 static doublereal bkmin;
1504 static integer nleft;
1505 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
1506 doublereal *, integer *), daxpy_(integer *, doublereal *,
1507 doublereal *, integer *, doublereal *, integer *);
1508 static doublereal f2_org__;
1509 static integer nbreak, ibkmin;
1510 extern /* Subroutine */ int hpsolb_(integer *, doublereal *, integer *,
1511 integer *);
1512 static integer pointr;
1513 static logical xlower, xupper;
1514
1515 /* ************ */
1516
1517 /* Subroutine cauchy */
1518
1519 /* For given x, l, u, g (with sbgnrm > 0), and a limited memory */
1520 /* BFGS matrix B defined in terms of matrices WY, WS, WT, and */
1521 /* scalars head, col, and theta, this subroutine computes the */
1522 /* generalized Cauchy point (GCP), defined as the first local */
1523 /* minimizer of the quadratic */
1524
1525 /* Q(x + s) = g's + 1/2 s'Bs */
1526
1527 /* along the projected gradient direction P(x-tg,l,u). */
1528 /* The routine returns the GCP in xcp. */
1529
1530 /* n is an integer variable. */
1531 /* On entry n is the dimension of the problem. */
1532 /* On exit n is unchanged. */
1533
1534 /* x is a double precision array of dimension n. */
1535 /* On entry x is the starting point for the GCP computation. */
1536 /* On exit x is unchanged. */
1537
1538 /* l is a double precision array of dimension n. */
1539 /* On entry l is the lower bound of x. */
1540 /* On exit l is unchanged. */
1541
1542 /* u is a double precision array of dimension n. */
1543 /* On entry u is the upper bound of x. */
1544 /* On exit u is unchanged. */
1545
1546 /* nbd is an integer array of dimension n. */
1547 /* On entry nbd represents the type of bounds imposed on the */
1548 /* variables, and must be specified as follows: */
1549 /* nbd(i)=0 if x(i) is unbounded, */
1550 /* 1 if x(i) has only a lower bound, */
1551 /* 2 if x(i) has both lower and upper bounds, and */
1552 /* 3 if x(i) has only an upper bound. */
1553 /* On exit nbd is unchanged. */
1554
1555 /* g is a double precision array of dimension n. */
1556 /* On entry g is the gradient of f(x). g must be a nonzero vector. */
1557 /* On exit g is unchanged. */
1558
1559 /* iorder is an integer working array of dimension n. */
1560 /* iorder will be used to store the breakpoints in the piecewise */
1561 /* linear path and free variables encountered. On exit, */
1562 /* iorder(1),...,iorder(nleft) are indices of breakpoints */
1563 /* which have not been encountered; */
1564 /* iorder(nleft+1),...,iorder(nbreak) are indices of */
1565 /* encountered breakpoints; and */
1566 /* iorder(nfree),...,iorder(n) are indices of variables which */
1567 /* have no bound constraits along the search direction. */
1568
1569 /* iwhere is an integer array of dimension n. */
1570 /* On entry iwhere indicates only the permanently fixed (iwhere=3) */
1571 /* or free (iwhere= -1) components of x. */
1572 /* On exit iwhere records the status of the current x variables. */
1573 /* iwhere(i)=-3 if x(i) is free and has bounds, but is not moved */
1574 /* 0 if x(i) is free and has bounds, and is moved */
1575 /* 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) */
1576 /* 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) */
1577 /* 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) */
1578 /* -1 if x(i) is always free, i.e., it has no bounds. */
1579
1580 /* t is a double precision working array of dimension n. */
1581 /* t will be used to store the break points. */
1582
1583 /* d is a double precision array of dimension n used to store */
1584 /* the Cauchy direction P(x-tg)-x. */
1585
1586 /* xcp is a double precision array of dimension n used to return the */
1587 /* GCP on exit. */
1588
1589 /* m is an integer variable. */
1590 /* On entry m is the maximum number of variable metric corrections */
1591 /* used to define the limited memory matrix. */
1592 /* On exit m is unchanged. */
1593
1594 /* ws, wy, sy, and wt are double precision arrays. */
1595 /* On entry they store information that defines the */
1596 /* limited memory BFGS matrix: */
1597 /* ws(n,m) stores S, a set of s-vectors; */
1598 /* wy(n,m) stores Y, a set of y-vectors; */
1599 /* sy(m,m) stores S'Y; */
1600 /* wt(m,m) stores the */
1601 /* Cholesky factorization of (theta*S'S+LD^(-1)L'). */
1602 /* On exit these arrays are unchanged. */
1603
1604 /* theta is a double precision variable. */
1605 /* On entry theta is the scaling factor specifying B_0 = theta I. */
1606 /* On exit theta is unchanged. */
1607
1608 /* col is an integer variable. */
1609 /* On entry col is the actual number of variable metric */
1610 /* corrections stored so far. */
1611 /* On exit col is unchanged. */
1612
1613 /* head is an integer variable. */
1614 /* On entry head is the location of the first s-vector (or y-vector) */
1615 /* in S (or Y). */
1616 /* On exit col is unchanged. */
1617
1618 /* p is a double precision working array of dimension 2m. */
1619 /* p will be used to store the vector p = W^(T)d. */
1620
1621 /* c is a double precision working array of dimension 2m. */
1622 /* c will be used to store the vector c = W^(T)(xcp-x). */
1623
1624 /* wbp is a double precision working array of dimension 2m. */
1625 /* wbp will be used to store the row of W corresponding */
1626 /* to a breakpoint. */
1627
1628 /* v is a double precision working array of dimension 2m. */
1629
1630 /* nseg is an integer variable. */
1631 /* On exit nseg records the number of quadratic segments explored */
1632 /* in searching for the GCP. */
1633
1634 /* sg and yg are double precision arrays of dimension m. */
1635 /* On entry sg and yg store S'g and Y'g correspondingly. */
1636 /* On exit they are unchanged. */
1637
1638 /* iprint is an INTEGER variable that must be set by the user. */
1639 /* It controls the frequency and type of output generated: */
1640 /* iprint<0 no output is generated; */
1641 /* iprint=0 print only one line at the last iteration; */
1642 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
1643 /* iprint=99 print details of every iteration except n-vectors; */
1644 /* iprint=100 print also the changes of active set and final x; */
1645 /* iprint>100 print details of every iteration including x and g; */
1646 /* When iprint > 0, the file iterate.dat will be created to */
1647 /* summarize the iteration. */
1648
1649 /* sbgnrm is a double precision variable. */
1650 /* On entry sbgnrm is the norm of the projected gradient at x. */
1651 /* On exit sbgnrm is unchanged. */
1652
1653 /* info is an integer variable. */
1654 /* On entry info is 0. */
1655 /* On exit info = 0 for normal return, */
1656 /* = nonzero for abnormal return when the the system */
1657 /* used in routine bmv is singular. */
1658
1659 /* Subprograms called: */
1660
1661 /* L-BFGS-B Library ... hpsolb, bmv. */
1662
1663 /* Linpack ... dscal dcopy, daxpy. */
1664
1665
1666 /* References: */
1667
1668 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
1669 /* memory algorithm for bound constrained optimization'', */
1670 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
1671
1672 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
1673 /* Subroutines for Large Scale Bound Constrained Optimization'' */
1674 /* Tech. Report, NAM-11, EECS Department, Northwestern University, */
1675 /* 1994. */
1676
1677 /* (Postscript files of these papers are available via anonymous */
1678 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
1679
1680 /* * * * */
1681
1682 /* NEOS, November 1994. (Latest revision June 1996.) */
1683 /* Optimization Technology Center. */
1684 /* Argonne National Laboratory and Northwestern University. */
1685 /* Written by */
1686 /* Ciyou Zhu */
1687 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
1688
1689
1690 /* ************ */
1691 /* Check the status of the variables, reset iwhere(i) if necessary; */
1692 /* compute the Cauchy direction d and the breakpoints t; initialize */
1693 /* the derivative f1 and the vector p = W'd (for theta = 1). */
1694 /* Parameter adjustments */
1695 --xcp;
1696 --d__;
1697 --t;
1698 --iwhere;
1699 --iorder;
1700 --g;
1701 --nbd;
1702 --u;
1703 --l;
1704 --x;
1705 --v;
1706 --wbp;
1707 --c__;
1708 --p;
1709 wt_dim1 = *m;
1710 wt_offset = 1 + wt_dim1;
1711 wt -= wt_offset;
1712 sy_dim1 = *m;
1713 sy_offset = 1 + sy_dim1;
1714 sy -= sy_offset;
1715 ws_dim1 = *n;
1716 ws_offset = 1 + ws_dim1;
1717 ws -= ws_offset;
1718 wy_dim1 = *n;
1719 wy_offset = 1 + wy_dim1;
1720 wy -= wy_offset;
1721
1722 /* Function Body */
1723 if (*sbgnrm <= 0.) {
1724 if (*iprint >= 0) {
1725 printf("Subgnorm = 0. GCP = X.\n");
1726 }
1727 dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
1728 return 0;
1729 }
1730 bnded = TRUE_;
1731 nfree = *n + 1;
1732 nbreak = 0;
1733 ibkmin = 0;
1734 bkmin = 0.;
1735 col2 = *col << 1;
1736 f1 = 0.;
1737 if (*iprint >= 99) {
1738 printf("---------------- CAUCHY entered-------------------\n");
1739 }
1740 /* We set p to zero and build it up as we determine d. */
1741 i__1 = col2;
1742 for (i__ = 1; i__ <= i__1; ++i__) {
1743 p[i__] = 0.;
1744 /* L20: */
1745 }
1746 /* In the following loop we determine for each variable its bound */
1747 /* status and its breakpoint, and update p accordingly. */
1748 /* Smallest breakpoint is identified. */
1749 i__1 = *n;
1750 for (i__ = 1; i__ <= i__1; ++i__) {
1751 neggi = -g[i__];
1752 if (iwhere[i__] != 3 && iwhere[i__] != -1) {
1753 /* if x(i) is not a constant and has bounds, */
1754 /* compute the difference between x(i) and its bounds. */
1755 if (nbd[i__] <= 2) {
1756 tl = x[i__] - l[i__];
1757 }
1758 if (nbd[i__] >= 2) {
1759 tu = u[i__] - x[i__];
1760 }
1761 /* If a variable is close enough to a bound */
1762 /* we treat it as at bound. */
1763 xlower = nbd[i__] <= 2 && tl <= 0.;
1764 xupper = nbd[i__] >= 2 && tu <= 0.;
1765 /* reset iwhere(i). */
1766 iwhere[i__] = 0;
1767 if (xlower) {
1768 if (neggi <= 0.) {
1769 iwhere[i__] = 1;
1770 }
1771 } else if (xupper) {
1772 if (neggi >= 0.) {
1773 iwhere[i__] = 2;
1774 }
1775 } else {
1776 if (abs(neggi) <= 0.) {
1777 iwhere[i__] = -3;
1778 }
1779 }
1780 }
1781 pointr = *head;
1782 if (iwhere[i__] != 0 && iwhere[i__] != -1) {
1783 d__[i__] = 0.;
1784 } else {
1785 d__[i__] = neggi;
1786 f1 -= neggi * neggi;
1787 /* calculate p := p - W'e_i* (g_i). */
1788 i__2 = *col;
1789 for (j = 1; j <= i__2; ++j) {
1790 p[j] += wy[i__ + pointr * wy_dim1] * neggi;
1791 p[*col + j] += ws[i__ + pointr * ws_dim1] * neggi;
1792 pointr = pointr % *m + 1;
1793 /* L40: */
1794 }
1795 if (nbd[i__] <= 2 && nbd[i__] != 0 && neggi < 0.) {
1796 /* x(i) + d(i) is bounded; compute t(i). */
1797 ++nbreak;
1798 iorder[nbreak] = i__;
1799 t[nbreak] = tl / (-neggi);
1800 if (nbreak == 1 || t[nbreak] < bkmin) {
1801 bkmin = t[nbreak];
1802 ibkmin = nbreak;
1803 }
1804 } else if (nbd[i__] >= 2 && neggi > 0.) {
1805 /* x(i) + d(i) is bounded; compute t(i). */
1806 ++nbreak;
1807 iorder[nbreak] = i__;
1808 t[nbreak] = tu / neggi;
1809 if (nbreak == 1 || t[nbreak] < bkmin) {
1810 bkmin = t[nbreak];
1811 ibkmin = nbreak;
1812 }
1813 } else {
1814 /* x(i) + d(i) is not bounded. */
1815 --nfree;
1816 iorder[nfree] = i__;
1817 if (abs(neggi) > 0.) {
1818 bnded = FALSE_;
1819 }
1820 }
1821 }
1822 /* L50: */
1823 }
1824 /* The indices of the nonzero components of d are now stored */
1825 /* in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */
1826 /* The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */
1827 if (*theta != 1.) {
1828 /* complete the initialization of p for theta not= one. */
1829 dscal_(col, theta, &p[*col + 1], &c__1);
1830 }
1831 /* Initialize GCP xcp = x. */
1832 dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
1833 if (nbreak == 0 && nfree == *n + 1) {
1834 /* is a zero vector, return with the initial xcp as GCP. */
1835 if (*iprint > 100) {
1836 i__1 = *n;
1837 lbfgsb_printf_vec("Cauchy X", xcp, i__1);
1838 }
1839 return 0;
1840 }
1841 /* Initialize c = W'(xcp - x) = 0. */
1842 i__1 = col2;
1843 for (j = 1; j <= i__1; ++j) {
1844 c__[j] = 0.;
1845 /* L60: */
1846 }
1847 /* Initialize derivative f2. */
1848 f2 = -(*theta) * f1;
1849 f2_org__ = f2;
1850 if (*col > 0) {
1851 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info);
1852 if (*info != 0) {
1853 return 0;
1854 }
1855 f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
1856 }
1857 dtm = -f1 / f2;
1858 tsum = 0.;
1859 *nseg = 1;
1860 if (*iprint >= 99) {
1861 printf("There are %ld breakpoints.\n", nbreak);
1862 }
1863 /* If there are no breakpoints, locate the GCP and return. */
1864 if (nbreak == 0) {
1865 goto L888;
1866 }
1867 nleft = nbreak;
1868 iter = 1;
1869 tj = 0.;
1870 /* ------------------- the beginning of the loop ------------------------- */
1871 L777:
1872 /* Find the next smallest breakpoint; */
1873 /* compute dt = t(nleft) - t(nleft + 1). */
1874 tj0 = tj;
1875 if (iter == 1) {
1876 /* Since we already have the smallest breakpoint we need not do */
1877 /* heapsort yet. Often only one breakpoint is used and the */
1878 /* cost of heapsort is avoided. */
1879 tj = bkmin;
1880 ibp = iorder[ibkmin];
1881 } else {
1882 if (iter == 2) {
1883 /* Replace the already used smallest breakpoint with the */
1884 /* breakpoint numbered nbreak > nlast, before heapsort call. */
1885 if (ibkmin != nbreak) {
1886 t[ibkmin] = t[nbreak];
1887 iorder[ibkmin] = iorder[nbreak];
1888 }
1889 /* Update heap structure of breakpoints */
1890 /* (if iter=2, initialize heap). */
1891 }
1892 i__1 = iter - 2;
1893 hpsolb_(&nleft, &t[1], &iorder[1], &i__1);
1894 tj = t[nleft];
1895 ibp = iorder[nleft];
1896 }
1897 dt = tj - tj0;
1898 if (dt != 0. && *iprint >= 100) {
1899 /*
1900 4010 format ('Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4))
1901 */
1902 printf("Piece %3ld --f1, f2 at start point %11.4g %11.5g\n",
1903 *nseg, f1, f2);
1904 /*
1905 5010 format ('Distance to the next break point = ',1p,d11.4)
1906 */
1907 printf("Distance to the next break point = %11.4g", dt);
1908 /*
1909 6010 format ('Distance to the stationary point = ',1p,d11.4)
1910 */
1911 printf("Distance to the stationary point = %11.4g", dtm);
1912 }
1913 /* If a minimizer is within this interval, locate the GCP and return. */
1914 if (dtm < dt) {
1915 goto L888;
1916 }
1917 /* Otherwise fix one variable and */
1918 /* reset the corresponding component of d to zero. */
1919 tsum += dt;
1920 --nleft;
1921 ++iter;
1922 dibp = d__[ibp];
1923 d__[ibp] = 0.;
1924 if (dibp > 0.) {
1925 zibp = u[ibp] - x[ibp];
1926 xcp[ibp] = u[ibp];
1927 iwhere[ibp] = 2;
1928 } else {
1929 zibp = l[ibp] - x[ibp];
1930 xcp[ibp] = l[ibp];
1931 iwhere[ibp] = 1;
1932 }
1933 if (*iprint >= 100) {
1934 printf("Variable %ld is fixed.\n", ibp);
1935 }
1936 if (nleft == 0 && nbreak == *n) {
1937 /* all n variables are fixed, */
1938 /* return with xcp as GCP. */
1939 dtm = dt;
1940 goto L999;
1941 }
1942 /* Update the derivative information. */
1943 ++(*nseg);
1944 /* Computing 2nd power */
1945 d__1 = dibp;
1946 dibp2 = d__1 * d__1;
1947 /* Update f1 and f2. */
1948 /* temporarily set f1 and f2 for col=0. */
1949 f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp;
1950 f2 -= *theta * dibp2;
1951 if (*col > 0) {
1952 /* update c = c + dt*p. */
1953 daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
1954 /* choose wbp, */
1955 /* the row of W corresponding to the breakpoint encountered. */
1956 pointr = *head;
1957 i__1 = *col;
1958 for (j = 1; j <= i__1; ++j) {
1959 wbp[j] = wy[ibp + pointr * wy_dim1];
1960 wbp[*col + j] = *theta * ws[ibp + pointr * ws_dim1];
1961 pointr = pointr % *m + 1;
1962 /* L70: */
1963 }
1964 /* compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */
1965 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info);
1966 if (*info != 0) {
1967 return 0;
1968 }
1969 wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
1970 wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
1971 wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
1972 /* update p = p - dibp*wbp. */
1973 d__1 = -dibp;
1974 daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
1975 /* complete updating f1 and f2 while col > 0. */
1976 f1 += dibp * wmc;
1977 f2 = f2 + dibp * 2. * wmp - dibp2 * wmw;
1978 }
1979 /* Computing MAX */
1980 d__1 = *epsmch * f2_org__;
1981 f2 = max(d__1,f2);
1982 if (nleft > 0) {
1983 dtm = -f1 / f2;
1984 goto L777;
1985 /* to repeat the loop for unsearched intervals. */
1986 } else if (bnded) {
1987 f1 = 0.;
1988 f2 = 0.;
1989 dtm = 0.;
1990 } else {
1991 dtm = -f1 / f2;
1992 }
1993 /* ------------------- the end of the loop ------------------------------- */
1994 L888:
1995 if (*iprint >= 99) {
1996 printf("\n");
1997 printf("GCP found in this segment\n");
1998 /*
1999 4010 format ('Piece ',i3,' --f1, f2 at start point ',1p,2(1x,d11.4))
2000 */
2001 printf("Piece %3ld --f1, f2 at start point %11.4g %11.4g\n",
2002 *nseg, f1, f2);
2003 /*
2004 6010 format ('Distance to the stationary point = ',1p,d11.4)
2005 */
2006 printf("Distance to the stationary point = %11.4g\n", dtm);
2007 }
2008 if (dtm <= 0.) {
2009 dtm = 0.;
2010 }
2011 tsum += dtm;
2012 /* Move free variables (i.e., the ones w/o breakpoints) and */
2013 /* the variables whose breakpoints haven't been reached. */
2014 daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
2015 L999:
2016 /* Update c = c + dtm*p = W'(x^c - x) */
2017 /* which will be used in computing r = Z'(B(x^c - x) + g). */
2018 if (*col > 0) {
2019 daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
2020 }
2021 if (*iprint > 100) {
2022 i__1 = *n;
2023 lbfgsb_printf_vec("Cauchy X", xcp, i__1);
2024 }
2025 if (*iprint >= 99) {
2026 printf("---------------- exit CAUCHY----------------------\n");
2027 }
2028 return 0;
2029 } /* cauchy_ */
2030
2031 /* ====================== The end of cauchy ============================== */
cmprlb_(integer * n,integer * m,doublereal * x,doublereal * g,doublereal * ws,doublereal * wy,doublereal * sy,doublereal * wt,doublereal * z__,doublereal * r__,doublereal * wa,integer * index,doublereal * theta,integer * col,integer * head,integer * nfree,logical * cnstnd,integer * info)2032 /* Subroutine */ int cmprlb_(integer *n, integer *m, doublereal *x,
2033 doublereal *g, doublereal *ws, doublereal *wy, doublereal *sy,
2034 doublereal *wt, doublereal *z__, doublereal *r__, doublereal *wa,
2035 integer *index, doublereal *theta, integer *col, integer *head,
2036 integer *nfree, logical *cnstnd, integer *info)
2037 {
2038 /* System generated locals */
2039 integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
2040 wt_dim1, wt_offset, i__1, i__2;
2041
2042 /* Local variables */
2043 static integer i__, j, k;
2044 static doublereal a1, a2;
2045 extern /* Subroutine */ int bmv_(integer *, doublereal *, doublereal *,
2046 integer *, doublereal *, doublereal *, integer *);
2047 static integer pointr;
2048
2049 /* ************ */
2050
2051 /* Subroutine cmprlb */
2052
2053 /* This subroutine computes r=-Z'B(xcp-xk)-Z'g by using */
2054 /* wa(2m+1)=W'(xcp-x) from subroutine cauchy. */
2055
2056 /* Subprograms called: */
2057
2058 /* L-BFGS-B Library ... bmv. */
2059
2060
2061 /* * * * */
2062
2063 /* NEOS, November 1994. (Latest revision June 1996.) */
2064 /* Optimization Technology Center. */
2065 /* Argonne National Laboratory and Northwestern University. */
2066 /* Written by */
2067 /* Ciyou Zhu */
2068 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2069
2070
2071 /* ************ */
2072 /* Parameter adjustments */
2073 --index;
2074 --r__;
2075 --z__;
2076 --g;
2077 --x;
2078 --wa;
2079 wt_dim1 = *m;
2080 wt_offset = 1 + wt_dim1;
2081 wt -= wt_offset;
2082 sy_dim1 = *m;
2083 sy_offset = 1 + sy_dim1;
2084 sy -= sy_offset;
2085 wy_dim1 = *n;
2086 wy_offset = 1 + wy_dim1;
2087 wy -= wy_offset;
2088 ws_dim1 = *n;
2089 ws_offset = 1 + ws_dim1;
2090 ws -= ws_offset;
2091
2092 /* Function Body */
2093 if (! (*cnstnd) && *col > 0) {
2094 i__1 = *n;
2095 for (i__ = 1; i__ <= i__1; ++i__) {
2096 r__[i__] = -g[i__];
2097 /* L26: */
2098 }
2099 } else {
2100 i__1 = *nfree;
2101 for (i__ = 1; i__ <= i__1; ++i__) {
2102 k = index[i__];
2103 r__[i__] = -(*theta) * (z__[k] - x[k]) - g[k];
2104 /* L30: */
2105 }
2106 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wa[(*m << 1) + 1], &wa[
2107 1], info);
2108 if (*info != 0) {
2109 *info = -8;
2110 return 0;
2111 }
2112 pointr = *head;
2113 i__1 = *col;
2114 for (j = 1; j <= i__1; ++j) {
2115 a1 = wa[j];
2116 a2 = *theta * wa[*col + j];
2117 i__2 = *nfree;
2118 for (i__ = 1; i__ <= i__2; ++i__) {
2119 k = index[i__];
2120 r__[i__] = r__[i__] + wy[k + pointr * wy_dim1] * a1 + ws[k +
2121 pointr * ws_dim1] * a2;
2122 /* L32: */
2123 }
2124 pointr = pointr % *m + 1;
2125 /* L34: */
2126 }
2127 }
2128 return 0;
2129 } /* cmprlb_ */
2130
2131 /* ======================= The end of cmprlb ============================= */
errclb_(integer * n,integer * m,doublereal * factr,doublereal * l,doublereal * u,integer * nbd,char * task,integer * info,integer * k,ftnlen task_len)2132 /* Subroutine */ int errclb_(integer *n, integer *m, doublereal *factr,
2133 doublereal *l, doublereal *u, integer *nbd, char *task, integer *info,
2134 integer *k, ftnlen task_len)
2135 {
2136 /* System generated locals */
2137 integer i__1;
2138
2139 /* Builtin functions */
2140 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
2141
2142 /* Local variables */
2143 static integer i__;
2144
2145 /* ************ */
2146
2147 /* Subroutine errclb */
2148
2149 /* This subroutine checks the validity of the input data. */
2150
2151
2152 /* * * * */
2153
2154 /* NEOS, November 1994. (Latest revision June 1996.) */
2155 /* Optimization Technology Center. */
2156 /* Argonne National Laboratory and Northwestern University. */
2157 /* Written by */
2158 /* Ciyou Zhu */
2159 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2160
2161
2162 /* ************ */
2163 /* Check the input arguments for errors. */
2164 /* Parameter adjustments */
2165 --nbd;
2166 --u;
2167 --l;
2168
2169 /* Function Body */
2170 if (*n <= 0) {
2171 s_copy(task, "ERROR: N .LE. 0", (ftnlen)60, (ftnlen)(15+1));
2172 }
2173 if (*m <= 0) {
2174 s_copy(task, "ERROR: M .LE. 0", (ftnlen)60, (ftnlen)(15+1));
2175 }
2176 if (*factr < 0.) {
2177 s_copy(task, "ERROR: FACTR .LT. 0", (ftnlen)60, (ftnlen)(19+1));
2178 }
2179 /* Check the validity of the arrays nbd(i), u(i), and l(i). */
2180 i__1 = *n;
2181 for (i__ = 1; i__ <= i__1; ++i__) {
2182 if (nbd[i__] < 0 || nbd[i__] > 3) {
2183 /* return */
2184 s_copy(task, "ERROR: INVALID NBD", (ftnlen)60, (ftnlen)(18+1));
2185 *info = -6;
2186 *k = i__;
2187 }
2188 if (nbd[i__] == 2) {
2189 if (l[i__] > u[i__]) {
2190 /* return */
2191 s_copy(task, "ERROR: NO FEASIBLE SOLUTION", (ftnlen)60, (
2192 ftnlen)(27+1));
2193 *info = -7;
2194 *k = i__;
2195 }
2196 }
2197 /* L10: */
2198 }
2199 return 0;
2200 } /* errclb_ */
2201
2202 /* ======================= The end of errclb ============================= */
formk_(integer * n,integer * nsub,integer * ind,integer * nenter,integer * ileave,integer * indx2,integer * iupdat,logical * updatd,doublereal * wn,doublereal * wn1,integer * m,doublereal * ws,doublereal * wy,doublereal * sy,doublereal * theta,integer * col,integer * head,integer * info)2203 /* Subroutine */ int formk_(integer *n, integer *nsub, integer *ind, integer *
2204 nenter, integer *ileave, integer *indx2, integer *iupdat, logical *
2205 updatd, doublereal *wn, doublereal *wn1, integer *m, doublereal *ws,
2206 doublereal *wy, doublereal *sy, doublereal *theta, integer *col,
2207 integer *head, integer *info)
2208 {
2209 /* System generated locals */
2210 integer wn_dim1, wn_offset, wn1_dim1, wn1_offset, ws_dim1, ws_offset,
2211 wy_dim1, wy_offset, sy_dim1, sy_offset, i__1, i__2, i__3;
2212
2213 /* Local variables */
2214 static integer i__, k, k1, m2, is, js, iy, jy, is1, js1, col2, dend, pend;
2215 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2216 integer *);
2217 static integer upcl;
2218 static doublereal temp1, temp2, temp3, temp4;
2219 /* Subroutine */ int lbfgsb_local_dpofa_(doublereal *, integer *, integer *,
2220 integer *), dcopy_(integer *, doublereal *, integer *, doublereal
2221 *, integer *), dtrsl_(doublereal *, integer *, integer *,
2222 doublereal *, integer *, integer *);
2223 static integer ipntr, jpntr, dbegin, pbegin;
2224
2225 /* ************ */
2226
2227 /* Subroutine formk */
2228
2229 /* This subroutine forms the LEL^T factorization of the indefinite */
2230
2231 /* matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
2232 /* [L_a -R_z theta*S'AA'S ] */
2233 /* where E = [-I 0] */
2234 /* [ 0 I] */
2235 /* The matrix K can be shown to be equal to the matrix M^[-1]N */
2236 /* occurring in section 5.1 of [1], as well as to the matrix */
2237 /* Mbar^[-1] Nbar in section 5.3. */
2238
2239 /* n is an integer variable. */
2240 /* On entry n is the dimension of the problem. */
2241 /* On exit n is unchanged. */
2242
2243 /* nsub is an integer variable */
2244 /* On entry nsub is the number of subspace variables in free set. */
2245 /* On exit nsub is not changed. */
2246
2247 /* ind is an integer array of dimension nsub. */
2248 /* On entry ind specifies the indices of subspace variables. */
2249 /* On exit ind is unchanged. */
2250
2251 /* nenter is an integer variable. */
2252 /* On entry nenter is the number of variables entering the */
2253 /* free set. */
2254 /* On exit nenter is unchanged. */
2255
2256 /* ileave is an integer variable. */
2257 /* On entry indx2(ileave),...,indx2(n) are the variables leaving */
2258 /* the free set. */
2259 /* On exit ileave is unchanged. */
2260
2261 /* indx2 is an integer array of dimension n. */
2262 /* On entry indx2(1),...,indx2(nenter) are the variables entering */
2263 /* the free set, while indx2(ileave),...,indx2(n) are the */
2264 /* variables leaving the free set. */
2265 /* On exit indx2 is unchanged. */
2266
2267 /* iupdat is an integer variable. */
2268 /* On entry iupdat is the total number of BFGS updates made so far. */
2269 /* On exit iupdat is unchanged. */
2270
2271 /* updatd is a logical variable. */
2272 /* On entry 'updatd' is true if the L-BFGS matrix is updatd. */
2273 /* On exit 'updatd' is unchanged. */
2274
2275 /* wn is a double precision array of dimension 2m x 2m. */
2276 /* On entry wn is unspecified. */
2277 /* On exit the upper triangle of wn stores the LEL^T factorization */
2278 /* of the 2*col x 2*col indefinite matrix */
2279 /* [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
2280 /* [L_a -R_z theta*S'AA'S ] */
2281
2282 /* wn1 is a double precision array of dimension 2m x 2m. */
2283 /* On entry wn1 stores the lower triangular part of */
2284 /* [Y' ZZ'Y L_a'+R_z'] */
2285 /* [L_a+R_z S'AA'S ] */
2286 /* in the previous iteration. */
2287 /* On exit wn1 stores the corresponding updated matrices. */
2288 /* The purpose of wn1 is just to store these inner products */
2289 /* so they can be easily updated and inserted into wn. */
2290
2291 /* m is an integer variable. */
2292 /* On entry m is the maximum number of variable metric corrections */
2293 /* used to define the limited memory matrix. */
2294 /* On exit m is unchanged. */
2295
2296 /* ws, wy, sy, and wtyy are double precision arrays; */
2297 /* theta is a double precision variable; */
2298 /* col is an integer variable; */
2299 /* head is an integer variable. */
2300 /* On entry they store the information defining the */
2301 /* limited memory BFGS matrix: */
2302 /* ws(n,m) stores S, a set of s-vectors; */
2303 /* wy(n,m) stores Y, a set of y-vectors; */
2304 /* sy(m,m) stores S'Y; */
2305 /* wtyy(m,m) stores the Cholesky factorization */
2306 /* of (theta*S'S+LD^(-1)L') */
2307 /* theta is the scaling factor specifying B_0 = theta I; */
2308 /* col is the number of variable metric corrections stored; */
2309 /* head is the location of the 1st s- (or y-) vector in S (or Y). */
2310 /* On exit they are unchanged. */
2311
2312 /* info is an integer variable. */
2313 /* On entry info is unspecified. */
2314 /* On exit info = 0 for normal return; */
2315 /* = -1 when the 1st Cholesky factorization failed; */
2316 /* = -2 when the 2st Cholesky factorization failed. */
2317
2318 /* Subprograms called: */
2319
2320 /* Linpack ... dcopy, dpofa, dtrsl. */
2321
2322
2323 /* References: */
2324 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
2325 /* memory algorithm for bound constrained optimization'', */
2326 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
2327
2328 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
2329 /* limited memory FORTRAN code for solving bound constrained */
2330 /* optimization problems'', Tech. Report, NAM-11, EECS Department, */
2331 /* Northwestern University, 1994. */
2332
2333 /* (Postscript files of these papers are available via anonymous */
2334 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
2335
2336 /* * * * */
2337
2338 /* NEOS, November 1994. (Latest revision June 1996.) */
2339 /* Optimization Technology Center. */
2340 /* Argonne National Laboratory and Northwestern University. */
2341 /* Written by */
2342 /* Ciyou Zhu */
2343 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2344
2345
2346 /* ************ */
2347 /* Form the lower triangular part of */
2348 /* WN1 = [Y' ZZ'Y L_a'+R_z'] */
2349 /* [L_a+R_z S'AA'S ] */
2350 /* where L_a is the strictly lower triangular part of S'AA'Y */
2351 /* R_z is the upper triangular part of S'ZZ'Y. */
2352 /* Parameter adjustments */
2353 --indx2;
2354 --ind;
2355 sy_dim1 = *m;
2356 sy_offset = 1 + sy_dim1;
2357 sy -= sy_offset;
2358 wy_dim1 = *n;
2359 wy_offset = 1 + wy_dim1;
2360 wy -= wy_offset;
2361 ws_dim1 = *n;
2362 ws_offset = 1 + ws_dim1;
2363 ws -= ws_offset;
2364 wn1_dim1 = 2 * *m;
2365 wn1_offset = 1 + wn1_dim1;
2366 wn1 -= wn1_offset;
2367 wn_dim1 = 2 * *m;
2368 wn_offset = 1 + wn_dim1;
2369 wn -= wn_offset;
2370
2371 /* Function Body */
2372 if (*updatd) {
2373 if (*iupdat > *m) {
2374 /* shift old part of WN1. */
2375 i__1 = *m - 1;
2376 for (jy = 1; jy <= i__1; ++jy) {
2377 js = *m + jy;
2378 i__2 = *m - jy;
2379 dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[
2380 jy + jy * wn1_dim1], &c__1);
2381 i__2 = *m - jy;
2382 dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[
2383 js + js * wn1_dim1], &c__1);
2384 i__2 = *m - 1;
2385 dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[
2386 *m + 1 + jy * wn1_dim1], &c__1);
2387 /* L10: */
2388 }
2389 }
2390 /* put new rows in blocks (1,1), (2,1) and (2,2). */
2391 pbegin = 1;
2392 pend = *nsub;
2393 dbegin = *nsub + 1;
2394 dend = *n;
2395 iy = *col;
2396 is = *m + *col;
2397 ipntr = *head + *col - 1;
2398 if (ipntr > *m) {
2399 ipntr -= *m;
2400 }
2401 jpntr = *head;
2402 i__1 = *col;
2403 for (jy = 1; jy <= i__1; ++jy) {
2404 js = *m + jy;
2405 temp1 = 0.;
2406 temp2 = 0.;
2407 temp3 = 0.;
2408 /* compute element jy of row 'col' of Y'ZZ'Y */
2409 i__2 = pend;
2410 for (k = pbegin; k <= i__2; ++k) {
2411 k1 = ind[k];
2412 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
2413 /* L15: */
2414 }
2415 /* compute elements jy of row 'col' of L_a and S'AA'S */
2416 i__2 = dend;
2417 for (k = dbegin; k <= i__2; ++k) {
2418 k1 = ind[k];
2419 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
2420 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
2421 /* L16: */
2422 }
2423 wn1[iy + jy * wn1_dim1] = temp1;
2424 wn1[is + js * wn1_dim1] = temp2;
2425 wn1[is + jy * wn1_dim1] = temp3;
2426 jpntr = jpntr % *m + 1;
2427 /* L20: */
2428 }
2429 /* put new column in block (2,1). */
2430 jy = *col;
2431 jpntr = *head + *col - 1;
2432 if (jpntr > *m) {
2433 jpntr -= *m;
2434 }
2435 ipntr = *head;
2436 i__1 = *col;
2437 for (i__ = 1; i__ <= i__1; ++i__) {
2438 is = *m + i__;
2439 temp3 = 0.;
2440 /* compute element i of column 'col' of R_z */
2441 i__2 = pend;
2442 for (k = pbegin; k <= i__2; ++k) {
2443 k1 = ind[k];
2444 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
2445 /* L25: */
2446 }
2447 ipntr = ipntr % *m + 1;
2448 wn1[is + jy * wn1_dim1] = temp3;
2449 /* L30: */
2450 }
2451 upcl = *col - 1;
2452 } else {
2453 upcl = *col;
2454 }
2455 /* modify the old parts in blocks (1,1) and (2,2) due to changes */
2456 /* in the set of free variables. */
2457 ipntr = *head;
2458 i__1 = upcl;
2459 for (iy = 1; iy <= i__1; ++iy) {
2460 is = *m + iy;
2461 jpntr = *head;
2462 i__2 = iy;
2463 for (jy = 1; jy <= i__2; ++jy) {
2464 js = *m + jy;
2465 temp1 = 0.;
2466 temp2 = 0.;
2467 temp3 = 0.;
2468 temp4 = 0.;
2469 i__3 = *nenter;
2470 for (k = 1; k <= i__3; ++k) {
2471 k1 = indx2[k];
2472 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
2473 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
2474 /* L35: */
2475 }
2476 i__3 = *n;
2477 for (k = *ileave; k <= i__3; ++k) {
2478 k1 = indx2[k];
2479 temp3 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
2480 temp4 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
2481 /* L36: */
2482 }
2483 wn1[iy + jy * wn1_dim1] = wn1[iy + jy * wn1_dim1] + temp1 - temp3;
2484 wn1[is + js * wn1_dim1] = wn1[is + js * wn1_dim1] - temp2 + temp4;
2485 jpntr = jpntr % *m + 1;
2486 /* L40: */
2487 }
2488 ipntr = ipntr % *m + 1;
2489 /* L45: */
2490 }
2491 /* modify the old parts in block (2,1). */
2492 ipntr = *head;
2493 i__1 = *m + upcl;
2494 for (is = *m + 1; is <= i__1; ++is) {
2495 jpntr = *head;
2496 i__2 = upcl;
2497 for (jy = 1; jy <= i__2; ++jy) {
2498 temp1 = 0.;
2499 temp3 = 0.;
2500 i__3 = *nenter;
2501 for (k = 1; k <= i__3; ++k) {
2502 k1 = indx2[k];
2503 temp1 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
2504 /* L50: */
2505 }
2506 i__3 = *n;
2507 for (k = *ileave; k <= i__3; ++k) {
2508 k1 = indx2[k];
2509 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
2510 /* L51: */
2511 }
2512 if (is <= jy + *m) {
2513 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] + temp1 -
2514 temp3;
2515 } else {
2516 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] - temp1 +
2517 temp3;
2518 }
2519 jpntr = jpntr % *m + 1;
2520 /* L55: */
2521 }
2522 ipntr = ipntr % *m + 1;
2523 /* L60: */
2524 }
2525 /* Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] */
2526 /* [-L_a +R_z S'AA'S*theta] */
2527 m2 = *m << 1;
2528 i__1 = *col;
2529 for (iy = 1; iy <= i__1; ++iy) {
2530 is = *col + iy;
2531 is1 = *m + iy;
2532 i__2 = iy;
2533 for (jy = 1; jy <= i__2; ++jy) {
2534 js = *col + jy;
2535 js1 = *m + jy;
2536 wn[jy + iy * wn_dim1] = wn1[iy + jy * wn1_dim1] / *theta;
2537 wn[js + is * wn_dim1] = wn1[is1 + js1 * wn1_dim1] * *theta;
2538 /* L65: */
2539 }
2540 i__2 = iy - 1;
2541 for (jy = 1; jy <= i__2; ++jy) {
2542 wn[jy + is * wn_dim1] = -wn1[is1 + jy * wn1_dim1];
2543 /* L66: */
2544 }
2545 i__2 = *col;
2546 for (jy = iy; jy <= i__2; ++jy) {
2547 wn[jy + is * wn_dim1] = wn1[is1 + jy * wn1_dim1];
2548 /* L67: */
2549 }
2550 wn[iy + iy * wn_dim1] += sy[iy + iy * sy_dim1];
2551 /* L70: */
2552 }
2553 /* Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] */
2554 /* [(-L_a +R_z)L'^-1 S'AA'S*theta ] */
2555 /* first Cholesky factor (1,1) block of wn to get LL' */
2556 /* with L' stored in the upper triangle of wn. */
2557 lbfgsb_local_dpofa_(&wn[wn_offset], &m2, col, info);
2558 if (*info != 0) {
2559 *info = -1;
2560 return 0;
2561 }
2562 /* then form L^-1(-L_a'+R_z') in the (1,2) block. */
2563 col2 = *col << 1;
2564 i__1 = col2;
2565 for (js = *col + 1; js <= i__1; ++js) {
2566 dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info);
2567 /* L71: */
2568 }
2569 /* Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the */
2570 /* upper triangle of (2,2) block of wn. */
2571 i__1 = col2;
2572 for (is = *col + 1; is <= i__1; ++is) {
2573 i__2 = col2;
2574 for (js = is; js <= i__2; ++js) {
2575 wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], &c__1,
2576 &wn[js * wn_dim1 + 1], &c__1);
2577 /* L74: */
2578 }
2579 /* L72: */
2580 }
2581 /* Cholesky factorization of (2,2) block of wn. */
2582 lbfgsb_local_dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info);
2583 if (*info != 0) {
2584 *info = -2;
2585 return 0;
2586 }
2587 return 0;
2588 } /* formk_ */
2589
2590 /* ======================= The end of formk ============================== */
formt_(integer * m,doublereal * wt,doublereal * sy,doublereal * ss,integer * col,doublereal * theta,integer * info)2591 /* Subroutine */ int formt_(integer *m, doublereal *wt, doublereal *sy,
2592 doublereal *ss, integer *col, doublereal *theta, integer *info)
2593 {
2594 /* System generated locals */
2595 integer wt_dim1, wt_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, i__1,
2596 i__2, i__3;
2597
2598 /* Local variables */
2599 static integer i__, j, k, k1;
2600 static doublereal ddum;
2601 /* Subroutine */ int lbfgsb_local_dpofa_(doublereal *, integer *, integer *,
2602 integer *);
2603
2604 /* ************ */
2605
2606 /* Subroutine formt */
2607
2608 /* This subroutine forms the upper half of the pos. def. and symm. */
2609 /* T = theta*SS + L*D^(-1)*L', stores T in the upper triangle */
2610 /* of the array wt, and performs the Cholesky factorization of T */
2611 /* to produce J*J', with J' stored in the upper triangle of wt. */
2612
2613 /* Subprograms called: */
2614
2615 /* Linpack ... dpofa. */
2616
2617
2618 /* * * * */
2619
2620 /* NEOS, November 1994. (Latest revision June 1996.) */
2621 /* Optimization Technology Center. */
2622 /* Argonne National Laboratory and Northwestern University. */
2623 /* Written by */
2624 /* Ciyou Zhu */
2625 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2626
2627
2628 /* ************ */
2629 /* Form the upper half of T = theta*SS + L*D^(-1)*L', */
2630 /* store T in the upper triangle of the array wt. */
2631 /* Parameter adjustments */
2632 ss_dim1 = *m;
2633 ss_offset = 1 + ss_dim1;
2634 ss -= ss_offset;
2635 sy_dim1 = *m;
2636 sy_offset = 1 + sy_dim1;
2637 sy -= sy_offset;
2638 wt_dim1 = *m;
2639 wt_offset = 1 + wt_dim1;
2640 wt -= wt_offset;
2641
2642 /* Function Body */
2643 i__1 = *col;
2644 for (j = 1; j <= i__1; ++j) {
2645 wt[j * wt_dim1 + 1] = *theta * ss[j * ss_dim1 + 1];
2646 /* L52: */
2647 }
2648 i__1 = *col;
2649 for (i__ = 2; i__ <= i__1; ++i__) {
2650 i__2 = *col;
2651 for (j = i__; j <= i__2; ++j) {
2652 k1 = min(i__,j) - 1;
2653 ddum = 0.;
2654 i__3 = k1;
2655 for (k = 1; k <= i__3; ++k) {
2656 ddum += sy[i__ + k * sy_dim1] * sy[j + k * sy_dim1] / sy[k +
2657 k * sy_dim1];
2658 /* L53: */
2659 }
2660 wt[i__ + j * wt_dim1] = ddum + *theta * ss[i__ + j * ss_dim1];
2661 /* L54: */
2662 }
2663 /* L55: */
2664 }
2665 /* Cholesky factorize T to J*J' with */
2666 /* J' stored in the upper triangle of wt. */
2667 lbfgsb_local_dpofa_(&wt[wt_offset], m, col, info);
2668 if (*info != 0) {
2669 *info = -3;
2670 }
2671 return 0;
2672 } /* formt_ */
2673
2674 /* ======================= The end of formt ============================== */
freev_(integer * n,integer * nfree,integer * index,integer * nenter,integer * ileave,integer * indx2,integer * iwhere,logical * wrk,logical * updatd,logical * cnstnd,integer * iprint,integer * iter)2675 /* Subroutine */ int freev_(integer *n, integer *nfree, integer *index,
2676 integer *nenter, integer *ileave, integer *indx2, integer *iwhere,
2677 logical *wrk, logical *updatd, logical *cnstnd, integer *iprint,
2678 integer *iter)
2679 {
2680 /* System generated locals */
2681 integer i__1;
2682
2683 /* Local variables */
2684 static integer i__, k, iact;
2685
2686 /* ************ */
2687
2688 /* Subroutine freev */
2689
2690 /* This subroutine counts the entering and leaving variables when */
2691 /* iter > 0, and finds the index set of free and active variables */
2692 /* at the GCP. */
2693
2694 /* cnstnd is a logical variable indicating whether bounds are present */
2695
2696 /* index is an integer array of dimension n */
2697 /* for i=1,...,nfree, index(i) are the indices of free variables */
2698 /* for i=nfree+1,...,n, index(i) are the indices of bound variables */
2699 /* On entry after the first iteration, index gives */
2700 /* the free variables at the previous iteration. */
2701 /* On exit it gives the free variables based on the determination */
2702 /* in cauchy using the array iwhere. */
2703
2704 /* indx2 is an integer array of dimension n */
2705 /* On entry indx2 is unspecified. */
2706 /* On exit with iter>0, indx2 indicates which variables */
2707 /* have changed status since the previous iteration. */
2708 /* For i= 1,...,nenter, indx2(i) have changed from bound to free. */
2709 /* For i= ileave+1,...,n, indx2(i) have changed from free to bound. */
2710
2711
2712 /* * * * */
2713
2714 /* NEOS, November 1994. (Latest revision June 1996.) */
2715 /* Optimization Technology Center. */
2716 /* Argonne National Laboratory and Northwestern University. */
2717 /* Written by */
2718 /* Ciyou Zhu */
2719 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2720
2721
2722 /* ************ */
2723 /* Parameter adjustments */
2724 --iwhere;
2725 --indx2;
2726 --index;
2727
2728 /* Function Body */
2729 *nenter = 0;
2730 *ileave = *n + 1;
2731 if (*iter > 0 && *cnstnd) {
2732 /* count the entering and leaving variables. */
2733 i__1 = *nfree;
2734 for (i__ = 1; i__ <= i__1; ++i__) {
2735 k = index[i__];
2736 /* write(6,*) ' k = index(i) ', k */
2737 /* write(6,*) ' index = ', i */
2738 if (iwhere[k] > 0) {
2739 --(*ileave);
2740 indx2[*ileave] = k;
2741 if (*iprint >= 100) {
2742 printf("Variable %ld leaves the set of free variables\n", k);
2743 }
2744 }
2745 /* L20: */
2746 }
2747 i__1 = *n;
2748 for (i__ = *nfree + 1; i__ <= i__1; ++i__) {
2749 k = index[i__];
2750 if (iwhere[k] <= 0) {
2751 ++(*nenter);
2752 indx2[*nenter] = k;
2753 if (*iprint >= 100) {
2754 printf("Variable %ld enters the set of free variables\n", k);
2755 }
2756 }
2757 /* L22: */
2758 }
2759 if (*iprint >= 99) {
2760 i__1 = *n + 1 - *ileave;
2761 printf("%ld variables leave; %ld variables enter\n", i__1, *nenter);
2762 }
2763 }
2764 *wrk = *ileave < *n + 1 || *nenter > 0 || *updatd;
2765 /* Find the index set of free and active variables at the GCP. */
2766 *nfree = 0;
2767 iact = *n + 1;
2768 i__1 = *n;
2769 for (i__ = 1; i__ <= i__1; ++i__) {
2770 if (iwhere[i__] <= 0) {
2771 ++(*nfree);
2772 index[*nfree] = i__;
2773 } else {
2774 --iact;
2775 index[iact] = i__;
2776 }
2777 /* L24: */
2778 }
2779 if (*iprint >= 99) {
2780 i__1 = *iter + 1;
2781 printf("%ld variables are free at GCP %ld\n", *nfree, i__1);
2782 }
2783 return 0;
2784 } /* freev_ */
2785
2786 /* ======================= The end of freev ============================== */
hpsolb_(integer * n,doublereal * t,integer * iorder,integer * iheap)2787 /* Subroutine */ int hpsolb_(integer *n, doublereal *t, integer *iorder,
2788 integer *iheap)
2789 {
2790 /* System generated locals */
2791 integer i__1;
2792
2793 /* Local variables */
2794 static integer i__, j, k;
2795 static doublereal out, ddum;
2796 static integer indxin, indxou;
2797
2798 /* ************ */
2799
2800 /* Subroutine hpsolb */
2801
2802 /* This subroutine sorts out the least element of t, and puts the */
2803 /* remaining elements of t in a heap. */
2804
2805 /* n is an integer variable. */
2806 /* On entry n is the dimension of the arrays t and iorder. */
2807 /* On exit n is unchanged. */
2808
2809 /* t is a double precision array of dimension n. */
2810 /* On entry t stores the elements to be sorted, */
2811 /* On exit t(n) stores the least elements of t, and t(1) to t(n-1) */
2812 /* stores the remaining elements in the form of a heap. */
2813
2814 /* iorder is an integer array of dimension n. */
2815 /* On entry iorder(i) is the index of t(i). */
2816 /* On exit iorder(i) is still the index of t(i), but iorder may be */
2817 /* permuted in accordance with t. */
2818
2819 /* iheap is an integer variable specifying the task. */
2820 /* On entry iheap should be set as follows: */
2821 /* iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, */
2822 /* iheap .ne. 0 if otherwise. */
2823 /* On exit iheap is unchanged. */
2824
2825
2826 /* References: */
2827 /* Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. */
2828
2829 /* * * * */
2830
2831 /* NEOS, November 1994. (Latest revision June 1996.) */
2832 /* Optimization Technology Center. */
2833 /* Argonne National Laboratory and Northwestern University. */
2834 /* Written by */
2835 /* Ciyou Zhu */
2836 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2837
2838 /* ************ */
2839 /* Parameter adjustments */
2840 --iorder;
2841 --t;
2842
2843 /* Function Body */
2844 if (*iheap == 0) {
2845 /* Rearrange the elements t(1) to t(n) to form a heap. */
2846 i__1 = *n;
2847 for (k = 2; k <= i__1; ++k) {
2848 ddum = t[k];
2849 indxin = iorder[k];
2850 /* Add ddum to the heap. */
2851 i__ = k;
2852 L10:
2853 if (i__ > 1) {
2854 j = i__ / 2;
2855 if (ddum < t[j]) {
2856 t[i__] = t[j];
2857 iorder[i__] = iorder[j];
2858 i__ = j;
2859 goto L10;
2860 }
2861 }
2862 t[i__] = ddum;
2863 iorder[i__] = indxin;
2864 /* L20: */
2865 }
2866 }
2867 /* Assign to 'out' the value of t(1), the least member of the heap, */
2868 /* and rearrange the remaining members to form a heap as */
2869 /* elements 1 to n-1 of t. */
2870 if (*n > 1) {
2871 i__ = 1;
2872 out = t[1];
2873 indxou = iorder[1];
2874 ddum = t[*n];
2875 indxin = iorder[*n];
2876 /* Restore the heap */
2877 L30:
2878 j = i__ + i__;
2879 if (j <= *n - 1) {
2880 if (t[j + 1] < t[j]) {
2881 ++j;
2882 }
2883 if (t[j] < ddum) {
2884 t[i__] = t[j];
2885 iorder[i__] = iorder[j];
2886 i__ = j;
2887 goto L30;
2888 }
2889 }
2890 t[i__] = ddum;
2891 iorder[i__] = indxin;
2892 /* Put the least member in t(n). */
2893 t[*n] = out;
2894 iorder[*n] = indxou;
2895 }
2896 return 0;
2897 } /* hpsolb_ */
2898
2899 /* ====================== The end of hpsolb ============================== */
lnsrlb_(integer * n,doublereal * l,doublereal * u,integer * nbd,doublereal * x,doublereal * f,doublereal * fold,doublereal * gd,doublereal * gdold,doublereal * g,doublereal * d__,doublereal * r__,doublereal * t,doublereal * z__,doublereal * stp,doublereal * dnorm,doublereal * dtd,doublereal * xstep,doublereal * stpmx,integer * iter,integer * ifun,integer * iback,integer * nfgv,integer * info,char * task,logical * boxed,logical * cnstnd,char * csave,integer * isave,doublereal * dsave,ftnlen task_len,ftnlen csave_len)2900 /* Subroutine */ int lnsrlb_(integer *n, doublereal *l, doublereal *u,
2901 integer *nbd, doublereal *x, doublereal *f, doublereal *fold,
2902 doublereal *gd, doublereal *gdold, doublereal *g, doublereal *d__,
2903 doublereal *r__, doublereal *t, doublereal *z__, doublereal *stp,
2904 doublereal *dnorm, doublereal *dtd, doublereal *xstep, doublereal *
2905 stpmx, integer *iter, integer *ifun, integer *iback, integer *nfgv,
2906 integer *info, char *task, logical *boxed, logical *cnstnd, char *
2907 csave, integer *isave, doublereal *dsave, ftnlen task_len, ftnlen
2908 csave_len)
2909 {
2910 /* System generated locals */
2911 integer i__1;
2912 doublereal d__1;
2913
2914 /* Builtin functions */
2915 integer s_cmp(char *, char *, ftnlen, ftnlen);
2916 double sqrt(doublereal);
2917 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
2918
2919 /* Local variables */
2920 static integer i__;
2921 static doublereal a1, a2;
2922 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2923 integer *);
2924 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
2925 doublereal *, integer *), dcsrch_(doublereal *, doublereal *,
2926 doublereal *, doublereal *, doublereal *, doublereal *,
2927 doublereal *, doublereal *, char *, integer *, doublereal *,
2928 ftnlen);
2929
2930 /* ********** */
2931
2932 /* Subroutine lnsrlb */
2933
2934 /* This subroutine calls subroutine dcsrch from the Minpack2 library */
2935 /* to perform the line search. Subroutine dscrch is safeguarded so */
2936 /* that all trial points lie within the feasible region. */
2937
2938 /* Subprograms called: */
2939
2940 /* Minpack2 Library ... dcsrch. */
2941
2942 /* Linpack ... dtrsl, ddot. */
2943
2944
2945 /* * * * */
2946
2947 /* NEOS, November 1994. (Latest revision June 1996.) */
2948 /* Optimization Technology Center. */
2949 /* Argonne National Laboratory and Northwestern University. */
2950 /* Written by */
2951 /* Ciyou Zhu */
2952 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
2953
2954
2955 /* ********** */
2956 /* Parameter adjustments */
2957 --z__;
2958 --t;
2959 --r__;
2960 --d__;
2961 --g;
2962 --x;
2963 --nbd;
2964 --u;
2965 --l;
2966 --isave;
2967 --dsave;
2968
2969 /* Function Body */
2970 if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
2971 goto L556;
2972 }
2973 *dtd = ddot_(n, &d__[1], &c__1, &d__[1], &c__1);
2974 *dnorm = sqrt(*dtd);
2975 /* Determine the maximum step length. */
2976 *stpmx = 1e10;
2977 if (*cnstnd) {
2978 if (*iter == 0) {
2979 *stpmx = 1.;
2980 } else {
2981 i__1 = *n;
2982 for (i__ = 1; i__ <= i__1; ++i__) {
2983 a1 = d__[i__];
2984 if (nbd[i__] != 0) {
2985 if (a1 < 0. && nbd[i__] <= 2) {
2986 a2 = l[i__] - x[i__];
2987 if (a2 >= 0.) {
2988 *stpmx = 0.;
2989 } else if (a1 * *stpmx < a2) {
2990 *stpmx = a2 / a1;
2991 }
2992 } else if (a1 > 0. && nbd[i__] >= 2) {
2993 a2 = u[i__] - x[i__];
2994 if (a2 <= 0.) {
2995 *stpmx = 0.;
2996 } else if (a1 * *stpmx > a2) {
2997 *stpmx = a2 / a1;
2998 }
2999 }
3000 }
3001 /* L43: */
3002 }
3003 }
3004 }
3005 if (*iter == 0 && ! (*boxed)) {
3006 /* Computing MIN */
3007 d__1 = 1. / *dnorm;
3008 *stp = min(d__1,*stpmx);
3009 } else {
3010 *stp = 1.;
3011 }
3012 dcopy_(n, &x[1], &c__1, &t[1], &c__1);
3013 dcopy_(n, &g[1], &c__1, &r__[1], &c__1);
3014 *fold = *f;
3015 *ifun = 0;
3016 *iback = 0;
3017 s_copy(csave, "START", (ftnlen)60, (ftnlen)(5+1));
3018 L556:
3019 *gd = ddot_(n, &g[1], &c__1, &d__[1], &c__1);
3020 if (*ifun == 0) {
3021 *gdold = *gd;
3022 if (*gd >= 0.) {
3023 /* the directional derivative >=0. */
3024 /* Line search is impossible. */
3025 /*printf(" ascent direction in projection gd = %lf", *gd );*/
3026 *info = -4;
3027 return 0;
3028 }
3029 }
3030 dcsrch_(f, gd, stp, &c_b280, &c_b281, &c_b282, &c_b9, stpmx, csave, &
3031 isave[1], &dsave[1], (ftnlen)60);
3032 *xstep = *stp * *dnorm;
3033 if (s_cmp(csave, "CONV", (ftnlen)4, (ftnlen)4) != 0 && s_cmp(csave, "WARN"
3034 , (ftnlen)4, (ftnlen)4) != 0) {
3035 s_copy(task, "FG_LNSRCH", (ftnlen)60, (ftnlen)(9+1));
3036 ++(*ifun);
3037 ++(*nfgv);
3038 *iback = *ifun - 1;
3039 if (*stp == 1.) {
3040 dcopy_(n, &z__[1], &c__1, &x[1], &c__1);
3041 } else {
3042 i__1 = *n;
3043 for (i__ = 1; i__ <= i__1; ++i__) {
3044 x[i__] = *stp * d__[i__] + t[i__];
3045 /* L41: */
3046 }
3047 }
3048 } else {
3049 s_copy(task, "NEW_X", (ftnlen)60, (ftnlen)(5+1));
3050 }
3051 return 0;
3052 } /* lnsrlb_ */
3053
3054 /* ======================= The end of lnsrlb ============================= */
matupd_(integer * n,integer * m,doublereal * ws,doublereal * wy,doublereal * sy,doublereal * ss,doublereal * d__,doublereal * r__,integer * itail,integer * iupdat,integer * col,integer * head,doublereal * theta,doublereal * rr,doublereal * dr,doublereal * stp,doublereal * dtd)3055 /* Subroutine */ int matupd_(integer *n, integer *m, doublereal *ws,
3056 doublereal *wy, doublereal *sy, doublereal *ss, doublereal *d__,
3057 doublereal *r__, integer *itail, integer *iupdat, integer *col,
3058 integer *head, doublereal *theta, doublereal *rr, doublereal *dr,
3059 doublereal *stp, doublereal *dtd)
3060 {
3061 /* System generated locals */
3062 integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
3063 ss_dim1, ss_offset, i__1, i__2;
3064
3065 /* Local variables */
3066 static integer j;
3067 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
3068 integer *);
3069 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
3070 doublereal *, integer *);
3071 static integer pointr;
3072
3073 /* ************ */
3074
3075 /* Subroutine matupd */
3076
3077 /* This subroutine updates matrices WS and WY, and forms the */
3078 /* middle matrix in B. */
3079
3080 /* Subprograms called: */
3081
3082 /* Linpack ... dcopy, ddot. */
3083
3084
3085 /* * * * */
3086
3087 /* NEOS, November 1994. (Latest revision June 1996.) */
3088 /* Optimization Technology Center. */
3089 /* Argonne National Laboratory and Northwestern University. */
3090 /* Written by */
3091 /* Ciyou Zhu */
3092 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3093
3094
3095 /* ************ */
3096 /* Set pointers for matrices WS and WY. */
3097 /* Parameter adjustments */
3098 --r__;
3099 --d__;
3100 ss_dim1 = *m;
3101 ss_offset = 1 + ss_dim1;
3102 ss -= ss_offset;
3103 sy_dim1 = *m;
3104 sy_offset = 1 + sy_dim1;
3105 sy -= sy_offset;
3106 wy_dim1 = *n;
3107 wy_offset = 1 + wy_dim1;
3108 wy -= wy_offset;
3109 ws_dim1 = *n;
3110 ws_offset = 1 + ws_dim1;
3111 ws -= ws_offset;
3112
3113 /* Function Body */
3114 if (*iupdat <= *m) {
3115 *col = *iupdat;
3116 *itail = (*head + *iupdat - 2) % *m + 1;
3117 } else {
3118 *itail = *itail % *m + 1;
3119 *head = *head % *m + 1;
3120 }
3121 /* Update matrices WS and WY. */
3122 dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1);
3123 dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1);
3124 /* Set theta=yy/ys. */
3125 *theta = *rr / *dr;
3126 /* Form the middle matrix in B. */
3127 /* update the upper triangle of SS, */
3128 /* and the lower triangle of SY: */
3129 if (*iupdat > *m) {
3130 /* move old information */
3131 i__1 = *col - 1;
3132 for (j = 1; j <= i__1; ++j) {
3133 dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1]
3134 , &c__1);
3135 i__2 = *col - j;
3136 dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j *
3137 sy_dim1], &c__1);
3138 /* L50: */
3139 }
3140 }
3141 /* add new information: the last row of SY */
3142 /* and the last column of SS: */
3143 pointr = *head;
3144 i__1 = *col - 1;
3145 for (j = 1; j <= i__1; ++j) {
3146 sy[*col + j * sy_dim1] = ddot_(n, &d__[1], &c__1, &wy[pointr *
3147 wy_dim1 + 1], &c__1);
3148 ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, &
3149 d__[1], &c__1);
3150 pointr = pointr % *m + 1;
3151 /* L51: */
3152 }
3153 if (*stp == 1.) {
3154 ss[*col + *col * ss_dim1] = *dtd;
3155 } else {
3156 ss[*col + *col * ss_dim1] = *stp * *stp * *dtd;
3157 }
3158 sy[*col + *col * sy_dim1] = *dr;
3159 return 0;
3160 } /* matupd_ */
3161
3162 /* ======================= The end of matupd ============================= */
prn1lb_(integer * n,integer * m,doublereal * l,doublereal * u,doublereal * x,integer * iprint,integer * itfile,doublereal * epsmch)3163 /* Subroutine */ int prn1lb_(integer *n, integer *m, doublereal *l,
3164 doublereal *u, doublereal *x, integer *iprint, integer *itfile,
3165 doublereal *epsmch)
3166 {
3167 /* ************ */
3168
3169 /* Subroutine prn1lb */
3170
3171 /* This subroutine prints the input data, initial point, upper and */
3172 /* lower bounds of each variable, machine precision, as well as */
3173 /* the headings of the output. */
3174
3175
3176 /* * * * */
3177
3178 /* NEOS, November 1994. (Latest revision June 1996.) */
3179 /* Optimization Technology Center. */
3180 /* Argonne National Laboratory and Northwestern University. */
3181 /* Written by */
3182 /* Ciyou Zhu */
3183 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3184
3185
3186 /* ************ */
3187 /* Parameter adjustments */
3188 --x;
3189 --u;
3190 --l;
3191
3192 /* Function Body */
3193 if (*iprint >= 0) {
3194 /*
3195 7001 format ('RUNNING THE L-BFGS-B CODE',/,/,
3196 + ' * * *',/,/,
3197 + 'Machine precision =',1p,d10.3)
3198 */
3199 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3200 fprintf(0,
3201 "RUNNING THE L-BFGS-B CODE\n"
3202 "\n"
3203 " * * *\n"
3204 "\n"
3205 "Machine precision = %10.3g\n", *epsmch);
3206 #endif
3207 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3208 fprintf(0, "N = %ld M = %ld\n", *n, *m);
3209 #endif
3210 if (*iprint >= 1) {
3211 /*
3212 2001 format ('RUNNING THE L-BFGS-B CODE',/,/,
3213 + 'it = iteration number',/,
3214 + 'nf = number of function evaluations',/,
3215 + 'nint = number of segments explored during the Cauchy search',/,
3216 + 'nact = number of active bounds at the generalized Cauchy point'
3217 + ,/,
3218 + 'sub = manner in which the subspace minimization terminated:'
3219 + ,/,' con = converged, bnd = a bound was reached',/,
3220 + 'itls = number of iterations performed in the line search',/,
3221 + 'stepl = step length used',/,
3222 + 'tstep = norm of the displacement (total step)',/,
3223 + 'projg = norm of the projected gradient',/,
3224 + 'f = function value',/,/,
3225 + ' * * *',/,/,
3226 + 'Machine precision =',1p,d10.3)
3227 */
3228 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3229 fprintf(0,
3230 "RUNNING THE L-BFGS-B CODE\n"
3231 "\n"
3232 "it = iteration number\n"
3233 "nf = number of function evaluations\n"
3234 "nint = number of segments explored during the Cauchy search\n"
3235 "nact = number of active bounds at the generalized Cauchy point\n"
3236 "sub = manner in which the subspace minimization terminated:\n"
3237 " con = converged, bnd = a bound was reached\n"
3238 "itls = number of iterations performed in the line search\n"
3239 "stepl = step length used\n"
3240 "tstep = norm of the displacement (total step)\n"
3241 "projg = norm of the projected gradient\n"
3242 "f = function value\n"
3243 "\n"
3244 " * * *\n"
3245 "\n"
3246 "Machine precision = %10.3g\n", *epsmch);
3247 #endif
3248 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3249 fprintf(0, "N = %ld M = %ld\n", *n, *m);
3250 #endif
3251 /*
3252 9001 format (/,3x,'it',3x,'nf',2x,'nint',2x,'nact',2x,'sub',2x,'itls',
3253 + 2x,'stepl',4x,'tstep',5x,'projg',8x,'f')
3254 */
3255 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3256 fprintf(0, " it nf nint nact sub itls stepl tstep projg f\n");
3257 #endif
3258 if (*iprint > 100) {
3259 lbfgsb_printf_vec("L", l, *n);
3260 lbfgsb_printf_vec("X0", x, *n);
3261 lbfgsb_printf_vec("U", u, *n);
3262 }
3263 }
3264 }
3265 return 0;
3266 } /* prn1lb_ */
3267
3268 /* ======================= The end of prn1lb ============================= */
prn2lb_(integer * n,doublereal * x,doublereal * f,doublereal * g,integer * iprint,integer * itfile,integer * iter,integer * nfgv,integer * nact,doublereal * sbgnrm,integer * nseg,char * word,integer * iword,integer * iback,doublereal * stp,doublereal * xstep,ftnlen word_len)3269 /* Subroutine */ int prn2lb_(integer *n, doublereal *x, doublereal *f,
3270 doublereal *g, integer *iprint, integer *itfile, integer *iter,
3271 integer *nfgv, integer *nact, doublereal *sbgnrm, integer *nseg, char
3272 *word, integer *iword, integer *iback, doublereal *stp, doublereal *
3273 xstep, ftnlen word_len)
3274 {
3275 /* Builtin functions */
3276 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
3277
3278 /* Local variables */
3279 static integer imod;
3280
3281 /* ************ */
3282
3283 /* Subroutine prn2lb */
3284
3285 /* This subroutine prints out new information after a successful */
3286 /* line search. */
3287
3288
3289 /* * * * */
3290
3291 /* NEOS, November 1994. (Latest revision June 1996.) */
3292 /* Optimization Technology Center. */
3293 /* Argonne National Laboratory and Northwestern University. */
3294 /* Written by */
3295 /* Ciyou Zhu */
3296 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3297
3298
3299 /* ************ */
3300 /* 'word' records the status of subspace solutions. */
3301 /* Parameter adjustments */
3302 --g;
3303 --x;
3304
3305 /* Function Body */
3306 if (*iword == 0) {
3307 /* the subspace minimization converged. */
3308 s_copy(word, "con", (ftnlen)3, (ftnlen)3);
3309 } else if (*iword == 1) {
3310 /* the subspace minimization stopped at a bound. */
3311 s_copy(word, "bnd", (ftnlen)3, (ftnlen)3);
3312 } else if (*iword == 5) {
3313 /* the truncated Newton step has been used. */
3314 s_copy(word, "TNT", (ftnlen)3, (ftnlen)3);
3315 } else {
3316 s_copy(word, "---", (ftnlen)3, (ftnlen)3);
3317 }
3318 if (*iprint >= 99) {
3319 printf("LINE SEARCH %ld times; norm of step = %g\n", *iback, *xstep);
3320 /*
3321 2001 format
3322 + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
3323 */
3324 printf("At iterate %5ld f= %12.5g |proj g|= %12.5g\n",
3325 *iter, *f, *sbgnrm);
3326 if (*iprint > 100) {
3327 lbfgsb_printf_vec("X", x, *n);
3328 lbfgsb_printf_vec("G", g, *n);
3329 }
3330 } else if (*iprint > 0) {
3331 imod = *iter % *iprint;
3332 if (imod == 0) {
3333 /*
3334 2001 format
3335 + (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
3336 */
3337 printf("At iterate %5ld f= %12.5g |proj g|= %12.5g\n",
3338 *iter, *f, *sbgnrm);
3339 }
3340 }
3341 if (*iprint >= 1) {
3342 /*
3343 < 3001 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1p,2(1x,d10.3)) >
3344 */
3345 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3346 fprintf(0, " %4ld %4ld %5ld %5ld %3s %4ld %7.1g %7.1g<1p> %10.3g %10.3g\n",
3347 *iter, *nfgv, *nint, *nact, word, *iback, *stp, *xstep,
3348 *sbgnrm, *f);
3349 #endif
3350 }
3351 return 0;
3352 } /* prn2lb_ */
3353
3354 /* ======================= The end of prn2lb ============================= */
prn3lb_(integer * n,doublereal * x,doublereal * f,char * task,integer * iprint,integer * info,integer * itfile,integer * iter,integer * nfgv,integer * nintol,integer * nskip,integer * nact,doublereal * sbgnrm,doublereal * time,integer * nseg,char * word,integer * iback,doublereal * stp,doublereal * xstep,integer * k,doublereal * cachyt,doublereal * sbtime,doublereal * lnscht,ftnlen task_len,ftnlen word_len)3355 /* Subroutine */ int prn3lb_(integer *n, doublereal *x, doublereal *f, char *
3356 task, integer *iprint, integer *info, integer *itfile, integer *iter,
3357 integer *nfgv, integer *nintol, integer *nskip, integer *nact,
3358 doublereal *sbgnrm, doublereal *time, integer *nseg, char *word,
3359 integer *iback, doublereal *stp, doublereal *xstep, integer *k,
3360 doublereal *cachyt, doublereal *sbtime, doublereal *lnscht, ftnlen
3361 task_len, ftnlen word_len)
3362 {
3363 /* Builtin functions */
3364 integer s_cmp(char *, char *, ftnlen, ftnlen);
3365
3366 /* ************ */
3367
3368 /* Subroutine prn3lb */
3369
3370 /* This subroutine prints out information when either a built-in */
3371 /* convergence test is satisfied or when an error message is */
3372 /* generated. */
3373
3374
3375 /* * * * */
3376
3377 /* NEOS, November 1994. (Latest revision June 1996.) */
3378 /* Optimization Technology Center. */
3379 /* Argonne National Laboratory and Northwestern University. */
3380 /* Written by */
3381 /* Ciyou Zhu */
3382 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3383
3384
3385 /* ************ */
3386 /* Parameter adjustments */
3387 --x;
3388
3389 /* Function Body */
3390 if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
3391 goto L999;
3392 }
3393 if (*iprint >= 0) {
3394 /*
3395 3003 format (/,
3396 + ' * * *',/,/,
3397 + 'Tit = total number of iterations',/,
3398 + 'Tnf = total number of function evaluations',/,
3399 + 'Tnint = total number of segments explored during',
3400 + ' Cauchy searches',/,
3401 + 'Skip = number of BFGS updates skipped',/,
3402 + 'Nact = number of active bounds at final generalized',
3403 + ' Cauchy point',/,
3404 + 'Projg = norm of the final projected gradient',/,
3405 + 'F = final function value',/,/,
3406 + ' * * *')
3407 */
3408 printf(
3409 " * * *\n"
3410 "\n"
3411 "Tit = total number of iterations\n"
3412 "Tnf = total number of function evaluations\n"
3413 "Tnint = total number of segments explored during Cauchy searches\n"
3414 "Skip = number of BFGS updates skipped\n"
3415 "Nact = number of active bounds at final generalized Cauchy point\n"
3416 "Projg = norm of the final projected gradient\n"
3417 "F = final function value\n"
3418 "\n"
3419 " * * *\n");
3420 /*
3421 3004 format (/,3x,'N',3x,'Tit',2x,'Tnf',2x,'Tnint',2x,
3422 + 'Skip',2x,'Nact',5x,'Projg',8x,'F')
3423 */
3424 printf(" N Tit Tnf Tnint Skip Nact Projg F\n");
3425 /*
3426 3005 format (i5,2(1x,i4),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d10.3))
3427 */
3428 printf(" %4ld %4ld %4ld %6ld %4ld %5ld %10.3g %10.3g\n",
3429 *n, *iter, *nfgv, *nintol, *nskip, *nact, *sbgnrm, *f);
3430 if (*iprint >= 100) {
3431 lbfgsb_printf_vec("X", x, *n);
3432 }
3433 if (*iprint >= 1) {
3434 printf("F = %g\n", *f);
3435 }
3436 }
3437 L999:
3438 if (*iprint >= 0) {
3439 printf("%60s\n", task);
3440 if (*info != 0) {
3441 if (*info == -1) {
3442 /*
3443 9011 format (/,
3444 +' Matrix in 1st Cholesky factorization in formk is not Pos. Def.')
3445 */
3446 printf(" Matrix in 1st Cholesky factorization in formk is not Pos. Def.\n");
3447 }
3448 if (*info == -2) {
3449 /*
3450 9012 format (/,
3451 +' Matrix in 2st Cholesky factorization in formk is not Pos. Def.')
3452 */
3453 printf(" Matrix in 2st Cholesky factorization in formk is not Pos. Def.\n");
3454 }
3455 if (*info == -3) {
3456 /*
3457 9013 format (/,
3458 +' Matrix in the Cholesky factorization in formt is not Pos. Def.')
3459 */
3460 printf(" Matrix in the Cholesky factorization in formk is not Pos. Def.\n");
3461 }
3462 if (*info == -4) {
3463 /*
3464 9014 format (/,
3465 +' Derivative >= 0, backtracking line search impossible.',/,
3466 +' Previous x, f and g restored.',/,
3467 +' Possible causes: 1 error in function or gradient evaluation;',/,
3468 +' 2 rounding errors dominate computation.')
3469 */
3470 printf(" Derivative >= 0, backtracking line search impossible.\n"
3471 " Previous x, f and g restored.\n"
3472 " Possible causes: 1 error in function or gradient evaluation;\n"
3473 " 2 rounding errors dominate computation.\n");
3474 }
3475 if (*info == -5) {
3476 /*
3477 9015 format (/,
3478 +' Warning: more than 10 function and gradient',/,
3479 +' evaluations in the last line search. Termination',/,
3480 +' may possibly be caused by a bad search direction.')
3481 */
3482 printf(" Warning: more than 10 function and gradient\n"
3483 " evaluations in the last line search. Termination\n"
3484 " may possibly be caused by a bad search direction.");
3485 }
3486 if (*info == -6) {
3487 printf(" Input nbd(%ld) is invalid.\n", *k);
3488 }
3489 if (*info == -7) {
3490 printf(" l(%ld) > u(%ld). No feasible solution.\n", *k, *k);
3491 }
3492 if (*info == -8) {
3493 /*
3494 9018 format (/,' The triangular system is singular.')
3495 */
3496 printf(" The triangular system is singular.\n");
3497 }
3498 if (*info == -9) {
3499 /*
3500 9019 format (/,
3501 +' Line search cannot locate an adequate point after 20 function',/
3502 +,' and gradient evaluations. Previous x, f and g restored.',/,
3503 +' Possible causes: 1 error in function or gradient evaluation;',/,
3504 +' 2 rounding error dominate computation.')
3505
3506 */
3507 printf(" Line search cannot locate an adequate point after 20 function\n"
3508 " and gradient evaluations. Previous x, f and g restored.\n"
3509 " Possible causes: 1 error in function or gradient evaluation;\n"
3510 " 2 rounding error dominate computation.\n");
3511 }
3512 }
3513 if (*iprint >= 1) {
3514 printf(" Cauchy time %10.3g seconds.\n"
3515 " Subspace minimization time %10.3g seconds.\n"
3516 " Line search time %10.3g seconds.\n",
3517 *cachyt, *sbtime, *lnscht);
3518 }
3519 printf(" Total User time %10.3g seconds.\n", *time);
3520 if (*iprint >= 1) {
3521 if (*info == -4 || *info == -9) {
3522 /*
3523 3002 format(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6x,'-',10x,'-')
3524 */
3525 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3526 fprintf(0,
3527 " %4ld %4ld %5ld %5ld %3s %4ld %7.1g %7.1g - -\n",
3528 *iter, *nfgv, *nint, *nact, word, *iback, *stp, *xstep);
3529 #endif
3530 }
3531 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3532 fprintf(0, "%60s\n", task);
3533 #endif
3534 if (*info != 0) {
3535 if (*info == -1) {
3536 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3537 fprintf(0, " Matrix in 1st Cholesky factorization in formk is not Pos. Def.\n");
3538 #endif
3539 }
3540 if (*info == -2) {
3541 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3542 fprintf(0, " Matrix in 2st Cholesky factorization in formk is not Pos. Def.\n");
3543 #endif
3544 }
3545 if (*info == -3) {
3546 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3547 fprintf(0, " Matrix in the Cholesky factorization in formk is not Pos. Def.\n");
3548 #endif
3549 }
3550 if (*info == -4) {
3551 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3552 fprintf(0, " Derivative >= 0, backtracking line search impossible.\n"
3553 " Previous x, f and g restored.\n"
3554 " Possible causes: 1 error in function or gradient evaluation;\n"
3555 " 2 rounding errors dominate computation.\n");
3556 #endif
3557 }
3558 if (*info == -5) {
3559 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3560 fprintf(0, " Warning: more than 10 function and gradient\n"
3561 " evaluations in the last line search. Termination\n"
3562 " may possibly be caused by a bad search direction.");
3563 #endif
3564 }
3565 if (*info == -8) {
3566 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3567 fprintf(0, " The triangular system is singular.\n");
3568 #endif
3569 }
3570 if (*info == -9) {
3571 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3572 fprintf(0, " Line search cannot locate an adequate point after 20 function\n"
3573 " and gradient evaluations. Previous x, f and g restored.\n"
3574 " Possible causes: 1 error in function or gradient evaluation;\n"
3575 " 2 rounding error dominate computation.\n");
3576 #endif
3577 }
3578 }
3579 #ifdef LBFGSB_ENABLE_ITERATE_FILE
3580 fprintf(0, " Total User time %10.3g seconds.\n", *time);
3581 #endif
3582 }
3583 }
3584 return 0;
3585 } /* prn3lb_ */
3586
3587 /* ======================= The end of prn3lb ============================= */
projgr_(integer * n,doublereal * l,doublereal * u,integer * nbd,doublereal * x,doublereal * g,doublereal * sbgnrm)3588 /* Subroutine */ int projgr_(integer *n, doublereal *l, doublereal *u,
3589 integer *nbd, doublereal *x, doublereal *g, doublereal *sbgnrm)
3590 {
3591 /* System generated locals */
3592 integer i__1;
3593 doublereal d__1, d__2;
3594
3595 /* Local variables */
3596 static integer i__;
3597 static doublereal gi;
3598
3599 /* ************ */
3600
3601 /* Subroutine projgr */
3602
3603 /* This subroutine computes the infinity norm of the projected */
3604 /* gradient. */
3605
3606
3607 /* * * * */
3608
3609 /* NEOS, November 1994. (Latest revision June 1996.) */
3610 /* Optimization Technology Center. */
3611 /* Argonne National Laboratory and Northwestern University. */
3612 /* Written by */
3613 /* Ciyou Zhu */
3614 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3615
3616
3617 /* ************ */
3618 /* Parameter adjustments */
3619 --g;
3620 --x;
3621 --nbd;
3622 --u;
3623 --l;
3624
3625 /* Function Body */
3626 *sbgnrm = 0.;
3627 i__1 = *n;
3628 for (i__ = 1; i__ <= i__1; ++i__) {
3629 gi = g[i__];
3630 if (nbd[i__] != 0) {
3631 if (gi < 0.) {
3632 if (nbd[i__] >= 2) {
3633 /* Computing MAX */
3634 d__1 = x[i__] - u[i__];
3635 gi = max(d__1,gi);
3636 }
3637 } else {
3638 if (nbd[i__] <= 2) {
3639 /* Computing MIN */
3640 d__1 = x[i__] - l[i__];
3641 gi = min(d__1,gi);
3642 }
3643 }
3644 }
3645 /* Computing MAX */
3646 d__1 = *sbgnrm, d__2 = abs(gi);
3647 *sbgnrm = max(d__1,d__2);
3648 /* L15: */
3649 }
3650 return 0;
3651 } /* projgr_ */
3652
3653 /* ======================= The end of projgr ============================= */
subsm_(integer * n,integer * m,integer * nsub,integer * ind,doublereal * l,doublereal * u,integer * nbd,doublereal * x,doublereal * d__,doublereal * xp,doublereal * ws,doublereal * wy,doublereal * theta,doublereal * xx,doublereal * gg,integer * col,integer * head,integer * iword,doublereal * wv,doublereal * wn,integer * iprint,integer * info)3654 /* Subroutine */ int subsm_(integer *n, integer *m, integer *nsub, integer *
3655 ind, doublereal *l, doublereal *u, integer *nbd, doublereal *x,
3656 doublereal *d__, doublereal *xp, doublereal *ws, doublereal *wy,
3657 doublereal *theta, doublereal *xx, doublereal *gg, integer *col,
3658 integer *head, integer *iword, doublereal *wv, doublereal *wn,
3659 integer *iprint, integer *info)
3660 {
3661 /* System generated locals */
3662 integer ws_dim1, ws_offset, wy_dim1, wy_offset, wn_dim1, wn_offset, i__1,
3663 i__2;
3664 doublereal d__1, d__2;
3665
3666 /* Local variables */
3667 static integer i__, j, k, m2;
3668 static doublereal dk;
3669 static integer js, jy;
3670 static doublereal xk;
3671 static integer ibd, col2;
3672 static doublereal dd_p__, temp1, temp2, alpha;
3673 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
3674 integer *), dcopy_(integer *, doublereal *, integer *, doublereal
3675 *, integer *), dtrsl_(doublereal *, integer *, integer *,
3676 doublereal *, integer *, integer *);
3677 static integer pointr;
3678
3679 /* ********************************************************************** */
3680
3681 /* This routine contains the major changes in the updated version. */
3682 /* The changes are described in the accompanying paper */
3683
3684 /* Jose Luis Morales, Jorge Nocedal */
3685 /* "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale */
3686 /* Bound Constrained Optimization". Decemmber 27, 2010. */
3687
3688 /* J.L. Morales Departamento de Matematicas, */
3689 /* Instituto Tecnologico Autonomo de Mexico */
3690 /* Mexico D.F. */
3691
3692 /* J, Nocedal Department of Electrical Engineering and */
3693 /* Computer Science. */
3694 /* Northwestern University. Evanston, IL. USA */
3695
3696 /* January 17, 2011 */
3697
3698 /* ********************************************************************** */
3699
3700
3701 /* Subroutine subsm */
3702
3703 /* Given xcp, l, u, r, an index set that specifies */
3704 /* the active set at xcp, and an l-BFGS matrix B */
3705 /* (in terms of WY, WS, SY, WT, head, col, and theta), */
3706 /* this subroutine computes an approximate solution */
3707 /* of the subspace problem */
3708
3709 /* (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) */
3710
3711 /* subject to l<=x<=u */
3712 /* x_i=xcp_i for all i in A(xcp) */
3713
3714 /* along the subspace unconstrained Newton direction */
3715
3716 /* d = -(Z'BZ)^(-1) r. */
3717
3718 /* The formula for the Newton direction, given the L-BFGS matrix */
3719 /* and the Sherman-Morrison formula, is */
3720
3721 /* d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. */
3722
3723 /* where */
3724 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
3725 /* [L_a -R_z theta*S'AA'S ] */
3726
3727 /* Note that this procedure for computing d differs */
3728 /* from that described in [1]. One can show that the matrix K is */
3729 /* equal to the matrix M^[-1]N in that paper. */
3730
3731 /* n is an integer variable. */
3732 /* On entry n is the dimension of the problem. */
3733 /* On exit n is unchanged. */
3734
3735 /* m is an integer variable. */
3736 /* On entry m is the maximum number of variable metric corrections */
3737 /* used to define the limited memory matrix. */
3738 /* On exit m is unchanged. */
3739
3740 /* nsub is an integer variable. */
3741 /* On entry nsub is the number of free variables. */
3742 /* On exit nsub is unchanged. */
3743
3744 /* ind is an integer array of dimension nsub. */
3745 /* On entry ind specifies the coordinate indices of free variables. */
3746 /* On exit ind is unchanged. */
3747
3748 /* l is a double precision array of dimension n. */
3749 /* On entry l is the lower bound of x. */
3750 /* On exit l is unchanged. */
3751
3752 /* u is a double precision array of dimension n. */
3753 /* On entry u is the upper bound of x. */
3754 /* On exit u is unchanged. */
3755
3756 /* nbd is a integer array of dimension n. */
3757 /* On entry nbd represents the type of bounds imposed on the */
3758 /* variables, and must be specified as follows: */
3759 /* nbd(i)=0 if x(i) is unbounded, */
3760 /* 1 if x(i) has only a lower bound, */
3761 /* 2 if x(i) has both lower and upper bounds, and */
3762 /* 3 if x(i) has only an upper bound. */
3763 /* On exit nbd is unchanged. */
3764
3765 /* x is a double precision array of dimension n. */
3766 /* On entry x specifies the Cauchy point xcp. */
3767 /* On exit x(i) is the minimizer of Q over the subspace of */
3768 /* free variables. */
3769
3770 /* d is a double precision array of dimension n. */
3771 /* On entry d is the reduced gradient of Q at xcp. */
3772 /* On exit d is the Newton direction of Q. */
3773
3774 /* xp is a double precision array of dimension n. */
3775 /* used to safeguard the projected Newton direction */
3776
3777 /* xx is a double precision array of dimension n */
3778 /* On entry it holds the current iterate */
3779 /* On output it is unchanged */
3780 /* gg is a double precision array of dimension n */
3781 /* On entry it holds the gradient at the current iterate */
3782 /* On output it is unchanged */
3783
3784 /* ws and wy are double precision arrays; */
3785 /* theta is a double precision variable; */
3786 /* col is an integer variable; */
3787 /* head is an integer variable. */
3788 /* On entry they store the information defining the */
3789 /* limited memory BFGS matrix: */
3790 /* ws(n,m) stores S, a set of s-vectors; */
3791 /* wy(n,m) stores Y, a set of y-vectors; */
3792 /* theta is the scaling factor specifying B_0 = theta I; */
3793 /* col is the number of variable metric corrections stored; */
3794 /* head is the location of the 1st s- (or y-) vector in S (or Y). */
3795 /* On exit they are unchanged. */
3796
3797 /* iword is an integer variable. */
3798 /* On entry iword is unspecified. */
3799 /* On exit iword specifies the status of the subspace solution. */
3800 /* iword = 0 if the solution is in the box, */
3801 /* 1 if some bound is encountered. */
3802
3803 /* wv is a double precision working array of dimension 2m. */
3804
3805 /* wn is a double precision array of dimension 2m x 2m. */
3806 /* On entry the upper triangle of wn stores the LEL^T factorization */
3807 /* of the indefinite matrix */
3808
3809 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
3810 /* [L_a -R_z theta*S'AA'S ] */
3811 /* where E = [-I 0] */
3812 /* [ 0 I] */
3813 /* On exit wn is unchanged. */
3814
3815 /* iprint is an INTEGER variable that must be set by the user. */
3816 /* It controls the frequency and type of output generated: */
3817 /* iprint<0 no output is generated; */
3818 /* iprint=0 print only one line at the last iteration; */
3819 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
3820 /* iprint=99 print details of every iteration except n-vectors; */
3821 /* iprint=100 print also the changes of active set and final x; */
3822 /* iprint>100 print details of every iteration including x and g; */
3823 /* When iprint > 0, the file iterate.dat will be created to */
3824 /* summarize the iteration. */
3825
3826 /* info is an integer variable. */
3827 /* On entry info is unspecified. */
3828 /* On exit info = 0 for normal return, */
3829 /* = nonzero for abnormal return */
3830 /* when the matrix K is ill-conditioned. */
3831
3832 /* Subprograms called: */
3833
3834 /* Linpack dtrsl. */
3835
3836
3837 /* References: */
3838
3839 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
3840 /* memory algorithm for bound constrained optimization'', */
3841 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
3842
3843
3844
3845 /* * * * */
3846
3847 /* NEOS, November 1994. (Latest revision June 1996.) */
3848 /* Optimization Technology Center. */
3849 /* Argonne National Laboratory and Northwestern University. */
3850 /* Written by */
3851 /* Ciyou Zhu */
3852 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
3853
3854
3855 /* ************ */
3856
3857 /* Parameter adjustments */
3858 --gg;
3859 --xx;
3860 --xp;
3861 --d__;
3862 --x;
3863 --nbd;
3864 --u;
3865 --l;
3866 wn_dim1 = 2 * *m;
3867 wn_offset = 1 + wn_dim1;
3868 wn -= wn_offset;
3869 --wv;
3870 wy_dim1 = *n;
3871 wy_offset = 1 + wy_dim1;
3872 wy -= wy_offset;
3873 ws_dim1 = *n;
3874 ws_offset = 1 + ws_dim1;
3875 ws -= ws_offset;
3876 --ind;
3877
3878 /* Function Body */
3879 if (*nsub <= 0) {
3880 return 0;
3881 }
3882 if (*iprint >= 99) {
3883 printf("----------------SUBSM entered-----------------");
3884 }
3885 /* Compute wv = W'Zd. */
3886 pointr = *head;
3887 i__1 = *col;
3888 for (i__ = 1; i__ <= i__1; ++i__) {
3889 temp1 = 0.;
3890 temp2 = 0.;
3891 i__2 = *nsub;
3892 for (j = 1; j <= i__2; ++j) {
3893 k = ind[j];
3894 temp1 += wy[k + pointr * wy_dim1] * d__[j];
3895 temp2 += ws[k + pointr * ws_dim1] * d__[j];
3896 /* L10: */
3897 }
3898 wv[i__] = temp1;
3899 wv[*col + i__] = *theta * temp2;
3900 pointr = pointr % *m + 1;
3901 /* L20: */
3902 }
3903 /* Compute wv:=K^(-1)wv. */
3904 m2 = *m << 1;
3905 col2 = *col << 1;
3906 dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info);
3907 if (*info != 0) {
3908 return 0;
3909 }
3910 i__1 = *col;
3911 for (i__ = 1; i__ <= i__1; ++i__) {
3912 wv[i__] = -wv[i__];
3913 /* L25: */
3914 }
3915 dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info);
3916 if (*info != 0) {
3917 return 0;
3918 }
3919 /* Compute d = (1/theta)d + (1/theta**2)Z'W wv. */
3920 pointr = *head;
3921 i__1 = *col;
3922 for (jy = 1; jy <= i__1; ++jy) {
3923 js = *col + jy;
3924 i__2 = *nsub;
3925 for (i__ = 1; i__ <= i__2; ++i__) {
3926 k = ind[i__];
3927 d__[i__] = d__[i__] + wy[k + pointr * wy_dim1] * wv[jy] / *theta
3928 + ws[k + pointr * ws_dim1] * wv[js];
3929 /* L30: */
3930 }
3931 pointr = pointr % *m + 1;
3932 /* L40: */
3933 }
3934 d__1 = 1. / *theta;
3935 dscal_(nsub, &d__1, &d__[1], &c__1);
3936
3937 /* ----------------------------------------------------------------- */
3938 /* Let us try the projection, d is the Newton direction */
3939 *iword = 0;
3940 dcopy_(n, &x[1], &c__1, &xp[1], &c__1);
3941
3942 i__1 = *nsub;
3943 for (i__ = 1; i__ <= i__1; ++i__) {
3944 k = ind[i__];
3945 dk = d__[i__];
3946 xk = x[k];
3947 if (nbd[k] != 0) {
3948
3949 if (nbd[k] == 1) {
3950 /* lower bounds only */
3951 /* Computing MAX */
3952 d__1 = l[k], d__2 = xk + dk;
3953 x[k] = max(d__1,d__2);
3954 if (x[k] == l[k]) {
3955 *iword = 1;
3956 }
3957 } else {
3958
3959 if (nbd[k] == 2) {
3960 /* upper and lower bounds */
3961 /* Computing MAX */
3962 d__1 = l[k], d__2 = xk + dk;
3963 xk = max(d__1,d__2);
3964 /* Computing MIN */
3965 d__1 = u[k];
3966 x[k] = min(d__1,xk);
3967 if (x[k] == l[k] || x[k] == u[k]) {
3968 *iword = 1;
3969 }
3970 } else {
3971
3972 if (nbd[k] == 3) {
3973 /* upper bounds only */
3974 /* Computing MIN */
3975 d__1 = u[k], d__2 = xk + dk;
3976 x[k] = min(d__1,d__2);
3977 if (x[k] == u[k]) {
3978 *iword = 1;
3979 }
3980 }
3981 }
3982 }
3983
3984 } else {
3985 /* free variables */
3986 x[k] = xk + dk;
3987 }
3988 /* L50: */
3989 }
3990
3991 if (*iword == 0) {
3992 goto L911;
3993 }
3994
3995 /* check sign of the directional derivative */
3996
3997 dd_p__ = 0.;
3998 i__1 = *n;
3999 for (i__ = 1; i__ <= i__1; ++i__) {
4000 dd_p__ += (x[i__] - xx[i__]) * gg[i__];
4001 /* L55: */
4002 }
4003 if (dd_p__ > 0.) {
4004 dcopy_(n, &xp[1], &c__1, &x[1], &c__1);
4005 if (*iprint > 0) {
4006 printf(" Positive dir derivative in projection ");
4007 printf(" Using the backtracking step ");
4008 }
4009 } else {
4010 goto L911;
4011 }
4012
4013 /* ----------------------------------------------------------------- */
4014
4015 alpha = 1.;
4016 temp1 = alpha;
4017 ibd = 0;
4018 i__1 = *nsub;
4019 for (i__ = 1; i__ <= i__1; ++i__) {
4020 k = ind[i__];
4021 dk = d__[i__];
4022 if (nbd[k] != 0) {
4023 if (dk < 0. && nbd[k] <= 2) {
4024 temp2 = l[k] - x[k];
4025 if (temp2 >= 0.) {
4026 temp1 = 0.;
4027 } else if (dk * alpha < temp2) {
4028 temp1 = temp2 / dk;
4029 }
4030 } else if (dk > 0. && nbd[k] >= 2) {
4031 temp2 = u[k] - x[k];
4032 if (temp2 <= 0.) {
4033 temp1 = 0.;
4034 } else if (dk * alpha > temp2) {
4035 temp1 = temp2 / dk;
4036 }
4037 }
4038 if (temp1 < alpha) {
4039 alpha = temp1;
4040 ibd = i__;
4041 }
4042 }
4043 /* L60: */
4044 }
4045 if (alpha < 1.) {
4046 dk = d__[ibd];
4047 k = ind[ibd];
4048 if (dk > 0.) {
4049 x[k] = u[k];
4050 d__[ibd] = 0.;
4051 } else if (dk < 0.) {
4052 x[k] = l[k];
4053 d__[ibd] = 0.;
4054 }
4055 }
4056 i__1 = *nsub;
4057 for (i__ = 1; i__ <= i__1; ++i__) {
4058 k = ind[i__];
4059 x[k] += alpha * d__[i__];
4060 /* L70: */
4061 }
4062 /* ccccc */
4063 L911:
4064 if (*iprint >= 99) {
4065 if (alpha < 1.) {
4066 printf("ALPHA = %7.5g backtrack to the BOX\n", alpha);
4067 } else {
4068 printf("SM solution inside the box\n");
4069 }
4070 /*
4071 1003 format ('Subspace solution X = ',/,(4x,1p,6(1x,d11.4)))
4072 */
4073 if (*iprint > 100) {
4074 i__1 = *n;
4075 lbfgsb_printf_vec("Subspace solution X", x, i__1);
4076 }
4077 }
4078 if (alpha < 1.) {
4079 *iword = 1;
4080 } else {
4081 *iword = 0;
4082 }
4083 /*
4084 1004 format (/,'----------------exit SUBSM --------------------',/)
4085 */
4086 if (*iprint >= 99) {
4087 printf("----------------exit SUBSM --------------------");
4088 }
4089 return 0;
4090 } /* subsm_ */
4091
4092 /* ====================== The end of subsm =============================== */
dcsrch_(doublereal * f,doublereal * g,doublereal * stp,doublereal * ftol,doublereal * gtol,doublereal * xtol,doublereal * stpmin,doublereal * stpmax,char * task,integer * isave,doublereal * dsave,ftnlen task_len)4093 /* Subroutine */ int dcsrch_(doublereal *f, doublereal *g, doublereal *stp,
4094 doublereal *ftol, doublereal *gtol, doublereal *xtol, doublereal *
4095 stpmin, doublereal *stpmax, char *task, integer *isave, doublereal *
4096 dsave, ftnlen task_len)
4097 {
4098 /* System generated locals */
4099 doublereal d__1;
4100
4101 /* Builtin functions */
4102 integer s_cmp(char *, char *, ftnlen, ftnlen);
4103 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
4104
4105 /* Local variables */
4106 static doublereal fm, gm, fx, fy, gx, gy, fxm, fym, gxm, gym, stx, sty;
4107 static integer stage;
4108 static doublereal finit, ginit, width, ftest, gtest, stmin, stmax, width1;
4109 static logical brackt;
4110 extern /* Subroutine */ int dcstep_(doublereal *, doublereal *,
4111 doublereal *, doublereal *, doublereal *, doublereal *,
4112 doublereal *, doublereal *, doublereal *, logical *, doublereal *,
4113 doublereal *);
4114
4115 /* ********** */
4116
4117 /* Subroutine dcsrch */
4118
4119 /* This subroutine finds a step that satisfies a sufficient */
4120 /* decrease condition and a curvature condition. */
4121
4122 /* Each call of the subroutine updates an interval with */
4123 /* endpoints stx and sty. The interval is initially chosen */
4124 /* so that it contains a minimizer of the modified function */
4125
4126 /* psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). */
4127
4128 /* If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
4129 /* interval is chosen so that it contains a minimizer of f. */
4130
4131 /* The algorithm is designed to find a step that satisfies */
4132 /* the sufficient decrease condition */
4133
4134 /* f(stp) <= f(0) + ftol*stp*f'(0), */
4135
4136 /* and the curvature condition */
4137
4138 /* abs(f'(stp)) <= gtol*abs(f'(0)). */
4139
4140 /* If ftol is less than gtol and if, for example, the function */
4141 /* is bounded below, then there is always a step which satisfies */
4142 /* both conditions. */
4143
4144 /* If no step can be found that satisfies both conditions, then */
4145 /* the algorithm stops with a warning. In this case stp only */
4146 /* satisfies the sufficient decrease condition. */
4147
4148 /* A typical invocation of dcsrch has the following outline: */
4149
4150 /* task = 'START' */
4151 /* 10 continue */
4152 /* call dcsrch( ... ) */
4153 /* if (task .eq. 'FG') then */
4154 /* Evaluate the function and the gradient at stp */
4155 /* goto 10 */
4156 /* end if */
4157
4158 /* NOTE: The user must no alter work arrays between calls. */
4159
4160 /* The subroutine statement is */
4161
4162 /* subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, */
4163 /* task,isave,dsave) */
4164 /* where */
4165
4166 /* f is a double precision variable. */
4167 /* On initial entry f is the value of the function at 0. */
4168 /* On subsequent entries f is the value of the */
4169 /* function at stp. */
4170 /* On exit f is the value of the function at stp. */
4171
4172 /* g is a double precision variable. */
4173 /* On initial entry g is the derivative of the function at 0. */
4174 /* On subsequent entries g is the derivative of the */
4175 /* function at stp. */
4176 /* On exit g is the derivative of the function at stp. */
4177
4178 /* stp is a double precision variable. */
4179 /* On entry stp is the current estimate of a satisfactory */
4180 /* step. On initial entry, a positive initial estimate */
4181 /* must be provided. */
4182 /* On exit stp is the current estimate of a satisfactory step */
4183 /* if task = 'FG'. If task = 'CONV' then stp satisfies */
4184 /* the sufficient decrease and curvature condition. */
4185
4186 /* ftol is a double precision variable. */
4187 /* On entry ftol specifies a nonnegative tolerance for the */
4188 /* sufficient decrease condition. */
4189 /* On exit ftol is unchanged. */
4190
4191 /* gtol is a double precision variable. */
4192 /* On entry gtol specifies a nonnegative tolerance for the */
4193 /* curvature condition. */
4194 /* On exit gtol is unchanged. */
4195
4196 /* xtol is a double precision variable. */
4197 /* On entry xtol specifies a nonnegative relative tolerance */
4198 /* for an acceptable step. The subroutine exits with a */
4199 /* warning if the relative difference between sty and stx */
4200 /* is less than xtol. */
4201 /* On exit xtol is unchanged. */
4202
4203 /* stpmin is a double precision variable. */
4204 /* On entry stpmin is a nonnegative lower bound for the step. */
4205 /* On exit stpmin is unchanged. */
4206
4207 /* stpmax is a double precision variable. */
4208 /* On entry stpmax is a nonnegative upper bound for the step. */
4209 /* On exit stpmax is unchanged. */
4210
4211 /* task is a character variable of length at least 60. */
4212 /* On initial entry task must be set to 'START'. */
4213 /* On exit task indicates the required action: */
4214
4215 /* If task(1:2) = 'FG' then evaluate the function and */
4216 /* derivative at stp and call dcsrch again. */
4217
4218 /* If task(1:4) = 'CONV' then the search is successful. */
4219
4220 /* If task(1:4) = 'WARN' then the subroutine is not able */
4221 /* to satisfy the convergence conditions. The exit value of */
4222 /* stp contains the best point found during the search. */
4223
4224 /* If task(1:5) = 'ERROR' then there is an error in the */
4225 /* input arguments. */
4226
4227 /* On exit with convergence, a warning or an error, the */
4228 /* variable task contains additional information. */
4229
4230 /* isave is an integer work array of dimension 2. */
4231
4232 /* dsave is a double precision work array of dimension 13. */
4233
4234 /* Subprograms called */
4235
4236 /* MINPACK-2 ... dcstep */
4237
4238 /* MINPACK-1 Project. June 1983. */
4239 /* Argonne National Laboratory. */
4240 /* Jorge J. More' and David J. Thuente. */
4241
4242 /* MINPACK-2 Project. October 1993. */
4243 /* Argonne National Laboratory and University of Minnesota. */
4244 /* Brett M. Averick, Richard G. Carter, and Jorge J. More'. */
4245
4246 /* ********** */
4247 /* Initialization block. */
4248 /* Parameter adjustments */
4249 --dsave;
4250 --isave;
4251
4252 /* Function Body */
4253 if (s_cmp(task, "START", (ftnlen)5, (ftnlen)5) == 0) {
4254 /* Check the input arguments for errors. */
4255 if (*stp < *stpmin) {
4256 s_copy(task, "ERROR: STP .LT. STPMIN", task_len, (ftnlen)(22+1));
4257 }
4258 if (*stp > *stpmax) {
4259 s_copy(task, "ERROR: STP .GT. STPMAX", task_len, (ftnlen)(22+1));
4260 }
4261 if (*g >= 0.) {
4262 s_copy(task, "ERROR: INITIAL G .GE. ZERO", task_len, (ftnlen)(26+1));
4263 }
4264 if (*ftol < 0.) {
4265 s_copy(task, "ERROR: FTOL .LT. ZERO", task_len, (ftnlen)(21+1));
4266 }
4267 if (*gtol < 0.) {
4268 s_copy(task, "ERROR: GTOL .LT. ZERO", task_len, (ftnlen)(21+1));
4269 }
4270 if (*xtol < 0.) {
4271 s_copy(task, "ERROR: XTOL .LT. ZERO", task_len, (ftnlen)(21+1));
4272 }
4273 if (*stpmin < 0.) {
4274 s_copy(task, "ERROR: STPMIN .LT. ZERO", task_len, (ftnlen)(23+1));
4275 }
4276 if (*stpmax < *stpmin) {
4277 s_copy(task, "ERROR: STPMAX .LT. STPMIN", task_len, (ftnlen)(25+1));
4278 }
4279 /* Exit if there are errors on input. */
4280 if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
4281 return 0;
4282 }
4283 /* Initialize local variables. */
4284 brackt = FALSE_;
4285 stage = 1;
4286 finit = *f;
4287 ginit = *g;
4288 gtest = *ftol * ginit;
4289 width = *stpmax - *stpmin;
4290 width1 = width / .5;
4291 /* The variables stx, fx, gx contain the values of the step, */
4292 /* function, and derivative at the best step. */
4293 /* The variables sty, fy, gy contain the value of the step, */
4294 /* function, and derivative at sty. */
4295 /* The variables stp, f, g contain the values of the step, */
4296 /* function, and derivative at stp. */
4297 stx = 0.;
4298 fx = finit;
4299 gx = ginit;
4300 sty = 0.;
4301 fy = finit;
4302 gy = ginit;
4303 stmin = 0.;
4304 stmax = *stp + *stp * 4.;
4305 s_copy(task, "FG", task_len, (ftnlen)(2+1));
4306 goto L1000;
4307 } else {
4308 /* Restore local variables. */
4309 if (isave[1] == 1) {
4310 brackt = TRUE_;
4311 } else {
4312 brackt = FALSE_;
4313 }
4314 stage = isave[2];
4315 ginit = dsave[1];
4316 gtest = dsave[2];
4317 gx = dsave[3];
4318 gy = dsave[4];
4319 finit = dsave[5];
4320 fx = dsave[6];
4321 fy = dsave[7];
4322 stx = dsave[8];
4323 sty = dsave[9];
4324 stmin = dsave[10];
4325 stmax = dsave[11];
4326 width = dsave[12];
4327 width1 = dsave[13];
4328 }
4329 /* If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
4330 /* algorithm enters the second stage. */
4331 ftest = finit + *stp * gtest;
4332 if (stage == 1 && *f <= ftest && *g >= 0.) {
4333 stage = 2;
4334 }
4335 /* Test for warnings. */
4336 if (brackt && (*stp <= stmin || *stp >= stmax)) {
4337 s_copy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS", task_len, (
4338 ftnlen)(41+1));
4339 }
4340 if (brackt && stmax - stmin <= *xtol * stmax) {
4341 s_copy(task, "WARNING: XTOL TEST SATISFIED", task_len, (ftnlen)(28+1));
4342 }
4343 if (*stp == *stpmax && *f <= ftest && *g <= gtest) {
4344 s_copy(task, "WARNING: STP = STPMAX", task_len, (ftnlen)(21+1));
4345 }
4346 if (*stp == *stpmin && (*f > ftest || *g >= gtest)) {
4347 s_copy(task, "WARNING: STP = STPMIN", task_len, (ftnlen)(21+1));
4348 }
4349 /* Test for convergence. */
4350 if (*f <= ftest && abs(*g) <= *gtol * (-ginit)) {
4351 s_copy(task, "CONVERGENCE", task_len, (ftnlen)(11+1));
4352 }
4353 /* Test for termination. */
4354 if (s_cmp(task, "WARN", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(task, "CONV",
4355 (ftnlen)4, (ftnlen)4) == 0) {
4356 goto L1000;
4357 }
4358 /* A modified function is used to predict the step during the */
4359 /* first stage if a lower function value has been obtained but */
4360 /* the decrease is not sufficient. */
4361 if (stage == 1 && *f <= fx && *f > ftest) {
4362 /* Define the modified function and derivative values. */
4363 fm = *f - *stp * gtest;
4364 fxm = fx - stx * gtest;
4365 fym = fy - sty * gtest;
4366 gm = *g - gtest;
4367 gxm = gx - gtest;
4368 gym = gy - gtest;
4369 /* Call dcstep to update stx, sty, and to compute the new step. */
4370 dcstep_(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, &fm, &gm, &brackt, &
4371 stmin, &stmax);
4372 /* Reset the function and derivative values for f. */
4373 fx = fxm + stx * gtest;
4374 fy = fym + sty * gtest;
4375 gx = gxm + gtest;
4376 gy = gym + gtest;
4377 } else {
4378 /* Call dcstep to update stx, sty, and to compute the new step. */
4379 dcstep_(&stx, &fx, &gx, &sty, &fy, &gy, stp, f, g, &brackt, &stmin, &
4380 stmax);
4381 }
4382 /* Decide if a bisection step is needed. */
4383 if (brackt) {
4384 if ((d__1 = sty - stx, abs(d__1)) >= width1 * .66) {
4385 *stp = stx + (sty - stx) * .5;
4386 }
4387 width1 = width;
4388 width = (d__1 = sty - stx, abs(d__1));
4389 }
4390 /* Set the minimum and maximum steps allowed for stp. */
4391 if (brackt) {
4392 stmin = min(stx,sty);
4393 stmax = max(stx,sty);
4394 } else {
4395 stmin = *stp + (*stp - stx) * 1.1;
4396 stmax = *stp + (*stp - stx) * 4.;
4397 }
4398 /* Force the step to be within the bounds stpmax and stpmin. */
4399 *stp = max(*stp,*stpmin);
4400 *stp = min(*stp,*stpmax);
4401 /* If further progress is not possible, let stp be the best */
4402 /* point obtained during the search. */
4403 if ( ( brackt && ( ( *stp <= stmin ) || ( *stp >= stmax ) ) ) ||
4404 ( brackt && ( ( stmax - stmin ) <= ( *xtol * stmax ) ) ) ) {
4405 *stp = stx;
4406 }
4407 /* Obtain another function and derivative. */
4408 s_copy(task, "FG", task_len, (ftnlen)(2+1));
4409 L1000:
4410 /* Save local variables. */
4411 if (brackt) {
4412 isave[1] = 1;
4413 } else {
4414 isave[1] = 0;
4415 }
4416 isave[2] = stage;
4417 dsave[1] = ginit;
4418 dsave[2] = gtest;
4419 dsave[3] = gx;
4420 dsave[4] = gy;
4421 dsave[5] = finit;
4422 dsave[6] = fx;
4423 dsave[7] = fy;
4424 dsave[8] = stx;
4425 dsave[9] = sty;
4426 dsave[10] = stmin;
4427 dsave[11] = stmax;
4428 dsave[12] = width;
4429 dsave[13] = width1;
4430 return 0;
4431 } /* dcsrch_ */
4432
4433 /* ====================== The end of dcsrch ============================== */
dcstep_(doublereal * stx,doublereal * fx,doublereal * dx,doublereal * sty,doublereal * fy,doublereal * dy,doublereal * stp,doublereal * fp,doublereal * dp,logical * brackt,doublereal * stpmin,doublereal * stpmax)4434 /* Subroutine */ int dcstep_(doublereal *stx, doublereal *fx, doublereal *dx,
4435 doublereal *sty, doublereal *fy, doublereal *dy, doublereal *stp,
4436 doublereal *fp, doublereal *dp, logical *brackt, doublereal *stpmin,
4437 doublereal *stpmax)
4438 {
4439 /* System generated locals */
4440 doublereal d__1, d__2, d__3;
4441
4442 /* Builtin functions */
4443 double sqrt(doublereal);
4444
4445 /* Local variables */
4446 static doublereal p, q, r__, s, sgnd, stpc, stpf, stpq, gamma, theta;
4447
4448 /* ********** */
4449
4450 /* Subroutine dcstep */
4451
4452 /* This subroutine computes a safeguarded step for a search */
4453 /* procedure and updates an interval that contains a step that */
4454 /* satisfies a sufficient decrease and a curvature condition. */
4455
4456 /* The parameter stx contains the step with the least function */
4457 /* value. If brackt is set to .true. then a minimizer has */
4458 /* been bracketed in an interval with endpoints stx and sty. */
4459 /* The parameter stp contains the current step. */
4460 /* The subroutine assumes that if brackt is set to .true. then */
4461
4462 /* min(stx,sty) < stp < max(stx,sty), */
4463
4464 /* and that the derivative at stx is negative in the direction */
4465 /* of the step. */
4466
4467 /* The subroutine statement is */
4468
4469 /* subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, */
4470 /* stpmin,stpmax) */
4471
4472 /* where */
4473
4474 /* stx is a double precision variable. */
4475 /* On entry stx is the best step obtained so far and is an */
4476 /* endpoint of the interval that contains the minimizer. */
4477 /* On exit stx is the updated best step. */
4478
4479 /* fx is a double precision variable. */
4480 /* On entry fx is the function at stx. */
4481 /* On exit fx is the function at stx. */
4482
4483 /* dx is a double precision variable. */
4484 /* On entry dx is the derivative of the function at */
4485 /* stx. The derivative must be negative in the direction of */
4486 /* the step, that is, dx and stp - stx must have opposite */
4487 /* signs. */
4488 /* On exit dx is the derivative of the function at stx. */
4489
4490 /* sty is a double precision variable. */
4491 /* On entry sty is the second endpoint of the interval that */
4492 /* contains the minimizer. */
4493 /* On exit sty is the updated endpoint of the interval that */
4494 /* contains the minimizer. */
4495
4496 /* fy is a double precision variable. */
4497 /* On entry fy is the function at sty. */
4498 /* On exit fy is the function at sty. */
4499
4500 /* dy is a double precision variable. */
4501 /* On entry dy is the derivative of the function at sty. */
4502 /* On exit dy is the derivative of the function at the exit sty. */
4503
4504 /* stp is a double precision variable. */
4505 /* On entry stp is the current step. If brackt is set to .true. */
4506 /* then on input stp must be between stx and sty. */
4507 /* On exit stp is a new trial step. */
4508
4509 /* fp is a double precision variable. */
4510 /* On entry fp is the function at stp */
4511 /* On exit fp is unchanged. */
4512
4513 /* dp is a double precision variable. */
4514 /* On entry dp is the the derivative of the function at stp. */
4515 /* On exit dp is unchanged. */
4516
4517 /* brackt is an logical variable. */
4518 /* On entry brackt specifies if a minimizer has been bracketed. */
4519 /* Initially brackt must be set to .false. */
4520 /* On exit brackt specifies if a minimizer has been bracketed. */
4521 /* When a minimizer is bracketed brackt is set to .true. */
4522
4523 /* stpmin is a double precision variable. */
4524 /* On entry stpmin is a lower bound for the step. */
4525 /* On exit stpmin is unchanged. */
4526
4527 /* stpmax is a double precision variable. */
4528 /* On entry stpmax is an upper bound for the step. */
4529 /* On exit stpmax is unchanged. */
4530
4531 /* MINPACK-1 Project. June 1983 */
4532 /* Argonne National Laboratory. */
4533 /* Jorge J. More' and David J. Thuente. */
4534
4535 /* MINPACK-2 Project. October 1993. */
4536 /* Argonne National Laboratory and University of Minnesota. */
4537 /* Brett M. Averick and Jorge J. More'. */
4538
4539 /* ********** */
4540 sgnd = *dp * (*dx / abs(*dx));
4541 /* First case: A higher function value. The minimum is bracketed. */
4542 /* If the cubic step is closer to stx than the quadratic step, the */
4543 /* cubic step is taken, otherwise the average of the cubic and */
4544 /* quadratic steps is taken. */
4545 if (*fp > *fx) {
4546 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
4547 /* Computing MAX */
4548 d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
4549 *dp);
4550 s = max(d__1,d__2);
4551 /* Computing 2nd power */
4552 d__1 = theta / s;
4553 gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
4554 if (*stp < *stx) {
4555 gamma = -gamma;
4556 }
4557 p = gamma - *dx + theta;
4558 q = gamma - *dx + gamma + *dp;
4559 r__ = p / q;
4560 stpc = *stx + r__ * (*stp - *stx);
4561 stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2. * (*stp
4562 - *stx);
4563 if ((d__1 = stpc - *stx, abs(d__1)) < (d__2 = stpq - *stx, abs(d__2)))
4564 {
4565 stpf = stpc;
4566 } else {
4567 stpf = stpc + (stpq - stpc) / 2.;
4568 }
4569 *brackt = TRUE_;
4570 /* Second case: A lower function value and derivatives of opposite */
4571 /* sign. The minimum is bracketed. If the cubic step is farther from */
4572 /* stp than the secant step, the cubic step is taken, otherwise the */
4573 /* secant step is taken. */
4574 } else if (sgnd < 0.) {
4575 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
4576 /* Computing MAX */
4577 d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
4578 *dp);
4579 s = max(d__1,d__2);
4580 /* Computing 2nd power */
4581 d__1 = theta / s;
4582 gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
4583 if (*stp > *stx) {
4584 gamma = -gamma;
4585 }
4586 p = gamma - *dp + theta;
4587 q = gamma - *dp + gamma + *dx;
4588 r__ = p / q;
4589 stpc = *stp + r__ * (*stx - *stp);
4590 stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
4591 if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(d__2)))
4592 {
4593 stpf = stpc;
4594 } else {
4595 stpf = stpq;
4596 }
4597 *brackt = TRUE_;
4598 /* Third case: A lower function value, derivatives of the same sign, */
4599 /* and the magnitude of the derivative decreases. */
4600 } else if (abs(*dp) < abs(*dx)) {
4601 /* The cubic step is computed only if the cubic tends to infinity */
4602 /* in the direction of the step or if the minimum of the cubic */
4603 /* is beyond stp. Otherwise the cubic step is defined to be the */
4604 /* secant step. */
4605 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
4606 /* Computing MAX */
4607 d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
4608 *dp);
4609 s = max(d__1,d__2);
4610 /* The case gamma = 0 only arises if the cubic does not tend */
4611 /* to infinity in the direction of the step. */
4612 /* Computing MAX */
4613 /* Computing 2nd power */
4614 d__3 = theta / s;
4615 d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s);
4616 gamma = s * sqrt((max(d__1,d__2)));
4617 if (*stp > *stx) {
4618 gamma = -gamma;
4619 }
4620 p = gamma - *dp + theta;
4621 q = gamma + (*dx - *dp) + gamma;
4622 r__ = p / q;
4623 if (r__ < 0. && gamma != 0.) {
4624 stpc = *stp + r__ * (*stx - *stp);
4625 } else if (*stp > *stx) {
4626 stpc = *stpmax;
4627 } else {
4628 stpc = *stpmin;
4629 }
4630 stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
4631 if (*brackt) {
4632 /* A minimizer has been bracketed. If the cubic step is */
4633 /* closer to stp than the secant step, the cubic step is */
4634 /* taken, otherwise the secant step is taken. */
4635 if ((d__1 = stpc - *stp, abs(d__1)) < (d__2 = stpq - *stp, abs(
4636 d__2))) {
4637 stpf = stpc;
4638 } else {
4639 stpf = stpq;
4640 }
4641 if (*stp > *stx) {
4642 /* Computing MIN */
4643 d__1 = *stp + (*sty - *stp) * .66;
4644 stpf = min(d__1,stpf);
4645 } else {
4646 /* Computing MAX */
4647 d__1 = *stp + (*sty - *stp) * .66;
4648 stpf = max(d__1,stpf);
4649 }
4650 } else {
4651 /* A minimizer has not been bracketed. If the cubic step is */
4652 /* farther from stp than the secant step, the cubic step is */
4653 /* taken, otherwise the secant step is taken. */
4654 if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(
4655 d__2))) {
4656 stpf = stpc;
4657 } else {
4658 stpf = stpq;
4659 }
4660 stpf = min(*stpmax,stpf);
4661 stpf = max(*stpmin,stpf);
4662 }
4663 /* Fourth case: A lower function value, derivatives of the same sign, */
4664 /* and the magnitude of the derivative does not decrease. If the */
4665 /* minimum is not bracketed, the step is either stpmin or stpmax, */
4666 /* otherwise the cubic step is taken. */
4667 } else {
4668 if (*brackt) {
4669 theta = (*fp - *fy) * 3. / (*sty - *stp) + *dy + *dp;
4670 /* Computing MAX */
4671 d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 =
4672 abs(*dp);
4673 s = max(d__1,d__2);
4674 /* Computing 2nd power */
4675 d__1 = theta / s;
4676 gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s));
4677 if (*stp > *sty) {
4678 gamma = -gamma;
4679 }
4680 p = gamma - *dp + theta;
4681 q = gamma - *dp + gamma + *dy;
4682 r__ = p / q;
4683 stpc = *stp + r__ * (*sty - *stp);
4684 stpf = stpc;
4685 } else if (*stp > *stx) {
4686 stpf = *stpmax;
4687 } else {
4688 stpf = *stpmin;
4689 }
4690 }
4691 /* Update the interval which contains a minimizer. */
4692 if (*fp > *fx) {
4693 *sty = *stp;
4694 *fy = *fp;
4695 *dy = *dp;
4696 } else {
4697 if (sgnd < 0.) {
4698 *sty = *stx;
4699 *fy = *fx;
4700 *dy = *dx;
4701 }
4702 *stx = *stp;
4703 *fx = *fp;
4704 *dx = *dp;
4705 }
4706 /* Compute the new step. */
4707 *stp = stpf;
4708 return 0;
4709 } /* dcstep_ */
4710
4711 /* timer.f -- translated by f2c (version 20100827). */
timer_(doublereal * ttime)4712 /* Subroutine */ int timer_(doublereal *ttime)
4713 {
4714 //struct timeval t;
4715 //gettimeofday(&t, 0);
4716 //*ttime=1.0*double(t.tv_sec)+ 0.000001*double(t.tv_usec);
4717 *ttime = 0.0;
4718 return 0;
4719 } /* timer_ */
4720
4721
4722 /* linpack.f -- translated by f2c (version 20100827).
4723 You must link the resulting object file with libf2c:
4724 on Microsoft Windows system, link with libf2c.lib;
4725 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
4726 or, if you install libf2c.a in a standard place, with -lf2c -lm
4727 -- in that order, at the end of the command line, as in
4728 cc *.o -lf2c -lm
4729 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
4730
4731 http://www.netlib.org/f2c/libf2c.zip
4732 */
4733
4734
4735 /* Table of constant values */
lbfgsb_local_dpofa_(doublereal * a,integer * lda,integer * n,integer * info)4736 /* Subroutine */ int lbfgsb_local_dpofa_(doublereal *a, integer *lda, integer *n, integer *
4737 info)
4738 {
4739 /* System generated locals */
4740 integer a_dim1, a_offset, i__1, i__2, i__3;
4741
4742 /* Builtin functions */
4743 double sqrt(doublereal);
4744
4745 /* Local variables */
4746 static integer j, k;
4747 static doublereal s, t;
4748 static integer jm1;
4749 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
4750 integer *);
4751
4752
4753 /* dpofa factors a double precision symmetric positive definite */
4754 /* matrix. */
4755
4756 /* dpofa is usually called by dpoco, but it can be called */
4757 /* directly with a saving in time if rcond is not needed. */
4758 /* (time for dpoco) = (1 + 18/n)*(time for dpofa) . */
4759
4760 /* on entry */
4761
4762 /* a double precision(lda, n) */
4763 /* the symmetric matrix to be factored. only the */
4764 /* diagonal and upper triangle are used. */
4765
4766 /* lda integer */
4767 /* the leading dimension of the array a . */
4768
4769 /* n integer */
4770 /* the order of the matrix a . */
4771
4772 /* on return */
4773
4774 /* a an upper triangular matrix r so that a = trans(r)*r */
4775 /* where trans(r) is the transpose. */
4776 /* the strict lower triangle is unaltered. */
4777 /* if info .ne. 0 , the factorization is not complete. */
4778
4779 /* info integer */
4780 /* = 0 for normal return. */
4781 /* = k signals an error condition. the leading minor */
4782 /* of order k is not positive definite. */
4783
4784 /* linpack. this version dated 08/14/78 . */
4785 /* cleve moler, university of new mexico, argonne national lab. */
4786
4787 /* subroutines and functions */
4788
4789 /* blas ddot */
4790 /* fortran sqrt */
4791
4792 /* internal variables */
4793
4794 /* begin block with ...exits to 40 */
4795
4796
4797 /* Parameter adjustments */
4798 a_dim1 = *lda;
4799 a_offset = 1 + a_dim1;
4800 a -= a_offset;
4801
4802 /* Function Body */
4803 i__1 = *n;
4804 for (j = 1; j <= i__1; ++j) {
4805 *info = j;
4806 s = 0.;
4807 jm1 = j - 1;
4808 if (jm1 < 1) {
4809 goto L20;
4810 }
4811 i__2 = jm1;
4812 for (k = 1; k <= i__2; ++k) {
4813 i__3 = k - 1;
4814 t = a[k + j * a_dim1] - ddot_(&i__3, &a[k * a_dim1 + 1], &c__1, &
4815 a[j * a_dim1 + 1], &c__1);
4816 t /= a[k + k * a_dim1];
4817 a[k + j * a_dim1] = t;
4818 s += t * t;
4819 /* L10: */
4820 }
4821 L20:
4822 s = a[j + j * a_dim1] - s;
4823 /* ......exit */
4824 if (s <= 0.) {
4825 goto L40;
4826 }
4827 a[j + j * a_dim1] = sqrt(s);
4828 /* L30: */
4829 }
4830 *info = 0;
4831 L40:
4832 return 0;
4833 } /* lbfgsb_local_dpofa_ */
4834
4835 /* ====================== The end of dpofa =============================== */
dtrsl_(doublereal * t,integer * ldt,integer * n,doublereal * b,integer * job,integer * info)4836 /* Subroutine */ int dtrsl_(doublereal *t, integer *ldt, integer *n,
4837 doublereal *b, integer *job, integer *info)
4838 {
4839 /* System generated locals */
4840 integer t_dim1, t_offset, i__1, i__2;
4841
4842 /* Local variables */
4843 static integer j, jj, case__;
4844 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
4845 integer *);
4846 static doublereal temp;
4847 extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
4848 integer *, doublereal *, integer *);
4849
4850
4851
4852 /* dtrsl solves systems of the form */
4853
4854 /* t * x = b */
4855 /* or */
4856 /* trans(t) * x = b */
4857
4858 /* where t is a triangular matrix of order n. here trans(t) */
4859 /* denotes the transpose of the matrix t. */
4860
4861 /* on entry */
4862
4863 /* t double precision(ldt,n) */
4864 /* t contains the matrix of the system. the zero */
4865 /* elements of the matrix are not referenced, and */
4866 /* the corresponding elements of the array can be */
4867 /* used to store other information. */
4868
4869 /* ldt integer */
4870 /* ldt is the leading dimension of the array t. */
4871
4872 /* n integer */
4873 /* n is the order of the system. */
4874
4875 /* b double precision(n). */
4876 /* b contains the right hand side of the system. */
4877
4878 /* job integer */
4879 /* job specifies what kind of system is to be solved. */
4880 /* if job is */
4881
4882 /* 00 solve t*x=b, t lower triangular, */
4883 /* 01 solve t*x=b, t upper triangular, */
4884 /* 10 solve trans(t)*x=b, t lower triangular, */
4885 /* 11 solve trans(t)*x=b, t upper triangular. */
4886
4887 /* on return */
4888
4889 /* b b contains the solution, if info .eq. 0. */
4890 /* otherwise b is unaltered. */
4891
4892 /* info integer */
4893 /* info contains zero if the system is nonsingular. */
4894 /* otherwise info contains the index of */
4895 /* the first zero diagonal element of t. */
4896
4897 /* linpack. this version dated 08/14/78 . */
4898 /* g. w. stewart, university of maryland, argonne national lab. */
4899
4900 /* subroutines and functions */
4901
4902 /* blas daxpy,ddot */
4903 /* fortran mod */
4904
4905 /* internal variables */
4906
4907
4908 /* begin block permitting ...exits to 150 */
4909
4910 /* check for zero diagonal elements. */
4911
4912 /* Parameter adjustments */
4913 t_dim1 = *ldt;
4914 t_offset = 1 + t_dim1;
4915 t -= t_offset;
4916 --b;
4917
4918 /* Function Body */
4919 i__1 = *n;
4920 for (*info = 1; *info <= i__1; ++(*info)) {
4921 /* ......exit */
4922 if (t[*info + *info * t_dim1] == 0.) {
4923 goto L150;
4924 }
4925 /* L10: */
4926 }
4927 *info = 0;
4928
4929 /* determine the task and go to it. */
4930
4931 case__ = 1;
4932 if (*job % 10 != 0) {
4933 case__ = 2;
4934 }
4935 if (*job % 100 / 10 != 0) {
4936 case__ += 2;
4937 }
4938 switch (case__) {
4939 case 1: goto L20;
4940 case 2: goto L50;
4941 case 3: goto L80;
4942 case 4: goto L110;
4943 }
4944
4945 /* solve t*x=b for t lower triangular */
4946
4947 L20:
4948 b[1] /= t[t_dim1 + 1];
4949 if (*n < 2) {
4950 goto L40;
4951 }
4952 i__1 = *n;
4953 for (j = 2; j <= i__1; ++j) {
4954 temp = -b[j - 1];
4955 i__2 = *n - j + 1;
4956 daxpy_(&i__2, &temp, &t[j + (j - 1) * t_dim1], &c__1, &b[j], &c__1);
4957 b[j] /= t[j + j * t_dim1];
4958 /* L30: */
4959 }
4960 L40:
4961 goto L140;
4962
4963 /* solve t*x=b for t upper triangular. */
4964
4965 L50:
4966 b[*n] /= t[*n + *n * t_dim1];
4967 if (*n < 2) {
4968 goto L70;
4969 }
4970 i__1 = *n;
4971 for (jj = 2; jj <= i__1; ++jj) {
4972 j = *n - jj + 1;
4973 temp = -b[j + 1];
4974 daxpy_(&j, &temp, &t[(j + 1) * t_dim1 + 1], &c__1, &b[1], &c__1);
4975 b[j] /= t[j + j * t_dim1];
4976 /* L60: */
4977 }
4978 L70:
4979 goto L140;
4980
4981 /* solve trans(t)*x=b for t lower triangular. */
4982
4983 L80:
4984 b[*n] /= t[*n + *n * t_dim1];
4985 if (*n < 2) {
4986 goto L100;
4987 }
4988 i__1 = *n;
4989 for (jj = 2; jj <= i__1; ++jj) {
4990 j = *n - jj + 1;
4991 i__2 = jj - 1;
4992 b[j] -= ddot_(&i__2, &t[j + 1 + j * t_dim1], &c__1, &b[j + 1], &c__1);
4993 b[j] /= t[j + j * t_dim1];
4994 /* L90: */
4995 }
4996 L100:
4997 goto L140;
4998
4999 /* solve trans(t)*x=b for t upper triangular. */
5000
5001 L110:
5002 b[1] /= t[t_dim1 + 1];
5003 if (*n < 2) {
5004 goto L130;
5005 }
5006 i__1 = *n;
5007 for (j = 2; j <= i__1; ++j) {
5008 i__2 = j - 1;
5009 b[j] -= ddot_(&i__2, &t[j * t_dim1 + 1], &c__1, &b[1], &c__1);
5010 b[j] /= t[j + j * t_dim1];
5011 /* L120: */
5012 }
5013 L130:
5014 L140:
5015 L150:
5016 return 0;
5017 } /* dtrsl_ */
5018
5019 #ifdef __cplusplus
5020 }
5021 #endif
5022