ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
lbfgs.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: lbfgs.cpp 1768 2014-03-07 01:27:19Z johnoel $
00003  *
00004  * Author: Jorge Nodecal
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  *
00007  * License: Public domain
00008  * Translated from fortran and extensively modified for compliance with
00009  * ADMB.
00010  */
00015 /* lbfgs.f -- translated by f2c (version 19950110).
00016    You must link the resulting object file with the libraries:
00017         -lf2c -lm   (in that order)
00018 */
00019 
00020 #ifdef __cplusplus
00021 extern "C" {
00022 #endif
00023 
00024 /* f2c.h  --  Standard Fortran to C header file */
00025 
00030 #ifndef F2C_INCLUDE
00031 #define F2C_INCLUDE
00032 
00033 typedef long int integer;
00034 typedef char *address;
00035 typedef short int shortint;
00036 typedef float real;
00037 typedef double doublereal;
00038 typedef struct { real r, i; } complex;
00039 typedef struct { doublereal r, i; } doublecomplex;
00040 typedef long int logical;
00041 typedef short int shortlogical;
00042 typedef char logical1;
00043 typedef char integer1;
00044 /* typedef long long longint; */ /* system-dependent */
00045 
00046 #define TRUE_ (1)
00047 #define FALSE_ (0)
00048 
00049 /* Extern is for use with -E */
00050 #ifndef Extern
00051 #define Extern extern
00052 #endif
00053 
00054 /* I/O stuff */
00055 
00056 #ifdef f2c_i2
00057 /* for -i2 */
00058 typedef short flag;
00059 typedef short ftnlen;
00060 typedef short ftnint;
00061 #else
00062 typedef long int flag;
00063 typedef long int ftnlen;
00064 typedef long int ftnint;
00065 #endif
00066 
00067 /*external read, write*/
00068 typedef struct
00069 {
00070   flag cierr;
00071   ftnint ciunit;
00072   flag ciend;
00073   char *cifmt;
00074   ftnint cirec;
00075 } cilist;
00076 
00077 /*internal read, write*/
00078 typedef struct
00079 {
00080   flag icierr;
00081   char *iciunit;
00082   flag iciend;
00083   char *icifmt;
00084   ftnint icirlen;
00085   ftnint icirnum;
00086 } icilist;
00087 
00088 /*open*/
00089 typedef struct
00090 {
00091   flag oerr;
00092   ftnint ounit;
00093   char *ofnm;
00094   ftnlen ofnmlen;
00095   char *osta;
00096   char *oacc;
00097   char *ofm;
00098   ftnint orl;
00099   char *oblnk;
00100 } olist;
00101 
00102 /*close*/
00103 typedef struct
00104 {
00105   flag cerr;
00106   ftnint cunit;
00107   char *csta;
00108 } cllist;
00109 
00110 /*rewind, backspace, endfile*/
00111 typedef struct
00112 {
00113   flag aerr;
00114   ftnint aunit;
00115 } alist;
00116 
00117 /* inquire */
00118 typedef struct
00119 {
00120   flag inerr;
00121   ftnint inunit;
00122   char *infile;
00123   ftnlen infilen;
00124   ftnint *inex;/*parameters in standard's order*/
00125   ftnint *inopen;
00126   ftnint *innum;
00127   ftnint *innamed;
00128   char *inname;
00129   ftnlen innamlen;
00130   char *inacc;
00131   ftnlen inacclen;
00132   char *inseq;
00133   ftnlen inseqlen;
00134   char *indir;
00135   ftnlen indirlen;
00136   char *infmt;
00137   ftnlen infmtlen;
00138   char *inform;
00139   ftnint informlen;
00140   char *inunf;
00141   ftnlen inunflen;
00142   ftnint *inrecl;
00143   ftnint *innrec;
00144   char *inblank;
00145   ftnlen inblanklen;
00146 } inlist;
00147 
00148 #define VOID void
00149 
00150 union Multitype {/* for multiple entry points */
00151   integer1 g;
00152   shortint h;
00153   integer i;
00154   /* longint j; */
00155   real r;
00156   doublereal d;
00157   complex c;
00158   doublecomplex z;
00159 };
00160 
00161 typedef union Multitype Multitype;
00162 
00163 /*typedef long int Long;*//* No longer used; formerly in Namelist */
00164 
00165 struct Vardesc {/* for Namelist */
00166   char *name;
00167   char *addr;
00168   ftnlen *dims;
00169   int  type;
00170 };
00171 typedef struct Vardesc Vardesc;
00172 
00173 struct Namelist {
00174   char *name;
00175   Vardesc **vars;
00176   int nvars;
00177 };
00178 typedef struct Namelist Namelist;
00179 
00180 #define abs(x) ((x) >= 0 ? (x) : -(x))
00181 #define dabs(x) (doublereal)abs(x)
00182 #define min(a,b) ((a) <= (b) ? (a) : (b))
00183 #define max(a,b) ((a) >= (b) ? (a) : (b))
00184 #define dmin(a,b) (doublereal)min(a,b)
00185 #define dmax(a,b) (doublereal)max(a,b)
00186 
00187 /* procedure parameter types for -A and -C++ */
00188 
00189 #define F2C_proc_par_types 1
00190 #ifdef __cplusplus
00191 typedef int /* Unknown procedure type */ (*U_fp)(...);
00192 typedef shortint (*J_fp)(...);
00193 typedef integer (*I_fp)(...);
00194 typedef real (*R_fp)(...);
00195 typedef doublereal (*D_fp)(...), (*E_fp)(...);
00196 typedef /* Complex */ VOID (*C_fp)(...);
00197 typedef /* Double Complex */ VOID (*Z_fp)(...);
00198 typedef logical (*L_fp)(...);
00199 typedef shortlogical (*K_fp)(...);
00200 typedef /* Character */ VOID (*H_fp)(...);
00201 typedef /* Subroutine */ int (*S_fp)(...);
00202 #else
00203 typedef int /* Unknown procedure type */ (*U_fp)();
00204 typedef shortint (*J_fp)();
00205 typedef integer (*I_fp)();
00206 typedef real (*R_fp)();
00207 typedef doublereal (*D_fp)(), (*E_fp)();
00208 typedef /* Complex */ VOID (*C_fp)();
00209 typedef /* Double Complex */ VOID (*Z_fp)();
00210 typedef logical (*L_fp)();
00211 typedef shortlogical (*K_fp)();
00212 typedef /* Character */ VOID (*H_fp)();
00213 typedef /* Subroutine */ int (*S_fp)();
00214 #endif
00215 /* E_fp is for real functions when -R is not specified */
00216 typedef VOID C_f;/* complex function */
00217 typedef VOID H_f;/* character function */
00218 typedef VOID Z_f;/* double complex function */
00219 typedef doublereal E_f;/* real function with -R not specified */
00220 
00221 /* undef any lower-case symbols that your C compiler predefines, e.g.: */
00222 
00223 #ifndef Skip_f2c_Undefs
00224 #undef cray
00225 #undef gcos
00226 #undef mc68010
00227 #undef mc68020
00228 #undef mips
00229 #undef pdp11
00230 #undef sgi
00231 #undef sparc
00232 #undef sun
00233 #undef sun2
00234 #undef sun3
00235 #undef sun4
00236 #undef u370
00237 #undef u3b
00238 #undef u3b2
00239 #undef u3b5
00240 #undef unix
00241 #undef vax
00242 #endif
00243 #endif
00244 /* Common Block Declarations */
00245 
00246 struct lb3_1_ {
00247     integer mp, lp;
00248     doublereal gtol, stpmin, stpmax;
00249 };
00250 
00251 #define lb3_1 (*(struct lb3_1_ *) &lb3_)
00252 
00253 /* Initialized data */
00254 
00255 struct {
00256     integer e_1[2];
00257     doublereal e_2[3];
00258     } lb3_ = { {6, 6}, {.9, 1e-20, 1e20}};
00259 
00260 
00261 /* Table of constant values */
00262 
00263 static integer c__1 = 1;
00264 
00265 /*    ----------------------------------------------------------------------*/
00266 /*     This file contains the LBFGS algorithm and supporting routines */
00267 
00268 /*     **************** */
00269 /*     LBFGS SUBROUTINE */
00270 /*     **************** */
00271 
00272 /* Subroutine */ int lbfgs_(integer *n, integer *m, doublereal *x, doublereal
00273         *f, doublereal *g, logical *diagco, doublereal *diag, integer *iprint,
00274          doublereal *eps, doublereal *xtol, doublereal *w, integer *iflag,
00275    integer* iter, integer * info)
00276 {
00277     /* Initialized data */
00278 
00279     static doublereal one = 1.;
00280     static doublereal zero = 0.;
00281 
00282     /* Format strings */
00283 /*
00284     static char fmt_245[] = "(/\002  GTOL IS LESS THAN OR EQUAL TO 1.D-04\
00285 \002,/\002 IT HAS BEEN RESET TO 9.D-01\002)";
00286     static char fmt_200[] = "(/\002 IFLAG= -1 \002,/\002 LINE SEARCH FAILED.\
00287  SEE DOCUMENTATION OF ROUTINE XXXXXX\002,/\002 ERROR RETURN OF LINE SEARCH: \
00288 INFO= \002,i2,/\002 POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT\002,/\
00289 \002 OR INCORRECT TOLERANCES\002)";
00290     static char fmt_235[] = "(/\002 IFLAG= -2\002,/\002 THE\002,i5,\002-TH D\
00291 IAGONAL ELEMENT OF THE\002,/,\002 INVERSE HESSIAN APPROXIMATION IS NOT POSIT\
00292 IVE\002)";
00293     static char fmt_240[] = "(/\002 IFLAG= -3\002,/\002 IMPROPER INPUT PARAM\
00294 ETERS (N OR M\002,\002 ARE NOT POSITIVE)\002)";
00295 */
00296 
00297     /* System generated locals */
00298     integer i__1;
00299     doublereal d__1;
00300 
00301     /* Builtin functions */
00302     //integer //s_wsfe(cilist *), e_wsfe();
00303     double sqrt(doublereal);
00304     //integer // do_fio(integer *, char *, ftnlen);
00305 
00306     /* Local variables */
00307     static doublereal beta;
00308     static integer inmc;
00309     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
00310             integer *);
00311     static integer iscn, nfev, iycn;
00312     static doublereal ftol;
00313     static integer nfun, ispt, iypt, i, bound;
00314     static doublereal gnorm;
00315     extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *,
00316             integer *, doublereal *, integer *);
00317     static integer point;
00318     static doublereal xnorm;
00319     static integer cp;
00320     static doublereal sq, yr, ys;
00321     extern /* Subroutine */ int mcsrch_(integer *, doublereal *, doublereal *,
00322              doublereal *, doublereal *, doublereal *, doublereal *,
00323             doublereal *, integer *, integer *, integer *, doublereal *);
00324     static logical finish;
00325     static doublereal yy;
00326     static integer maxfev;
00327 /*
00328 //Subroutine 
00329     extern int lb1_(integer *, integer *, integer *,
00330             doublereal *, integer *, integer *, doublereal *, doublereal *,
00331             doublereal *, doublereal *, logical *);
00332 */
00333     static integer npt;
00334     static doublereal stp, stp1;
00335 
00336     /* Fortran I/O blocks */
00337 /*
00338     static cilist io___4 = { 0, 0, 0, fmt_245, 0 };
00339     static cilist io___30 = { 0, 0, 0, fmt_200, 0 };
00340     static cilist io___31 = { 0, 0, 0, fmt_235, 0 };
00341     static cilist io___32 = { 0, 0, 0, fmt_240, 0 };
00342 */
00343 
00344 
00345 /*        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION */
00346 /*                          JORGE NOCEDAL */
00347 /*                        *** July 1990 *** */
00348 
00349 
00350 /*     This subroutine solves the unconstrained minimization problem */
00351 
00352 /*                      min F(x),    x= (x1,x2,...,xN), */
00353 
00354 /*      using the limited memory BFGS method. The routine is especially */
00355 /*      effective on problems involving a large number of variables. In */
00356 /*      a typical iteration of this method an approximation Hk to the */
00357 /*      inverse of the Hessian is obtained by applying M BFGS updates to
00358 */
00359 /*     a diagonal matrix Hk0, using information from the previous M steps.
00360 */
00361 /*      The user specifies the number M, which determines the amount of */
00362 /*      storage required by the routine. The user may also provide the */
00363 /*      diagonal matrices Hk0 if not satisfied with the default choice. */
00364 /*      The algorithm is described in "On the limited memory BFGS method
00365 */
00366 /*      for large scale optimization", by D. Liu and J. Nocedal, */
00367 /*      Mathematical Programming B 45 (1989) 503-528. */
00368 
00369 /*      The user is required to calculate the function value F and its */
00370 /*      gradient G. In order to allow the user complete control over */
00371 /*      these computations, reverse  communication is used. The routine */
00372 /*      must be called repeatedly under the control of the parameter */
00373 /*      IFLAG. */
00374 
00375 /*      The steplength is determined at each iteration by means of the */
00376 /*      line search routine MCVSRCH, which is a slight modification of */
00377 /*      the routine CSRCH written by More' and Thuente. */
00378 
00379 /*      The calling statement is */
00380 
00381 /*          CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) */
00382 
00383 /*      where */
00384 
00385 /*     N       is an INTEGER variable that must be set by the user to the
00386 */
00387 /*             number of variables. It is not altered by the routine. */
00388 /*             Restriction: N>0. */
00389 
00390 /*     M       is an INTEGER variable that must be set by the user to */
00391 /*             the number of corrections used in the BFGS update. It */
00392 /*             is not altered by the routine. Values of M less than 3 are
00393 */
00394 /*             not recommended; large values of M will result in excessive
00395  */
00396 /*             computing time. 3<= M <=7 is recommended. Restriction: M>0.
00397  */
00398 
00399 /*     X       is a DOUBLE PRECISION array of length N. On initial entry
00400 */
00401 /*             it must be set by the user to the values of the initial */
00402 /*             estimate of the solution vector. On exit with IFLAG=0, it
00403 */
00404 /*             contains the values of the variables at the best point */
00405 /*             found (usually a solution). */
00406 
00407 /*     F       is a DOUBLE PRECISION variable. Before initial entry and on
00408  */
00409 /*             a re-entry with IFLAG=1, it must be set by the user to */
00410 /*             contain the value of the function F at the point X. */
00411 
00412 /*     G       is a DOUBLE PRECISION array of length N. Before initial */
00413 /*             entry and on a re-entry with IFLAG=1, it must be set by */
00414 /*             the user to contain the components of the gradient G at */
00415 /*             the point X. */
00416 
00417 /*     DIAGCO  is a LOGICAL variable that must be set to .TRUE. if the */
00418 /*             user  wishes to provide the diagonal matrix Hk0 at each */
00419 /*             iteration. Otherwise it should be set to .FALSE., in which
00420 */
00421 /*             case  LBFGS will use a default value described below. If */
00422 /*             DIAGCO is set to .TRUE. the routine will return at each */
00423 /*             iteration of the algorithm with IFLAG=2, and the diagonal
00424 */
00425 /*              matrix Hk0  must be provided in the array DIAG. */
00426 
00427 
00428 /*     DIAG    is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE.,
00429 */
00430 /*             then on initial entry or on re-entry with IFLAG=2, DIAG */
00431 /*             it must be set by the user to contain the values of the */
00432 /*             diagonal matrix Hk0.  Restriction: all elements of DIAG */
00433 /*             must be positive. */
00434 
00435 /*     IPRINT  is an INTEGER array of length two which must be set by the
00436 */
00437 /*             user. */
00438 
00439 /*             IPRINT(1) specifies the frequency of the output: */
00440 /*                IPRINT(1) < 0 : no output is generated, */
00441 /*                IPRINT(1) = 0 : output only at first and last iteration,
00442  */
00443 /*                IPRINT(1) > 0 : output every IPRINT(1) iterations. */
00444 
00445 /*             IPRINT(2) specifies the type of output generated: */
00446 /*                IPRINT(2) = 0 : iteration count, number of function */
00447 /*                                evaluations, function value, norm of the
00448  */
00449 /*                                gradient, and steplength, */
00450 /*                IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of */
00451 /*                                variables and  gradient vector at the */
00452 /*                                initial point, */
00453 /*                IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of */
00454 /*                                variables, */
00455 /*               IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector
00456 .*/
00457 
00458 
00459 /*     EPS     is a positive DOUBLE PRECISION variable that must be set by
00460  */
00461 /*            the user, and determines the accuracy with which the solutio
00462 n*/
00463 /*             is to be found. The subroutine terminates when */
00464 
00465 /*                         ||G|| < EPS max(1,||X||), */
00466 
00467 /*             where ||.|| denotes the Euclidean norm. */
00468 
00469 /*    XTOL    is a  positive DOUBLE PRECISION variable that must be set by
00470 */
00471 /*             the user to an estimate of the machine precision (e.g. */
00472 /*            10**(-16) on a SUN station 3/60). The line search routine wi
00473 ll*/
00474 /*            terminate if the relative width of the interval of uncertain
00475 ty*/
00476 /*             is less than XTOL. */
00477 
00478 /*     W       is a DOUBLE PRECISION array of length N(2M+1)+2M used as */
00479 /*             workspace for LBFGS. This array must not be altered by the
00480 */
00481 /*             user. */
00482 
00483 /*    IFLAG   is an INTEGER variable that must be set to 0 on initial entr
00484 y*/
00485 /*            to the subroutine. A return with IFLAG<0 indicates an error,
00486 */
00487 /*            and IFLAG=0 indicates that the routine has terminated withou
00488 t*/
00489 /*             detecting errors. On a return with IFLAG=1, the user must
00490 */
00491 /*             evaluate the function F and gradient G. On a return with */
00492 /*             IFLAG=2, the user must provide the diagonal matrix Hk0. */
00493 
00494 /*             The following negative values of IFLAG, detecting an error,
00495  */
00496 /*             are possible: */
00497 
00498 /*              IFLAG=-1  The line search routine MCSRCH failed. The */
00499 /*                       parameter INFO provides more detailed information
00500 */
00501 /*                        (see also the documentation of MCSRCH): */
00502 
00503 /*                       INFO = 0  IMPROPER INPUT PARAMETERS. */
00504 
00505 /*                       INFO = 2  RELATIVE WIDTH OF THE INTERVAL OF */
00506 /*                                 UNCERTAINTY IS AT MOST XTOL. */
00507 
00508 /*                       INFO = 3  MORE THAN 20 FUNCTION EVALUATIONS WERE
00509 */
00510 /*                                 REQUIRED AT THE PRESENT ITERATION. */
00511 
00512 /*                       INFO = 4  THE STEP IS TOO SMALL. */
00513 
00514 /*                       INFO = 5  THE STEP IS TOO LARGE. */
00515 
00516 /*                      INFO = 6  ROUNDING ERRORS PREVENT FURTHER PROGRESS
00517 .*/
00518 /*                                 THERE MAY NOT BE A STEP WHICH SATISFIES
00519  */
00520 /*                                 THE SUFFICIENT DECREASE AND CURVATURE
00521 */
00522 /*                                CONDITIONS. TOLERANCES MAY BE TOO SMALL.
00523 */
00524 
00525 
00526 /*             IFLAG=-2  The i-th diagonal element of the diagonal inverse
00527 */
00528 /*                        Hessian approximation, given in DIAG, is not */
00529 /*                        positive. */
00530 
00531 /*              IFLAG=-3  Improper input parameters for LBFGS (N or M are
00532 */
00533 /*                        not positive). */
00534 
00535 
00536 
00537 /*    ON THE DRIVER: */
00538 
00539 /*    The program that calls LBFGS must contain the declaration: */
00540 
00541 /*                       EXTERNAL LB2 */
00542 
00543 /*    LB2 is a BLOCK DATA that defines the default values of several */
00544 /*    parameters described in the COMMON section. */
00545 
00546 
00547 
00548 /*    COMMON: */
00549 
00550 /*     The subroutine contains one common area, which the user may wish to
00551  */
00552 /*    reference: */
00553 
00554 
00555 /*    MP  is an INTEGER variable with default value 6. It is used as the
00556 */
00557 /*        unit number for the printing of the monitoring information */
00558 /*        controlled by IPRINT. */
00559 
00560 /*    LP  is an INTEGER variable with default value 6. It is used as the
00561 */
00562 /*        unit number for the printing of error messages. This printing */
00563 /*        may be suppressed by setting LP to a non-positive value. */
00564 
00565 /*    GTOL is a DOUBLE PRECISION variable with default value 0.9, which */
00566 /*        controls the accuracy of the line search routine MCSRCH. If the
00567 */
00568 /*        function and gradient evaluations are inexpensive with respect
00569 */
00570 /*        to the cost of the iteration (which is sometimes the case when
00571 */
00572 /*        solving very large problems) it may be advantageous to set GTOL
00573 */
00574 /*        to a small value. A typical small value is 0.1.  Restriction: */
00575 /*        GTOL should be greater than 1.D-04. */
00576 
00577 /*    STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
00578 */
00579 /*        specify lower and uper bounds for the step in the line search.
00580 */
00581 /*        Their default values are 1.D-20 and 1.D+20, respectively. These
00582 */
00583 /*        values need not be modified unless the exponents are too large
00584 */
00585 /*        for the machine being used, or unless the problem is extremely
00586 */
00587 /*        badly scaled (in which case the exponents should be increased).
00588 */
00589 
00590 
00591 /*  MACHINE DEPENDENCIES */
00592 
00593 /*        The only variables that are machine-dependent are XTOL, */
00594 /*        STPMIN and STPMAX. */
00595 
00596 
00597 /*  GENERAL INFORMATION */
00598 
00599 /*    Other routines called directly:  DAXPY, DDOT, LB1, MCSRCH */
00600 
00601 /*    Input/Output  :  No input; diagnostic messages on unit MP and */
00602 /*                     error messages on unit LP. */
00603 
00604 
00605 /*    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00606 -*/
00607 
00608 
00609     /* Parameter adjustments */
00610     --diag;
00611     --g;
00612     --x;
00613     --w;
00614     --iprint;
00615 
00616     /* Function Body */
00617 
00618 /*     INITIALIZE */
00619 /*     ---------- */
00620 
00621     if (*iflag == 0) {
00622         goto L10;
00623     }
00624     switch (*iflag) {
00625         case 1:  goto L172;
00626         case 2:  goto L100;
00627     }
00628 L10:
00629     //iter = 0;
00630     if (*n <= 0 || *m <= 0) {
00631         goto L196;
00632     }
00633     if (lb3_1.gtol <= 1e-4) {
00634         if (lb3_1.lp > 0) {
00635             //io___4.ciunit = lb3_1.lp;
00636             //s_wsfe(&io___4);
00637             //e_wsfe();
00638         }
00639         lb3_1.gtol = .9;
00640     }
00641     nfun = 1;
00642     point = 0;
00643     finish = FALSE_;
00644     if (*diagco) {
00645         i__1 = *n;
00646         for (i = 1; i <= i__1; ++i) {
00647 /* L30: */
00648             if (diag[i] <= zero) {
00649                 goto L195;
00650             }
00651         }
00652     } else {
00653         i__1 = *n;
00654         for (i = 1; i <= i__1; ++i) {
00655 /* L40: */
00656             diag[i] = 1.;
00657         }
00658     }
00659 
00660 /*     THE WORK VECTOR W IS DIVIDED AS FOLLOWS: */
00661 /*     --------------------------------------- */
00662 /*     THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND */
00663 /*         OTHER TEMPORARY INFORMATION. */
00664 /*     LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO. */
00665 /*     LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED */
00666 /*         IN THE FORMULA THAT COMPUTES H*G. */
00667 /*     LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH */
00668 /*         STEPS. */
00669 /*     LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M */
00670 /*         GRADIENT DIFFERENCES. */
00671 
00672 /*     THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A */
00673 /*     CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT. */
00674 
00675     ispt = *n + (*m << 1);
00676     iypt = ispt + *n * *m;
00677     i__1 = *n;
00678     for (i = 1; i <= i__1; ++i) {
00679 /* L50: */
00680         w[ispt + i] = -g[i] * diag[i];
00681     }
00682     gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
00683     stp1 = one / gnorm;
00684 
00685 /*     PARAMETERS FOR LINE SEARCH ROUTINE */
00686 
00687     ftol = 1e-4;
00688     maxfev = 20;
00689 
00690    /*
00691     if (iprint[1] >= 0) {
00692         lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &
00693                 finish);
00694     }
00695     */
00696 
00697 /*    -------------------- */
00698 /*     MAIN ITERATION LOOP */
00699 /*    -------------------- */
00700 
00701 L80:
00702     ++(*iter);
00703     *info = 0;
00704     bound = *iter - 1;
00705     if (*iter == 1) {
00706         goto L165;
00707     }
00708     if (*iter > *m) {
00709         bound = *m;
00710     }
00711 
00712     ys = ddot_(n, &w[iypt + npt + 1], &c__1, &w[ispt + npt + 1], &c__1);
00713     if (ys==0.0)
00714         ys=1.e-10;
00715     if (! (*diagco)) {
00716         yy = ddot_(n, &w[iypt + npt + 1], &c__1, &w[iypt + npt + 1], &c__1);
00717         i__1 = *n;
00718         for (i = 1; i <= i__1; ++i) {
00719 /* L90: */
00720             diag[i] = ys / yy;
00721         }
00722     } else {
00723         *iflag = 2;
00724         return 0;
00725     }
00726 L100:
00727     if (*diagco) {
00728         i__1 = *n;
00729         for (i = 1; i <= i__1; ++i) {
00730 /* L110: */
00731             if (diag[i] <= zero) {
00732                 goto L195;
00733             }
00734         }
00735     }
00736 
00737 /*     COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, */
00738 /*     "Updating quasi-Newton matrices with limited storage", */
00739 /*     Mathematics of Computation, Vol.24, No.151, pp. 773-782. */
00740 /*     --------------------------------------------------------- */
00741 
00742     cp = point;
00743     if (point == 0) {
00744         cp = *m;
00745     }
00746     w[*n + cp] = one / ys;
00747     i__1 = *n;
00748     for (i = 1; i <= i__1; ++i) {
00749 /* L112: */
00750         w[i] = -g[i];
00751     }
00752     cp = point;
00753     i__1 = bound;
00754     for (i = 1; i <= i__1; ++i) {
00755         --cp;
00756         if (cp == -1) {
00757             cp = *m - 1;
00758         }
00759         sq = ddot_(n, &w[ispt + cp * *n + 1], &c__1, &w[1], &c__1);
00760         inmc = *n + *m + cp + 1;
00761         iycn = iypt + cp * *n;
00762         w[inmc] = w[*n + cp + 1] * sq;
00763         d__1 = -w[inmc];
00764         daxpy_(n, &d__1, &w[iycn + 1], &c__1, &w[1], &c__1);
00765 /* L125: */
00766     }
00767 
00768     i__1 = *n;
00769     for (i = 1; i <= i__1; ++i) {
00770 /* L130: */
00771         w[i] = diag[i] * w[i];
00772     }
00773 
00774     i__1 = bound;
00775     for (i = 1; i <= i__1; ++i) {
00776         yr = ddot_(n, &w[iypt + cp * *n + 1], &c__1, &w[1], &c__1);
00777         beta = w[*n + cp + 1] * yr;
00778         inmc = *n + *m + cp + 1;
00779         beta = w[inmc] - beta;
00780         iscn = ispt + cp * *n;
00781         daxpy_(n, &beta, &w[iscn + 1], &c__1, &w[1], &c__1);
00782         ++cp;
00783         if (cp == *m) {
00784             cp = 0;
00785         }
00786 /* L145: */
00787     }
00788 
00789 /*     STORE THE NEW SEARCH DIRECTION */
00790 /*     ------------------------------ */
00791 
00792     i__1 = *n;
00793     for (i = 1; i <= i__1; ++i) {
00794 /* L160: */
00795         w[ispt + point * *n + i] = w[i];
00796     }
00797 
00798 /*     OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION */
00799 /*     BY USING THE LINE SEARCH ROUTINE MCSRCH */
00800 /*     ---------------------------------------------------- */
00801 L165:
00802     nfev = 0;
00803     stp = one;
00804     if (*iter == 1) {
00805         stp = stp1;
00806     }
00807     i__1 = *n;
00808     for (i = 1; i <= i__1; ++i) {
00809 /* L170: */
00810         w[i] = g[i];
00811     }
00812 L172:
00813     mcsrch_(n, &x[1], f, &g[1], &w[ispt + point * *n + 1], &stp, &ftol, xtol,
00814             &maxfev, info, &nfev, &diag[1]);
00815     if (*info == -1) {
00816         *iflag = 1;
00817         return 0;
00818     }
00819     if (*info != 1) {
00820         goto L190;
00821     }
00822     nfun += nfev;
00823 
00824 /*     COMPUTE THE NEW STEP AND GRADIENT CHANGE */
00825 /*     ----------------------------------------- */
00826 
00827     npt = point * *n;
00828     i__1 = *n;
00829     for (i = 1; i <= i__1; ++i) {
00830         w[ispt + npt + i] = stp * w[ispt + npt + i];
00831 /* L175: */
00832         w[iypt + npt + i] = g[i] - w[i];
00833     }
00834     ++point;
00835     if (point == *m) {
00836         point = 0;
00837     }
00838 
00839 /*     TERMINATION TEST */
00840 /*     ---------------- */
00841 
00842     gnorm = sqrt(ddot_(n, &g[1], &c__1, &g[1], &c__1));
00843     xnorm = sqrt(ddot_(n, &x[1], &c__1, &x[1], &c__1));
00844     xnorm = max(1.,xnorm);
00845     if (gnorm / xnorm <= *eps) {
00846         finish = TRUE_;
00847     }
00848 
00849    /*
00850     if (iprint[1] >= 0) {
00851         lb1_(&iprint[1], &iter, &nfun, &gnorm, n, m, &x[1], f, &g[1], &stp, &
00852                 finish);
00853     }
00854     */
00855     if (finish) {
00856         *iflag = 0;
00857         return 0;
00858     }
00859     goto L80;
00860 
00861 /*     ------------------------------------------------------------ */
00862 /*     END OF MAIN ITERATION LOOP. ERROR EXITS. */
00863 /*     ------------------------------------------------------------ */
00864 
00865 L190:
00866     *iflag = -1;
00867     if (lb3_1.lp > 0) {
00868         //io___30.ciunit = lb3_1.lp;
00869         //s_wsfe(&io___30);
00871         //e_wsfe();
00872     }
00873     return 0;
00874 L195:
00875     *iflag = -2;
00876     if (lb3_1.lp > 0) {
00877         //io___31.ciunit = lb3_1.lp;
00878         //s_wsfe(&io___31);
00880         //e_wsfe();
00881     }
00882     return 0;
00883 L196:
00884     *iflag = -3;
00885     if (lb3_1.lp > 0) {
00886         //io___32.ciunit = lb3_1.lp;
00887         //s_wsfe(&io___32);
00888         //e_wsfe();
00889     }
00890 
00891 /*     FORMATS */
00892 /*     ------- */
00893 
00894     return 0;
00895 } /* lbfgs_ */
00896 
00897 
00898 /*     LAST LINE OF SUBROUTINE LBFGS */
00899 
00900 
00901 /* 
00902 //Subroutine 
00903 int lb1_(integer *iprint, integer *iter, integer *nfun,
00904         doublereal *gnorm, integer *n, integer *m, doublereal *x, doublereal *
00905         f, doublereal *g, doublereal *stp, logical *finish)
00906 {
00907     // Format strings
00908     static char fmt_10[] = "(\002*******************************************\
00909 ******\002)";
00910     static char fmt_20[] = "(\002  N=\002,i5,\002   NUMBER OF CORRECTIONS\
00911 =\002,i2,/,\002       INITIAL VALUES\002)";
00912     static char fmt_30[] = "(\002 F= \002,1pd10.3,\002   GNORM= \002,1pd10.3)"
00913             ;
00914     static char fmt_40[] = "(\002 VECTOR X= \002)";
00915     static char fmt_50[] = "(6(2x,1pd10.3))";
00916     static char fmt_60[] = "(\002 GRADIENT VECTOR G= \002)";
00917     static char fmt_70[] = "(/\002   I   NFN\002,4x,\002FUNC\002,8x,\002GN\
00918 ORM\002,7x,\002STEPLENGTH\002/)";
00919     static char fmt_80[] = "(2(i4,1x),3x,3(1pd10.3,2x))";
00920     static char fmt_90[] = "(\002 FINAL POINT X= \002)";
00921     static char fmt_100[] = "(/\002 THE MINIMIZATION TERMINATED WITHOUT DETE\
00922 CTING ERRORS.\002,/\002 IFLAG = 0\002)";
00923 
00924     // System generated locals
00925     integer i__1;
00926 
00927     // Builtin functions
00928    //integer //s_wsfe(cilist *), e_wsfe(), // do_fio(integer *, char *, ftnlen);
00929 
00930     // Local variables
00931     static integer i;
00932 
00933     // Fortran I/O blocks
00934     static cilist io___33 = { 0, 0, 0, fmt_10, 0 };
00935     static cilist io___34 = { 0, 0, 0, fmt_20, 0 };
00936     static cilist io___35 = { 0, 0, 0, fmt_30, 0 };
00937     static cilist io___36 = { 0, 0, 0, fmt_40, 0 };
00938     static cilist io___37 = { 0, 0, 0, fmt_50, 0 };
00939     static cilist io___39 = { 0, 0, 0, fmt_60, 0 };
00940     static cilist io___40 = { 0, 0, 0, fmt_50, 0 };
00941     static cilist io___41 = { 0, 0, 0, fmt_10, 0 };
00942     static cilist io___42 = { 0, 0, 0, fmt_70, 0 };
00943     static cilist io___43 = { 0, 0, 0, fmt_70, 0 };
00944     static cilist io___44 = { 0, 0, 0, fmt_80, 0 };
00945     static cilist io___45 = { 0, 0, 0, fmt_70, 0 };
00946     static cilist io___46 = { 0, 0, 0, fmt_80, 0 };
00947     static cilist io___47 = { 0, 0, 0, fmt_90, 0 };
00948     static cilist io___48 = { 0, 0, 0, fmt_40, 0 };
00949     static cilist io___49 = { 0, 0, 0, fmt_50, 0 };
00950     static cilist io___50 = { 0, 0, 0, fmt_60, 0 };
00951     static cilist io___51 = { 0, 0, 0, fmt_50, 0 };
00952     static cilist io___52 = { 0, 0, 0, fmt_100, 0 };
00953 
00954 
00955 
00956 //     -------------------------------------------------------------
00957 //     THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND
00958 //     AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT.
00959 //     -------------------------------------------------------------
00960 
00961 
00962     // Parameter adjustments
00963     --iprint;
00964     --g;
00965     --x;
00966 
00967     // Function Body
00968     if (*iter == 0) {
00969         io___33.ciunit = lb3_1.mp;
00970         // s_wsfe(&io___33);
00971         // e_wsfe();
00972         io___34.ciunit = lb3_1.mp;
00973         // s_wsfe(&io___34);
00974         // do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
00975         // do_fio(&c__1, (char *)&(*m), (ftnlen)sizeof(integer));
00976         // e_wsfe();
00977         io___35.ciunit = lb3_1.mp;
00978         // s_wsfe(&io___35);
00979         // do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(doublereal));
00980         // do_fio(&c__1, (char *)&(*gnorm), (ftnlen)sizeof(doublereal));
00981         // e_wsfe();
00982         if (iprint[2] >= 1) {
00983             io___36.ciunit = lb3_1.mp;
00984             //s_wsfe(&io___36);
00985             // e_wsfe();
00986             io___37.ciunit = lb3_1.mp;
00987             //s_wsfe(&io___37);
00988             i__1 = *n;
00989             for (i = 1; i <= i__1; ++i) {
00990                 // do_fio(&c__1, (char *)&x[i], (ftnlen)sizeof(doublereal));
00991             }
00992             // e_wsfe();
00993             io___39.ciunit = lb3_1.mp;
00994             //s_wsfe(&io___39);
00995             // e_wsfe();
00996             io___40.ciunit = lb3_1.mp;
00997             //s_wsfe(&io___40);
00998             i__1 = *n;
00999             for (i = 1; i <= i__1; ++i) {
01000                 // do_fio(&c__1, (char *)&g[i], (ftnlen)sizeof(doublereal));
01001             }
01002             // e_wsfe();
01003         }
01004         io___41.ciunit = lb3_1.mp;
01005         // s_wsfe(&io___41);
01006         // e_wsfe();
01007         io___42.ciunit = lb3_1.mp;
01008         // s_wsfe(&io___42);
01009         // e_wsfe();
01010     } else {
01011         if (iprint[1] == 0 && (*iter != 1 && ! (*finish))) {
01012             return 0;
01013         }
01014         if (iprint[1] != 0) {
01015             if ((*iter - 1) % iprint[1] == 0 || *finish) {
01016                 if (iprint[2] > 1 && *iter > 1) {
01017                     io___43.ciunit = lb3_1.mp;
01018                     //s_wsfe(&io___43);
01019                     // e_wsfe();
01020                 }
01021                 io___44.ciunit = lb3_1.mp;
01022                 // s_wsfe(&io___44);
01023                 // do_fio(&c__1, (char *)&(*iter), (ftnlen)sizeof(integer));
01024                 // do_fio(&c__1, (char *)&(*nfun), (ftnlen)sizeof(integer));
01025                 // do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(doublereal));
01026                 // do_fio(&c__1, (char *)&(*gnorm), (ftnlen)sizeof(doublereal));
01027                 // do_fio(&c__1, (char *)&(*stp), (ftnlen)sizeof(doublereal));
01028                 // e_wsfe();
01029             } else {
01030                 return 0;
01031             }
01032         } else {
01033             if (iprint[2] > 1 && *finish) {
01034                 io___45.ciunit = lb3_1.mp;
01035                 // s_wsfe(&io___45);
01036                 // e_wsfe();
01037             }
01038             io___46.ciunit = lb3_1.mp;
01039             //s_wsfe(&io___46);
01040             // do_fio(&c__1, (char *)&(*iter), (ftnlen)sizeof(integer));
01041             // do_fio(&c__1, (char *)&(*nfun), (ftnlen)sizeof(integer));
01042             // do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(doublereal));
01043             // do_fio(&c__1, (char *)&(*gnorm), (ftnlen)sizeof(doublereal));
01044             // do_fio(&c__1, (char *)&(*stp), (ftnlen)sizeof(doublereal));
01045             // e_wsfe();
01046         }
01047         if (iprint[2] == 2 || iprint[2] == 3) {
01048             if (*finish) {
01049                 io___47.ciunit = lb3_1.mp;
01050                 // s_wsfe(&io___47);
01051                 // e_wsfe();
01052             } else {
01053                 io___48.ciunit = lb3_1.mp;
01054                 // s_wsfe(&io___48);
01055                 // e_wsfe();
01056             }
01057             io___49.ciunit = lb3_1.mp;
01058             //s_wsfe(&io___49);
01059             i__1 = *n;
01060             for (i = 1; i <= i__1; ++i) {
01061                 // do_fio(&c__1, (char *)&x[i], (ftnlen)sizeof(doublereal));
01062             }
01063             // e_wsfe();
01064             if (iprint[2] == 3) {
01065                 io___50.ciunit = lb3_1.mp;
01066                 // s_wsfe(&io___50);
01067                 // e_wsfe();
01068                 io___51.ciunit = lb3_1.mp;
01069                 // s_wsfe(&io___51);
01070                 i__1 = *n;
01071                 for (i = 1; i <= i__1; ++i) {
01072                     // do_fio(&c__1, (char *)&g[i], (ftnlen)sizeof(doublereal));
01073                 }
01074                 // e_wsfe();
01075             }
01076         }
01077         if (*finish) {
01078             io___52.ciunit = lb3_1.mp;
01079             //s_wsfe(&io___52);
01080             // e_wsfe();
01081         }
01082     }
01083 
01084 
01085     return 0;
01086 }
01087 */
01088 
01089 /*     ****** */
01090 
01091 
01092 /*   ---------------------------------------------------------- */
01093 /*     DATA */
01094 /*   ---------------------------------------------------------- */
01095 
01096 
01097 
01098 
01099 /*   ---------------------------------------------------------- */
01100 
01101 /* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal *dx,
01102         integer *incx, doublereal *dy, integer *incy)
01103 {
01104     /* System generated locals */
01105     integer i__1;
01106 
01107     /* Local variables */
01108     static integer i, m, ix, iy, mp1;
01109 
01110 
01111 /*     constant times a vector plus a vector. */
01112 /*     uses unrolled loops for increments equal to one. */
01113 /*     jack dongarra, linpack, 3/11/78. */
01114 
01115 
01116     /* Parameter adjustments */
01117     --dy;
01118     --dx;
01119 
01120     /* Function Body */
01121     if (*n <= 0) {
01122         return 0;
01123     }
01124     if (*da == 0.) {
01125         return 0;
01126     }
01127     if (*incx == 1 && *incy == 1) {
01128         goto L20;
01129     }
01130 
01131 /*        code for unequal increments or equal increments */
01132 /*          not equal to 1 */
01133 
01134     ix = 1;
01135     iy = 1;
01136     if (*incx < 0) {
01137         ix = (-(*n) + 1) * *incx + 1;
01138     }
01139     if (*incy < 0) {
01140         iy = (-(*n) + 1) * *incy + 1;
01141     }
01142     i__1 = *n;
01143     for (i = 1; i <= i__1; ++i) {
01144         dy[iy] += *da * dx[ix];
01145         ix += *incx;
01146         iy += *incy;
01147 /* L10: */
01148     }
01149     return 0;
01150 
01151 /*        code for both increments equal to 1 */
01152 
01153 
01154 /*        clean-up loop */
01155 
01156 L20:
01157     m = *n % 4;
01158     if (m == 0) {
01159         goto L40;
01160     }
01161     i__1 = m;
01162     for (i = 1; i <= i__1; ++i) {
01163         dy[i] += *da * dx[i];
01164 /* L30: */
01165     }
01166     if (*n < 4) {
01167         return 0;
01168     }
01169 L40:
01170     mp1 = m + 1;
01171     i__1 = *n;
01172     for (i = mp1; i <= i__1; i += 4) {
01173         dy[i] += *da * dx[i];
01174         dy[i + 1] += *da * dx[i + 1];
01175         dy[i + 2] += *da * dx[i + 2];
01176         dy[i + 3] += *da * dx[i + 3];
01177 /* L50: */
01178     }
01179     return 0;
01180 } /* daxpy_ */
01181 
01182 
01183 
01184 /*   ---------------------------------------------------------- */
01185 
01186 doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
01187         integer *incy)
01188 {
01189     /* System generated locals */
01190     integer i__1;
01191     doublereal ret_val;
01192 
01193     /* Local variables */
01194     static integer i, m;
01195     static doublereal dtemp;
01196     static integer ix, iy, mp1;
01197 
01198 
01199 /*     forms the dot product of two vectors. */
01200 /*     uses unrolled loops for increments equal to one. */
01201 /*     jack dongarra, linpack, 3/11/78. */
01202 
01203 
01204     /* Parameter adjustments */
01205     --dy;
01206     --dx;
01207 
01208     /* Function Body */
01209     ret_val = 0.;
01210     dtemp = 0.;
01211     if (*n <= 0) {
01212         return ret_val;
01213     }
01214     if (*incx == 1 && *incy == 1) {
01215         goto L20;
01216     }
01217 
01218 /*        code for unequal increments or equal increments */
01219 /*          not equal to 1 */
01220 
01221     ix = 1;
01222     iy = 1;
01223     if (*incx < 0) {
01224         ix = (-(*n) + 1) * *incx + 1;
01225     }
01226     if (*incy < 0) {
01227         iy = (-(*n) + 1) * *incy + 1;
01228     }
01229     i__1 = *n;
01230     for (i = 1; i <= i__1; ++i) {
01231         dtemp += dx[ix] * dy[iy];
01232         ix += *incx;
01233         iy += *incy;
01234 /* L10: */
01235     }
01236     ret_val = dtemp;
01237     return ret_val;
01238 
01239 /*        code for both increments equal to 1 */
01240 
01241 
01242 /*        clean-up loop */
01243 
01244 L20:
01245     m = *n % 5;
01246     if (m == 0) {
01247         goto L40;
01248     }
01249     i__1 = m;
01250     for (i = 1; i <= i__1; ++i) {
01251         dtemp += dx[i] * dy[i];
01252 /* L30: */
01253     }
01254     if (*n < 5) {
01255         goto L60;
01256     }
01257 L40:
01258     mp1 = m + 1;
01259     i__1 = *n;
01260     for (i = mp1; i <= i__1; i += 5) {
01261         dtemp = dtemp + dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] *
01262                 dy[i + 2] + dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4];
01263 /* L50: */
01264     }
01265 L60:
01266     ret_val = dtemp;
01267     return ret_val;
01268 } /* ddot_ */
01269 
01270 /*    ------------------------------------------------------------------ */
01271 
01272 /*     ************************** */
01273 /*     LINE SEARCH ROUTINE MCSRCH */
01274 /*     ************************** */
01275 
01276 /* Subroutine */ int mcsrch_(integer *n, doublereal *x, doublereal *f,
01277         doublereal *g, doublereal *s, doublereal *stp, doublereal *ftol,
01278         doublereal *xtol, integer *maxfev, integer *info, integer *nfev,
01279         doublereal *wa)
01280 {
01281     /* Initialized data */
01282 
01283     static doublereal p5 = .5;
01284     static doublereal p66 = .66;
01285     static doublereal xtrapf = 4.;
01286     static doublereal zero = 0.;
01287 
01288     /* Format strings */
01289 /*
01290     static char fmt_15[] = "(/\002  THE SEARCH DIRECTION IS NOT A DESCENT DI\
01291 RECTION\002)";
01292 */
01293 
01294     /* System generated locals */
01295     integer i__1;
01296     doublereal d__1;
01297 
01298     /* Builtin functions */
01299     //integer //s_wsfe(cilist *), e_wsfe();
01300 
01301     /* Local variables */
01302     static doublereal dgxm, dgym;
01303     static integer j, infoc;
01304     static doublereal finit, width, stmin, stmax;
01305     static logical stage1;
01306     static doublereal width1, ftest1, dg, fm, fx, fy;
01307     static logical brackt;
01308     static doublereal dginit, dgtest;
01309     extern /* Subroutine */ int mcstep_(doublereal *, doublereal *,
01310             doublereal *, doublereal *, doublereal *, doublereal *,
01311             doublereal *, doublereal *, doublereal *, logical *, doublereal *,
01312              doublereal *, integer *);
01313     static doublereal dgm, dgx, dgy, fxm, fym, stx, sty;
01314 
01315     /* Fortran I/O blocks */
01316     //static cilist io___71 = { 0, 0, 0, fmt_15, 0 };
01317 
01318 
01319 
01320 /*                     SUBROUTINE MCSRCH */
01321 
01322 /*     A slight modification of the subroutine CSRCH of More' and Thuente.
01323  */
01324 /*     The changes are to allow reverse communication, and do not affect
01325 */
01326 /*       NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF */
01327 /*         CALLS TO FCN. */
01328 
01329 /*       WA IS A WORK ARRAY OF LENGTH N. */
01330 
01331 /*     SUBPROGRAMS CALLED */
01332 
01333 /*       MCSTEP */
01334 
01335 /*       FORTRAN-SUPPLIED...ABS,MAX,MIN */
01336 
01337 /*     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 */
01338 /*     JORGE J. MORE', DAVID J. THUENTE */
01339 
01340 /*     ********** */
01341     /* Parameter adjustments */
01342     --wa;
01343     --s;
01344     --g;
01345     --x;
01346 
01347     /* Function Body */
01348     if (*info == -1) {
01349         goto L45;
01350     }
01351     infoc = 1;
01352 
01353 /*     CHECK THE INPUT PARAMETERS FOR ERRORS. */
01354 
01355     if (*n <= 0 || *stp <= zero || *ftol < zero || lb3_1.gtol < zero || *xtol
01356             < zero || lb3_1.stpmin < zero || lb3_1.stpmax < lb3_1.stpmin || *
01357             maxfev <= 0) {
01358         return 0;
01359     }
01360 
01361 /*     COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION */
01362 /*     AND CHECK THAT S IS A DESCENT DIRECTION. */
01363 
01364     dginit = zero;
01365     i__1 = *n;
01366     for (j = 1; j <= i__1; ++j) {
01367         dginit += g[j] * s[j];
01368 /* L10: */
01369     }
01370     if (dginit >= zero) {
01371         //io___71.ciunit = lb3_1.lp;
01372         // s_wsfe(&io___71);
01373         // e_wsfe();
01374         return 0;
01375     }
01376 
01377 /*     INITIALIZE LOCAL VARIABLES. */
01378 
01379     brackt = FALSE_;
01380     stage1 = TRUE_;
01381     *nfev = 0;
01382     finit = *f;
01383     dgtest = *ftol * dginit;
01384     width = lb3_1.stpmax - lb3_1.stpmin;
01385     width1 = width / p5;
01386     i__1 = *n;
01387     for (j = 1; j <= i__1; ++j) {
01388         wa[j] = x[j];
01389 /* L20: */
01390     }
01391 
01392 /*     THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP, */
01393 /*     FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP. */
01394 /*     THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP, */
01395 /*     FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF */
01396 /*     THE INTERVAL OF UNCERTAINTY. */
01397 /*     THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP, */
01398 /*     FUNCTION, AND DERIVATIVE AT THE CURRENT STEP. */
01399 
01400     stx = zero;
01401     fx = finit;
01402     dgx = dginit;
01403     sty = zero;
01404     fy = finit;
01405     dgy = dginit;
01406 
01407 /*     START OF ITERATION. */
01408 
01409 L30:
01410 
01411 /*        SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND */
01412 /*        TO THE PRESENT INTERVAL OF UNCERTAINTY. */
01413 
01414     if (brackt) {
01415         stmin = min(stx,sty);
01416         stmax = max(stx,sty);
01417     } else {
01418         stmin = stx;
01419         stmax = *stp + xtrapf * (*stp - stx);
01420     }
01421 
01422 /*        FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN. */
01423 
01424     *stp = max(*stp,lb3_1.stpmin);
01425     *stp = min(*stp,lb3_1.stpmax);
01426 
01427 /*        IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET */
01428 /*        STP BE THE LOWEST POINT OBTAINED SO FAR. */
01429 
01430   if ( ( brackt && (*stp <= stmin || *stp >= stmax) ) || *nfev >= *maxfev - 1 ||
01431             infoc == 0 || (brackt && stmax - stmin <= *xtol * stmax) ) {
01432         *stp = stx;
01433     }
01434 
01435 /*        EVALUATE THE FUNCTION AND GRADIENT AT STP */
01436 /*        AND COMPUTE THE DIRECTIONAL DERIVATIVE. */
01437 /*        We return to main program to obtain F and G. */
01438 
01439     i__1 = *n;
01440     for (j = 1; j <= i__1; ++j) {
01441         x[j] = wa[j] + *stp * s[j];
01442 /* L40: */
01443     }
01444     *info = -1;
01445     return 0;
01446 
01447 L45:
01448     *info = 0;
01449     ++(*nfev);
01450     dg = zero;
01451     i__1 = *n;
01452     for (j = 1; j <= i__1; ++j) {
01453         dg += g[j] * s[j];
01454 /* L50: */
01455     }
01456     ftest1 = finit + *stp * dgtest;
01457 
01458 /*        TEST FOR CONVERGENCE. */
01459 
01460     if ( (brackt && (*stp <= stmin || *stp >= stmax)) || infoc == 0) {
01461         *info = 6;
01462     }
01463     if (*stp == lb3_1.stpmax && *f <= ftest1 && dg <= dgtest) {
01464         *info = 5;
01465     }
01466     if (*stp == lb3_1.stpmin && (*f > ftest1 || dg >= dgtest)) {
01467         *info = 4;
01468     }
01469     if (*nfev >= *maxfev) {
01470         *info = 3;
01471     }
01472     if (brackt && stmax - stmin <= *xtol * stmax) {
01473         *info = 2;
01474     }
01475     if (*f <= ftest1 && abs(dg) <= lb3_1.gtol * (-dginit)) {
01476         *info = 1;
01477     }
01478 
01479 /*        CHECK FOR TERMINATION. */
01480 
01481     if (*info != 0) {
01482         return 0;
01483     }
01484 
01485 /*        IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED */
01486 /*        FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE. */
01487 
01488     if (stage1 && *f <= ftest1 && dg >= min(*ftol,lb3_1.gtol) * dginit) {
01489         stage1 = FALSE_;
01490     }
01491 
01492 /*        A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF */
01493 /*        WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED */
01494 /*        FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE */
01495 /*        DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN */
01496 /*        OBTAINED BUT THE DECREASE IS NOT SUFFICIENT. */
01497 
01498     if (stage1 && *f <= fx && *f > ftest1) {
01499 /*           DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES. */
01500 
01501         fm = *f - *stp * dgtest;
01502         fxm = fx - stx * dgtest;
01503         fym = fy - sty * dgtest;
01504         dgm = dg - dgtest;
01505         dgxm = dgx - dgtest;
01506         dgym = dgy - dgtest;
01507 
01508 /*           CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY */
01509 /*           AND TO COMPUTE THE NEW STEP. */
01510 
01511         mcstep_(&stx, &fxm, &dgxm, &sty, &fym, &dgym, stp, &fm, &dgm, &brackt,
01512                  &stmin, &stmax, &infoc);
01513 
01514 /*           RESET THE FUNCTION AND GRADIENT VALUES FOR F. */
01515 
01516         fx = fxm + stx * dgtest;
01517         fy = fym + sty * dgtest;
01518         dgx = dgxm + dgtest;
01519         dgy = dgym + dgtest;
01520     } else {
01521 /*           CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY */
01522 /*           AND TO COMPUTE THE NEW STEP. */
01523 
01524         mcstep_(&stx, &fx, &dgx, &sty, &fy, &dgy, stp, f, &dg, &brackt, &
01525                 stmin, &stmax, &infoc);
01526     }
01527 
01528 /*        FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE */
01529 /*        INTERVAL OF UNCERTAINTY. */
01530 
01531     if (brackt) {
01532         if ((d__1 = sty - stx, abs(d__1)) >= p66 * width1) {
01533             *stp = stx + p5 * (sty - stx);
01534         }
01535         width1 = width;
01536         width = (d__1 = sty - stx, abs(d__1));
01537     }
01538 
01539 /*        END OF ITERATION. */
01540 
01541     goto L30;
01542 
01543 /*     LAST LINE OF SUBROUTINE MCSRCH. */
01544 } /* mcsrch_ */
01545 
01546 /* Subroutine */ int mcstep_(doublereal *stx, doublereal *fx, doublereal *dx,
01547         doublereal *sty, doublereal *fy, doublereal *dy, doublereal *stp,
01548         doublereal *fp, doublereal *dp, logical *brackt, doublereal *stpmin,
01549         doublereal *stpmax, integer *info)
01550 {
01551     /* System generated locals */
01552     doublereal d__1, d__2, d__3;
01553 
01554     /* Builtin functions */
01555     double sqrt(doublereal);
01556 
01557     /* Local variables */
01558     static doublereal sgnd, stpc, stpf, stpq, p, q, gamma, r, s, theta;
01559     static logical bound;
01560 
01561 
01562 /*     SUBROUTINE MCSTEP */
01563 
01564 /*     THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR */
01565 /*     A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR */
01566 /*     A MINIMIZER OF THE FUNCTION. */
01567 
01568 /*     THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION */
01569 /*     VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS */
01570 /*     ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE */
01571 /*     DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A */
01572 /*     MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY */
01573 /*     WITH ENDPOINTS STX AND STY. */
01574 
01575 /*     THE SUBROUTINE STATEMENT IS */
01576 
01577 /*       SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT, */
01578 /*                        STPMIN,STPMAX,INFO) */
01579 
01580 /*     WHERE */
01581 
01582 /*       STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP, */
01583 /*         THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED */
01584 /*         SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION */
01585 /*         OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE */
01586 /*         SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY. */
01587 
01588 /*       STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP, */
01589 /*         THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF */
01590 /*         THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE */
01591 /*         UPDATED APPROPRIATELY. */
01592 
01593 /*       STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP, */
01594 /*         THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP. */
01595 /*         IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE */
01596 /*         BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP. */
01597 
01598 /*       BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER */
01599 /*         HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED */
01600 /*         THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER */
01601 /*         IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE. */
01602 
01603 /*       STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER */
01604 /*         AND UPPER BOUNDS FOR THE STEP. */
01605 
01606 /*       INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS: */
01607 /*         IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED */
01608 /*         ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE */
01609 /*         INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS. */
01610 
01611 /*     SUBPROGRAMS CALLED */
01612 
01613 /*       FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT */
01614 
01615 /*     ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983 */
01616 /*     JORGE J. MORE', DAVID J. THUENTE */
01617 
01618     *info = 0;
01619 
01620 /*     CHECK THE INPUT PARAMETERS FOR ERRORS. */
01621 
01622     if ((*brackt && (*stp <= min(*stx,*sty) || *stp >= max(*stx,*sty))) || *dx *
01623              (*stp - *stx) >= (float)0. || *stpmax < *stpmin) {
01624         return 0;
01625     }
01626 
01627 /*     DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN. */
01628 
01629     sgnd = *dp * (*dx / abs(*dx));
01630 
01631 /*     FIRST CASE. A HIGHER FUNCTION VALUE. */
01632 /*     THE MINIMUM IS BRACKETED. IF THE CUBIC STEP IS CLOSER */
01633 /*     TO STX THAN THE QUADRATIC STEP, THE CUBIC STEP IS TAKEN, */
01634 /*     ELSE THE AVERAGE OF THE CUBIC AND QUADRATIC STEPS IS TAKEN. */
01635 
01636     if (*fp > *fx) {
01637         *info = 1;
01638         bound = TRUE_;
01639         theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
01640 /* Computing MAX */
01641         d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
01642                 *dp);
01643         s = max(d__1,d__2);
01644 /* Computing 2nd power */
01645         d__1 = theta / s;
01646         gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
01647         if (*stp < *stx) {
01648             gamma = -gamma;
01649         }
01650         p = gamma - *dx + theta;
01651         q = gamma - *dx + gamma + *dp;
01652         r = p / q;
01653         stpc = *stx + r * (*stp - *stx);
01654         stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2 * (*stp -
01655                 *stx);
01656         if ((d__1 = stpc - *stx, abs(d__1)) < (d__2 = stpq - *stx, abs(d__2)))
01657                  {
01658             stpf = stpc;
01659         } else {
01660             stpf = stpc + (stpq - stpc) / 2;
01661         }
01662         *brackt = TRUE_;
01663 
01664 /*     SECOND CASE. A LOWER FUNCTION VALUE AND DERIVATIVES OF */
01665 /*     OPPOSITE SIGN. THE MINIMUM IS BRACKETED. IF THE CUBIC */
01666 /*     STEP IS CLOSER TO STX THAN THE QUADRATIC (SECANT) STEP, */
01667 /*     THE CUBIC STEP IS TAKEN, ELSE THE QUADRATIC STEP IS TAKEN. */
01668 
01669     } else if (sgnd < (float)0.) {
01670         *info = 2;
01671         bound = FALSE_;
01672         theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
01673 /* Computing MAX */
01674         d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
01675                 *dp);
01676         s = max(d__1,d__2);
01677 /* Computing 2nd power */
01678         d__1 = theta / s;
01679         gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
01680         if (*stp > *stx) {
01681             gamma = -gamma;
01682         }
01683         p = gamma - *dp + theta;
01684         q = gamma - *dp + gamma + *dx;
01685         r = p / q;
01686         stpc = *stp + r * (*stx - *stp);
01687         stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
01688         if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(d__2)))
01689                  {
01690             stpf = stpc;
01691         } else {
01692             stpf = stpq;
01693         }
01694         *brackt = TRUE_;
01695 
01696 /*     THIRD CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
01697 /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DECREASES. */
01698 /*     THE CUBIC STEP IS ONLY USED IF THE CUBIC TENDS TO INFINITY */
01699 /*     IN THE DIRECTION OF THE STEP OR IF THE MINIMUM OF THE CUBIC */
01700 /*     IS BEYOND STP. OTHERWISE THE CUBIC STEP IS DEFINED TO BE */
01701 /*     EITHER STPMIN OR STPMAX. THE QUADRATIC (SECANT) STEP IS ALSO */
01702 /*     COMPUTED AND IF THE MINIMUM IS BRACKETED THEN THE THE STEP */
01703 /*     CLOSEST TO STX IS TAKEN, ELSE THE STEP FARTHEST AWAY IS TAKEN.
01704 */
01705 
01706     } else if (abs(*dp) < abs(*dx)) {
01707         *info = 3;
01708         bound = TRUE_;
01709         theta = (*fx - *fp) * 3 / (*stp - *stx) + *dx + *dp;
01710 /* Computing MAX */
01711         d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs(
01712                 *dp);
01713         s = max(d__1,d__2);
01714 
01715 /*        THE CASE GAMMA = 0 ONLY ARISES IF THE CUBIC DOES NOT TEND */
01716 /*        TO INFINITY IN THE DIRECTION OF THE STEP. */
01717 
01718 /* Computing MAX */
01719 /* Computing 2nd power */
01720         d__3 = theta / s;
01721         d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s);
01722         gamma = s * sqrt((max(d__1,d__2)));
01723         if (*stp > *stx) {
01724             gamma = -gamma;
01725         }
01726         p = gamma - *dp + theta;
01727         q = gamma + (*dx - *dp) + gamma;
01728         r = p / q;
01729         if (r < (float)0. && gamma != (float)0.) {
01730             stpc = *stp + r * (*stx - *stp);
01731         } else if (*stp > *stx) {
01732             stpc = *stpmax;
01733         } else {
01734             stpc = *stpmin;
01735         }
01736         stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
01737         if (*brackt) {
01738             if ((d__1 = *stp - stpc, abs(d__1)) < (d__2 = *stp - stpq, abs(
01739                     d__2))) {
01740                 stpf = stpc;
01741             } else {
01742                 stpf = stpq;
01743             }
01744         } else {
01745             if ((d__1 = *stp - stpc, abs(d__1)) > (d__2 = *stp - stpq, abs(
01746                     d__2))) {
01747                 stpf = stpc;
01748             } else {
01749                 stpf = stpq;
01750             }
01751         }
01752 
01753 /*     FOURTH CASE. A LOWER FUNCTION VALUE, DERIVATIVES OF THE */
01754 /*     SAME SIGN, AND THE MAGNITUDE OF THE DERIVATIVE DOES */
01755 /*     NOT DECREASE. IF THE MINIMUM IS NOT BRACKETED, THE STEP */
01756 /*     IS EITHER STPMIN OR STPMAX, ELSE THE CUBIC STEP IS TAKEN. */
01757 
01758     } else {
01759         *info = 4;
01760         bound = FALSE_;
01761         if (*brackt) {
01762             theta = (*fp - *fy) * 3 / (*sty - *stp) + *dy + *dp;
01763 /* Computing MAX */
01764             d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 =
01765                     abs(*dp);
01766             s = max(d__1,d__2);
01767 /* Computing 2nd power */
01768             d__1 = theta / s;
01769             gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s));
01770             if (*stp > *sty) {
01771                 gamma = -gamma;
01772             }
01773             p = gamma - *dp + theta;
01774             q = gamma - *dp + gamma + *dy;
01775             r = p / q;
01776             stpc = *stp + r * (*sty - *stp);
01777             stpf = stpc;
01778         } else if (*stp > *stx) {
01779             stpf = *stpmax;
01780         } else {
01781             stpf = *stpmin;
01782         }
01783     }
01784 
01785 /*     UPDATE THE INTERVAL OF UNCERTAINTY. THIS UPDATE DOES NOT */
01786 /*     DEPEND ON THE NEW STEP OR THE CASE ANALYSIS ABOVE. */
01787 
01788     if (*fp > *fx) {
01789         *sty = *stp;
01790         *fy = *fp;
01791         *dy = *dp;
01792     } else {
01793         if (sgnd < (float)0.) {
01794             *sty = *stx;
01795             *fy = *fx;
01796             *dy = *dx;
01797         }
01798         *stx = *stp;
01799         *fx = *fp;
01800         *dx = *dp;
01801     }
01802 
01803 /*     COMPUTE THE NEW STEP AND SAFEGUARD IT. */
01804 
01805     stpf = min(*stpmax,stpf);
01806     stpf = max(*stpmin,stpf);
01807     *stp = stpf;
01808     if (*brackt && bound) {
01809         if (*sty > *stx) {
01810 /* Computing MIN */
01811             d__1 = *stx + (*sty - *stx) * (float).66;
01812             *stp = min(d__1,*stp);
01813         } else {
01814 /* Computing MAX */
01815             d__1 = *stx + (*sty - *stx) * (float).66;
01816             *stp = max(d__1,*stp);
01817         }
01818     }
01819     return 0;
01820 
01821 /*     LAST LINE OF SUBROUTINE MCSTEP. */
01822 } /* mcstep_ */
01823 
01824 #ifdef __cplusplus
01825         }
01826 #endif