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