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