--- /dev/null
+/* tensor.f -- translated by f2c (version 20050501).
+ You must link the resulting object file with libf2c:
+ on Microsoft Windows system, link with libf2c.lib;
+ on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+ or, if you install libf2c.a in a standard place, with -lf2c -lm
+ -- in that order, at the end of the command line, as in
+ cc *.o -lf2c -lm
+ Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+ http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#include "f2c.h"
+
+/* Table of constant values */
+
+static integer c__5 = 5;
+static integer c__1 = 1;
+static integer c__9 = 9;
+static integer c__3 = 3;
+static doublereal c_b111 = 2.;
+static doublereal c_b134 = 10.;
+static doublereal c_b250 = .33333333333333331;
+static doublereal c_b324 = 1.;
+static doublereal c_b384 = 3.;
+
+/* ALGORITHM 739, COLLECTED ALGORITHMS FROM ACM. */
+/* THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
+/* VOL. 20, NO. 4, DECEMBER, 1994, P.518-530. */
+/* *** driver.f */
+/* Main program */ int MAIN__(void)
+{
+ /* Format strings */
+ static char fmt_211[] = "(\002 G=\002,999e20.13)";
+
+ /* System generated locals */
+ integer i__1;
+ olist o__1;
+ cllist cl__1;
+ alist al__1, al__2;
+
+ /* Builtin functions */
+ integer f_open(olist *), s_rsle(cilist *), do_lio(integer *, integer *,
+ char *, ftnlen), e_rsle(void), s_wsle(cilist *), e_wsle(void),
+ s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
+ f_rew(alist *), f_end(alist *), f_clos(cllist *);
+ /* Subroutine */ int s_stop(char *, ftnlen);
+
+ /* Local variables */
+ doublereal h__[200] /* was [50][4] */;
+ integer i__, n;
+ doublereal x[50];
+ integer nr;
+ extern /* Subroutine */ int fcn_(), grd_();
+ integer msg;
+ extern /* Subroutine */ int hsn_();
+ integer ipr;
+ doublereal wrk[400] /* was [50][8] */, fpls, gpls[50];
+ integer iwrk[4];
+ doublereal xpls[50], typx[50];
+ integer itnno;
+ doublereal fscale;
+ integer grdflg, hesflg;
+ doublereal gradtl;
+ integer ndigit;
+ extern /* Subroutine */ int dfault_(integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *, doublereal *, integer *,
+ integer *, integer *, integer *, integer *, integer *);
+ integer method;
+ extern /* Subroutine */ int prtfcn_(integer *);
+ integer itnlim;
+ extern /* Subroutine */ int tensrd_(integer *, integer *, doublereal *,
+ U_fp, integer *, doublereal *, doublereal *, doublereal *,
+ doublereal *, integer *, doublereal *, integer *), tensor_(
+ integer *, integer *, doublereal *, U_fp, U_fp, U_fp, doublereal *
+ , doublereal *, doublereal *, doublereal *, integer *, doublereal
+ *, integer *, integer *, integer *, integer *, integer *, integer
+ *, doublereal *, doublereal *, doublereal *, doublereal *,
+ integer *, doublereal *, integer *);
+ doublereal steptl, stepmx;
+
+ /* Fortran I/O blocks */
+ static cilist io___3 = { 0, 5, 0, 0, 0 };
+ static cilist io___6 = { 0, 6, 0, 0, 0 };
+ static cilist io___7 = { 0, 6, 0, 0, 0 };
+ static cilist io___8 = { 0, 6, 0, 0, 0 };
+ static cilist io___9 = { 0, 6, 0, 0, 0 };
+ static cilist io___10 = { 0, 6, 0, 0, 0 };
+ static cilist io___11 = { 0, 6, 0, 0, 0 };
+ static cilist io___12 = { 0, 6, 0, 0, 0 };
+ static cilist io___21 = { 0, 6, 0, 0, 0 };
+ static cilist io___22 = { 0, 6, 0, fmt_211, 0 };
+ static cilist io___23 = { 0, 6, 0, 0, 0 };
+ static cilist io___24 = { 0, 6, 0, 0, 0 };
+ static cilist io___25 = { 0, 5, 0, 0, 0 };
+ static cilist io___26 = { 0, 6, 0, 0, 0 };
+ static cilist io___27 = { 0, 6, 0, 0, 0 };
+ static cilist io___28 = { 0, 6, 0, 0, 0 };
+ static cilist io___29 = { 0, 6, 0, 0, 0 };
+ static cilist io___30 = { 0, 6, 0, 0, 0 };
+ static cilist io___31 = { 0, 6, 0, 0, 0 };
+ static cilist io___32 = { 0, 6, 0, 0, 0 };
+ static cilist io___33 = { 0, 6, 0, 0, 0 };
+ static cilist io___45 = { 0, 6, 0, 0, 0 };
+ static cilist io___46 = { 0, 6, 0, fmt_211, 0 };
+ static cilist io___47 = { 0, 6, 0, 0, 0 };
+ static cilist io___48 = { 0, 6, 0, 0, 0 };
+
+
+ nr = 50;
+ o__1.oerr = 0;
+ o__1.ounit = 5;
+ o__1.ofnmlen = 8;
+ o__1.ofnm = "wood.inp";
+ o__1.orl = 0;
+ o__1.osta = 0;
+ o__1.oacc = 0;
+ o__1.ofm = 0;
+ o__1.oblnk = 0;
+ f_open(&o__1);
+ prtfcn_(&n);
+ s_rsle(&io___3);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_rsle();
+ s_wsle(&io___6);
+ do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
+ 41);
+ e_wsle();
+ s_wsle(&io___7);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_wsle();
+/* SHORT FORM */
+ s_wsle(&io___8);
+ do_lio(&c__9, &c__1, " ", (ftnlen)1);
+ e_wsle();
+ s_wsle(&io___9);
+ do_lio(&c__9, &c__1, "CALLING TENSRD - ALL INPUT PARAMETERS ARE SET", (
+ ftnlen)46);
+ e_wsle();
+ s_wsle(&io___10);
+ do_lio(&c__9, &c__1, " TO THEIR DEFAULT VALUES.", (ftnlen)
+ 41);
+ e_wsle();
+ s_wsle(&io___11);
+ do_lio(&c__9, &c__1, " OUTPUT WILL BE WRITTEN TO THE STANDARD OUTPUT.", (
+ ftnlen)47);
+ e_wsle();
+ s_wsle(&io___12);
+ do_lio(&c__9, &c__1, " ", (ftnlen)1);
+ e_wsle();
+ tensrd_(&nr, &n, x, (U_fp)fcn_, &msg, xpls, &fpls, gpls, h__, &itnno, wrk,
+ iwrk);
+ s_wsle(&io___21);
+ do_lio(&c__9, &c__1, "RESULTS OF TENSRD:", (ftnlen)18);
+ e_wsle();
+ s_wsfe(&io___22);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ s_wsle(&io___23);
+ do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
+ doublereal));
+ }
+ e_wsle();
+ s_wsle(&io___24);
+ do_lio(&c__9, &c__1, "FPLS=", (ftnlen)5);
+ do_lio(&c__5, &c__1, (char *)&fpls, (ftnlen)sizeof(doublereal));
+ do_lio(&c__9, &c__1, " ITNNO=", (ftnlen)8);
+ do_lio(&c__3, &c__1, (char *)&itnno, (ftnlen)sizeof(integer));
+ do_lio(&c__9, &c__1, " MSG=", (ftnlen)6);
+ do_lio(&c__3, &c__1, (char *)&msg, (ftnlen)sizeof(integer));
+ e_wsle();
+/* LONG FORM */
+ prtfcn_(&n);
+ al__1.aerr = 0;
+ al__1.aunit = 5;
+ f_rew(&al__1);
+ s_rsle(&io___25);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_rsle();
+ s_wsle(&io___26);
+ do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
+ 41);
+ e_wsle();
+ s_wsle(&io___27);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_wsle();
+ s_wsle(&io___28);
+ do_lio(&c__9, &c__1, " ", (ftnlen)1);
+ e_wsle();
+ s_wsle(&io___29);
+ do_lio(&c__9, &c__1, "CALLING TENSOR AFTER RESETTING INPUT PARAMETERS", (
+ ftnlen)47);
+ e_wsle();
+ s_wsle(&io___30);
+ do_lio(&c__9, &c__1, "IPR, MSG AND METHOD.", (ftnlen)20);
+ e_wsle();
+ s_wsle(&io___31);
+ do_lio(&c__9, &c__1, "THE STANDARD METHOD IS USED IN THIS EXAMPLE.", (
+ ftnlen)44);
+ e_wsle();
+ s_wsle(&io___32);
+ do_lio(&c__9, &c__1, "ADDITIONAL OUTPUT WILL BE WRITTEN TO FILE 'FORT8'.",
+ (ftnlen)50);
+ e_wsle();
+ s_wsle(&io___33);
+ do_lio(&c__9, &c__1, " ", (ftnlen)1);
+ e_wsle();
+/* L997: */
+ dfault_(&n, typx, &fscale, &gradtl, &steptl, &itnlim, &stepmx, &ipr, &
+ method, &grdflg, &hesflg, &ndigit, &msg);
+ o__1.oerr = 0;
+ o__1.ounit = 8;
+ o__1.ofnmlen = 5;
+ o__1.ofnm = "FORT8";
+ o__1.orl = 0;
+ o__1.osta = 0;
+ o__1.oacc = 0;
+ o__1.ofm = 0;
+ o__1.oblnk = 0;
+ f_open(&o__1);
+ ipr = 8;
+ msg = 2;
+ method = 0;
+ tensor_(&nr, &n, x, (U_fp)fcn_, (U_fp)grd_, (U_fp)hsn_, typx, &fscale, &
+ gradtl, &steptl, &itnlim, &stepmx, &ipr, &method, &grdflg, &
+ hesflg, &ndigit, &msg, xpls, &fpls, gpls, h__, &itnno, wrk, iwrk);
+ s_wsle(&io___45);
+ do_lio(&c__9, &c__1, "RESULTS OF TENSOR:", (ftnlen)18);
+ e_wsle();
+ s_wsfe(&io___46);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ s_wsle(&io___47);
+ do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
+ i__1 = n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
+ doublereal));
+ }
+ e_wsle();
+ s_wsle(&io___48);
+ do_lio(&c__9, &c__1, "FPLS=", (ftnlen)5);
+ do_lio(&c__5, &c__1, (char *)&fpls, (ftnlen)sizeof(doublereal));
+ do_lio(&c__9, &c__1, " ITNNO=", (ftnlen)8);
+ do_lio(&c__3, &c__1, (char *)&itnno, (ftnlen)sizeof(integer));
+ do_lio(&c__9, &c__1, " MSG=", (ftnlen)6);
+ do_lio(&c__3, &c__1, (char *)&msg, (ftnlen)sizeof(integer));
+ e_wsle();
+ al__2.aerr = 0;
+ al__2.aunit = 8;
+ f_end(&al__2);
+ cl__1.cerr = 0;
+ cl__1.cunit = 8;
+ cl__1.csta = 0;
+ f_clos(&cl__1);
+ cl__1.cerr = 0;
+ cl__1.cunit = 5;
+ cl__1.csta = 0;
+ f_clos(&cl__1);
+ s_stop("", (ftnlen)0);
+ return 0;
+} /* MAIN__ */
+
+/* Subroutine */ int fcn_(integer *n, doublereal *x, doublereal *f)
+{
+ /* System generated locals */
+ doublereal d__1, d__2, d__3, d__4, d__5, d__6;
+
+ /* Builtin functions */
+ double pow_dd(doublereal *, doublereal *);
+
+ /* Parameter adjustments */
+ --x;
+
+ /* Function Body */
+ if (*n <= 0) {
+ *n = 4;
+ } else {
+ d__1 = x[1] * x[1] - x[2];
+ d__2 = 1. - x[1];
+ d__3 = x[3] * x[3] - x[4];
+ d__4 = 1. - x[3];
+ d__5 = 1. - x[2];
+ d__6 = 1. - x[4];
+ *f = pow_dd(&d__1, &c_b111) * 100. + pow_dd(&d__2, &c_b111) + pow_dd(&
+ d__3, &c_b111) * 90. + pow_dd(&d__4, &c_b111) + (pow_dd(&d__5,
+ &c_b111) + pow_dd(&d__6, &c_b111)) * 10.1 + (1. - x[2]) *
+ 19.8 * (1. - x[4]);
+ }
+ return 0;
+} /* fcn_ */
+
+/* Subroutine */ int prtfcn_(integer *n)
+{
+ /* Builtin functions */
+ integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
+ e_wsle(void);
+
+ /* Fortran I/O blocks */
+ static cilist io___49 = { 0, 6, 0, 0, 0 };
+ static cilist io___50 = { 0, 6, 0, 0, 0 };
+ static cilist io___51 = { 0, 6, 0, 0, 0 };
+
+
+ *n = 4;
+ s_wsle(&io___49);
+ do_lio(&c__9, &c__1, "__________________________________________", (
+ ftnlen)42);
+ e_wsle();
+ s_wsle(&io___50);
+ do_lio(&c__9, &c__1, "f(x)= Wood function", (ftnlen)19);
+ e_wsle();
+ s_wsle(&io___51);
+ do_lio(&c__9, &c__1, "__________________________________________", (
+ ftnlen)42);
+ e_wsle();
+ return 0;
+} /* prtfcn_ */
+
+/* ---------------------- */
+/* | G R D | */
+/* ---------------------- */
+/* Subroutine */ int grd_(integer *n, real *x, real *g)
+{
+/* DUMMY ROUTINE */
+ return 0;
+} /* grd_ */
+
+/* ---------------------- */
+/* | H S N | */
+/* ---------------------- */
+/* Subroutine */ int hsn_(integer *nr, integer *n, real *x, real *h__)
+{
+/* DUMMY ROUTINE */
+ return 0;
+} /* hsn_ */
+
+/* *** tensrd.f */
+/* ---------------------- */
+/* | T E N S O R | */
+/* ---------------------- */
+/* Subroutine */ int tensor_(integer *nr, integer *n, doublereal *x, U_fp fcn,
+ U_fp grd, U_fp hsn, doublereal *typx, doublereal *fscale, doublereal
+ *gradtl, doublereal *steptl, integer *itnlim, doublereal *stepmx,
+ integer *ipr, integer *method, integer *grdflg, integer *hesflg,
+ integer *ndigit, integer *msg, doublereal *xpls, doublereal *fpls,
+ doublereal *gpls, doublereal *h__, integer *itnno, doublereal *wrk,
+ integer *iwrk)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, wrk_dim1, wrk_offset;
+
+ /* Local variables */
+ extern /* Subroutine */ int opt_(integer *, integer *, doublereal *, U_fp,
+ U_fp, U_fp, doublereal *, doublereal *, doublereal *, doublereal
+ *, integer *, doublereal *, integer *, integer *, integer *,
+ integer *, integer *, integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, integer *);
+
+
+/* AUTHORS: TA-TUNG CHOW, ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
+
+/* DATE: MARCH 29, 1991 */
+
+/* PURPOSE: */
+/* DERIVATIVE TENSOR METHOD FOR UNCONSTRAINED OPTIMIZATION */
+/* THE METHOD BASES EACH ITERATION ON A SPECIALLY CONSTRUCTED */
+/* FOURTH ORDER MODEL OF THE OBJECTIVE FUNCTION. THE MODEL */
+/* INTERPOLATES THE FUNCTION VALUE AND GRADIENT FROM THE PREVIOUS */
+/* ITERATE AND THE CURRENT FUNCTION VALUE, GRADIENT AND HESSIAN */
+/* MATRIX. */
+
+/* BLAS SUBROUTINES: DCOPY,DDOT,DSCAL */
+/* UNCMIN SUBROUTINES: DFAULT, OPTCHK, GRDCHK, HESCHK, LNSRCH, FSTOFD, */
+/* SNDOFD, BAKSLV, FORSLV, OPTSTP */
+/* MODCHL SUBROUTINES: MODCHL, INIT, GERCH, FIN2X2 */
+
+/* ----------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
+/* GRD --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT */
+/* HSN --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN */
+/* HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
+/* AND DIAGONAL OF H */
+/* TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X */
+/* FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
+/* GRADTL --> GRADIENT TOLERANCE */
+/* STEPTL --> STEP TOLERANCE */
+/* ITNLIM --> ITERATION LIMIT */
+/* STEPMX --> MAXIMUM STEP LENGTH ALLOWED */
+/* IPR --> OUTPUT UNIT NUMBER */
+/* METHOD --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT */
+/* EACH ITERATION */
+/* IF VALUE IS 1 THEN TRY BOTH TENSOR AND NEWTON */
+/* STEPS AT EACH ITERATION */
+/* GRDFLG --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED */
+/* HESFLG --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED */
+/* NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
+/* MSG --> OUTPUT MESSAGE CONTROL */
+/* XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) */
+/* FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
+/* GPLS(N) <-- CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) */
+/* H(N,N) --> HESSIAN */
+/* ITNNO <-- NUMBER OF ITERATIONS */
+/* WRK(N,8)--> WORK SPACE */
+/* IWRK(N) --> WORK SPACE */
+
+
+/* EQUIVALENCE WRK(N,1) = G(N) */
+/* WRK(N,2) = S(N) */
+/* WRK(N,3) = D(N) */
+/* WRK(N,4) = DN(N) */
+/* WRK(N,6) = E(N) */
+/* WRK(N,7) = WK1(N) */
+/* WRK(N,8) = WK2(N) */
+
+ /* Parameter adjustments */
+ wrk_dim1 = *nr;
+ wrk_offset = 1 + wrk_dim1;
+ wrk -= wrk_offset;
+ --iwrk;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+ --gpls;
+ --xpls;
+ --typx;
+ --x;
+
+ /* Function Body */
+ opt_(nr, n, &x[1], (U_fp)fcn, (U_fp)grd, (U_fp)hsn, &typx[1], fscale,
+ gradtl, steptl, itnlim, stepmx, ipr, method, grdflg, hesflg,
+ ndigit, msg, &xpls[1], fpls, &gpls[1], &h__[h_offset], itnno, &
+ wrk[wrk_dim1 + 1], &wrk[(wrk_dim1 << 1) + 1], &wrk[wrk_dim1 * 3 +
+ 1], &wrk[(wrk_dim1 << 2) + 1], &wrk[wrk_dim1 * 6 + 1], &wrk[
+ wrk_dim1 * 7 + 1], &wrk[(wrk_dim1 << 3) + 1], &iwrk[1]);
+ return 0;
+} /* tensor_ */
+
+/* ---------------------- */
+/* | T E N S R D | */
+/* ---------------------- */
+/* Subroutine */ int tensrd_(integer *nr, integer *n, doublereal *x, U_fp fcn,
+ integer *msg, doublereal *xpls, doublereal *fpls, doublereal *gpls,
+ doublereal *h__, integer *itnno, doublereal *wrk, integer *iwrk)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, wrk_dim1, wrk_offset;
+
+ /* Local variables */
+ extern /* Subroutine */ int grd_(integer *, real *, real *), hsn_(integer
+ *, integer *, real *, real *);
+ integer ipr;
+ doublereal fscale;
+ integer grdflg, hesflg;
+ doublereal gradtl;
+ integer ndigit;
+ extern /* Subroutine */ int dfault_(integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *, doublereal *, integer *,
+ integer *, integer *, integer *, integer *, integer *);
+ integer method, itnlim;
+ extern /* Subroutine */ int tensor_(integer *, integer *, doublereal *,
+ U_fp, S_fp, S_fp, doublereal *, doublereal *, doublereal *,
+ doublereal *, integer *, doublereal *, integer *, integer *,
+ integer *, integer *, integer *, integer *, doublereal *,
+ doublereal *, doublereal *, doublereal *, integer *, doublereal *,
+ integer *);
+ doublereal steptl, stepmx;
+
+
+/* PURPOSE: */
+/* SHORT CALLING SEQUENCE SUBROUTINE */
+
+/* ----------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
+/* MSG --> OUTPUT MESSAGE CONTROL */
+/* XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) */
+/* FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
+/* GPLS(N) <-- GRADIENT AT FINAL POINT (OUTPUT) */
+/* H(N,N) --> HESSIAN */
+/* ITNNO <-- NUMBER OF ITERATIONS */
+/* WRK(N,8) --> WORK SPACE */
+/* IWRK(N) --> WORK SPACE */
+
+
+/* EQUIVALENCE WRK(N,1) = G(N) */
+/* WRK(N,2) = S(N) */
+/* WRK(N,3) = D(N) */
+/* WRK(N,4) = DN(N) */
+/* WRK(N,5) = TYPX(N) */
+/* WRK(N,6) = E(N) */
+/* WRK(N,7) = WK1(N) */
+/* WRK(N,8) = WK2(N) */
+
+ /* Parameter adjustments */
+ wrk_dim1 = *nr;
+ wrk_offset = 1 + wrk_dim1;
+ wrk -= wrk_offset;
+ --iwrk;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+ --gpls;
+ --xpls;
+ --x;
+
+ /* Function Body */
+ dfault_(n, &wrk[wrk_dim1 * 5 + 1], &fscale, &gradtl, &steptl, &itnlim, &
+ stepmx, &ipr, &method, &grdflg, &hesflg, &ndigit, msg);
+ tensor_(nr, n, &x[1], (U_fp)fcn, (S_fp)grd_, (S_fp)hsn_, &wrk[wrk_dim1 *
+ 5 + 1], &fscale, &gradtl, &steptl, &itnlim, &stepmx, &ipr, &
+ method, &grdflg, &hesflg, &ndigit, msg, &xpls[1], fpls, &gpls[1],
+ &h__[h_offset], itnno, &wrk[wrk_offset], &iwrk[1]);
+ return 0;
+} /* tensrd_ */
+
+/* ---------------------- */
+/* | O P T | */
+/* ---------------------- */
+/* Subroutine */ int opt_(integer *nr, integer *n, doublereal *x, S_fp fcn,
+ S_fp grd, S_fp hsn, doublereal *typx, doublereal *fscale, doublereal *
+ gradtl, doublereal *steptl, integer *itnlim, doublereal *stepmx,
+ integer *ipr, integer *method, integer *grdflg, integer *hesflg,
+ integer *ndigit, integer *msg, doublereal *xpls, doublereal *fpls,
+ doublereal *gpls, doublereal *h__, integer *itnno, doublereal *g,
+ doublereal *s, doublereal *d__, doublereal *dn, doublereal *e,
+ doublereal *wk1, doublereal *wk2, integer *pivot)
+{
+ /* Format strings */
+ static char fmt_25[] = "(\002 INITIAL FUNCTION VALUE F=\002,e20.13)";
+ static char fmt_30[] = "(\002 INITIAL GRADIENT G=\002,999e20.13)";
+ static char fmt_103[] = "(\002 ----------- ITERATION \002,i4,\002 ---"
+ "-------------\002)";
+ static char fmt_104[] = "(\002 X=\002,999e20.13)";
+ static char fmt_105[] = "(\002 F(X)=\002,e20.13)";
+ static char fmt_106[] = "(\002 G(X)=\002,999e20.13)";
+ static char fmt_108[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
+ ". MAX:\002,e20.13)";
+ static char fmt_110[] = "(\002REL. STEP MAX :\002,e20.13)";
+ static char fmt_370[] = "(\002 FINAL X=\002,999e20.13)";
+ static char fmt_380[] = "(\002 GRADIENT G=\002,999e20.13)";
+ static char fmt_390[] = "(\002 FUNCTION VALUE F(X)=\002,e20.13,\002 AT I"
+ "TERATION \002,i4)";
+ static char fmt_400[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
+ ". MAX:\002,e20.13)";
+ static char fmt_410[] = "(\002REL. STEP MAX :\002,e20.13)";
+ static char fmt_560[] = "(\002 ----------- ITERATION \002,i4,\002 ---"
+ "-------------\002)";
+ static char fmt_570[] = "(\002 X=\002,999e20.13)";
+ static char fmt_580[] = "(\002 F(X)=\002,e20.13)";
+ static char fmt_590[] = "(\002 G(X)=\002,999e20.13)";
+ static char fmt_600[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
+ ". MAX:\002,e20.13)";
+ static char fmt_610[] = "(\002REL. STEP MAX :\002,e20.13)";
+
+ /* System generated locals */
+ integer h_dim1, h_offset, i__1, i__2;
+ doublereal d__1, d__2;
+
+ /* Builtin functions */
+ double pow_di(doublereal *, integer *), sqrt(doublereal);
+ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
+
+ /* Local variables */
+ doublereal f;
+ integer i__;
+ doublereal gd, fn, fp, rnf, eps, rgx, rsx, beta;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ integer imsg;
+ doublereal temp, alpha;
+ extern /* Subroutine */ int mkmdl_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *);
+ integer nfcnt;
+ doublereal finit;
+ extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
+ doublereal *, integer *);
+ logical nomin;
+ doublereal gnorm, addmax;
+ extern /* Subroutine */ int grdchk_(integer *, doublereal *, S_fp,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *, integer *,
+ integer *, integer *), heschk_(integer *, integer *, doublereal *
+ , S_fp, S_fp, S_fp, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *, integer *,
+ doublereal *, doublereal *, doublereal *, integer *, integer *,
+ integer *);
+ integer iretcd;
+ doublereal analtl;
+ extern /* Subroutine */ int choldr_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, integer *, doublereal *, doublereal *,
+ doublereal *), mcheps_(doublereal *), sndofd_(integer *, integer
+ *, doublereal *, S_fp, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, integer *);
+ integer itrmcd;
+ extern /* Subroutine */ int bakslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *), fstofd_(integer *, integer *,
+ integer *, doublereal *, S_fp, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, integer *, integer *);
+ integer icscmx;
+ extern /* Subroutine */ int optchk_(integer *, integer *, integer *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, integer *, integer *, doublereal *, integer *,
+ integer *, integer *, doublereal *, integer *);
+ logical mxtake;
+ extern /* Subroutine */ int lnsrch_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, logical *, integer *, doublereal *, doublereal *,
+ doublereal *, S_fp, doublereal *, integer *), slvmdl_(integer *,
+ integer *, doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, integer *, doublereal *
+ , doublereal *, doublereal *, doublereal *, doublereal *, logical
+ *, doublereal *), forslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *);
+ extern doublereal twonrm_(integer *, doublereal *);
+ extern /* Subroutine */ int optstp_(integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *, integer *, integer *,
+ doublereal *, doublereal *, doublereal *, integer *, integer *,
+ logical *, integer *, integer *, doublereal *, doublereal *);
+
+ /* Fortran I/O blocks */
+ static cilist io___70 = { 0, 0, 0, fmt_25, 0 };
+ static cilist io___71 = { 0, 0, 0, fmt_30, 0 };
+ static cilist io___81 = { 0, 0, 0, fmt_103, 0 };
+ static cilist io___82 = { 0, 0, 0, fmt_104, 0 };
+ static cilist io___83 = { 0, 0, 0, fmt_105, 0 };
+ static cilist io___84 = { 0, 0, 0, fmt_106, 0 };
+ static cilist io___85 = { 0, 0, 0, fmt_108, 0 };
+ static cilist io___86 = { 0, 0, 0, fmt_110, 0 };
+ static cilist io___93 = { 0, 0, 0, fmt_370, 0 };
+ static cilist io___94 = { 0, 0, 0, fmt_380, 0 };
+ static cilist io___95 = { 0, 0, 0, fmt_390, 0 };
+ static cilist io___96 = { 0, 0, 0, fmt_400, 0 };
+ static cilist io___97 = { 0, 0, 0, fmt_410, 0 };
+ static cilist io___98 = { 0, 0, 0, fmt_560, 0 };
+ static cilist io___99 = { 0, 0, 0, fmt_570, 0 };
+ static cilist io___100 = { 0, 0, 0, fmt_580, 0 };
+ static cilist io___101 = { 0, 0, 0, fmt_590, 0 };
+ static cilist io___102 = { 0, 0, 0, fmt_600, 0 };
+ static cilist io___103 = { 0, 0, 0, fmt_610, 0 };
+
+
+
+/* PURPOSE: */
+/* DERIVATIVE TENSOR METHODS FOR UNCONSTRAINED OPTIMIZATION */
+
+/* ----------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> INITIAL GUESS (INPUT) AND CURRENT POINT */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
+/* FCN: R(N) --> R(1) */
+/* GRD --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT */
+/* FCN: R(N) --> R(N) */
+/* HSN --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN */
+/* FCN: R(N) --> R(N X N) */
+/* HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
+/* AND DIAGONAL OF H */
+/* TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X */
+/* FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
+/* GRADTL --> GRADIENT TOLERANCE */
+/* STEPTL --> STEP TOLERANCE */
+/* ITNLIM --> ITERATION LIMIT */
+/* STEPMX --> MAXIMUM STEP LENGTH ALLOWED */
+/* IPR --> OUTPUT UNIT NUMBER */
+/* METHOD --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT */
+/* EACH ITERATION */
+/* GRDFLG --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED */
+/* HESFLG --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED */
+/* NDIGIT --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
+/* MSG --> OUTPUT MESSAGE CONTRO */
+/* XPLS(N) <-- NEW POINT AND FINAL POINT (OUTPUT) */
+/* FPLS <-- NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
+/* GPLS(N) <-- CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) */
+/* H(N,N) --> HESSIAN */
+/* ITNNO <-- NUMBER OF ITERATIONS */
+/* G(N) --> PREVIOUS GRADIENT */
+/* S --> STEP TO PREVIOUS ITERATE (FOR TENSOR MODEL) */
+/* D --> CURRENT STEP (TENSOR OR NEWTON) */
+/* DN --> NEWTON STEP */
+/* E(N) --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION */
+/* WK1(N) --> WORKSPACE */
+/* WK2(N) --> WORKSPACE */
+/* PIVOT(N)--> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION */
+
+/* -------------------------------- */
+/* INITIALIZATION | */
+/* -------------------------------- */
+
+ /* Parameter adjustments */
+ --pivot;
+ --wk2;
+ --wk1;
+ --e;
+ --dn;
+ --d__;
+ --s;
+ --g;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+ --gpls;
+ --xpls;
+ --typx;
+ --x;
+
+ /* Function Body */
+ *itnno = 0;
+ icscmx = 0;
+ nfcnt = 0;
+ mcheps_(&eps);
+ optchk_(nr, n, msg, &x[1], &typx[1], fscale, gradtl, steptl, itnlim,
+ ndigit, &eps, method, grdflg, hesflg, stepmx, ipr);
+ if (*msg < 0) {
+ return 0;
+ }
+
+/* SCALE X */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] /= typx[i__];
+/* L10: */
+ }
+
+/* -------------------------------- */
+/* INITIAL ITERATION | */
+/* -------------------------------- */
+
+/* COMPUTE TYPICAL SIZE OF X */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dn[i__] = 1. / typx[i__];
+/* L15: */
+ }
+/* Computing MAX */
+ i__1 = -(*ndigit);
+ d__1 = pow_di(&c_b134, &i__1);
+ rnf = max(d__1,eps);
+/* Computing MAX */
+ d__1 = .01, d__2 = sqrt(rnf);
+ analtl = max(d__1,d__2);
+/* UNSCALE X AND COMPUTE F AND G */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ wk1[i__] = x[i__] * typx[i__];
+/* L20: */
+ }
+ (*fcn)(n, &wk1[1], &f);
+ ++nfcnt;
+ if (*grdflg >= 1) {
+ (*grd)(n, &wk1[1], &g[1]);
+ if (*grdflg == 1) {
+ *fscale = 1.;
+ grdchk_(n, &wk1[1], (S_fp)fcn, &f, &g[1], &dn[1], &typx[1],
+ fscale, &rnf, &analtl, &wk2[1], msg, ipr, &nfcnt);
+ }
+ } else {
+ fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, &f, &g[1], &typx[1], &
+ rnf, &wk2[1], &c__1, &nfcnt);
+ }
+ if (*msg < -20) {
+ return 0;
+ }
+
+ gnorm = twonrm_(n, &g[1]);
+
+/* PRINT OUT 1ST ITERATION? */
+ if (*msg >= 1) {
+ io___70.ciunit = *ipr;
+ s_wsfe(&io___70);
+ do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___71.ciunit = *ipr;
+ s_wsfe(&io___71);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ }
+
+/* TEST WHETHER THE INITIAL GUESS SATISFIES THE STOPPING CRITERIA */
+ if (gnorm <= *gradtl) {
+ *fpls = f;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ xpls[i__] = x[i__];
+ gpls[i__] = g[i__];
+/* L40: */
+ }
+ optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
+ gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &
+ rgx, &rsx);
+ goto L350;
+ }
+ finit = f;
+
+/* ------------------------ */
+/* ITERATION 1 | */
+/* ------------------------ */
+
+ ++(*itnno);
+
+/* COMPUTE HESSIAN */
+ if (*hesflg == 0) {
+ if (*grdflg == 1) {
+ fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
+ typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
+ } else {
+ sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
+ rnf, &wk2[1], &d__[1], &nfcnt);
+ }
+ } else {
+ if (*hesflg == 2) {
+ (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
+ } else {
+ if (*hesflg == 1) {
+/* IN HESCHK GPLS,XPLS AND E ARE USED AS WORK SPACE */
+ heschk_(nr, n, &wk1[1], (S_fp)fcn, (S_fp)grd, (S_fp)hsn, &f, &
+ g[1], &h__[h_offset], &dn[1], &typx[1], &rnf, &analtl,
+ grdflg, &gpls[1], &xpls[1], &e[1], msg, ipr, &nfcnt);
+ }
+ }
+ }
+ if (*msg < -20) {
+ return 0;
+ }
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ i__2 = i__ - 1;
+ dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
+/* L50: */
+ }
+
+/* CHOLESKY DECOMPOSITION FOR H (H=LLT) */
+ choldr_(nr, n, &h__[h_offset], &d__[1], &eps, &pivot[1], &e[1], &wk1[1], &
+ addmax);
+
+/* SOLVE FOR NEWTON STEP D */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* L60: */
+ wk2[i__] = -g[i__];
+ }
+ forslv_(nr, n, &h__[h_offset], &wk1[1], &wk2[1]);
+ bakslv_(nr, n, &h__[h_offset], &d__[1], &wk1[1]);
+
+/* APPLY LINESEARCH TO THE NEWTON STEP */
+ lnsrch_(nr, n, &x[1], &f, &g[1], &d__[1], &xpls[1], fpls, &mxtake, &
+ iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
+
+/* UPDATE G */
+/* CALL DCOPY(N,GPLS(1),1,GP(1),1) */
+
+/* UNSCALE XPLS AND COMPUTE GPLS */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ wk1[i__] = xpls[i__] * typx[i__];
+/* L80: */
+ }
+ if (*grdflg >= 1) {
+ (*grd)(n, &wk1[1], &gpls[1]);
+ } else {
+ fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
+ &rnf, &wk2[1], &c__1, &nfcnt);
+ }
+
+/* CHECK STOPPING CONDITIONS */
+ optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
+ gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &rgx,
+ &rsx);
+
+/* IF ITRMCD > 0 THEN STOPPING CONDITIONS SATISFIED */
+ if (itrmcd > 0) {
+ goto L350;
+ }
+
+/* UPDATE X,F AND S FOR TENSOR MODEL */
+ fp = f;
+ f = *fpls;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = xpls[i__];
+ s[i__] = x[i__] - temp;
+ x[i__] = temp;
+/* L90: */
+ }
+
+/* IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
+ if (*msg >= 2) {
+ io___81.ciunit = *ipr;
+ s_wsfe(&io___81);
+ do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
+ e_wsfe();
+ io___82.ciunit = *ipr;
+ s_wsfe(&io___82);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ io___83.ciunit = *ipr;
+ s_wsfe(&io___83);
+ do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___84.ciunit = *ipr;
+ s_wsfe(&io___84);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ }
+ if (*msg >= 3) {
+ io___85.ciunit = *ipr;
+ s_wsfe(&io___85);
+ do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___86.ciunit = *ipr;
+ s_wsfe(&io___86);
+ do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ }
+
+/* ------------------------ */
+/* ITERATION > 1 | */
+/* ------------------------ */
+
+/* UNSCALE X AND COMPUTE H */
+L200:
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ wk1[i__] = x[i__] * typx[i__];
+/* L210: */
+ }
+ if (*hesflg == 0) {
+ if (*grdflg == 1) {
+ fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
+ typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
+ } else {
+ sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
+ rnf, &wk2[1], &d__[1], &nfcnt);
+ }
+ } else {
+ (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
+ }
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ i__2 = i__ - 1;
+ dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
+/* L230: */
+ }
+
+/* IF METHOD = 0 THEN USE NEWTON STEP ONLY */
+ if (*method == 0) {
+
+/* CHOLESKY DECOMPOSITION FOR H */
+ choldr_(nr, n, &h__[h_offset], &wk2[1], &eps, &pivot[1], &e[1], &wk1[
+ 1], &addmax);
+
+/* COMPUTE NEWTON STEP */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ wk1[i__] = -gpls[i__];
+/* L240: */
+ }
+ forslv_(nr, n, &h__[h_offset], &wk2[1], &wk1[1]);
+ bakslv_(nr, n, &h__[h_offset], &d__[1], &wk2[1]);
+
+/* NO TENSOR STEP */
+ nomin = TRUE_;
+ goto L300;
+
+ }
+
+/* FORM TENSOR MODEL */
+ mkmdl_(nr, n, &f, &fp, &gpls[1], &g[1], &s[1], &h__[h_offset], &alpha, &
+ beta, &wk1[1], &d__[1]);
+
+/* SOLVE TENSOR MODEL AND COMPUTE NEWTON STEP */
+/* ON INPUT : SH IS STORED IN WK1 */
+/* A=(G-GPLS-SH-S*BETA/(6*STS)) IS STORED IN D */
+/* ON OUTPUT: NEWTON STEP IS STORED IN DN */
+/* TENSOR STEP IS STORED IN D */
+ slvmdl_(nr, n, &h__[h_offset], &xpls[1], &wk2[1], &e[1], &g[1], &s[1], &
+ gpls[1], &pivot[1], &d__[1], &wk1[1], &dn[1], &alpha, &beta, &
+ nomin, &eps);
+
+/* IF TENSOR MODEL HAS NO MINIMIZER THEN USE NEWTON STEP */
+ if (nomin) {
+ dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
+ goto L300;
+ }
+
+/* IF TENSOR STEP IS NOT IN DESCENT DIRECTION THEN USE NEWTON STEP */
+ gd = ddot_(n, &gpls[1], &c__1, &d__[1], &c__1);
+ if (gd > 0.) {
+ dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
+ nomin = TRUE_;
+ }
+
+L300:
+ ++(*itnno);
+ dcopy_(n, &gpls[1], &c__1, &g[1], &c__1);
+
+/* APPLY LINESEARCH TO TENSOR (OR NEWTON) STEP */
+ lnsrch_(nr, n, &x[1], &f, &g[1], &d__[1], &xpls[1], fpls, &mxtake, &
+ iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
+
+ if (! nomin) {
+/* TENSOR STEP IS FOUND AND IN DESCENT DIRECTION, */
+/* APPLY LINESEARCH TO NEWTON STEP */
+/* NEW NEWTON POINT IN WK2 */
+ lnsrch_(nr, n, &x[1], &f, &g[1], &dn[1], &wk2[1], &fn, &mxtake, &
+ iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
+
+/* COMPARE TENSOR STEP TO NEWTON STEP */
+/* IF NEWTON STEP IS BETTER, SET NEXT ITERATE TO NEW NEWTON POINT */
+ if (fn < *fpls) {
+ *fpls = fn;
+ dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
+ dcopy_(n, &wk2[1], &c__1, &xpls[1], &c__1);
+ }
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ d__[i__] = xpls[i__] - x[i__];
+/* L320: */
+ }
+
+/* UNSCALE XPLS, AND COMPUTE FPLS AND GPLS */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ wk1[i__] = xpls[i__] * typx[i__];
+/* L330: */
+ }
+ (*fcn)(n, &wk1[1], fpls);
+ ++nfcnt;
+ if (*grdflg >= 1) {
+ (*grd)(n, &wk1[1], &gpls[1]);
+ } else {
+ fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
+ &rnf, &wk2[1], &c__1, &nfcnt);
+ }
+
+/* CHECK STOPPING CONDITIONS */
+ imsg = *msg;
+ optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
+ gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &rgx,
+ &rsx);
+
+/* IF ITRMCD = 0 THEN NOT OVER YET */
+ if (itrmcd == 0) {
+ goto L500;
+ }
+
+/* IF MSG >= 1 THEN PRINT OUT FINAL ITERATION */
+L350:
+ if (imsg >= 1) {
+
+/* TRANSFORM X BACK TO ORIGINAL SPACE */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ xpls[i__] *= typx[i__];
+/* L360: */
+ }
+ io___93.ciunit = *ipr;
+ s_wsfe(&io___93);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&xpls[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ io___94.ciunit = *ipr;
+ s_wsfe(&io___94);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ io___95.ciunit = *ipr;
+ s_wsfe(&io___95);
+ do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
+ do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
+ e_wsfe();
+ if (imsg >= 3) {
+ io___96.ciunit = *ipr;
+ s_wsfe(&io___96);
+ do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___97.ciunit = *ipr;
+ s_wsfe(&io___97);
+ do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ }
+/* UPDATE THE HESSIAN */
+ if (*hesflg == 0) {
+ if (*grdflg == 1) {
+ fstofd_(nr, n, n, &xpls[1], (S_fp)grd, &gpls[1], &h__[
+ h_offset], &typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
+ } else {
+ sndofd_(nr, n, &xpls[1], (S_fp)fcn, fpls, &h__[h_offset], &
+ typx[1], &rnf, &wk2[1], &d__[1], &nfcnt);
+ }
+ } else {
+ (*hsn)(nr, n, &xpls[1], &h__[h_offset]);
+ }
+ }
+ return 0;
+
+/* UPDATE INFORMATION AT THE CURRENT POINT */
+L500:
+ dcopy_(n, &xpls[1], &c__1, &x[1], &c__1);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ s[i__] = -d__[i__];
+/* L550: */
+ }
+
+/* IF TOO MANY ITERATIONS THEN RETURN */
+ if (*itnno > *itnlim) {
+ goto L350;
+ }
+
+/* IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
+ if (*msg >= 2) {
+ io___98.ciunit = *ipr;
+ s_wsfe(&io___98);
+ do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
+ e_wsfe();
+ io___99.ciunit = *ipr;
+ s_wsfe(&io___99);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ io___100.ciunit = *ipr;
+ s_wsfe(&io___100);
+ do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___101.ciunit = *ipr;
+ s_wsfe(&io___101);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ }
+ if (*msg >= 3) {
+ io___102.ciunit = *ipr;
+ s_wsfe(&io___102);
+ do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ io___103.ciunit = *ipr;
+ s_wsfe(&io___103);
+ do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
+ e_wsfe();
+ }
+/* UPDATE F */
+ fp = f;
+ f = *fpls;
+
+/* PERFORM NEXT ITERATION */
+ goto L200;
+
+/* END OF ITERATION > 1 */
+
+} /* opt_ */
+
+/* ---------------------- */
+/* | D F A U L T | */
+/* ---------------------- */
+/* Subroutine */ int dfault_(integer *n, doublereal *typx, doublereal *fscale,
+ doublereal *gradtl, doublereal *steptl, integer *itnlim, doublereal *
+ stepmx, integer *ipr, integer *method, integer *grdflg, integer *
+ hesflg, integer *ndigit, integer *msg)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Builtin functions */
+ double pow_dd(doublereal *, doublereal *), d_lg10(doublereal *);
+
+ /* Local variables */
+ integer i__;
+ doublereal eps, temp;
+ extern /* Subroutine */ int mcheps_(doublereal *);
+
+ /* Parameter adjustments */
+ --typx;
+
+ /* Function Body */
+ mcheps_(&eps);
+ *method = 1;
+ *fscale = 1.;
+ *grdflg = 0;
+ *hesflg = 0;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ typx[i__] = 1.;
+/* L1: */
+ }
+ temp = pow_dd(&eps, &c_b250);
+ *gradtl = temp;
+ *steptl = temp * temp;
+ *ndigit = (integer) (-d_lg10(&eps));
+/* SET ACTUAL DFAULT VALUE OF STEPMX IN OPTCHK */
+ *stepmx = 0.;
+ *itnlim = 100;
+ *ipr = 6;
+ *msg = 1;
+ return 0;
+} /* dfault_ */
+
+/* ---------------------- */
+/* | O P T C H K | */
+/* ---------------------- */
+/* Subroutine */ int optchk_(integer *nr, integer *n, integer *msg,
+ doublereal *x, doublereal *typx, doublereal *fscale, doublereal *
+ gradtl, doublereal *steptl, integer *itnlim, integer *ndigit,
+ doublereal *eps, integer *method, integer *grdflg, integer *hesflg,
+ doublereal *stepmx, integer *ipr)
+{
+ /* Format strings */
+ static char fmt_901[] = "(\0020OPTCHK ILLEGAL DIMENSION, N=\002,i5)";
+ static char fmt_902[] = "(\002OPTCHK ILLEGAL INPUT VALUES OF: NR=\002"
+ ",i5,\002, N=\002,i5,\002, NR MUST BE <= N.\002)";
+
+ /* System generated locals */
+ integer i__1;
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal), pow_dd(doublereal *, doublereal *), d_lg10(
+ doublereal *), pow_di(doublereal *, integer *);
+ integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
+
+ /* Local variables */
+ integer i__;
+ doublereal temp, stpsiz;
+
+ /* Fortran I/O blocks */
+ static cilist io___110 = { 0, 0, 0, fmt_901, 0 };
+ static cilist io___111 = { 0, 0, 0, fmt_902, 0 };
+
+
+
+/* PURPOSE */
+/* ------- */
+/* CHECK INPUT FOR REASONABLENESS */
+
+/* PARAMETERS */
+/* ---------- */
+/* NR <--> ROW DIMENSION OF H AND WRK */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> ON ENTRY, ESTIMATE TO ROOT OF FCN */
+/* TYPX(N) <--> TYPICAL SIZE OF EACH COMPONENT OF X */
+/* FSCALE <--> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
+/* GRADTL <--> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE */
+/* ENOUGH TO ZERO TO TERMINATE ALGORITHM */
+/* STEPTL <--> TOLERANCE AT WHICH STEP LENGTH CONSIDERED CLOSE */
+/* ENOUGH TO ZERO TO TERMINATE ALGORITHM */
+/* ITNLIM <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS */
+/* NDIGIT <--> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
+/* EPS --> MACHINE PRECISION */
+/* METHOD <--> ALGORITHM INDICATOR */
+/* GRDFLG <--> =1 IF ANALYTIC GRADIENT SUPPLIED */
+/* HESFLG <--> =1 IF ANALYTIC HESSIAN SUPPLIED */
+/* STEPMX <--> MAXIMUM STEP SIZE */
+/* MSG <--> MESSAGE AND ERROR CODE */
+/* IPR --> DEVICE TO WHICH TO SEND OUTPUT */
+
+
+/* CHECK DIMENSION OF PROBLEM */
+
+ /* Parameter adjustments */
+ --typx;
+ --x;
+
+ /* Function Body */
+ if (*n <= 0) {
+ goto L805;
+ }
+ if (*nr < *n) {
+ goto L806;
+ }
+
+/* CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES. */
+/* IF NOT, SET THEM TO DEFAULT VALUES. */
+ if (*method != 0) {
+ *method = 1;
+ }
+ if (*grdflg != 2 && *grdflg != 1) {
+ *grdflg = 0;
+ }
+ if (*hesflg != 2 && *hesflg != 1) {
+ *hesflg = 0;
+ }
+ if (*msg > 3 || (real) (*msg) < 0.f) {
+ *msg = 1;
+ }
+
+/* COMPUTE SCALE MATRIX */
+
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (typx[i__] == 0.f) {
+ typx[i__] = 1.;
+ }
+ if (typx[i__] < 0.f) {
+ typx[i__] = -typx[i__];
+ }
+/* L10: */
+ }
+
+/* CHECK MAXIMUM STEP SIZE */
+
+ if (*stepmx > 0.) {
+ goto L20;
+ }
+ stpsiz = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ stpsiz += x[i__] * x[i__] / typx[i__] * typx[i__];
+/* L15: */
+ }
+ stpsiz = sqrt(stpsiz);
+/* Computing MAX */
+ d__1 = stpsiz * 1e3;
+ *stepmx = max(d__1,1e3);
+L20:
+/* CHECK FUNCTION SCALE */
+ if (*fscale == 0.f) {
+ *fscale = 1.;
+ }
+ if (*fscale < 0.f) {
+ *fscale = -(*fscale);
+ }
+/* CHECK GRADIENT TOLERANCE */
+ if (*gradtl < 0.f) {
+ *gradtl = pow_dd(eps, &c_b250);
+ }
+/* CHECK STEP TOLERANCE */
+ if (*steptl < 0.f) {
+ temp = pow_dd(eps, &c_b250);
+ *steptl = temp * temp;
+ }
+
+/* CHECK ITERATION LIMIT */
+ if (*itnlim <= 0) {
+ *itnlim = 100;
+ }
+
+/* CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN */
+ if (*ndigit <= 0) {
+ *ndigit = (integer) (-d_lg10(eps));
+ }
+ i__1 = -(*ndigit);
+ if (pow_di(&c_b134, &i__1) <= *eps) {
+ *ndigit = (integer) (-d_lg10(eps));
+ }
+ return 0;
+
+/* ERROR EXITS */
+
+L805:
+ io___110.ciunit = *ipr;
+ s_wsfe(&io___110);
+ do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
+ e_wsfe();
+ *msg = -20;
+ return 0;
+L806:
+ io___111.ciunit = *ipr;
+ s_wsfe(&io___111);
+ do_fio(&c__1, (char *)&(*nr), (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
+ e_wsfe();
+ *msg = -20;
+ return 0;
+} /* optchk_ */
+
+/* ---------------------- */
+/* | G R D C H K | */
+/* ---------------------- */
+/* Subroutine */ int grdchk_(integer *n, doublereal *x, S_fp fcn, doublereal *
+ f, doublereal *g, doublereal *typsiz, doublereal *typx, doublereal *
+ fscale, doublereal *rnf, doublereal *analtl, doublereal *wrk1,
+ integer *msg, integer *ipr, integer *ifn)
+{
+ /* Format strings */
+ static char fmt_901[] = "(\0020GRDCHK PROBABLE ERROR IN CODING OF ANA"
+ "LYTIC\002,\002 GRADIENT FUNCTION.\002/\002 GRDCHK COMP\002,1"
+ "2x,\002ANALYTIC\002,12x,\002ESTIMATE\002)";
+ static char fmt_902[] = "(\002 GRDCHK \002,i5,3x,e20.13,3x,e20.13)";
+
+ /* System generated locals */
+ integer i__1;
+ doublereal d__1, d__2, d__3, d__4;
+
+ /* Builtin functions */
+ integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
+
+ /* Local variables */
+ integer i__;
+ doublereal gs;
+ integer ker;
+ doublereal wrk;
+ extern /* Subroutine */ int fstofd_(integer *, integer *, integer *,
+ doublereal *, S_fp, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *, integer *);
+
+ /* Fortran I/O blocks */
+ static cilist io___116 = { 0, 0, 0, fmt_901, 0 };
+ static cilist io___117 = { 0, 0, 0, fmt_902, 0 };
+
+
+
+/* PURPOSE */
+/* ------- */
+/* CHECK ANALYTIC GRADIENT AGAINST ESTIMATED GRADIENT */
+
+/* PARAMETERS */
+/* ---------- */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> ESTIMATE TO A ROOT OF FCN */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION */
+/* MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
+/* FCN: R(N) --> R(1) */
+/* F --> FUNCTION VALUE: FCN(X) */
+/* G(N) --> GRADIENT: G(X) */
+/* TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X */
+/* FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
+/* RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN */
+/* ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND */
+/* ANALYTICAL GRADIENTS */
+/* WRK1(N) --> WORKSPACE */
+/* MSG <-- MESSAGE OR ERROR CODE */
+/* ON OUTPUT: =-21, PROBABLE CODING ERROR OF GRADIENT */
+/* IPR --> DEVICE TO WHICH TO SEND OUTPUT */
+/* IFN <--> NUMBER OF FUNCTION EVALUATIONS */
+
+/* COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO */
+/* ANALYTIC GRADIENT. */
+
+ /* Parameter adjustments */
+ --wrk1;
+ --typx;
+ --typsiz;
+ --g;
+ --x;
+
+ /* Function Body */
+ fstofd_(&c__1, &c__1, n, &x[1], (S_fp)fcn, f, &wrk1[1], &typx[1], rnf, &
+ wrk, &c__1, ifn);
+ ker = 0;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__2 = abs(*f);
+/* Computing MAX */
+ d__3 = (d__1 = x[i__], abs(d__1)), d__4 = typsiz[i__];
+ gs = max(d__2,*fscale) / max(d__3,d__4);
+/* Computing MAX */
+ d__3 = (d__1 = g[i__], abs(d__1));
+ if ((d__2 = g[i__] - wrk1[i__], abs(d__2)) > max(d__3,gs) * *analtl) {
+ ker = 1;
+ }
+/* L5: */
+ }
+ if (ker == 0) {
+ goto L20;
+ }
+ io___116.ciunit = *ipr;
+ s_wsfe(&io___116);
+ e_wsfe();
+ io___117.ciunit = *ipr;
+ s_wsfe(&io___117);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
+ do_fio(&c__1, (char *)&wrk1[i__], (ftnlen)sizeof(doublereal));
+ }
+ e_wsfe();
+ *msg = -21;
+L20:
+ return 0;
+} /* grdchk_ */
+
+/* ---------------------- */
+/* | H E S C H K | */
+/* ---------------------- */
+/* Subroutine */ int heschk_(integer *nr, integer *n, doublereal *x, S_fp fcn,
+ S_fp grd, S_fp hsn, doublereal *f, doublereal *g, doublereal *a,
+ doublereal *typsiz, doublereal *typx, doublereal *rnf, doublereal *
+ analtl, integer *iagflg, doublereal *udiag, doublereal *wrk1,
+ doublereal *wrk2, integer *msg, integer *ipr, integer *ifn)
+{
+ /* Format strings */
+ static char fmt_901[] = "(\002 HESCHK PROBABLE ERROR IN CODING OF ANA"
+ "LYTIC\002,\002 HESSIAN FUNCTION.\002/\002 HESCHK ROW CO"
+ "L\002,14x,\002ANALYTIC\002,14x,\002(ESTIMATE)\002)";
+ static char fmt_902[] = "(\002 HESCHK \002,2i5,2x,e20.13,2x,\002(\002"
+ ",e20.13,\002)\002)";
+
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+ doublereal d__1, d__2, d__3, d__4, d__5;
+
+ /* Builtin functions */
+ integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
+
+ /* Local variables */
+ integer i__, j;
+ doublereal hs;
+ integer im1, jp1, ker;
+ extern /* Subroutine */ int sndofd_(integer *, integer *, doublereal *,
+ S_fp, doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, integer *), fstofd_(integer *,
+ integer *, integer *, doublereal *, S_fp, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *, integer *,
+ integer *);
+
+ /* Fortran I/O blocks */
+ static cilist io___123 = { 0, 0, 0, fmt_901, 0 };
+ static cilist io___125 = { 0, 0, 0, fmt_902, 0 };
+ static cilist io___126 = { 0, 0, 0, fmt_902, 0 };
+
+
+
+/* PURPOSE */
+/* ------- */
+/* CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN */
+/* (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN */
+/* HSN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF A) */
+
+/* PARAMETERS */
+/* ---------- */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> ESTIMATE TO A ROOT OF FCN */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION */
+/* MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
+/* FCN: R(N) --> R(1) */
+/* GRD --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN. */
+/* MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
+/* HSN --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN. */
+/* MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
+/* HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
+/* AND DIAGONAL OF A */
+/* F --> FUNCTION VALUE: FCN(X) */
+/* G(N) <-- GRADIENT: G(X) */
+/* A(N,N) <-- ON EXIT: HESSIAN IN LOWER TRIANGULAR PART AND DIAG */
+/* TYPSIZ(N) --> TYPICAL SIZE FOR EACH COMPONENT OF X */
+/* RNF --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN */
+/* ANALTL --> TOLERANCE FOR COMPARISON OF ESTIMATED AND */
+/* ANALYTICAL GRADIENTS */
+/* IAGFLG --> =1 IF ANALYTIC GRADIENT SUPPLIED */
+/* UDIAG(N) --> WORKSPACE */
+/* WRK1(N) --> WORKSPACE */
+/* WRK2(N) --> WORKSPACE */
+/* MSG <--> MESSAGE OR ERROR CODE */
+/* ON INPUT : IF =1XX DO NOT COMPARE ANAL + EST HESS */
+/* ON OUTPUT: =-22, PROBABLE CODING ERROR OF HESSIAN */
+/* IPR --> DEVICE TO WHICH TO SEND OUTPUT */
+/* IFN <--> NUMBER OF FUNCTION EVALUTATIONS */
+
+/* COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN. */
+
+ /* Parameter adjustments */
+ a_dim1 = *nr;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --wrk2;
+ --wrk1;
+ --udiag;
+ --typx;
+ --typsiz;
+ --g;
+ --x;
+
+ /* Function Body */
+ if (*iagflg == 1) {
+ fstofd_(nr, n, n, &x[1], (S_fp)grd, &g[1], &a[a_offset], &typx[1],
+ rnf, &wrk1[1], &c__3, ifn);
+ }
+ if (*iagflg != 1) {
+ sndofd_(nr, n, &x[1], (S_fp)fcn, f, &a[a_offset], &typx[1], rnf, &
+ wrk1[1], &wrk2[1], ifn);
+ }
+
+ ker = 0;
+
+/* COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART */
+/* AND DIAGONAL OF "A" TO UDIAG */
+
+ i__1 = *n;
+ for (j = 1; j <= i__1; ++j) {
+ udiag[j] = a[j + j * a_dim1];
+ if (j == *n) {
+ goto L30;
+ }
+ jp1 = j + 1;
+ i__2 = *n;
+ for (i__ = jp1; i__ <= i__2; ++i__) {
+ a[j + i__ * a_dim1] = a[i__ + j * a_dim1];
+/* L25: */
+ }
+L30:
+ ;
+ }
+
+/* COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE */
+/* APPROXIMATION. */
+
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ i__2 = *n;
+ for (j = i__; j <= i__2; ++j) {
+ a[j + i__ * a_dim1] = 0.;
+/* L32: */
+ }
+ }
+
+ (*hsn)(nr, n, &x[1], &a[a_offset]);
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+/* Computing MAX */
+ d__3 = (d__1 = g[j], abs(d__1));
+/* Computing MAX */
+ d__4 = (d__2 = x[j], abs(d__2)), d__5 = typsiz[j];
+ hs = max(d__3,1.) / max(d__4,d__5);
+/* Computing MAX */
+ d__3 = (d__1 = udiag[j], abs(d__1));
+ if ((d__2 = a[j + j * a_dim1] - udiag[j], abs(d__2)) > max(d__3,hs) *
+ *analtl) {
+ ker = 1;
+ }
+ if (j == *n) {
+ goto L40;
+ }
+ jp1 = j + 1;
+ i__1 = *n;
+ for (i__ = jp1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__3 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
+ if ((d__2 = a[i__ + j * a_dim1] - a[j + i__ * a_dim1], abs(d__2))
+ > max(d__3,hs) * *analtl) {
+ ker = 1;
+ }
+/* L35: */
+ }
+L40:
+ ;
+ }
+
+ if (ker == 0) {
+ goto L90;
+ }
+ io___123.ciunit = *ipr;
+ s_wsfe(&io___123);
+ e_wsfe();
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ if (i__ == 1) {
+ goto L45;
+ }
+ im1 = i__ - 1;
+ i__1 = im1;
+ for (j = 1; j <= i__1; ++j) {
+ io___125.ciunit = *ipr;
+ s_wsfe(&io___125);
+ do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&a[i__ + j * a_dim1], (ftnlen)sizeof(
+ doublereal));
+ do_fio(&c__1, (char *)&a[j + i__ * a_dim1], (ftnlen)sizeof(
+ doublereal));
+ e_wsfe();
+/* L43: */
+ }
+L45:
+ io___126.ciunit = *ipr;
+ s_wsfe(&io___126);
+ do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
+ do_fio(&c__1, (char *)&a[i__ + i__ * a_dim1], (ftnlen)sizeof(
+ doublereal));
+ do_fio(&c__1, (char *)&udiag[i__], (ftnlen)sizeof(doublereal));
+ e_wsfe();
+/* L50: */
+ }
+ *msg = -22;
+/* ENDIF */
+L90:
+ return 0;
+} /* heschk_ */
+
+/* ----------------------- */
+/* | M C H E P S | */
+/* ----------------------- */
+/* Subroutine */ int mcheps_(doublereal *eps)
+{
+ doublereal temp, temp1;
+
+
+/* PURPOSE: */
+/* COMPUTE MACHINE PRECISION */
+
+/* ------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* EPS <-- MACHINE PRECISION */
+
+ temp = 1.;
+L20:
+ temp /= 2.;
+ temp1 = temp + 1.;
+ if (temp1 != 1.) {
+ goto L20;
+ }
+ *eps = temp * 2.;
+ return 0;
+} /* mcheps_ */
+
+
+doublereal twonrm_(integer *n, doublereal *v)
+{
+ /* System generated locals */
+ doublereal ret_val;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ doublereal temp;
+
+
+/* PURPOSE: */
+/* COMPUTER L-2 NORM */
+
+/* -------------------------------------------------------------------------- */
+
+/* PARAMETER: */
+
+/* N --> DIMENSION OF PROBLEM */
+/* V(N) --> VECTOR WHICH L-2 NORM IS EVALUATED */
+
+ /* Parameter adjustments */
+ --v;
+
+ /* Function Body */
+ temp = ddot_(n, &v[1], &c__1, &v[1], &c__1);
+ ret_val = sqrt(temp);
+ return ret_val;
+} /* twonrm_ */
+
+/* ---------------------- */
+/* | L N S R C H | */
+/* ---------------------- */
+/* Subroutine */ int lnsrch_(integer *nr, integer *n, doublereal *x,
+ doublereal *f, doublereal *g, doublereal *p, doublereal *xpls,
+ doublereal *fpls, logical *mxtake, integer *iretcd, doublereal *
+ stepmx, doublereal *steptl, doublereal *typx, S_fp fcn, doublereal *
+ w2, integer *nfcnt)
+{
+ /* System generated locals */
+ integer i__1;
+ doublereal d__1, d__2;
+
+ /* Builtin functions */
+ double sqrt(doublereal), d_sign(doublereal *, doublereal *);
+
+ /* Local variables */
+ doublereal a, b;
+ integer i__, k;
+ doublereal t1, t2, t3, scl, rln, sln, slp, tmp, disc;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ doublereal temp, temp1, temp2, alpha;
+ extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
+ integer *);
+ doublereal pfpls, almbda, plmbda, tlmbda, almbmn;
+
+
+/* THE ALPHA CONDITION ONLY LINE SEARCH */
+
+/* PURPOSE */
+/* ------- */
+/* FIND A NEXT NEWTON ITERATE BY LINE SEARCH. */
+
+/* PARAMETERS */
+/* ---------- */
+/* N --> DIMENSION OF PROBLEM */
+/* X(N) --> OLD ITERATE: X[K-1] */
+/* F --> FUNCTION VALUE AT OLD ITERATE, F(X) */
+/* G(N) --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE */
+/* P(N) --> NON-ZERO NEWTON STEP */
+/* XPLS(N) <-- NEW ITERATE X[K] */
+/* FPLS <-- FUNCTION VALUE AT NEW ITERATE, F(XPLS) */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
+/* IRETCD <-- RETURN CODE */
+/* MXTAKE <-- BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED */
+/* STEPMX --> MAXIMUM ALLOWABLE STEP SIZE */
+/* STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES */
+/* CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM */
+/* TYPX(N) --> DIAGONAL SCALING MATRIX FOR X (NOT IN UNCMIN) */
+/* IPR --> DEVICE TO WHICH TO SEND OUTPUT */
+/* W2 --> WORKING SPACE */
+/* NFCNT <--> NUMBER OF FUNCTION EVALUTIONS */
+
+/* INTERNAL VARIABLES */
+/* ------------------ */
+/* SLN NEWTON LENGTH */
+/* RLN RELATIVE LENGTH OF NEWTON STEP */
+
+
+ /* Parameter adjustments */
+ --w2;
+ --typx;
+ --xpls;
+ --p;
+ --g;
+ --x;
+
+ /* Function Body */
+ *mxtake = FALSE_;
+ *iretcd = 2;
+ alpha = 1e-4;
+/* $ WRITE(IPR,954) */
+/* $ WRITE(IPR,955) (P(I),I=1,N) */
+ tmp = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ tmp += p[i__] * p[i__];
+/* L5: */
+ }
+ sln = sqrt(tmp);
+ if (sln <= *stepmx) {
+ goto L10;
+ }
+
+/* NEWTON STEP LONGER THAN MAXIMUM ALLOWED */
+ scl = *stepmx / sln;
+ dscal_(n, &scl, &p[1], &c__1);
+ sln = *stepmx;
+/* $ WRITE(IPR,954) */
+/* $ WRITE(IPR,955) (P(I),I=1,N) */
+L10:
+ slp = ddot_(n, &g[1], &c__1, &p[1], &c__1);
+ rln = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = 1.;
+ temp1 = (d__1 = x[i__], abs(d__1));
+ temp2 = max(temp1,temp);
+ temp1 = (d__1 = p[i__], abs(d__1));
+/* Computing MAX */
+ d__1 = rln, d__2 = temp1 / temp2;
+ rln = max(d__1,d__2);
+/* L15: */
+ }
+ almbmn = *steptl / rln;
+ almbda = 1.;
+/* $ WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL */
+
+/* LOOP */
+/* CHECK IF NEW ITERATE SATISFACTORY. GENERATE NEW LAMBDA IF NECESSARY. */
+
+L100:
+ if (*iretcd < 2) {
+ return 0;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ xpls[i__] = x[i__] + almbda * p[i__];
+/* L105: */
+ }
+ i__1 = *n;
+ for (k = 1; k <= i__1; ++k) {
+ w2[k] = xpls[k] * typx[k];
+/* L101: */
+ }
+ (*fcn)(n, &w2[1], fpls);
+ ++(*nfcnt);
+/* $ WRITE(IPR,956) ALMBDA */
+/* $ WRITE(IPR,951) */
+/* $ WRITE(IPR,955) (XPLS(I),I=1,N) */
+/* $ WRITE(IPR,953) FPLS */
+ if (*fpls > *f + slp * alpha * almbda) {
+ goto L130;
+ }
+/* IF(FPLS.LE. F+SLP*1.E-4*ALMBDA) */
+/* THEN */
+
+/* SOLUTION FOUND */
+
+ *iretcd = 0;
+ if (almbda == 1. && sln > *stepmx * .99) {
+ *mxtake = TRUE_;
+ }
+ goto L100;
+
+/* SOLUTION NOT (YET) FOUND */
+
+/* ELSE */
+L130:
+ if (almbda >= almbmn) {
+ goto L140;
+ }
+/* IF(ALMBDA .LT. ALMBMN) */
+/* THEN */
+
+/* NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X */
+
+ *iretcd = 1;
+ goto L100;
+/* ELSE */
+
+/* CALCULATE NEW LAMBDA */
+
+L140:
+ if (almbda != 1.) {
+ goto L150;
+ }
+/* IF(ALMBDA.EQ.1.0) */
+/* THEN */
+
+/* FIRST BACKTRACK: QUADRATIC FIT */
+
+ tlmbda = -slp / ((*fpls - *f - slp) * 2.);
+ goto L170;
+/* ELSE */
+
+/* ALL SUBSEQUENT BACKTRACKS: CUBIC FIT */
+
+L150:
+ t1 = *fpls - *f - almbda * slp;
+ t2 = pfpls - *f - plmbda * slp;
+ t3 = 1. / (almbda - plmbda);
+ a = t3 * (t1 / (almbda * almbda) - t2 / (plmbda * plmbda));
+ b = t3 * (t2 * almbda / (plmbda * plmbda) - t1 * plmbda / (almbda *
+ almbda));
+ disc = b * b - a * 3.f * slp;
+ if (disc <= b * b) {
+ goto L160;
+ }
+/* IF(DISC.GT. B*B) */
+/* THEN */
+
+/* ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM */
+
+ tlmbda = (-b + d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
+ goto L165;
+/* ELSE */
+
+/* BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM */
+
+L160:
+ tlmbda = (-b - d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
+/* ENDIF */
+L165:
+ if (tlmbda > almbda * .5) {
+ tlmbda = almbda * .5;
+ }
+/* ENDIF */
+L170:
+ plmbda = almbda;
+ pfpls = *fpls;
+ if (tlmbda >= almbda * .1) {
+ goto L180;
+ }
+/* IF(TLMBDA.LT.ALMBDA/10.) */
+/* THEN */
+ almbda *= .1f;
+ goto L190;
+/* ELSE */
+L180:
+ almbda = tlmbda;
+/* ENDIF */
+/* ENDIF */
+/* ENDIF */
+L190:
+ goto L100;
+/* L956: */
+/* L951: */
+/* L952: */
+/* L953: */
+/* L954: */
+/* L955: */
+} /* lnsrch_ */
+
+/* ---------------------- */
+/* | Z H Z | */
+/* ---------------------- */
+/* Subroutine */ int zhz_(integer *nr, integer *n, doublereal *y, doublereal *
+ h__, doublereal *u, doublereal *t)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, i__1, i__2;
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ doublereal d__;
+ integer i__, j;
+ doublereal s, sgn;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ doublereal temp1, temp2;
+ extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
+ doublereal *, integer *);
+ doublereal ynorm;
+
+
+/* PURPOSE: */
+/* COMPUTE QTHQ(N,N) AND ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND */
+/* FIRST N-1 COLUMNS OF QTHQ */
+
+/* --------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* Y(N) --> FIRST BASIS IN Q */
+/* H(N,N) <--> ON INPUT : HESSIAN */
+/* ON OUTPUT: QTHQ (ZTHZ) */
+/* U(N) <-- VECTOR TO FORM Q AND Z */
+/* T(N) --> WORKSPACE */
+
+
+/* U=Y+SGN(Y(N))||Y||E(N) */
+ /* Parameter adjustments */
+ --t;
+ --u;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+ --y;
+
+ /* Function Body */
+ if (y[*n] != 0.) {
+ sgn = y[*n] / (d__1 = y[*n], abs(d__1));
+ } else {
+ sgn = 1.;
+ }
+ ynorm = ddot_(n, &y[1], &c__1, &y[1], &c__1);
+ ynorm = sqrt(ynorm);
+ u[*n] = y[*n] + sgn * ynorm;
+ i__1 = *n - 1;
+ dcopy_(&i__1, &y[1], &c__1, &u[1], &c__1);
+
+/* D=UTU/2 */
+ d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
+ d__ /= 2.;
+
+/* T=2HU/UTU */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ t[i__] = 0.;
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ t[i__] = h__[i__ + j * h_dim1] * u[j] + t[i__];
+/* L30: */
+ }
+ t[i__] /= d__;
+/* L40: */
+ }
+
+/* S=4UHU/(UTU)**2 */
+ s = ddot_(n, &u[1], &c__1, &t[1], &c__1);
+ s /= d__;
+
+/* COMPUTE QTHQ (ZTHZ) */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp1 = u[i__];
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ temp2 = u[j];
+ h__[i__ + j * h_dim1] = h__[i__ + j * h_dim1] - u[i__] * t[j] - t[
+ i__] * u[j] + u[i__] * u[j] * s;
+/* L60: */
+ }
+/* L70: */
+ }
+ return 0;
+} /* zhz_ */
+
+/* ---------------------- */
+/* | S O L V E W | */
+/* ---------------------- */
+/* Subroutine */ int solvew_(integer *nr, integer *n, doublereal *al,
+ doublereal *u, doublereal *w, doublereal *b)
+{
+ /* System generated locals */
+ integer al_dim1, al_offset, i__1, i__2;
+
+ /* Local variables */
+ doublereal d__;
+ integer i__, j;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *);
+
+
+/* PURPOSE: */
+/* SOLVE L*W=ZT*V */
+
+/* ---------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* AL(N-1,N-1) --> LOWER TRIAGULAR MATRIX */
+/* U(N) --> VECTOR TO FORM Z */
+/* W(N) --> ON INPUT : VECTOR V IN SYSTEM OF LINEAR EQUATIONS */
+/* ON OUTPUT: SOLUTION OF SYSTEM OF LINEAR EQUATIONS */
+/* B(N) --> WORKSPACE TO STORE ZT*V */
+
+/* FORM ZT*V (STORED IN B) */
+ /* Parameter adjustments */
+ al_dim1 = *nr;
+ al_offset = 1 + al_dim1;
+ al -= al_offset;
+ --b;
+ --w;
+ --u;
+
+ /* Function Body */
+ d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
+ d__ /= 2.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ b[i__] = 0.;
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ b[i__] += u[j] * u[i__] * w[j] / d__;
+/* L15: */
+ }
+ b[i__] = w[i__] - b[i__];
+/* L20: */
+ }
+
+/* SOLVE LW=ZT*V */
+ i__1 = *n - 1;
+ forslv_(nr, &i__1, &al[al_offset], &w[1], &b[1]);
+ return 0;
+} /* solvew_ */
+
+/* ---------------------- */
+/* | D S T A R | */
+/* ---------------------- */
+/* Subroutine */ int dstar_(integer *nr, integer *n, doublereal *u,
+ doublereal *s, doublereal *w1, doublereal *w2, doublereal *w3,
+ doublereal *sigma, doublereal *al, doublereal *d__)
+{
+ /* System generated locals */
+ integer al_dim1, al_offset, i__1;
+
+ /* Local variables */
+ integer i__;
+ doublereal utt, utu;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+ doublereal temp;
+ extern /* Subroutine */ int bakslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *);
+
+
+/* PURPOSE: */
+/* COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) */
+
+/* ------------------------------------------------------------------------ */
+
+/* PARAMETERS: */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* U(N) --> VECTOR TO FORM Z */
+/* S(N) --> PREVIOUS STEP */
+/* W1(N) --> L**-1*ZT*A, WHERE A IS DESCRIBED IN SUBROUTINE SLVMDL */
+/* W2(N) --> L**-1*ZT*SH, WHERE H IS CURRENT HESSIAN */
+/* W3(N) --> L**-1*ZT*G, WHERE G IS CURRENT GRADIENT */
+/* SIGMA --> SOLUTION FOR REDUCED ONE VARIABLE MODEL */
+/* AL(N-1,N-1) --> LOWER TRIANGULAR MATRIX L */
+/* D(N) --> TENSOR STEP */
+
+ /* Parameter adjustments */
+ al_dim1 = *nr;
+ al_offset = 1 + al_dim1;
+ al -= al_offset;
+ --d__;
+ --w3;
+ --w2;
+ --w1;
+ --s;
+ --u;
+
+ /* Function Body */
+ if (*n == 1) {
+ d__[1] = *sigma * s[1];
+ } else {
+
+/* COMPUTE T(SIGMA)=-(ZTHZ)*ZT*(G+SIGMA*SH+SIGMA**2*A/2) (STORED IN D) */
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w2[i__] = w3[i__] + *sigma * w2[i__] + w1[i__] * .5 * *sigma * *
+ sigma;
+/* L10: */
+ }
+ i__1 = *n - 1;
+ bakslv_(nr, &i__1, &al[al_offset], &d__[1], &w2[1]);
+ d__[*n] = 0.;
+
+/* COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) */
+ utu = ddot_(n, &u[1], &c__1, &u[1], &c__1);
+ utt = ddot_(n, &u[1], &c__1, &d__[1], &c__1);
+ temp = utt / utu;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ d__[i__] = *sigma * s[i__] - (d__[i__] - u[i__] * 2. * temp);
+/* L50: */
+ }
+ }
+ return 0;
+} /* dstar_ */
+
+/* ---------------------- */
+/* | M K M D L | */
+/* ---------------------- */
+/* Subroutine */ int mkmdl_(integer *nr, integer *n, doublereal *f,
+ doublereal *fp, doublereal *g, doublereal *gp, doublereal *s,
+ doublereal *h__, doublereal *alpha, doublereal *beta, doublereal *sh,
+ doublereal *a)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, i__1, i__2;
+
+ /* Local variables */
+ integer i__, j;
+ doublereal b1, b2, gs, gps, shs, sts;
+ extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
+ integer *);
+
+
+/* PURPOSE: */
+/* FORM TENSOR MODEL */
+
+/* ----------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* F --> CURRENT FUNCTION VALUE */
+/* FP --> PREVIOUS FUNCTION VALUE */
+/* G(N) --> CURRENT GRADIENT */
+/* GP(N) --> PREVIOUS GRADIENT */
+/* S(N) --> STEP TO PREVIOUS POINT */
+/* H(N,N) --> HESSIAN */
+/* ALPHA <-- SCALAR TO FORM 3RD ORDER TERM OF TENSOR MODEL */
+/* BETA <-- SCALAR TO FORM 4TH ORDER TERM OF TENSOR MODEL */
+/* SH(N) <-- SH */
+/* A(N) <-- A=2*(GP-G-SH-S*BETA/(6*STS)) */
+
+
+/* COMPUTE SH */
+ /* Parameter adjustments */
+ --a;
+ --sh;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+ --s;
+ --gp;
+ --g;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ sh[i__] = 0.;
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ sh[i__] += s[j] * h__[j + i__ * h_dim1];
+/* L10: */
+ }
+/* L20: */
+ }
+ gs = ddot_(n, &g[1], &c__1, &s[1], &c__1);
+ gps = ddot_(n, &gp[1], &c__1, &s[1], &c__1);
+ shs = ddot_(n, &sh[1], &c__1, &s[1], &c__1);
+ b1 = gps - gs - shs;
+ b2 = *fp - *f - gs - shs * .5;
+ *alpha = b2 * 24. - b1 * 6.;
+ *beta = b1 * 24. - b2 * 72.;
+
+/* COMPUTE A */
+ sts = ddot_(n, &s[1], &c__1, &s[1], &c__1);
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ a[i__] = (gp[i__] - g[i__] - sh[i__] - s[i__] * *beta / (sts * 6.)) *
+ 2.;
+/* L50: */
+ }
+ return 0;
+} /* mkmdl_ */
+
+/* ---------------------- */
+/* | S I G M A | */
+/* ---------------------- */
+/* Subroutine */ int sigma_(doublereal *sgstar, doublereal *a, doublereal *b,
+ doublereal *c__, doublereal *d__)
+{
+ doublereal s1, s2, s3;
+ extern /* Subroutine */ int roots_(doublereal *, doublereal *, doublereal
+ *, doublereal *, doublereal *, doublereal *, doublereal *),
+ sortrt_(doublereal *, doublereal *, doublereal *);
+
+
+/* PURPOSE: */
+/* COMPUTE DESIRABLE ROOT OF REDUCED ONE VARIABLE EQUATION */
+
+/* ------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+/* SGSTAR --> DESIRABLE ROOT */
+/* A --> COEFFICIENT OF 3RD ORDER TERM */
+/* B --> COEFFICIENT OF 2ND ORDER TERM */
+/* C --> COEFFICIENT OF 1ST ORDER TERM */
+/* D --> COEFFICIENT OF CONSTANT TERM */
+
+
+/* COMPUTE ALL THREE ROOTS */
+ roots_(&s1, &s2, &s3, a, b, c__, d__);
+
+/* SORT ROOTS */
+ sortrt_(&s1, &s2, &s3);
+
+/* CHOOSE DESIRABLE ROOT */
+ if (*a > 0.) {
+ *sgstar = s3;
+ if (s2 >= 0.) {
+ *sgstar = s1;
+ }
+ } else {
+/* IF A < 0 THEN */
+ *sgstar = s2;
+ if (s1 > 0. || s3 < 0.) {
+ if (s1 > 0.) {
+ *sgstar = s1;
+ } else {
+ *sgstar = s3;
+ }
+ *a = 0.;
+ }
+ }
+ return 0;
+} /* sigma_ */
+
+/* ---------------------- */
+/* | R O O T S | */
+/* ---------------------- */
+/* Subroutine */ int roots_(doublereal *s1, doublereal *s2, doublereal *s3,
+ doublereal *a, doublereal *b, doublereal *c__, doublereal *d__)
+{
+ /* System generated locals */
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal), pow_dd(doublereal *, doublereal *), acos(
+ doublereal), cos(doublereal);
+
+ /* Local variables */
+ doublereal q, r__, s, t, v, a1, a2, a3, pi, temp, theta;
+
+
+/* PURPOSE: */
+/* COMPUTE ROOT(S) OF 3RD ORDER EQUATION */
+
+/* --------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+/* S1 <-- ROOT (IF THREE ROOTS ARE */
+/* S2 <-- ROOT EQUAL, THEN S1=S2=S3) */
+/* S3 <-- ROOT */
+/* A --> COEFFICIENT OF 3RD ORDER TERM */
+/* B --> COEFFICIENT OF 2ND ORDER TERM */
+/* C --> COEFFICIENT OF 1ST ORDER TERM */
+/* D --> COEFFICIENT OF CONSTANT TERM */
+
+/* SET VALUE OF PI */
+ pi = 3.141592653589793;
+ a1 = *b / *a;
+ a2 = *c__ / *a;
+ a3 = *d__ / *a;
+ q = (a2 * 3. - a1 * a1) / 9.;
+ r__ = (a1 * 9. * a2 - a3 * 27. - a1 * 2. * a1 * a1) / 54.;
+ v = q * q * q + r__ * r__;
+ if (v > 0.) {
+ s = r__ + sqrt(v);
+ t = r__ - sqrt(v);
+ if (t < 0.) {
+ d__1 = -t;
+ t = -pow_dd(&d__1, &c_b250);
+ } else {
+ t = pow_dd(&t, &c_b250);
+ }
+ if (s < 0.) {
+ d__1 = -s;
+ s = -pow_dd(&d__1, &c_b250);
+ } else {
+ s = pow_dd(&s, &c_b250);
+ }
+ *s1 = s + t - a1 / 3.;
+ *s3 = *s1;
+ *s2 = *s1;
+ } else {
+ temp = r__ / sqrt(-pow_dd(&q, &c_b384));
+ theta = acos(temp);
+ theta /= 3.;
+ temp = sqrt(-q) * 2.;
+ *s1 = temp * cos(theta) - a1 / 3.;
+ *s2 = temp * cos(theta + pi * 2. / 3.) - a1 / 3.;
+ *s3 = temp * cos(theta + pi * 4. / 3.) - a1 / 3.;
+ }
+ return 0;
+} /* roots_ */
+
+/* ---------------------- */
+/* | S O R T R T | */
+/* ---------------------- */
+/* Subroutine */ int sortrt_(doublereal *s1, doublereal *s2, doublereal *s3)
+{
+ doublereal t;
+
+
+/* PURPOSE: */
+/* SORT ROOTS INTO ASCENDING ORDER */
+
+/* ----------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+/* S1 <--> ROOT */
+/* S2 <--> ROOT */
+/* S3 <--> ROOT */
+
+ if (*s1 > *s2) {
+ t = *s1;
+ *s1 = *s2;
+ *s2 = t;
+ }
+ if (*s2 > *s3) {
+ t = *s2;
+ *s2 = *s3;
+ *s3 = t;
+ }
+ if (*s1 > *s2) {
+ t = *s1;
+ *s1 = *s2;
+ *s2 = t;
+ }
+ return 0;
+} /* sortrt_ */
+
+/* ---------------------- */
+/* | F S T O F D | */
+/* ---------------------- */
+/* Subroutine */ int fstofd_(integer *nr, integer *m, integer *n, doublereal *
+ xpls, S_fp fcn, doublereal *fpls, doublereal *a, doublereal *typx,
+ doublereal *rnoise, doublereal *fhat, integer *icase, integer *nfcnt)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+ doublereal d__1, d__2;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ integer i__, j, jp1, nm1;
+ doublereal xtmpj, stepsz;
+
+/* PURPOSE */
+/* ------- */
+/* FIND FIRST ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" TO THE */
+/* FIRST DERIVATIVE OF THE FUNCTION DEFINED BY THE SUBPROGRAM "FNAME" */
+/* EVALUATED AT THE NEW ITERATE "XPLS". */
+
+
+/* FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE: */
+/* 1) THE FIRST DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN */
+/* ANALYTIC USER ROUTINE HAS BEEN SUPPLIED; */
+/* 2) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION */
+/* IF NO ANALYTIC USER ROUTINE HAS BEEN SUPPLIED FOR THE HESSIAN BUT */
+/* ONE HAS BEEN SUPPLIED FOR THE GRADIENT ("FCN") AND IF THE */
+/* OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE */
+
+/* NOTE */
+/* ---- */
+/* _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE GRADIENT OF THE FUNCTION */
+/* (FCN). FCN(X) # F: R(N)-->R(1) */
+/* _M=N (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION */
+/* FCN(X) # F: R(N)-->R(N). */
+/* _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN OF THE OPTIMIZATIO */
+/* FUNCTION, WHERE THE HESSIAN IS THE FIRST DERIVATIVE OF "FCN" */
+
+/* PARAMETERS */
+/* ---------- */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* M --> NUMBER OF ROWS IN A */
+/* N --> NUMBER OF COLUMNS IN A; DIMENSION OF PROBLEM */
+/* XPLS(N) --> NEW ITERATE: X[K] */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
+/* FPLS(M) --> _M=1 (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE: */
+/* FCN(XPLS) */
+/* _M=N (OPTIMIZATION) VALUE OF FIRST DERIVATIVE */
+/* (GRADIENT) GIVEN BY USER FUNCTION FCN */
+/* _M=N (SYSTEMS) FUNCTION VALUE OF ASSOCIATED */
+/* MINIMIZATION FUNCTION */
+/* A(NR,N) <-- FINITE DIFFERENCE APPROXIMATION (SEE NOTE). ONLY */
+/* LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED */
+/* RNOISE --> RELATIVE NOISE IN FCN [F(X)] */
+/* FHAT(M) --> WORKSPACE */
+/* ICASE --> =1 OPTIMIZATION (GRADIENT) */
+/* =2 SYSTEMS */
+/* =3 OPTIMIZATION (HESSIAN) */
+/* NFCNT <--> NUMBER OF FUNCTION EVALUTIONS */
+
+/* INTERNAL VARIABLES */
+/* ------------------ */
+/* STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION */
+
+
+/* FIND J-TH COLUMN OF A */
+/* EACH COLUMN IS DERIVATIVE OF F(FCN) WITH RESPECT TO XPLS(J) */
+
+ /* Parameter adjustments */
+ --fhat;
+ --fpls;
+ --typx;
+ a_dim1 = *nr;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --xpls;
+
+ /* Function Body */
+ i__1 = *n;
+ for (j = 1; j <= i__1; ++j) {
+ xtmpj = xpls[j];
+/* Computing MAX */
+ d__2 = (d__1 = xpls[j], abs(d__1));
+ stepsz = sqrt(*rnoise) * max(d__2,1.);
+ xpls[j] = xtmpj + stepsz;
+ (*fcn)(n, &xpls[1], &fhat[1]);
+ ++(*nfcnt);
+ xpls[j] = xtmpj;
+ i__2 = *m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] = (fhat[i__] - fpls[i__]) / stepsz;
+ a[i__ + j * a_dim1] *= typx[j];
+/* L20: */
+ }
+/* L30: */
+ }
+ if (*icase != 3) {
+ return 0;
+ }
+
+/* IF COMPUTING HESSIAN, A MUST BE SYMMETRIC */
+
+ if (*n == 1) {
+ return 0;
+ }
+ nm1 = *n - 1;
+ i__1 = nm1;
+ for (j = 1; j <= i__1; ++j) {
+ jp1 = j + 1;
+ i__2 = *m;
+ for (i__ = jp1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] = (a[i__ + j * a_dim1] + a[j + i__ * a_dim1])
+ / 2.f;
+/* L40: */
+ }
+/* L50: */
+ }
+ return 0;
+} /* fstofd_ */
+
+/* ---------------------- */
+/* | S N D O F D | */
+/* ---------------------- */
+/* Subroutine */ int sndofd_(integer *nr, integer *n, doublereal *xpls, S_fp
+ fcn, doublereal *fpls, doublereal *a, doublereal *typx, doublereal *
+ rnoise, doublereal *stepsz, doublereal *anbr, integer *nfcnt)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+ doublereal d__1, d__2;
+
+ /* Builtin functions */
+ double pow_dd(doublereal *, doublereal *);
+
+ /* Local variables */
+ integer i__, j, ip1;
+ doublereal ov3, fhat, xtmpi, xtmpj;
+
+/* PURPOSE */
+/* ------- */
+/* FIND SECOND ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" */
+/* TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE SUBP */
+/* "FCN" EVALUATED AT THE NEW ITERATE "XPLS" */
+
+/* FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE */
+/* 1) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION */
+/* IF NO ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER */
+/* THE GRADIENT OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION */
+/* "FCN" IS INEXPENSIVE TO EVALUATE. */
+
+/* PARAMETERS */
+/* ---------- */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* XPLS(N) --> NEW ITERATE: X[K] */
+/* FCN --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
+/* FPLS --> FUNCTION VALUE AT NEW ITERATE, F(XPLS) */
+/* A(N,N) <-- FINITE DIFFERENCE APPROXIMATION TO HESSIAN */
+/* ONLY LOWER TRIANGULAR MATRIX AND DIAGONAL */
+/* ARE RETURNED */
+/* RNOISE --> RELATIVE NOISE IN FNAME [F(X)] */
+/* STEPSZ(N) --> WORKSPACE (STEPSIZE IN I-TH COMPONENT DIRECTION) */
+/* ANBR(N) --> WORKSPACE (NEIGHBOR IN I-TH DIRECTION) */
+/* NFCNT <--> NUMBER OF FUNCTION EVALUATIONS */
+
+
+
+/* FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION */
+/* OF I-TH UNIT VECTOR. */
+
+ /* Parameter adjustments */
+ --anbr;
+ --stepsz;
+ --typx;
+ a_dim1 = *nr;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ --xpls;
+
+ /* Function Body */
+ ov3 = .33333333333333331;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ xtmpi = xpls[i__];
+/* Computing MAX */
+ d__2 = (d__1 = xpls[i__], abs(d__1));
+ stepsz[i__] = pow_dd(rnoise, &ov3) * max(d__2,1.);
+ xpls[i__] = xtmpi + stepsz[i__];
+ (*fcn)(n, &xpls[1], &anbr[i__]);
+ ++(*nfcnt);
+ xpls[i__] = xtmpi;
+/* L10: */
+ }
+
+/* CALCULATE COLUMN I OF A */
+
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ xtmpi = xpls[i__];
+ xpls[i__] = xtmpi + stepsz[i__] * 2.f;
+ (*fcn)(n, &xpls[1], &fhat);
+ ++(*nfcnt);
+ a[i__ + i__ * a_dim1] = (*fpls - anbr[i__] + (fhat - anbr[i__])) / (
+ stepsz[i__] * stepsz[i__]);
+ a[i__ + i__ * a_dim1] *= typx[i__] * typx[i__];
+
+/* CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN */
+ if (i__ == *n) {
+ goto L25;
+ }
+ xpls[i__] = xtmpi + stepsz[i__];
+ ip1 = i__ + 1;
+ i__2 = *n;
+ for (j = ip1; j <= i__2; ++j) {
+ xtmpj = xpls[j];
+ xpls[j] = xtmpj + stepsz[j];
+ (*fcn)(n, &xpls[1], &fhat);
+ ++(*nfcnt);
+ a[j + i__ * a_dim1] = (*fpls - anbr[i__] + (fhat - anbr[j])) / (
+ stepsz[i__] * stepsz[j]);
+ a[j + i__ * a_dim1] *= typx[i__] * typx[j];
+ xpls[j] = xtmpj;
+/* L20: */
+ }
+L25:
+ xpls[i__] = xtmpi;
+/* L30: */
+ }
+ return 0;
+} /* sndofd_ */
+
+/* ---------------------- */
+/* | B A K S L V | */
+/* ---------------------- */
+/* Subroutine */ int bakslv_(integer *nr, integer *n, doublereal *a,
+ doublereal *x, doublereal *b)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1;
+
+ /* Local variables */
+ integer i__, j, ip1;
+ doublereal sum;
+
+
+/* PURPOSE */
+/* ------- */
+/* SOLVE AX=B WHERE A IS UPPER TRIANGULAR MATRIX. */
+/* NOTE THAT A IS INPUT AS A LOWER TRIANGULAR MATRIX AND */
+/* THAT THIS ROUTINE TAKES ITS TRANSPOSE IMPLICITLY. */
+
+/* PARAMETERS */
+/* ---------- */
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* A(N,N) --> LOWER TRIANGULAR MATRIX (PRESERVED) */
+/* X(N) <-- SOLUTION VECTOR */
+/* B(N) --> RIGHT-HAND SIDE VECTOR */
+
+/* NOTE */
+/* ---- */
+/* IF B IS NO LONGER REQUIRED BY CALLING ROUTINE, */
+/* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE. */
+
+
+/* SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE) */
+
+ /* Parameter adjustments */
+ --b;
+ --x;
+ a_dim1 = *nr;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ i__ = *n;
+ x[i__] = b[i__] / a[i__ + i__ * a_dim1];
+ if (*n == 1) {
+ return 0;
+ }
+L30:
+ ip1 = i__;
+ --i__;
+ sum = 0.f;
+ i__1 = *n;
+ for (j = ip1; j <= i__1; ++j) {
+ sum += a[j + i__ * a_dim1] * x[j];
+/* L40: */
+ }
+ x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
+ if (i__ > 1) {
+ goto L30;
+ }
+ return 0;
+} /* bakslv_ */
+
+/* ---------------------- */
+/* | F O R S L V | */
+/* ---------------------- */
+/* Subroutine */ int forslv_(integer *nr, integer *n, doublereal *a,
+ doublereal *x, doublereal *b)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+
+ /* Local variables */
+ integer i__, j, im1;
+ doublereal sum;
+
+
+/* PURPOSE */
+/* -------- */
+/* SOLVE AX=B WHERE A IS LOWER TRIANGULAR MATRIX */
+
+/* PARAMETERS */
+/* --------- */
+
+/* NR -----> ROW DIMENSION OF MATRIX */
+/* N -----> DIMENSION OF PROBLEM */
+/* A(N,N) -----> LOWER TRIANGULAR MATRIX (PRESERVED) */
+/* X(N) <---- SOLUTION VECTOR */
+/* B(N) ----> RIGHT-HAND SIDE VECTOR */
+
+/* NOTE */
+/* ----- */
+/* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE */
+
+
+/* SOLVE LX=B. (FOREWARD SOLVE) */
+
+ /* Parameter adjustments */
+ --b;
+ --x;
+ a_dim1 = *nr;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ x[1] = b[1] / a[a_dim1 + 1];
+ if (*n == 1) {
+ return 0;
+ }
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ sum = 0.f;
+ im1 = i__ - 1;
+ i__2 = im1;
+ for (j = 1; j <= i__2; ++j) {
+ sum += a[i__ + j * a_dim1] * x[j];
+/* L10: */
+ }
+ x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
+/* L20: */
+ }
+ return 0;
+} /* forslv_ */
+
+/* ---------------------- */
+/* | C H O L D R | */
+/* ---------------------- */
+/* Subroutine */ int choldr_(integer *nr, integer *n, doublereal *h__,
+ doublereal *g, doublereal *eps, integer *pivot, doublereal *e,
+ doublereal *diag, doublereal *addmax)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, i__1, i__2, i__3;
+
+ /* Builtin functions */
+ double pow_dd(doublereal *, doublereal *), sqrt(doublereal);
+
+ /* Local variables */
+ integer i__, j, k;
+ doublereal tau1, tau2;
+ logical redo;
+ doublereal temp;
+ extern /* Subroutine */ int modchl_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *, integer *,
+ doublereal *);
+
+
+/* PURPOSE: */
+/* DRIVER FOR CHOLESKY DECOMPOSITION */
+
+/* ---------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION */
+/* N --> DIMENSION OF PROBLEM */
+/* H(N,N) --> MATRIX */
+/* G(N) --> WORK SPACE */
+/* EPS --> MACHINE EPSILON */
+/* PIVOT(N) --> PIVOTING VECTOR */
+/* E(N) --> DIAGONAL MATRIX ADDED TO H FOR MAKING H P.D. */
+/* DIAG(N) --> DIAGONAL OF H */
+/* ADDMAX --> ADDMAX * I IS ADDED TO H */
+ /* Parameter adjustments */
+ --diag;
+ --e;
+ --pivot;
+ --g;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+
+ /* Function Body */
+ redo = FALSE_;
+
+/* SAVE DIAGONAL OF H */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ diag[i__] = h__[i__ + i__ * h_dim1];
+/* L10: */
+ }
+ tau1 = pow_dd(eps, &c_b250);
+ tau2 = tau1;
+ modchl_(nr, n, &h__[h_offset], &g[1], eps, &tau1, &tau2, &pivot[1], &e[1])
+ ;
+ *addmax = e[*n];
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (pivot[i__] != i__) {
+ redo = TRUE_;
+ }
+/* L22: */
+ }
+ if (*addmax > 0. || redo) {
+/* ******************************** */
+/* * */
+/* H IS NOT P.D. * */
+/* * */
+/* ******************************** */
+
+/* H=H+UI */
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ i__2 = i__ - 1;
+ for (j = 1; j <= i__2; ++j) {
+ h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
+/* L32: */
+ }
+/* L30: */
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ pivot[i__] = i__;
+ h__[i__ + i__ * h_dim1] = diag[i__] + *addmax;
+/* L34: */
+ }
+/* ******************************** */
+/* * */
+/* COMPUTE L * */
+/* * */
+/* ******************************** */
+ i__1 = *n;
+ for (j = 1; j <= i__1; ++j) {
+
+/* COMPUTE L(J,J) */
+ temp = 0.;
+ if (j > 1) {
+ i__2 = j - 1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp += h__[j + i__ * h_dim1] * h__[j + i__ * h_dim1];
+/* L41: */
+ }
+ }
+ h__[j + j * h_dim1] -= temp;
+ h__[j + j * h_dim1] = sqrt(h__[j + j * h_dim1]);
+
+/* COMPUTE L(I,J) */
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ temp = 0.;
+ if (j > 1) {
+ i__3 = j - 1;
+ for (k = 1; k <= i__3; ++k) {
+ temp += h__[i__ + k * h_dim1] * h__[j + k * h_dim1];
+/* L45: */
+ }
+ }
+ h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1] - temp;
+ h__[i__ + j * h_dim1] /= h__[j + j * h_dim1];
+/* L43: */
+ }
+/* L40: */
+ }
+ }
+ return 0;
+} /* choldr_ */
+
+/* ---------------------- */
+/* | M O D C H L | */
+/* ---------------------- */
+/* ********************************************************************* */
+
+/* SUBROUTINE NAME: MODCHL */
+
+/* AUTHORS : ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
+
+/* DATE : DECEMBER, 1988 */
+
+/* PURPOSE : PERFORM A MODIFIED CHOLESKY FACTORIZATION */
+/* OF THE FORM (PTRANSPOSE)AP + E = L(LTRANSPOSE), */
+/* WHERE L IS STORED IN THE LOWER TRIANGLE OF THE */
+/* ORIGINAL MATRIX A. */
+/* THE FACTORIZATION HAS 2 PHASES: */
+/* PHASE 1: PIVOT ON THE MAXIMUM DIAGONAL ELEMENT. */
+/* CHECK THAT THE NORMAL CHOLESKY UPDATE */
+/* WOULD RESULT IN A POSITIVE DIAGONAL */
+/* AT THE CURRENT ITERATION, AND */
+/* IF SO, DO THE NORMAL CHOLESKY UPDATE, */
+/* OTHERWISE SWITCH TO PHASE 2. */
+/* PHASE 2: PIVOT ON THE MINIMUM OF THE NEGATIVES */
+/* OF THE LOWER GERSCHGORIN BOUND */
+/* ESTIMATES. */
+/* COMPUTE THE AMOUNT TO ADD TO THE */
+/* PIVOT ELEMENT AND ADD THIS */
+/* TO THE PIVOT ELEMENT. */
+/* DO THE CHOLESKY UPDATE. */
+/* UPDATE THE ESTIMATES OF THE */
+/* GERSCHGORIN BOUNDS. */
+
+/* INPUT : NDIM - LARGEST DIMENSION OF MATRIX THAT */
+/* WILL BE USED */
+
+/* N - DIMENSION OF MATRIX A */
+
+/* A - N*N SYMMETRIC MATRIX (ONLY LOWER TRIANGULAR */
+/* PORTION OF A, INCLUDING THE MAIN DIAGONAL, IS USED) */
+
+/* G - N*1 WORK ARRAY */
+
+/* MCHEPS - MACHINE PRECISION */
+
+/* TAU1 - TOLERANCE USED FOR DETERMINING WHEN TO SWITCH TO */
+/* PHASE 2 */
+
+/* TAU2 - TOLERANCE USED FOR DETERMINING THE MAXIMUM */
+/* CONDITION NUMBER OF THE FINAL 2X2 SUBMATRIX. */
+
+
+/* OUTPUT : L - STORED IN THE MATRIX A (IN LOWER TRIANGULAR */
+/* PORTION OF A, INCLUDING THE MAIN DIAGONAL) */
+
+/* P - A RECORD OF HOW THE ROWS AND COLUMNS */
+/* OF THE MATRIX WERE PERMUTED WHILE */
+/* PERFORMING THE DECOMPOSITION */
+
+/* E - N*1 ARRAY, THE ITH ELEMENT IS THE */
+/* AMOUNT ADDED TO THE DIAGONAL OF A */
+/* AT THE ITH ITERATION */
+
+
+/* ************************************************************************ */
+/* Subroutine */ int modchl_(integer *ndim, integer *n, doublereal *a,
+ doublereal *g, doublereal *mcheps, doublereal *tau1, doublereal *tau2,
+ integer *p, doublereal *e)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2, i__3;
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ integer i__, j, k, jp1;
+ doublereal maxd, ming;
+ extern /* Subroutine */ int init_(integer *, integer *, doublereal *,
+ logical *, doublereal *, integer *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *);
+ doublereal temp, gamma, delta, jdmin;
+ integer imaxd, iming;
+ extern /* Subroutine */ int fin2x2_(integer *, integer *, doublereal *,
+ doublereal *, integer *, doublereal *, doublereal *, doublereal *)
+ ;
+ doublereal tdmin;
+ integer itemp;
+ doublereal normj, delta1;
+ logical phase1;
+ extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
+ integer *, doublereal *);
+ doublereal taugam, tempjj;
+
+
+/* J - CURRENT ITERATION NUMBER */
+/* IMING - INDEX OF THE ROW WITH THE MIN. OF THE */
+/* NEG. LOWER GERSCH. BOUNDS */
+/* IMAXD - INDEX OF THE ROW WITH THE MAXIMUM DIAG. */
+/* ELEMENT */
+/* I,ITEMP,JPL,K - TEMPORARY INTEGER VARIABLES */
+/* DELTA - AMOUNT TO ADD TO AJJ AT THE JTH ITERATION */
+/* GAMMA - THE MAXIMUM DIAGONAL ELEMENT OF THE ORIGINAL */
+/* MATRIX A. */
+/* NORMJ - THE 1 NORM OF A(COLJ), ROWS J+1 --> N. */
+/* MING - THE MINIMUM OF THE NEG. LOWER GERSCH. BOUNDS */
+/* MAXD - THE MAXIMUM DIAGONAL ELEMENT */
+/* TAUGAM - TAU1 * GAMMA */
+/* PHASE1 - LOGICAL, TRUE IF IN PHASE1, OTHERWISE FALSE */
+/* DELTA1,TEMP,JDMIN,TDMIN,TEMPJJ - TEMPORARY DOUBLE PRECISION VARS. */
+
+ /* Parameter adjustments */
+ --e;
+ --p;
+ --g;
+ a_dim1 = *ndim;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ init_(n, ndim, &a[a_offset], &phase1, &delta, &p[1], &g[1], &e[1], &ming,
+ tau1, &gamma, &taugam);
+/* CHECK FOR N=1 */
+ if (*n == 1) {
+ delta = *tau2 * (d__1 = a[a_dim1 + 1], abs(d__1)) - a[a_dim1 + 1];
+ if (delta > 0.) {
+ e[1] = delta;
+ }
+ if (a[a_dim1 + 1] == 0.) {
+ e[1] = *tau2;
+ }
+ a[a_dim1 + 1] = sqrt(a[a_dim1 + 1] + e[1]);
+ }
+
+
+ i__1 = *n - 1;
+ for (j = 1; j <= i__1; ++j) {
+
+/* PHASE 1 */
+
+ if (phase1) {
+
+/* FIND INDEX OF MAXIMUM DIAGONAL ELEMENT A(I,I) WHERE I>=J */
+
+ maxd = a[j + j * a_dim1];
+ imaxd = j;
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ if (maxd < a[i__ + i__ * a_dim1]) {
+ maxd = a[i__ + i__ * a_dim1];
+ imaxd = i__;
+ }
+/* L20: */
+ }
+
+/* PIVOT TO THE TOP THE ROW AND COLUMN WITH THE MAX DIAG */
+
+ if (imaxd != j) {
+
+/* SWAP ROW J WITH ROW OF MAX DIAG */
+
+ i__2 = j - 1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = a[imaxd + i__ * a_dim1];
+ a[imaxd + i__ * a_dim1] = temp;
+/* L30: */
+ }
+
+/* SWAP COLJ AND ROW MAXDIAG BETWEEN J AND MAXDIAG */
+
+ i__2 = imaxd - 1;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = a[imaxd + i__ * a_dim1];
+ a[imaxd + i__ * a_dim1] = temp;
+/* L35: */
+ }
+
+/* SWAP COLUMN J WITH COLUMN OF MAX DIAG */
+
+ i__2 = *n;
+ for (i__ = imaxd + 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = a[i__ + imaxd * a_dim1];
+ a[i__ + imaxd * a_dim1] = temp;
+/* L40: */
+ }
+
+/* SWAP DIAG ELEMENTS */
+
+ temp = a[j + j * a_dim1];
+ a[j + j * a_dim1] = a[imaxd + imaxd * a_dim1];
+ a[imaxd + imaxd * a_dim1] = temp;
+
+/* SWAP ELEMENTS OF THE PERMUTATION VECTOR */
+
+ itemp = p[j];
+ p[j] = p[imaxd];
+ p[imaxd] = itemp;
+ }
+/* CHECK TO SEE WHETHER THE NORMAL CHOLESKY UPDATE FOR THIS */
+/* ITERATION WOULD RESULT IN A POSITIVE DIAGONAL, */
+/* AND IF NOT THEN SWITCH TO PHASE 2. */
+ jp1 = j + 1;
+ tempjj = a[j + j * a_dim1];
+ if (tempjj > 0.) {
+ jdmin = a[jp1 + jp1 * a_dim1];
+ i__2 = *n;
+ for (i__ = jp1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1] * a[i__ + j * a_dim1] / tempjj;
+ tdmin = a[i__ + i__ * a_dim1] - temp;
+ jdmin = min(jdmin,tdmin);
+/* L60: */
+ }
+ if (jdmin < taugam) {
+ phase1 = FALSE_;
+ }
+ } else {
+ phase1 = FALSE_;
+ }
+ if (phase1) {
+
+/* DO THE NORMAL CHOLESKY UPDATE IF STILL IN PHASE 1 */
+
+ a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
+ tempjj = a[j + j * a_dim1];
+ i__2 = *n;
+ for (i__ = jp1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] /= tempjj;
+/* L70: */
+ }
+ i__2 = *n;
+ for (i__ = jp1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ i__3 = i__;
+ for (k = jp1; k <= i__3; ++k) {
+ a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
+/* L75: */
+ }
+/* L80: */
+ }
+ if (j == *n - 1) {
+ a[*n + *n * a_dim1] = sqrt(a[*n + *n * a_dim1]);
+ }
+ } else {
+
+/* CALCULATE THE NEGATIVES OF THE LOWER GERSCHGORIN BOUNDS */
+
+ gersch_(ndim, n, &a[a_offset], &j, &g[1]);
+ }
+ }
+
+/* PHASE 2 */
+
+ if (! phase1) {
+ if (j != *n - 1) {
+
+/* FIND THE MINIMUM NEGATIVE GERSHGORIN BOUND */
+
+ iming = j;
+ ming = g[j];
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ if (ming > g[i__]) {
+ ming = g[i__];
+ iming = i__;
+ }
+/* L90: */
+ }
+
+/* PIVOT TO THE TOP THE ROW AND COLUMN WITH THE */
+/* MINIMUM NEGATIVE GERSCHGORIN BOUND */
+
+ if (iming != j) {
+
+/* SWAP ROW J WITH ROW OF MIN GERSCH BOUND */
+
+ i__2 = j - 1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = a[iming + i__ * a_dim1];
+ a[iming + i__ * a_dim1] = temp;
+/* L100: */
+ }
+
+/* SWAP COLJ WITH ROW IMING FROM J TO IMING */
+
+ i__2 = iming - 1;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = a[iming + i__ * a_dim1];
+ a[iming + i__ * a_dim1] = temp;
+/* L105: */
+ }
+
+/* SWAP COLUMN J WITH COLUMN OF MIN GERSCH BOUND */
+
+ i__2 = *n;
+ for (i__ = iming + 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = a[i__ + iming * a_dim1];
+ a[i__ + iming * a_dim1] = temp;
+/* L110: */
+ }
+
+/* SWAP DIAGONAL ELEMENTS */
+
+ temp = a[j + j * a_dim1];
+ a[j + j * a_dim1] = a[iming + iming * a_dim1];
+ a[iming + iming * a_dim1] = temp;
+
+/* SWAP ELEMENTS OF THE PERMUTATION VECTOR */
+
+ itemp = p[j];
+ p[j] = p[iming];
+ p[iming] = itemp;
+
+/* SWAP ELEMENTS OF THE NEGATIVE GERSCHGORIN BOUNDS VECTOR */
+
+ temp = g[j];
+ g[j] = g[iming];
+ g[iming] = temp;
+ }
+
+/* CALCULATE DELTA AND ADD TO THE DIAGONAL. */
+/* DELTA=MAX{0,-A(J,J) + MAX{NORMJ,TAUGAM},DELTA_PREVIOUS} */
+/* WHERE NORMJ=SUM OF |A(I,J)|,FOR I=1,N, */
+/* DELTA_PREVIOUS IS THE DELTA COMPUTED AT THE PREVIOUS ITERATION, */
+/* AND TAUGAM IS TAU1*GAMMA. */
+
+ normj = 0.f;
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ normj += (d__1 = a[i__ + j * a_dim1], abs(d__1));
+/* L140: */
+ }
+ temp = max(normj,taugam);
+ delta1 = temp - a[j + j * a_dim1];
+ temp = 0.f;
+ delta1 = max(temp,delta1);
+ delta = max(delta1,delta);
+ e[j] = delta;
+ a[j + j * a_dim1] += e[j];
+
+/* UPDATE THE GERSCHGORIN BOUND ESTIMATES */
+/* (NOTE: G(I) IS THE NEGATIVE OF THE */
+/* GERSCHGORIN LOWER BOUND.) */
+
+ if (a[j + j * a_dim1] != normj) {
+ temp = normj / a[j + j * a_dim1] - 1.f;
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ g[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)) *
+ temp;
+/* L150: */
+ }
+ }
+
+/* DO THE CHOLESKY UPDATE */
+
+ a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
+ tempjj = a[j + j * a_dim1];
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ a[i__ + j * a_dim1] /= tempjj;
+/* L160: */
+ }
+ i__2 = *n;
+ for (i__ = j + 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ i__3 = i__;
+ for (k = j + 1; k <= i__3; ++k) {
+ a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
+/* L170: */
+ }
+/* L180: */
+ }
+ } else {
+ fin2x2_(ndim, n, &a[a_offset], &e[1], &j, tau2, &delta, &
+ gamma);
+ }
+ }
+/* L200: */
+ }
+ return 0;
+} /* modchl_ */
+
+/* ************************************************************************ */
+/* SUBROUTINE NAME : INIT */
+
+/* PURPOSE : SET UP FOR START OF CHOLESKY FACTORIZATION */
+
+/* INPUT : N, NDIM, A, TAU1 */
+
+/* OUTPUT : PHASE1 - BOOLEAN VALUE SET TO TRUE IF IN PHASE ONE, */
+/* OTHERWISE FALSE. */
+/* DELTA - AMOUNT TO ADD TO AJJ AT ITERATION J */
+/* P,G,E - DESCRIBED ABOVE IN MODCHL */
+/* MING - THE MINIMUM NEGATIVE GERSCHGORIN BOUND */
+/* GAMMA - THE MAXIMUM DIAGONAL ELEMENT OF A */
+/* TAUGAM - TAU1 * GAMMA */
+
+/* ************************************************************************ */
+/* Subroutine */ int init_(integer *n, integer *ndim, doublereal *a, logical *
+ phase1, doublereal *delta, integer *p, doublereal *g, doublereal *e,
+ doublereal *ming, doublereal *tau1, doublereal *gamma, doublereal *
+ taugam)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1;
+ doublereal d__1, d__2, d__3;
+
+ /* Local variables */
+ integer i__;
+ extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
+ integer *, doublereal *);
+
+ /* Parameter adjustments */
+ --e;
+ --g;
+ --p;
+ a_dim1 = *ndim;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ *phase1 = TRUE_;
+ *delta = 0.f;
+ *ming = 0.f;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ p[i__] = i__;
+ g[i__] = 0.f;
+ e[i__] = 0.f;
+/* L10: */
+ }
+
+/* FIND THE MAXIMUM MAGNITUDE OF THE DIAGONAL ELEMENTS. */
+/* IF ANY DIAGONAL ELEMENT IS NEGATIVE, THEN PHASE1 IS FALSE. */
+
+ *gamma = 0.f;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__2 = *gamma, d__3 = (d__1 = a[i__ + i__ * a_dim1], abs(d__1));
+ *gamma = max(d__2,d__3);
+ if (a[i__ + i__ * a_dim1] < 0.f) {
+ *phase1 = FALSE_;
+ }
+/* L20: */
+ }
+ *taugam = *tau1 * *gamma;
+
+/* IF NOT IN PHASE1, THEN CALCULATE THE INITIAL GERSCHGORIN BOUNDS */
+/* NEEDED FOR THE START OF PHASE2. */
+
+ if (! (*phase1)) {
+ gersch_(ndim, n, &a[a_offset], &c__1, &g[1]);
+ }
+ return 0;
+} /* init_ */
+
+/* ************************************************************************ */
+
+/* SUBROUTINE NAME : GERSCH */
+
+/* PURPOSE : CALCULATE THE NEGATIVE OF THE GERSCHGORIN BOUNDS */
+/* CALLED ONCE AT THE START OF PHASE II. */
+
+/* INPUT : NDIM, N, A, J */
+
+/* OUTPUT : G - AN N VECTOR CONTAINING THE NEGATIVES OF THE */
+/* GERSCHGORIN BOUNDS. */
+
+/* ************************************************************************ */
+/* Subroutine */ int gersch_(integer *ndim, integer *n, doublereal *a,
+ integer *j, doublereal *g)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+ doublereal d__1;
+
+ /* Local variables */
+ integer i__, k;
+ doublereal offrow;
+
+ /* Parameter adjustments */
+ --g;
+ a_dim1 = *ndim;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = *j; i__ <= i__1; ++i__) {
+ offrow = 0.f;
+ i__2 = i__ - 1;
+ for (k = *j; k <= i__2; ++k) {
+ offrow += (d__1 = a[i__ + k * a_dim1], abs(d__1));
+/* L10: */
+ }
+ i__2 = *n;
+ for (k = i__ + 1; k <= i__2; ++k) {
+ offrow += (d__1 = a[k + i__ * a_dim1], abs(d__1));
+/* L20: */
+ }
+ g[i__] = offrow - a[i__ + i__ * a_dim1];
+/* L30: */
+ }
+ return 0;
+} /* gersch_ */
+
+/* ************************************************************************ */
+
+/* SUBROUTINE NAME : FIN2X2 */
+
+/* PURPOSE : HANDLES FINAL 2X2 SUBMATRIX IN PHASE II. */
+/* FINDS EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX, */
+/* CALCULATES THE AMOUNT TO ADD TO THE DIAGONAL, */
+/* ADDS TO THE FINAL 2 DIAGONAL ELEMENTS, */
+/* AND DOES THE FINAL UPDATE. */
+
+/* INPUT : NDIM, N, A, E, J, TAU2, */
+/* DELTA - AMOUNT ADDED TO THE DIAGONAL IN THE */
+/* PREVIOUS ITERATION */
+
+/* OUTPUT : A - MATRIX WITH COMPLETE L FACTOR IN THE LOWER TRIANGLE, */
+/* E - N*1 VECTOR CONTAINING THE AMOUNT ADDED TO THE DIAGONAL */
+/* AT EACH ITERATION, */
+/* DELTA - AMOUNT ADDED TO DIAGONAL ELEMENTS N-1 AND N. */
+
+/* ************************************************************************ */
+/* Subroutine */ int fin2x2_(integer *ndim, integer *n, doublereal *a,
+ doublereal *e, integer *j, doublereal *tau2, doublereal *delta,
+ doublereal *gamma)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset;
+ doublereal d__1;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ doublereal t1, t2, t3, t1a, t2a, temp, lmbd1, lmbd2, delta1, lmbdhi,
+ lmbdlo;
+
+
+/* FIND EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX */
+
+ /* Parameter adjustments */
+ --e;
+ a_dim1 = *ndim;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ t1 = a[*n - 1 + (*n - 1) * a_dim1] + a[*n + *n * a_dim1];
+ t2 = a[*n - 1 + (*n - 1) * a_dim1] - a[*n + *n * a_dim1];
+ t1a = abs(t2);
+ t2a = (d__1 = a[*n + (*n - 1) * a_dim1], abs(d__1)) * 2.;
+ if (t1a >= t2a) {
+ if (t1a > 0.) {
+ t2a /= t1a;
+ }
+/* Computing 2nd power */
+ d__1 = t2a;
+ t3 = t1a * sqrt(d__1 * d__1 + 1.);
+ } else {
+ t1a /= t2a;
+/* Computing 2nd power */
+ d__1 = t1a;
+ t3 = t2a * sqrt(d__1 * d__1 + 1.);
+ }
+ lmbd1 = (t1 - t3) / 2.f;
+ lmbd2 = (t1 + t3) / 2.f;
+ lmbdhi = max(lmbd1,lmbd2);
+ lmbdlo = min(lmbd1,lmbd2);
+
+/* FIND DELTA SUCH THAT: */
+/* 1. THE L2 CONDITION NUMBER OF THE FINAL */
+/* 2X2 SUBMATRIX + DELTA*I <= TAU2 */
+/* 2. DELTA >= PREVIOUS DELTA, */
+/* 3. LMBDLO + DELTA >= TAU2 * GAMMA, */
+/* WHERE LMBDLO IS THE SMALLEST EIGENVALUE OF THE FINAL */
+/* 2X2 SUBMATRIX */
+
+ delta1 = (lmbdhi - lmbdlo) / (1.f - *tau2);
+ delta1 = max(delta1,*gamma);
+ delta1 = *tau2 * delta1 - lmbdlo;
+ temp = 0.f;
+ *delta = max(*delta,temp);
+ *delta = max(*delta,delta1);
+ if (*delta > 0.f) {
+ a[*n - 1 + (*n - 1) * a_dim1] += *delta;
+ a[*n + *n * a_dim1] += *delta;
+ e[*n - 1] = *delta;
+ e[*n] = *delta;
+ }
+
+/* FINAL UPDATE */
+
+ a[*n - 1 + (*n - 1) * a_dim1] = sqrt(a[*n - 1 + (*n - 1) * a_dim1]);
+ a[*n + (*n - 1) * a_dim1] /= a[*n - 1 + (*n - 1) * a_dim1];
+/* Computing 2nd power */
+ d__1 = a[*n + (*n - 1) * a_dim1];
+ a[*n + *n * a_dim1] -= d__1 * d__1;
+ a[*n + *n * a_dim1] = sqrt(a[*n + *n * a_dim1]);
+ return 0;
+} /* fin2x2_ */
+
+/* ---------------------- */
+/* | S L V M D L | */
+/* ---------------------- */
+/* Subroutine */ int slvmdl_(integer *nr, integer *n, doublereal *h__,
+ doublereal *u, doublereal *t, doublereal *e, doublereal *diag,
+ doublereal *s, doublereal *g, integer *pivot, doublereal *w1,
+ doublereal *w2, doublereal *w3, doublereal *alpha, doublereal *beta,
+ logical *nomin, doublereal *eps)
+{
+ /* System generated locals */
+ integer h_dim1, h_offset, i__1, i__2;
+
+ /* Builtin functions */
+ double sqrt(doublereal);
+
+ /* Local variables */
+ integer i__, j;
+ doublereal r__, r1, r2, ca, cb, cc, cd, w11, sg, w22, w12, w33, w13, w23,
+ ss, uu, shs;
+ extern /* Subroutine */ int zhz_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, doublereal *);
+ doublereal temp;
+ extern /* Subroutine */ int sigma_(doublereal *, doublereal *, doublereal
+ *, doublereal *, doublereal *), dstar_(integer *, integer *,
+ doublereal *, doublereal *, doublereal *, doublereal *,
+ doublereal *, doublereal *, doublereal *, doublereal *);
+ doublereal addmax;
+ extern /* Subroutine */ int choldr_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *, integer *, doublereal *, doublereal *,
+ doublereal *), bakslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *);
+ doublereal sgstar;
+ extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
+ doublereal *, doublereal *), solvew_(integer *, integer *,
+ doublereal *, doublereal *, doublereal *, doublereal *);
+
+
+/* PURPOSE: */
+/* COMPUTE TENSOR AND NEWTON STEPS */
+
+/* ---------------------------------------------------------------------------- */
+
+/* PARAMETERS: */
+
+/* NR --> ROW DIMENSION OF MATRIX */
+/* N --> DIMENSION OF PROBLEM */
+/* H(N,N) --> HESSIAN */
+/* U(N) --> VECTOR TO FORM Q IN QR */
+/* T(N) --> WORKSPACE */
+/* E(N) --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION */
+/* DIAG(N) --> DIAGONAL OF HESSIAN */
+/* S(N) --> STEP TO PREVIOUS POINT (FOR TENSOR MODEL) */
+/* G(N) --> CURRENT GRADIENT */
+/* PIVOT(N) --> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION */
+/* W1(N) <--> ON INPUT: A=2*(GP-G-HS-S*BETA/(6*STS)) */
+/* ON OUTPUT: TENSOR STEP */
+/* W2(N) --> SH */
+/* W3(N) <-- NEWTON STEP */
+/* ALPHA --> SCALAR FOR 3RD ORDER TERM OF TENSOR MODEL */
+/* BETA --> SCALAR FOR 4TH ORDER TERM OF TENSOR MODEL */
+/* NOMIN <-- =.TRUE. IF TENSOR MODEL HAS NO MINIMIZER */
+/* EPS --> MACHINE EPSILON */
+
+
+/* S O L V E M O D E L */
+
+ /* Parameter adjustments */
+ --w3;
+ --w2;
+ --w1;
+ --pivot;
+ --g;
+ --s;
+ --diag;
+ --e;
+ --t;
+ --u;
+ h_dim1 = *nr;
+ h_offset = 1 + h_dim1;
+ h__ -= h_offset;
+
+ /* Function Body */
+ *nomin = FALSE_;
+
+/* COMPUTE QTHQ(N,N), ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND N-1 */
+/* COLUMNS OF QTHQ */
+ if (*n > 1) {
+ zhz_(nr, n, &s[1], &h__[h_offset], &u[1], &t[1]);
+
+/* IN CHOLESKY DECOMPOSITION WILL STORE H(1,1) ... H(N-1,N-1) */
+/* IN DIAG(1) ... DIAG(N-1), STORE H(N,N) IN DIAG(N) FIRST */
+ diag[*n] = h__[*n + *n * h_dim1];
+
+/* COLESKY DECOMPOSITION FOR FIRST N-1 ROWS AND N-1 COLUMNS OF ZTHZ */
+/* ZTHZ(N-1,N-1)=LLT */
+ i__1 = *n - 1;
+ choldr_(nr, &i__1, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
+ diag[1], &addmax);
+ }
+
+/* ON INPUT: SH IS STORED IN W2 */
+ shs = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ shs += w2[i__] * s[i__];
+ w3[i__] = g[i__];
+/* L100: */
+ }
+
+/* COMPUTE W1,W2,W3 */
+/* W1=L**-1*ZT*A */
+/* W2=L**-1*ZT*SH */
+/* W3=L**-1*ZT*G */
+
+ if (*n > 1) {
+ solvew_(nr, n, &h__[h_offset], &u[1], &w1[1], &t[1]);
+ solvew_(nr, n, &h__[h_offset], &u[1], &w2[1], &t[1]);
+ solvew_(nr, n, &h__[h_offset], &u[1], &w3[1], &t[1]);
+ }
+
+/* COMPUTE COEFFICIENTS CA,CB,CC AND CD FOR REDUCED ONE VARIABLE */
+/* 3RD ORDER EQUATION */
+ w11 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w11 += w1[i__] * w1[i__];
+/* L110: */
+ }
+ ca = *beta / 6. - w11 / 2.;
+ w12 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w12 += w1[i__] * w2[i__];
+/* L120: */
+ }
+ cb = *alpha / 2. - w12 * 3. / 2.;
+ w13 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w13 += w1[i__] * w3[i__];
+/* L130: */
+ }
+ w22 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w22 += w2[i__] * w2[i__];
+/* L133: */
+ }
+ cc = shs - w22 - w13;
+ sg = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ sg += s[i__] * g[i__];
+/* L140: */
+ }
+ w23 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w23 += w2[i__] * w3[i__];
+/* L145: */
+ }
+ cd = sg - w23;
+ w33 = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w33 += w3[i__] * w3[i__];
+/* L147: */
+ }
+
+/* COMPUTE DESIRABLE ROOT, SGSTAR, OF 3RD ORDER EQUATION */
+ if (ca != 0.) {
+ sigma_(&sgstar, &ca, &cb, &cc, &cd);
+ if (ca == 0.) {
+ *nomin = TRUE_;
+ goto L200;
+ }
+ } else {
+/* 2ND ORDER ( CA=0 ) */
+ if (cb != 0.) {
+ r__ = cc * cc - cb * 4. * cd;
+ if (r__ < 0.) {
+ *nomin = TRUE_;
+ goto L200;
+ } else {
+ r1 = (-cc + sqrt(r__)) / (cb * 2.);
+ r2 = (-cc - sqrt(r__)) / (cb * 2.);
+ if (r2 < r1) {
+ temp = r1;
+ r1 = r2;
+ r2 = temp;
+ }
+ if (cb > 0.) {
+ sgstar = r2;
+ } else {
+ sgstar = r1;
+ }
+ if (r1 > 0. && sgstar == r2 || r2 < 0. && sgstar == r1) {
+ *nomin = TRUE_;
+ goto L200;
+ }
+ }
+ } else {
+/* 1ST ORDER (CA=0,CB=0) */
+ if (cc > 0.) {
+ sgstar = -cd / cc;
+ } else {
+ *nomin = TRUE_;
+ goto L200;
+ }
+ }
+ }
+
+/* FIND TENSOR STEP, W1 (FUNCTION OF SGSTAR) */
+ dstar_(nr, n, &u[1], &s[1], &w1[1], &w2[1], &w3[1], &sgstar, &h__[
+ h_offset], &w1[1]);
+
+/* COMPUTE DN */
+L200:
+ uu = 0.;
+ ss = 0.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ uu += u[i__] * u[i__];
+ ss += s[i__] * s[i__];
+/* L202: */
+ }
+ uu /= 2.;
+ ss = sqrt(ss);
+ if (*n == 1) {
+ choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &diag[1],
+ &addmax);
+ } else {
+
+/* COMPUTE LAST ROW OF L(N,N) */
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = 0.;
+ if (i__ > 1) {
+ i__2 = i__ - 1;
+ for (j = 1; j <= i__2; ++j) {
+ temp += h__[*n + j * h_dim1] * h__[i__ + j * h_dim1];
+/* L210: */
+ }
+ }
+ h__[*n + i__ * h_dim1] = (h__[i__ + *n * h_dim1] - temp) / h__[
+ i__ + i__ * h_dim1];
+/* L220: */
+ }
+ temp = 0.;
+ i__1 = *n - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp += h__[*n + i__ * h_dim1] * h__[*n + i__ * h_dim1];
+/* L224: */
+ }
+ h__[*n + *n * h_dim1] = diag[*n] - temp + addmax;
+ if (h__[*n + *n * h_dim1] > 0.) {
+ h__[*n + *n * h_dim1] = sqrt(h__[*n + *n * h_dim1]);
+ } else {
+/* AFTER ADDING THE LAST COLUMN AND ROW */
+/* QTHQ IS NOT P.D., NEED TO REDO CHOLESKY DECOMPOSITION */
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ i__2 = i__ - 1;
+ for (j = 1; j <= i__2; ++j) {
+ h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
+/* L230: */
+ }
+ h__[i__ + i__ * h_dim1] = diag[i__];
+/* L232: */
+ }
+ h__[h_dim1 + 1] = diag[1];
+ choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
+ diag[1], &addmax);
+ }
+ }
+
+/* SOLVE QTHQ*QT*W3=-QT*G, WHERE W3 IS NEWTON STEP */
+/* W2=-QT*G */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w2[i__] = 0.;
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ w2[i__] += u[j] * u[i__] * g[j] / uu;
+/* L300: */
+ }
+ w2[i__] -= g[i__];
+/* L302: */
+ }
+ forslv_(nr, n, &h__[h_offset], &w3[1], &w2[1]);
+ bakslv_(nr, n, &h__[h_offset], &w2[1], &w3[1]);
+/* W2=QT*W3 => W3=Q*W2 --- NEWTON STEP */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ w3[i__] = 0.;
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ w3[i__] += u[i__] * u[j] * w2[j] / uu;
+/* L310: */
+ }
+ w3[i__] = w2[i__] - w3[i__];
+/* L312: */
+ }
+ return 0;
+} /* slvmdl_ */
+
+/* ---------------------- */
+/* | O P T S T P | */
+/* ---------------------- */
+/* Subroutine */ int optstp_(integer *n, doublereal *xpls, doublereal *fpls,
+ doublereal *gpls, doublereal *x, integer *itncnt, integer *icscmx,
+ integer *itrmcd, doublereal *gradtl, doublereal *steptl, doublereal *
+ fscale, integer *itnlim, integer *iretcd, logical *mxtake, integer *
+ ipr, integer *msg, doublereal *rgx, doublereal *rsx)
+{
+ /* Format strings */
+ static char fmt_900[] = "(\002 OPTSTP STEP OF MAXIMUM LENGTH (STEPMX)"
+ " TAKEN\002)";
+ static char fmt_901[] = "(\002 OPTSTP RELATIVE GRADIENT CLOSE TO ZE"
+ "RO.\002/\002 OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.\002)"
+ ;
+ static char fmt_902[] = "(\002 OPTSTP SUCCESSIVE ITERATES WITHIN TOLE"
+ "RANCE.\002/\002 OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION"
+ ".\002)";
+ static char fmt_903[] = "(\002 OPTSTP LAST GLOBAL STEP FAILED TO LOCA"
+ "TE A POINT\002,\002 LOWER THAN X.\002/\002 OPTSTP EITHER X IS"
+ " AN APPROXIMATE LOCAL MINIMUM\002,\002 OF THE FUNCTION,\002/\002"
+ " OPTSTP THE FUNCTION IS TOO NON-LINEAR FOR THIS\002,\002 ALGO"
+ "RITHM,\002/\002 OPTSTP OR STEPTL IS TOO LARGE.\002)";
+ static char fmt_904[] = "(\002 OPTSTP ITERATION LIMIT EXCEEDED.\002"
+ "/\002 OPTSTP ALGORITHM FAILED.\002)";
+ static char fmt_905[] = "(\002 OPTSTP MAXIMUM STEP SIZE EXCEEDED 5"
+ "\002,\002 CONSECUTIVE TIMES.\002/\002 OPTSTP EITHER THE FUNCT"
+ "ION IS UNBOUNDED BELOW,\002/\002 OPTSTP BECOMES ASYMPTOTIC TO"
+ " A FINITE VALUE\002,\002 FROM ABOVE IN SOME DIRECTION,\002/\002 "
+ "OPTSTP OR STEPMX IS TOO SMALL\002)";
+
+ /* System generated locals */
+ integer i__1;
+ doublereal d__1, d__2, d__3;
+
+ /* Builtin functions */
+ integer s_wsfe(cilist *), e_wsfe(void);
+
+ /* Local variables */
+ doublereal d__;
+ integer i__, imsg;
+ doublereal relgrd;
+ integer jtrmcd;
+ doublereal relstp;
+
+ /* Fortran I/O blocks */
+ static cilist io___280 = { 0, 0, 0, fmt_900, 0 };
+ static cilist io___281 = { 0, 0, 0, fmt_901, 0 };
+ static cilist io___282 = { 0, 0, 0, fmt_902, 0 };
+ static cilist io___283 = { 0, 0, 0, fmt_903, 0 };
+ static cilist io___284 = { 0, 0, 0, fmt_904, 0 };
+ static cilist io___285 = { 0, 0, 0, fmt_905, 0 };
+
+
+
+/* UNCONSTRAINED MINIMIZATION STOPPING CRITERIA */
+/* -------------------------------------------- */
+/* FIND WHETHER THE ALGORITHM SHOULD TERMINATE, DUE TO ANY */
+/* OF THE FOLLOWING: */
+/* 1) PROBLEM SOLVED WITHIN USER TOLERANCE */
+/* 2) CONVERGENCE WITHIN USER TOLERANCE */
+/* 3) ITERATION LIMIT REACHED */
+/* 4) DIVERGENCE OR TOO RESTRICTIVE MAXIMUM STEP (STEPMX) SUSPECTED */
+
+
+/* PARAMETERS */
+/* ---------- */
+/* N --> DIMENSION OF PROBLEM */
+/* XPLS(N) --> NEW ITERATE X[K] */
+/* FPLS --> FUNCTION VALUE AT NEW ITERATE F(XPLS) */
+/* GPLS(N) --> GRADIENT AT NEW ITERATE, G(XPLS), OR APPROXIMATE */
+/* X(N) --> OLD ITERATE X[K-1] */
+/* ITNCNT --> CURRENT ITERATION K */
+/* ICSCMX <--> NUMBER CONSECUTIVE STEPS .GE. STEPMX */
+/* [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] */
+/* ITRMCD <-- TERMINATION CODE */
+/* GRADTL --> TOLERANCE AT WHICH RELATIVE GRADIENT CONSIDERED CLOSE */
+/* ENOUGH TO ZERO TO TERMINATE ALGORITHM */
+/* STEPTL --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES */
+/* CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM */
+/* FSCALE --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION */
+/* ITNLIM --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS */
+/* IRETCD --> RETURN CODE */
+/* MXTAKE --> BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED */
+/* IPR --> DEVICE TO WHICH TO SEND OUTPUT */
+/* MSG <--> CONTROL OUTPUT ON INPUT AND CONTAIN STOPPING */
+/* CONDITION ON OUTPUT */
+
+
+
+ /* Parameter adjustments */
+ --x;
+ --gpls;
+ --xpls;
+
+ /* Function Body */
+ *itrmcd = 0;
+ imsg = *msg;
+ *rgx = 0.;
+ *rsx = 0.;
+
+/* LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X */
+ if (*iretcd != 1) {
+ goto L50;
+ }
+/* IF(IRETCD.EQ.1) */
+/* THEN */
+ jtrmcd = 3;
+ goto L600;
+/* ENDIF */
+L50:
+
+/* FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM. */
+/* CHECK WHETHER WITHIN TOLERANCE */
+
+/* Computing MAX */
+ d__1 = abs(*fpls);
+ d__ = max(d__1,*fscale);
+/* D=1D0 */
+ *rgx = 0.f;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__3 = (d__2 = xpls[i__], abs(d__2));
+ relgrd = (d__1 = gpls[i__], abs(d__1)) * max(d__3,1.) / d__;
+ *rgx = max(*rgx,relgrd);
+/* L100: */
+ }
+ jtrmcd = 1;
+ if (*rgx <= *gradtl) {
+ goto L600;
+ }
+
+ if (*itncnt == 0) {
+ return 0;
+ }
+
+/* FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM */
+/* CHECK WHETHER WITHIN TOLERANCE. */
+
+ *rsx = 0.f;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* Computing MAX */
+ d__3 = (d__2 = xpls[i__], abs(d__2));
+ relstp = (d__1 = xpls[i__] - x[i__], abs(d__1)) / max(d__3,1.);
+ *rsx = max(*rsx,relstp);
+/* L120: */
+ }
+ jtrmcd = 2;
+ if (*rsx <= *steptl) {
+ goto L600;
+ }
+
+/* CHECK ITERATION LIMIT */
+
+ jtrmcd = 4;
+ if (*itncnt >= *itnlim) {
+ goto L600;
+ }
+
+/* CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX */
+
+ if (*mxtake) {
+ goto L140;
+ }
+/* IF(.NOT.MXTAKE) */
+/* THEN */
+ *icscmx = 0;
+ return 0;
+/* ELSE */
+L140:
+/* IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900) */
+ if (imsg >= 1) {
+ io___280.ciunit = *ipr;
+ s_wsfe(&io___280);
+ e_wsfe();
+ }
+ ++(*icscmx);
+ if (*icscmx < 5) {
+ return 0;
+ }
+ jtrmcd = 5;
+/* ENDIF */
+
+
+/* PRINT TERMINATION CODE */
+
+L600:
+ *itrmcd = jtrmcd;
+/* IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD */
+ if (imsg >= 1) {
+ switch (*itrmcd) {
+ case 1: goto L601;
+ case 2: goto L602;
+ case 3: goto L603;
+ case 4: goto L604;
+ case 5: goto L605;
+ }
+ }
+ goto L700;
+L601:
+ io___281.ciunit = *ipr;
+ s_wsfe(&io___281);
+ e_wsfe();
+ goto L700;
+L602:
+ io___282.ciunit = *ipr;
+ s_wsfe(&io___282);
+ e_wsfe();
+ goto L700;
+L603:
+ io___283.ciunit = *ipr;
+ s_wsfe(&io___283);
+ e_wsfe();
+ goto L700;
+L604:
+ io___284.ciunit = *ipr;
+ s_wsfe(&io___284);
+ e_wsfe();
+ goto L700;
+L605:
+ io___285.ciunit = *ipr;
+ s_wsfe(&io___285);
+ e_wsfe();
+
+L700:
+ *msg = -(*itrmcd);
+ return 0;
+
+} /* optstp_ */
+
+
+/* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,
+ doublereal *dy, integer *incy)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ integer i__, m, ix, iy, mp1;
+
+
+/* COPIES A VECTOR, X, TO A VECTOR, Y. */
+/* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
+/* JACK DONGARRA, LINPACK, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dy;
+ --dx;
+
+ /* Function Body */
+ if (*n <= 0) {
+ return 0;
+ }
+ if (*incx == 1 && *incy == 1) {
+ goto L20;
+ }
+
+/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
+/* NOT EQUAL TO 1 */
+
+ ix = 1;
+ iy = 1;
+ if (*incx < 0) {
+ ix = (-(*n) + 1) * *incx + 1;
+ }
+ if (*incy < 0) {
+ iy = (-(*n) + 1) * *incy + 1;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dy[iy] = dx[ix];
+ ix += *incx;
+ iy += *incy;
+/* L10: */
+ }
+ return 0;
+
+/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
+
+
+/* CLEAN-UP LOOP */
+
+L20:
+ m = *n % 7;
+ if (m == 0) {
+ goto L40;
+ }
+ i__1 = m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dy[i__] = dx[i__];
+/* L30: */
+ }
+ if (*n < 7) {
+ return 0;
+ }
+L40:
+ mp1 = m + 1;
+ i__1 = *n;
+ for (i__ = mp1; i__ <= i__1; i__ += 7) {
+ dy[i__] = dx[i__];
+ dy[i__ + 1] = dx[i__ + 1];
+ dy[i__ + 2] = dx[i__ + 2];
+ dy[i__ + 3] = dx[i__ + 3];
+ dy[i__ + 4] = dx[i__ + 4];
+ dy[i__ + 5] = dx[i__ + 5];
+ dy[i__ + 6] = dx[i__ + 6];
+/* L50: */
+ }
+ return 0;
+} /* dcopy_ */
+
+
+doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
+ integer *incy)
+{
+ /* System generated locals */
+ integer i__1;
+ doublereal ret_val;
+
+ /* Local variables */
+ integer i__, m, ix, iy, mp1;
+ doublereal dtemp;
+
+
+/* FORMS THE DOT PRODUCT OF TWO VECTORS. */
+/* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
+/* JACK DONGARRA, LINPACK, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dy;
+ --dx;
+
+ /* Function Body */
+ ret_val = 0.;
+ dtemp = 0.;
+ if (*n <= 0) {
+ return ret_val;
+ }
+ if (*incx == 1 && *incy == 1) {
+ goto L20;
+ }
+
+/* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
+/* NOT EQUAL TO 1 */
+
+ ix = 1;
+ iy = 1;
+ if (*incx < 0) {
+ ix = (-(*n) + 1) * *incx + 1;
+ }
+ if (*incy < 0) {
+ iy = (-(*n) + 1) * *incy + 1;
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dtemp += dx[ix] * dy[iy];
+ ix += *incx;
+ iy += *incy;
+/* L10: */
+ }
+ ret_val = dtemp;
+ return ret_val;
+
+/* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
+
+
+/* CLEAN-UP LOOP */
+
+L20:
+ m = *n % 5;
+ if (m == 0) {
+ goto L40;
+ }
+ i__1 = m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ dtemp += dx[i__] * dy[i__];
+/* L30: */
+ }
+ if (*n < 5) {
+ goto L60;
+ }
+L40:
+ mp1 = m + 1;
+ i__1 = *n;
+ for (i__ = mp1; i__ <= i__1; i__ += 5) {
+ dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
+ i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
+ 4] * dy[i__ + 4];
+/* L50: */
+ }
+L60:
+ ret_val = dtemp;
+ return ret_val;
+} /* ddot_ */
+
+
+/* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,
+ integer *incx)
+{
+ /* System generated locals */
+ integer i__1, i__2;
+
+ /* Local variables */
+ integer i__, m, mp1, nincx;
+
+
+/* SCALES A VECTOR BY A CONSTANT. */
+/* USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
+/* JACK DONGARRA, LINPACK, 3/11/78. */
+
+
+ /* Parameter adjustments */
+ --dx;
+
+ /* Function Body */
+ if (*n <= 0) {
+ return 0;
+ }
+ if (*incx == 1) {
+ goto L20;
+ }
+
+/* CODE FOR INCREMENT NOT EQUAL TO 1 */
+
+ nincx = *n * *incx;
+ i__1 = nincx;
+ i__2 = *incx;
+ for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
+ dx[i__] = *da * dx[i__];
+/* L10: */
+ }
+ return 0;
+
+/* CODE FOR INCREMENT EQUAL TO 1 */
+
+
+/* CLEAN-UP LOOP */
+
+L20:
+ m = *n % 5;
+ if (m == 0) {
+ goto L40;
+ }
+ i__2 = m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ dx[i__] = *da * dx[i__];
+/* L30: */
+ }
+ if (*n < 5) {
+ return 0;
+ }
+L40:
+ mp1 = m + 1;
+ i__2 = *n;
+ for (i__ = mp1; i__ <= i__2; i__ += 5) {
+ dx[i__] = *da * dx[i__];
+ dx[i__ + 1] = *da * dx[i__ + 1];
+ dx[i__ + 2] = *da * dx[i__ + 2];
+ dx[i__ + 3] = *da * dx[i__ + 3];
+ dx[i__ + 4] = *da * dx[i__ + 4];
+/* L50: */
+ }
+ return 0;
+} /* dscal_ */
+
+/* Main program alias */ int driver_ () { MAIN__ (); return 0; }
~ianmdlvl / nlopt.git / commitdiff