tensor/tensor.c

   1 /* tensor.f -- translated by f2c (version 20050501).
   2    You must link the resulting object file with libf2c:
   3         on Microsoft Windows system, link with libf2c.lib;
   4         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
   5         or, if you install libf2c.a in a standard place, with -lf2c -lm
   6         -- in that order, at the end of the command line, as in
   7                 cc *.o -lf2c -lm
   8         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
   9
  10                 http://www.netlib.org/f2c/libf2c.zip
  11 */
  12
  13 #include "f2c.h"
  14
  15 /* Table of constant values */
  16
  17 static integer c__5 = 5;
  18 static integer c__1 = 1;
  19 static integer c__9 = 9;
  20 static integer c__3 = 3;
  21 static doublereal c_b111 = 2.;
  22 static doublereal c_b134 = 10.;
  23 static doublereal c_b250 = .33333333333333331;
  24 static doublereal c_b324 = 1.;
  25 static doublereal c_b384 = 3.;
  26
  27 /*      ALGORITHM 739, COLLECTED ALGORITHMS FROM ACM. */
  28 /*      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, */
  29 /*      VOL. 20, NO. 4, DECEMBER, 1994, P.518-530. */
  30 /* *** driver.f */
  31 /* Main program */ int MAIN__(void)
  32 {
  33     /* Format strings */
  34     static char fmt_211[] = "(\002    G=\002,999e20.13)";
  35
  36     /* System generated locals */
  37     integer i__1;
  38     olist o__1;
  39     cllist cl__1;
  40     alist al__1, al__2;
  41
  42     /* Builtin functions */
  43     integer f_open(olist *), s_rsle(cilist *), do_lio(integer *, integer *,
  44             char *, ftnlen), e_rsle(void), s_wsle(cilist *), e_wsle(void),
  45             s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void),
  46              f_rew(alist *), f_end(alist *), f_clos(cllist *);
  47     /* Subroutine */ int s_stop(char *, ftnlen);
  48
  49     /* Local variables */
  50     doublereal h__[200] /* was [50][4] */;
  51     integer i__, n;
  52     doublereal x[50];
  53     integer nr;
  54     extern /* Subroutine */ int fcn_(), grd_();
  55     integer msg;
  56     extern /* Subroutine */ int hsn_();
  57     integer ipr;
  58     doublereal wrk[400] /* was [50][8] */, fpls, gpls[50];
  59     integer iwrk[4];
  60     doublereal xpls[50], typx[50];
  61     integer itnno;
  62     doublereal fscale;
  63     integer grdflg, hesflg;
  64     doublereal gradtl;
  65     integer ndigit;
  66     extern /* Subroutine */ int dfault_(integer *, doublereal *, doublereal *,
  67              doublereal *, doublereal *, integer *, doublereal *, integer *,
  68             integer *, integer *, integer *, integer *, integer *);
  69     integer method;
  70     extern /* Subroutine */ int prtfcn_(integer *);
  71     integer itnlim;
  72     extern /* Subroutine */ int tensrd_(integer *, integer *, doublereal *,
  73             U_fp, integer *, doublereal *, doublereal *, doublereal *,
  74             doublereal *, integer *, doublereal *, integer *), tensor_(
  75             integer *, integer *, doublereal *, U_fp, U_fp, U_fp, doublereal *
  76             , doublereal *, doublereal *, doublereal *, integer *, doublereal
  77             *, integer *, integer *, integer *, integer *, integer *, integer
  78             *, doublereal *, doublereal *, doublereal *, doublereal *,
  79             integer *, doublereal *, integer *);
  80     doublereal steptl, stepmx;
  81
  82     /* Fortran I/O blocks */
  83     static cilist io___3 = { 0, 5, 0, 0, 0 };
  84     static cilist io___6 = { 0, 6, 0, 0, 0 };
  85     static cilist io___7 = { 0, 6, 0, 0, 0 };
  86     static cilist io___8 = { 0, 6, 0, 0, 0 };
  87     static cilist io___9 = { 0, 6, 0, 0, 0 };
  88     static cilist io___10 = { 0, 6, 0, 0, 0 };
  89     static cilist io___11 = { 0, 6, 0, 0, 0 };
  90     static cilist io___12 = { 0, 6, 0, 0, 0 };
  91     static cilist io___21 = { 0, 6, 0, 0, 0 };
  92     static cilist io___22 = { 0, 6, 0, fmt_211, 0 };
  93     static cilist io___23 = { 0, 6, 0, 0, 0 };
  94     static cilist io___24 = { 0, 6, 0, 0, 0 };
  95     static cilist io___25 = { 0, 5, 0, 0, 0 };
  96     static cilist io___26 = { 0, 6, 0, 0, 0 };
  97     static cilist io___27 = { 0, 6, 0, 0, 0 };
  98     static cilist io___28 = { 0, 6, 0, 0, 0 };
  99     static cilist io___29 = { 0, 6, 0, 0, 0 };
 100     static cilist io___30 = { 0, 6, 0, 0, 0 };
 101     static cilist io___31 = { 0, 6, 0, 0, 0 };
 102     static cilist io___32 = { 0, 6, 0, 0, 0 };
 103     static cilist io___33 = { 0, 6, 0, 0, 0 };
 104     static cilist io___45 = { 0, 6, 0, 0, 0 };
 105     static cilist io___46 = { 0, 6, 0, fmt_211, 0 };
 106     static cilist io___47 = { 0, 6, 0, 0, 0 };
 107     static cilist io___48 = { 0, 6, 0, 0, 0 };
 108
 109
 110     nr = 50;
 111     o__1.oerr = 0;
 112     o__1.ounit = 5;
 113     o__1.ofnmlen = 8;
 114     o__1.ofnm = "wood.inp";
 115     o__1.orl = 0;
 116     o__1.osta = 0;
 117     o__1.oacc = 0;
 118     o__1.ofm = 0;
 119     o__1.oblnk = 0;
 120     f_open(&o__1);
 121     prtfcn_(&n);
 122     s_rsle(&io___3);
 123     i__1 = n;
 124     for (i__ = 1; i__ <= i__1; ++i__) {
 125         do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
 126     }
 127     e_rsle();
 128     s_wsle(&io___6);
 129     do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
 130             41);
 131     e_wsle();
 132     s_wsle(&io___7);
 133     i__1 = n;
 134     for (i__ = 1; i__ <= i__1; ++i__) {
 135         do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
 136     }
 137     e_wsle();
 138 /*  SHORT FORM */
 139     s_wsle(&io___8);
 140     do_lio(&c__9, &c__1, " ", (ftnlen)1);
 141     e_wsle();
 142     s_wsle(&io___9);
 143     do_lio(&c__9, &c__1, "CALLING TENSRD - ALL INPUT PARAMETERS  ARE SET", (
 144             ftnlen)46);
 145     e_wsle();
 146     s_wsle(&io___10);
 147     do_lio(&c__9, &c__1, "                 TO THEIR DEFAULT VALUES.", (ftnlen)
 148             41);
 149     e_wsle();
 150     s_wsle(&io___11);
 151     do_lio(&c__9, &c__1, " OUTPUT WILL BE WRITTEN TO THE STANDARD OUTPUT.", (
 152             ftnlen)47);
 153     e_wsle();
 154     s_wsle(&io___12);
 155     do_lio(&c__9, &c__1, " ", (ftnlen)1);
 156     e_wsle();
 157     tensrd_(&nr, &n, x, (U_fp)fcn_, &msg, xpls, &fpls, gpls, h__, &itnno, wrk,
 158              iwrk);
 159     s_wsle(&io___21);
 160     do_lio(&c__9, &c__1, "RESULTS OF TENSRD:", (ftnlen)18);
 161     e_wsle();
 162     s_wsfe(&io___22);
 163     i__1 = n;
 164     for (i__ = 1; i__ <= i__1; ++i__) {
 165         do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
 166     }
 167     e_wsfe();
 168     s_wsle(&io___23);
 169     do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
 170     i__1 = n;
 171     for (i__ = 1; i__ <= i__1; ++i__) {
 172         do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
 173                 doublereal));
 174     }
 175     e_wsle();
 176     s_wsle(&io___24);
 177     do_lio(&c__9, &c__1, "FPLS=", (ftnlen)5);
 178     do_lio(&c__5, &c__1, (char *)&fpls, (ftnlen)sizeof(doublereal));
 179     do_lio(&c__9, &c__1, "  ITNNO=", (ftnlen)8);
 180     do_lio(&c__3, &c__1, (char *)&itnno, (ftnlen)sizeof(integer));
 181     do_lio(&c__9, &c__1, "  MSG=", (ftnlen)6);
 182     do_lio(&c__3, &c__1, (char *)&msg, (ftnlen)sizeof(integer));
 183     e_wsle();
 184 /*  LONG FORM */
 185     prtfcn_(&n);
 186     al__1.aerr = 0;
 187     al__1.aunit = 5;
 188     f_rew(&al__1);
 189     s_rsle(&io___25);
 190     i__1 = n;
 191     for (i__ = 1; i__ <= i__1; ++i__) {
 192         do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
 193     }
 194     e_rsle();
 195     s_wsle(&io___26);
 196     do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
 197             41);
 198     e_wsle();
 199     s_wsle(&io___27);
 200     i__1 = n;
 201     for (i__ = 1; i__ <= i__1; ++i__) {
 202         do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
 203     }
 204     e_wsle();
 205     s_wsle(&io___28);
 206     do_lio(&c__9, &c__1, " ", (ftnlen)1);
 207     e_wsle();
 208     s_wsle(&io___29);
 209     do_lio(&c__9, &c__1, "CALLING TENSOR AFTER RESETTING INPUT PARAMETERS", (
 210             ftnlen)47);
 211     e_wsle();
 212     s_wsle(&io___30);
 213     do_lio(&c__9, &c__1, "IPR, MSG AND METHOD.", (ftnlen)20);
 214     e_wsle();
 215     s_wsle(&io___31);
 216     do_lio(&c__9, &c__1, "THE STANDARD METHOD IS USED IN THIS EXAMPLE.", (
 217             ftnlen)44);
 218     e_wsle();
 219     s_wsle(&io___32);
 220     do_lio(&c__9, &c__1, "ADDITIONAL OUTPUT WILL BE WRITTEN TO FILE 'FORT8'.",
 221              (ftnlen)50);
 222     e_wsle();
 223     s_wsle(&io___33);
 224     do_lio(&c__9, &c__1, " ", (ftnlen)1);
 225     e_wsle();
 226 /* L997: */
 227     dfault_(&n, typx, &fscale, &gradtl, &steptl, &itnlim, &stepmx, &ipr, &
 228             method, &grdflg, &hesflg, &ndigit, &msg);
 229     o__1.oerr = 0;
 230     o__1.ounit = 8;
 231     o__1.ofnmlen = 5;
 232     o__1.ofnm = "FORT8";
 233     o__1.orl = 0;
 234     o__1.osta = 0;
 235     o__1.oacc = 0;
 236     o__1.ofm = 0;
 237     o__1.oblnk = 0;
 238     f_open(&o__1);
 239     ipr = 8;
 240     msg = 2;
 241     method = 0;
 242     tensor_(&nr, &n, x, (U_fp)fcn_, (U_fp)grd_, (U_fp)hsn_, typx, &fscale, &
 243             gradtl, &steptl, &itnlim, &stepmx, &ipr, &method, &grdflg, &
 244             hesflg, &ndigit, &msg, xpls, &fpls, gpls, h__, &itnno, wrk, iwrk);
 245     s_wsle(&io___45);
 246     do_lio(&c__9, &c__1, "RESULTS OF TENSOR:", (ftnlen)18);
 247     e_wsle();
 248     s_wsfe(&io___46);
 249     i__1 = n;
 250     for (i__ = 1; i__ <= i__1; ++i__) {
 251         do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
 252     }
 253     e_wsfe();
 254     s_wsle(&io___47);
 255     do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
 256     i__1 = n;
 257     for (i__ = 1; i__ <= i__1; ++i__) {
 258         do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
 259                 doublereal));
 260     }
 261     e_wsle();
 262     s_wsle(&io___48);
 263     do_lio(&c__9, &c__1, "FPLS=", (ftnlen)5);
 264     do_lio(&c__5, &c__1, (char *)&fpls, (ftnlen)sizeof(doublereal));
 265     do_lio(&c__9, &c__1, "  ITNNO=", (ftnlen)8);
 266     do_lio(&c__3, &c__1, (char *)&itnno, (ftnlen)sizeof(integer));
 267     do_lio(&c__9, &c__1, "  MSG=", (ftnlen)6);
 268     do_lio(&c__3, &c__1, (char *)&msg, (ftnlen)sizeof(integer));
 269     e_wsle();
 270     al__2.aerr = 0;
 271     al__2.aunit = 8;
 272     f_end(&al__2);
 273     cl__1.cerr = 0;
 274     cl__1.cunit = 8;
 275     cl__1.csta = 0;
 276     f_clos(&cl__1);
 277     cl__1.cerr = 0;
 278     cl__1.cunit = 5;
 279     cl__1.csta = 0;
 280     f_clos(&cl__1);
 281     s_stop("", (ftnlen)0);
 282     return 0;
 283 } /* MAIN__ */
 284
 285 /* Subroutine */ int fcn_(integer *n, doublereal *x, doublereal *f)
 286 {
 287     /* System generated locals */
 288     doublereal d__1, d__2, d__3, d__4, d__5, d__6;
 289
 290     /* Builtin functions */
 291     double pow_dd(doublereal *, doublereal *);
 292
 293     /* Parameter adjustments */
 294     --x;
 295
 296     /* Function Body */
 297     if (*n <= 0) {
 298         *n = 4;
 299     } else {
 300         d__1 = x[1] * x[1] - x[2];
 301         d__2 = 1. - x[1];
 302         d__3 = x[3] * x[3] - x[4];
 303         d__4 = 1. - x[3];
 304         d__5 = 1. - x[2];
 305         d__6 = 1. - x[4];
 306         *f = pow_dd(&d__1, &c_b111) * 100. + pow_dd(&d__2, &c_b111) + pow_dd(&
 307                 d__3, &c_b111) * 90. + pow_dd(&d__4, &c_b111) + (pow_dd(&d__5,
 308                  &c_b111) + pow_dd(&d__6, &c_b111)) * 10.1 + (1. - x[2]) *
 309                 19.8 * (1. - x[4]);
 310     }
 311     return 0;
 312 } /* fcn_ */
 313
 314 /* Subroutine */ int prtfcn_(integer *n)
 315 {
 316     /* Builtin functions */
 317     integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
 318             e_wsle(void);
 319
 320     /* Fortran I/O blocks */
 321     static cilist io___49 = { 0, 6, 0, 0, 0 };
 322     static cilist io___50 = { 0, 6, 0, 0, 0 };
 323     static cilist io___51 = { 0, 6, 0, 0, 0 };
 324
 325
 326     *n = 4;
 327     s_wsle(&io___49);
 328     do_lio(&c__9, &c__1, "__________________________________________", (
 329             ftnlen)42);
 330     e_wsle();
 331     s_wsle(&io___50);
 332     do_lio(&c__9, &c__1, "f(x)= Wood function", (ftnlen)19);
 333     e_wsle();
 334     s_wsle(&io___51);
 335     do_lio(&c__9, &c__1, "__________________________________________", (
 336             ftnlen)42);
 337     e_wsle();
 338     return 0;
 339 } /* prtfcn_ */
 340
 341 /*  ---------------------- */
 342 /*  |  G R D              | */
 343 /*  ---------------------- */
 344 /* Subroutine */ int grd_(integer *n, real *x, real *g)
 345 {
 346 /*     DUMMY ROUTINE */
 347     return 0;
 348 } /* grd_ */
 349
 350 /*  ---------------------- */
 351 /*  |  H S N              | */
 352 /*  ---------------------- */
 353 /* Subroutine */ int hsn_(integer *nr, integer *n, real *x, real *h__)
 354 {
 355 /*     DUMMY ROUTINE */
 356     return 0;
 357 } /* hsn_ */
 358
 359 /* *** tensrd.f */
 360 /*  ---------------------- */
 361 /*  |  T E N S O R       | */
 362 /*  ---------------------- */
 363 /* Subroutine */ int tensor_(integer *nr, integer *n, doublereal *x, U_fp fcn,
 364          U_fp grd, U_fp hsn, doublereal *typx, doublereal *fscale, doublereal
 365         *gradtl, doublereal *steptl, integer *itnlim, doublereal *stepmx,
 366         integer *ipr, integer *method, integer *grdflg, integer *hesflg,
 367         integer *ndigit, integer *msg, doublereal *xpls, doublereal *fpls,
 368         doublereal *gpls, doublereal *h__, integer *itnno, doublereal *wrk,
 369         integer *iwrk)
 370 {
 371     /* System generated locals */
 372     integer h_dim1, h_offset, wrk_dim1, wrk_offset;
 373
 374     /* Local variables */
 375     extern /* Subroutine */ int opt_(integer *, integer *, doublereal *, U_fp,
 376              U_fp, U_fp, doublereal *, doublereal *, doublereal *, doublereal
 377             *, integer *, doublereal *, integer *, integer *, integer *,
 378             integer *, integer *, integer *, doublereal *, doublereal *,
 379             doublereal *, doublereal *, integer *, doublereal *, doublereal *,
 380              doublereal *, doublereal *, doublereal *, doublereal *,
 381             doublereal *, integer *);
 382
 383
 384 /* AUTHORS:   TA-TUNG CHOW, ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
 385
 386 /* DATE:      MARCH 29, 1991 */
 387
 388 /* PURPOSE: */
 389 /*   DERIVATIVE TENSOR METHOD FOR UNCONSTRAINED OPTIMIZATION */
 390 /*     THE METHOD BASES EACH ITERATION ON A SPECIALLY CONSTRUCTED */
 391 /*     FOURTH ORDER MODEL OF THE OBJECTIVE FUNCTION.  THE MODEL */
 392 /*     INTERPOLATES THE FUNCTION VALUE AND GRADIENT FROM THE PREVIOUS */
 393 /*     ITERATE AND THE CURRENT FUNCTION VALUE, GRADIENT AND HESSIAN */
 394 /*     MATRIX. */
 395
 396 /* BLAS SUBROUTINES: DCOPY,DDOT,DSCAL */
 397 /* UNCMIN SUBROUTINES: DFAULT, OPTCHK, GRDCHK, HESCHK, LNSRCH, FSTOFD, */
 398 /*                     SNDOFD, BAKSLV, FORSLV, OPTSTP */
 399 /* MODCHL SUBROUTINES: MODCHL, INIT, GERCH, FIN2X2 */
 400
 401 /* ----------------------------------------------------------------------- */
 402
 403 /* PARAMETERS: */
 404
 405 /*   NR      --> ROW DIMENSION OF MATRIX */
 406 /*   N       --> DIMENSION OF PROBLEM */
 407 /*   X(N)    --> INITIAL GUESS (INPUT) AND CURRENT POINT */
 408 /*   FCN     --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
 409 /*   GRD     --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT */
 410 /*   HSN     --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN */
 411 /*               HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
 412 /*               AND DIAGONAL OF H */
 413 /*   TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X */
 414 /*   FSCALE  --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
 415 /*   GRADTL  --> GRADIENT TOLERANCE */
 416 /*   STEPTL  --> STEP TOLERANCE */
 417 /*   ITNLIM  --> ITERATION LIMIT */
 418 /*   STEPMX  --> MAXIMUM STEP LENGTH ALLOWED */
 419 /*   IPR     --> OUTPUT UNIT NUMBER */
 420 /*   METHOD  --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT */
 421 /*               EACH ITERATION */
 422 /*               IF VALUE IS 1 THEN TRY BOTH TENSOR AND NEWTON */
 423 /*               STEPS AT EACH ITERATION */
 424 /*   GRDFLG  --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED */
 425 /*   HESFLG  --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED */
 426 /*   NDIGIT  --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
 427 /*   MSG     --> OUTPUT MESSAGE CONTROL */
 428 /*   XPLS(N) <--  NEW POINT AND FINAL POINT (OUTPUT) */
 429 /*   FPLS    <--  NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
 430 /*   GPLS(N) <--  CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) */
 431 /*   H(N,N)  --> HESSIAN */
 432 /*   ITNNO   <--  NUMBER OF ITERATIONS */
 433 /*   WRK(N,8)--> WORK SPACE */
 434 /*   IWRK(N) --> WORK SPACE */
 435
 436
 437 /*   EQUIVALENCE WRK(N,1) = G(N) */
 438 /*               WRK(N,2) = S(N) */
 439 /*               WRK(N,3) = D(N) */
 440 /*               WRK(N,4) = DN(N) */
 441 /*               WRK(N,6) = E(N) */
 442 /*               WRK(N,7) = WK1(N) */
 443 /*               WRK(N,8) = WK2(N) */
 444
 445     /* Parameter adjustments */
 446     wrk_dim1 = *nr;
 447     wrk_offset = 1 + wrk_dim1;
 448     wrk -= wrk_offset;
 449     --iwrk;
 450     h_dim1 = *nr;
 451     h_offset = 1 + h_dim1;
 452     h__ -= h_offset;
 453     --gpls;
 454     --xpls;
 455     --typx;
 456     --x;
 457
 458     /* Function Body */
 459     opt_(nr, n, &x[1], (U_fp)fcn, (U_fp)grd, (U_fp)hsn, &typx[1], fscale,
 460             gradtl, steptl, itnlim, stepmx, ipr, method, grdflg, hesflg,
 461             ndigit, msg, &xpls[1], fpls, &gpls[1], &h__[h_offset], itnno, &
 462             wrk[wrk_dim1 + 1], &wrk[(wrk_dim1 << 1) + 1], &wrk[wrk_dim1 * 3 +
 463             1], &wrk[(wrk_dim1 << 2) + 1], &wrk[wrk_dim1 * 6 + 1], &wrk[
 464             wrk_dim1 * 7 + 1], &wrk[(wrk_dim1 << 3) + 1], &iwrk[1]);
 465     return 0;
 466 } /* tensor_ */
 467
 468 /*  ---------------------- */
 469 /*  |  T E N S R D       | */
 470 /*  ---------------------- */
 471 /* Subroutine */ int tensrd_(integer *nr, integer *n, doublereal *x, U_fp fcn,
 472          integer *msg, doublereal *xpls, doublereal *fpls, doublereal *gpls,
 473         doublereal *h__, integer *itnno, doublereal *wrk, integer *iwrk)
 474 {
 475     /* System generated locals */
 476     integer h_dim1, h_offset, wrk_dim1, wrk_offset;
 477
 478     /* Local variables */
 479     extern /* Subroutine */ int grd_(integer *, real *, real *), hsn_(integer
 480             *, integer *, real *, real *);
 481     integer ipr;
 482     doublereal fscale;
 483     integer grdflg, hesflg;
 484     doublereal gradtl;
 485     integer ndigit;
 486     extern /* Subroutine */ int dfault_(integer *, doublereal *, doublereal *,
 487              doublereal *, doublereal *, integer *, doublereal *, integer *,
 488             integer *, integer *, integer *, integer *, integer *);
 489     integer method, itnlim;
 490     extern /* Subroutine */ int tensor_(integer *, integer *, doublereal *,
 491             U_fp, S_fp, S_fp, doublereal *, doublereal *, doublereal *,
 492             doublereal *, integer *, doublereal *, integer *, integer *,
 493             integer *, integer *, integer *, integer *, doublereal *,
 494             doublereal *, doublereal *, doublereal *, integer *, doublereal *,
 495              integer *);
 496     doublereal steptl, stepmx;
 497
 498
 499 /* PURPOSE: */
 500 /*   SHORT CALLING SEQUENCE SUBROUTINE */
 501
 502 /* ----------------------------------------------------------------------- */
 503
 504 /* PARAMETERS: */
 505
 506 /*   NR         --> ROW DIMENSION OF MATRIX */
 507 /*   N          --> DIMENSION OF PROBLEM */
 508 /*   X(N)       --> INITIAL GUESS (INPUT) AND CURRENT POINT */
 509 /*   FCN        --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
 510 /*   MSG        --> OUTPUT MESSAGE CONTROL */
 511 /*   XPLS(N)    <--  NEW POINT AND FINAL POINT (OUTPUT) */
 512 /*   FPLS       <--  NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
 513 /*   GPLS(N)    <--  GRADIENT AT FINAL POINT (OUTPUT) */
 514 /*   H(N,N)     --> HESSIAN */
 515 /*   ITNNO      <--  NUMBER OF ITERATIONS */
 516 /*   WRK(N,8)   --> WORK SPACE */
 517 /*   IWRK(N)    --> WORK SPACE */
 518
 519
 520 /*   EQUIVALENCE WRK(N,1) = G(N) */
 521 /*               WRK(N,2) = S(N) */
 522 /*               WRK(N,3) = D(N) */
 523 /*               WRK(N,4) = DN(N) */
 524 /*               WRK(N,5) = TYPX(N) */
 525 /*               WRK(N,6) = E(N) */
 526 /*               WRK(N,7) = WK1(N) */
 527 /*               WRK(N,8) = WK2(N) */
 528
 529     /* Parameter adjustments */
 530     wrk_dim1 = *nr;
 531     wrk_offset = 1 + wrk_dim1;
 532     wrk -= wrk_offset;
 533     --iwrk;
 534     h_dim1 = *nr;
 535     h_offset = 1 + h_dim1;
 536     h__ -= h_offset;
 537     --gpls;
 538     --xpls;
 539     --x;
 540
 541     /* Function Body */
 542     dfault_(n, &wrk[wrk_dim1 * 5 + 1], &fscale, &gradtl, &steptl, &itnlim, &
 543             stepmx, &ipr, &method, &grdflg, &hesflg, &ndigit, msg);
 544     tensor_(nr, n, &x[1], (U_fp)fcn, (S_fp)grd_, (S_fp)hsn_, &wrk[wrk_dim1 *
 545             5 + 1], &fscale, &gradtl, &steptl, &itnlim, &stepmx, &ipr, &
 546             method, &grdflg, &hesflg, &ndigit, msg, &xpls[1], fpls, &gpls[1],
 547             &h__[h_offset], itnno, &wrk[wrk_offset], &iwrk[1]);
 548     return 0;
 549 } /* tensrd_ */
 550
 551 /*  ---------------------- */
 552 /*  |       O P T        | */
 553 /*  ---------------------- */
 554 /* Subroutine */ int opt_(integer *nr, integer *n, doublereal *x, S_fp fcn,
 555         S_fp grd, S_fp hsn, doublereal *typx, doublereal *fscale, doublereal *
 556         gradtl, doublereal *steptl, integer *itnlim, doublereal *stepmx,
 557         integer *ipr, integer *method, integer *grdflg, integer *hesflg,
 558         integer *ndigit, integer *msg, doublereal *xpls, doublereal *fpls,
 559         doublereal *gpls, doublereal *h__, integer *itnno, doublereal *g,
 560         doublereal *s, doublereal *d__, doublereal *dn, doublereal *e,
 561         doublereal *wk1, doublereal *wk2, integer *pivot)
 562 {
 563     /* Format strings */
 564     static char fmt_25[] = "(\002 INITIAL FUNCTION VALUE F=\002,e20.13)";
 565     static char fmt_30[] = "(\002 INITIAL GRADIENT G=\002,999e20.13)";
 566     static char fmt_103[] = "(\002 -----------    ITERATION \002,i4,\002 ---"
 567             "-------------\002)";
 568     static char fmt_104[] = "(\002 X=\002,999e20.13)";
 569     static char fmt_105[] = "(\002 F(X)=\002,e20.13)";
 570     static char fmt_106[] = "(\002 G(X)=\002,999e20.13)";
 571     static char fmt_108[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
 572             ". MAX:\002,e20.13)";
 573     static char fmt_110[] = "(\002REL. STEP MAX :\002,e20.13)";
 574     static char fmt_370[] = "(\002 FINAL X=\002,999e20.13)";
 575     static char fmt_380[] = "(\002 GRADIENT G=\002,999e20.13)";
 576     static char fmt_390[] = "(\002 FUNCTION VALUE F(X)=\002,e20.13,\002 AT I"
 577             "TERATION \002,i4)";
 578     static char fmt_400[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
 579             ". MAX:\002,e20.13)";
 580     static char fmt_410[] = "(\002REL. STEP MAX :\002,e20.13)";
 581     static char fmt_560[] = "(\002 -----------    ITERATION \002,i4,\002 ---"
 582             "-------------\002)";
 583     static char fmt_570[] = "(\002 X=\002,999e20.13)";
 584     static char fmt_580[] = "(\002 F(X)=\002,e20.13)";
 585     static char fmt_590[] = "(\002 G(X)=\002,999e20.13)";
 586     static char fmt_600[] = "(\002FUNCTION EVAL COUNT:\002,i6,\002 REL. GRAD"
 587             ". MAX:\002,e20.13)";
 588     static char fmt_610[] = "(\002REL. STEP MAX :\002,e20.13)";
 589
 590     /* System generated locals */
 591     integer h_dim1, h_offset, i__1, i__2;
 592     doublereal d__1, d__2;
 593
 594     /* Builtin functions */
 595     double pow_di(doublereal *, integer *), sqrt(doublereal);
 596     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
 597
 598     /* Local variables */
 599     doublereal f;
 600     integer i__;
 601     doublereal gd, fn, fp, rnf, eps, rgx, rsx, beta;
 602     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
 603             integer *);
 604     integer imsg;
 605     doublereal temp, alpha;
 606     extern /* Subroutine */ int mkmdl_(integer *, integer *, doublereal *,
 607             doublereal *, doublereal *, doublereal *, doublereal *,
 608             doublereal *, doublereal *, doublereal *, doublereal *,
 609             doublereal *);
 610     integer nfcnt;
 611     doublereal finit;
 612     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
 613             doublereal *, integer *);
 614     logical nomin;
 615     doublereal gnorm, addmax;
 616     extern /* Subroutine */ int grdchk_(integer *, doublereal *, S_fp,
 617             doublereal *, doublereal *, doublereal *, doublereal *,
 618             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
 619              integer *, integer *), heschk_(integer *, integer *, doublereal *
 620             , S_fp, S_fp, S_fp, doublereal *, doublereal *, doublereal *,
 621             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
 622              doublereal *, doublereal *, doublereal *, integer *, integer *,
 623             integer *);
 624     integer iretcd;
 625     doublereal analtl;
 626     extern /* Subroutine */ int choldr_(integer *, integer *, doublereal *,
 627             doublereal *, doublereal *, integer *, doublereal *, doublereal *,
 628              doublereal *), mcheps_(doublereal *), sndofd_(integer *, integer
 629             *, doublereal *, S_fp, doublereal *, doublereal *, doublereal *,
 630             doublereal *, doublereal *, doublereal *, integer *);
 631     integer itrmcd;
 632     extern /* Subroutine */ int bakslv_(integer *, integer *, doublereal *,
 633             doublereal *, doublereal *), fstofd_(integer *, integer *,
 634             integer *, doublereal *, S_fp, doublereal *, doublereal *,
 635             doublereal *, doublereal *, doublereal *, integer *, integer *);
 636     integer icscmx;
 637     extern /* Subroutine */ int optchk_(integer *, integer *, integer *,
 638             doublereal *, doublereal *, doublereal *, doublereal *,
 639             doublereal *, integer *, integer *, doublereal *, integer *,
 640             integer *, integer *, doublereal *, integer *);
 641     logical mxtake;
 642     extern /* Subroutine */ int lnsrch_(integer *, integer *, doublereal *,
 643             doublereal *, doublereal *, doublereal *, doublereal *,
 644             doublereal *, logical *, integer *, doublereal *, doublereal *,
 645             doublereal *, S_fp, doublereal *, integer *), slvmdl_(integer *,
 646             integer *, doublereal *, doublereal *, doublereal *, doublereal *,
 647              doublereal *, doublereal *, doublereal *, integer *, doublereal *
 648             , doublereal *, doublereal *, doublereal *, doublereal *, logical
 649             *, doublereal *), forslv_(integer *, integer *, doublereal *,
 650             doublereal *, doublereal *);
 651     extern doublereal twonrm_(integer *, doublereal *);
 652     extern /* Subroutine */ int optstp_(integer *, doublereal *, doublereal *,
 653              doublereal *, doublereal *, integer *, integer *, integer *,
 654             doublereal *, doublereal *, doublereal *, integer *, integer *,
 655             logical *, integer *, integer *, doublereal *, doublereal *);
 656
 657     /* Fortran I/O blocks */
 658     static cilist io___70 = { 0, 0, 0, fmt_25, 0 };
 659     static cilist io___71 = { 0, 0, 0, fmt_30, 0 };
 660     static cilist io___81 = { 0, 0, 0, fmt_103, 0 };
 661     static cilist io___82 = { 0, 0, 0, fmt_104, 0 };
 662     static cilist io___83 = { 0, 0, 0, fmt_105, 0 };
 663     static cilist io___84 = { 0, 0, 0, fmt_106, 0 };
 664     static cilist io___85 = { 0, 0, 0, fmt_108, 0 };
 665     static cilist io___86 = { 0, 0, 0, fmt_110, 0 };
 666     static cilist io___93 = { 0, 0, 0, fmt_370, 0 };
 667     static cilist io___94 = { 0, 0, 0, fmt_380, 0 };
 668     static cilist io___95 = { 0, 0, 0, fmt_390, 0 };
 669     static cilist io___96 = { 0, 0, 0, fmt_400, 0 };
 670     static cilist io___97 = { 0, 0, 0, fmt_410, 0 };
 671     static cilist io___98 = { 0, 0, 0, fmt_560, 0 };
 672     static cilist io___99 = { 0, 0, 0, fmt_570, 0 };
 673     static cilist io___100 = { 0, 0, 0, fmt_580, 0 };
 674     static cilist io___101 = { 0, 0, 0, fmt_590, 0 };
 675     static cilist io___102 = { 0, 0, 0, fmt_600, 0 };
 676     static cilist io___103 = { 0, 0, 0, fmt_610, 0 };
 677
 678
 679
 680 /* PURPOSE: */
 681 /*   DERIVATIVE TENSOR METHODS FOR UNCONSTRAINED OPTIMIZATION */
 682
 683 /* ----------------------------------------------------------------------- */
 684
 685 /* PARAMETERS: */
 686
 687 /*   NR      --> ROW DIMENSION OF MATRIX */
 688 /*   N       --> DIMENSION OF PROBLEM */
 689 /*   X(N)    --> INITIAL GUESS (INPUT) AND CURRENT POINT */
 690 /*   FCN     --> NAME OF SUBROUTINE TO EVALUATE FUNCTION VALUE */
 691 /*               FCN: R(N) --> R(1) */
 692 /*   GRD     --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL GRADIENT */
 693 /*               FCN: R(N) --> R(N) */
 694 /*   HSN     --> NAME OF SUBROUTINE TO EVALUATE ANALYTICAL HESSIAN */
 695 /*               FCN: R(N) --> R(N X N) */
 696 /*               HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
 697 /*               AND DIAGONAL OF H */
 698 /*   TYPX(N) --> TYPICAL SIZE OF EACH COMPONENT OF X */
 699 /*   FSCALE  --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
 700 /*   GRADTL  --> GRADIENT TOLERANCE */
 701 /*   STEPTL  --> STEP TOLERANCE */
 702 /*   ITNLIM  --> ITERATION LIMIT */
 703 /*   STEPMX  --> MAXIMUM STEP LENGTH ALLOWED */
 704 /*   IPR     --> OUTPUT UNIT NUMBER */
 705 /*   METHOD  --> IF VALUE IS 0 THEN USE ONLY NEWTON STEP AT */
 706 /*               EACH ITERATION */
 707 /*   GRDFLG  --> = 1 OR 2 IF ANALYTICAL GRADIENT IS SUPPLIED */
 708 /*   HESFLG  --> = 1 OR 2 IF ANALYTICAL HESSIAN IS SUPPLIED */
 709 /*   NDIGIT  --> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
 710 /*   MSG     --> OUTPUT MESSAGE CONTRO */
 711 /*   XPLS(N) <--  NEW POINT AND FINAL POINT (OUTPUT) */
 712 /*   FPLS    <--  NEW FUNCTION VALUE AND FINAL FUNCTION VALUE (OUTPUT) */
 713 /*   GPLS(N) <--  CURRENT GRADIENT AND GRADIENT AT FINAL POINT (OUTPUT) */
 714 /*   H(N,N)  --> HESSIAN */
 715 /*   ITNNO   <--  NUMBER OF ITERATIONS */
 716 /*   G(N)    --> PREVIOUS GRADIENT */
 717 /*   S       --> STEP TO PREVIOUS ITERATE (FOR TENSOR MODEL) */
 718 /*   D       --> CURRENT STEP (TENSOR OR NEWTON) */
 719 /*   DN      --> NEWTON STEP */
 720 /*   E(N)    --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION */
 721 /*   WK1(N)  --> WORKSPACE */
 722 /*   WK2(N)  --> WORKSPACE */
 723 /*   PIVOT(N)--> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION */
 724
 725 /* -------------------------------- */
 726 /*     INITIALIZATION           | */
 727 /* -------------------------------- */
 728
 729     /* Parameter adjustments */
 730     --pivot;
 731     --wk2;
 732     --wk1;
 733     --e;
 734     --dn;
 735     --d__;
 736     --s;
 737     --g;
 738     h_dim1 = *nr;
 739     h_offset = 1 + h_dim1;
 740     h__ -= h_offset;
 741     --gpls;
 742     --xpls;
 743     --typx;
 744     --x;
 745
 746     /* Function Body */
 747     *itnno = 0;
 748     icscmx = 0;
 749     nfcnt = 0;
 750     mcheps_(&eps);
 751     optchk_(nr, n, msg, &x[1], &typx[1], fscale, gradtl, steptl, itnlim,
 752             ndigit, &eps, method, grdflg, hesflg, stepmx, ipr);
 753     if (*msg < 0) {
 754         return 0;
 755     }
 756
 757 /*     SCALE X */
 758     i__1 = *n;
 759     for (i__ = 1; i__ <= i__1; ++i__) {
 760         x[i__] /= typx[i__];
 761 /* L10: */
 762     }
 763
 764 /* -------------------------------- */
 765 /*     INITIAL ITERATION | */
 766 /* -------------------------------- */
 767
 768 /*     COMPUTE TYPICAL SIZE OF X */
 769     i__1 = *n;
 770     for (i__ = 1; i__ <= i__1; ++i__) {
 771         dn[i__] = 1. / typx[i__];
 772 /* L15: */
 773     }
 774 /* Computing MAX */
 775     i__1 = -(*ndigit);
 776     d__1 = pow_di(&c_b134, &i__1);
 777     rnf = max(d__1,eps);
 778 /* Computing MAX */
 779     d__1 = .01, d__2 = sqrt(rnf);
 780     analtl = max(d__1,d__2);
 781 /*     UNSCALE X AND COMPUTE F AND G */
 782     i__1 = *n;
 783     for (i__ = 1; i__ <= i__1; ++i__) {
 784         wk1[i__] = x[i__] * typx[i__];
 785 /* L20: */
 786     }
 787     (*fcn)(n, &wk1[1], &f);
 788     ++nfcnt;
 789     if (*grdflg >= 1) {
 790         (*grd)(n, &wk1[1], &g[1]);
 791         if (*grdflg == 1) {
 792             *fscale = 1.;
 793             grdchk_(n, &wk1[1], (S_fp)fcn, &f, &g[1], &dn[1], &typx[1],
 794                     fscale, &rnf, &analtl, &wk2[1], msg, ipr, &nfcnt);
 795         }
 796     } else {
 797         fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, &f, &g[1], &typx[1], &
 798                 rnf, &wk2[1], &c__1, &nfcnt);
 799     }
 800     if (*msg < -20) {
 801         return 0;
 802     }
 803
 804     gnorm = twonrm_(n, &g[1]);
 805
 806 /*     PRINT OUT 1ST ITERATION? */
 807     if (*msg >= 1) {
 808         io___70.ciunit = *ipr;
 809         s_wsfe(&io___70);
 810         do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
 811         e_wsfe();
 812         io___71.ciunit = *ipr;
 813         s_wsfe(&io___71);
 814         i__1 = *n;
 815         for (i__ = 1; i__ <= i__1; ++i__) {
 816             do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
 817         }
 818         e_wsfe();
 819     }
 820
 821 /*     TEST WHETHER THE INITIAL GUESS SATISFIES THE STOPPING CRITERIA */
 822     if (gnorm <= *gradtl) {
 823         *fpls = f;
 824         i__1 = *n;
 825         for (i__ = 1; i__ <= i__1; ++i__) {
 826             xpls[i__] = x[i__];
 827             gpls[i__] = g[i__];
 828 /* L40: */
 829         }
 830         optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
 831                 gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &
 832                 rgx, &rsx);
 833         goto L350;
 834     }
 835     finit = f;
 836
 837 /* ------------------------ */
 838 /*     ITERATION 1 | */
 839 /* ------------------------ */
 840
 841     ++(*itnno);
 842
 843 /*     COMPUTE HESSIAN */
 844     if (*hesflg == 0) {
 845         if (*grdflg == 1) {
 846             fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
 847                     typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
 848         } else {
 849             sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
 850                     rnf, &wk2[1], &d__[1], &nfcnt);
 851         }
 852     } else {
 853         if (*hesflg == 2) {
 854             (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
 855         } else {
 856             if (*hesflg == 1) {
 857 /*           IN HESCHK GPLS,XPLS AND E ARE USED AS WORK SPACE */
 858                 heschk_(nr, n, &wk1[1], (S_fp)fcn, (S_fp)grd, (S_fp)hsn, &f, &
 859                         g[1], &h__[h_offset], &dn[1], &typx[1], &rnf, &analtl,
 860                          grdflg, &gpls[1], &xpls[1], &e[1], msg, ipr, &nfcnt);
 861             }
 862         }
 863     }
 864     if (*msg < -20) {
 865         return 0;
 866     }
 867     i__1 = *n;
 868     for (i__ = 2; i__ <= i__1; ++i__) {
 869         i__2 = i__ - 1;
 870         dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
 871 /* L50: */
 872     }
 873
 874 /*     CHOLESKY DECOMPOSITION FOR H (H=LLT) */
 875     choldr_(nr, n, &h__[h_offset], &d__[1], &eps, &pivot[1], &e[1], &wk1[1], &
 876             addmax);
 877
 878 /*     SOLVE FOR NEWTON STEP D */
 879     i__1 = *n;
 880     for (i__ = 1; i__ <= i__1; ++i__) {
 881 /* L60: */
 882         wk2[i__] = -g[i__];
 883     }
 884     forslv_(nr, n, &h__[h_offset], &wk1[1], &wk2[1]);
 885     bakslv_(nr, n, &h__[h_offset], &d__[1], &wk1[1]);
 886
 887 /*     APPLY LINESEARCH TO THE NEWTON STEP */
 888     lnsrch_(nr, n, &x[1], &f, &g[1], &d__[1], &xpls[1], fpls, &mxtake, &
 889             iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
 890
 891 /*     UPDATE G */
 892 /*      CALL DCOPY(N,GPLS(1),1,GP(1),1) */
 893
 894 /*     UNSCALE XPLS AND COMPUTE GPLS */
 895     i__1 = *n;
 896     for (i__ = 1; i__ <= i__1; ++i__) {
 897         wk1[i__] = xpls[i__] * typx[i__];
 898 /* L80: */
 899     }
 900     if (*grdflg >= 1) {
 901         (*grd)(n, &wk1[1], &gpls[1]);
 902     } else {
 903         fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
 904                  &rnf, &wk2[1], &c__1, &nfcnt);
 905     }
 906
 907 /*     CHECK STOPPING CONDITIONS */
 908     optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
 909             gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &rgx,
 910             &rsx);
 911
 912 /*     IF ITRMCD > 0 THEN STOPPING CONDITIONS SATISFIED */
 913     if (itrmcd > 0) {
 914         goto L350;
 915     }
 916
 917 /*     UPDATE X,F AND S FOR TENSOR MODEL */
 918     fp = f;
 919     f = *fpls;
 920     i__1 = *n;
 921     for (i__ = 1; i__ <= i__1; ++i__) {
 922         temp = xpls[i__];
 923         s[i__] = x[i__] - temp;
 924         x[i__] = temp;
 925 /* L90: */
 926     }
 927
 928 /*     IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
 929     if (*msg >= 2) {
 930         io___81.ciunit = *ipr;
 931         s_wsfe(&io___81);
 932         do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
 933         e_wsfe();
 934         io___82.ciunit = *ipr;
 935         s_wsfe(&io___82);
 936         i__1 = *n;
 937         for (i__ = 1; i__ <= i__1; ++i__) {
 938             do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
 939         }
 940         e_wsfe();
 941         io___83.ciunit = *ipr;
 942         s_wsfe(&io___83);
 943         do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
 944         e_wsfe();
 945         io___84.ciunit = *ipr;
 946         s_wsfe(&io___84);
 947         i__1 = *n;
 948         for (i__ = 1; i__ <= i__1; ++i__) {
 949             do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
 950         }
 951         e_wsfe();
 952     }
 953     if (*msg >= 3) {
 954         io___85.ciunit = *ipr;
 955         s_wsfe(&io___85);
 956         do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
 957         do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
 958         e_wsfe();
 959         io___86.ciunit = *ipr;
 960         s_wsfe(&io___86);
 961         do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
 962         e_wsfe();
 963     }
 964
 965 /* ------------------------ */
 966 /*     ITERATION > 1     | */
 967 /* ------------------------ */
 968
 969 /*     UNSCALE X AND COMPUTE H */
 970 L200:
 971     i__1 = *n;
 972     for (i__ = 1; i__ <= i__1; ++i__) {
 973         wk1[i__] = x[i__] * typx[i__];
 974 /* L210: */
 975     }
 976     if (*hesflg == 0) {
 977         if (*grdflg == 1) {
 978             fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
 979                     typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
 980         } else {
 981             sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
 982                     rnf, &wk2[1], &d__[1], &nfcnt);
 983         }
 984     } else {
 985         (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
 986     }
 987     i__1 = *n;
 988     for (i__ = 2; i__ <= i__1; ++i__) {
 989         i__2 = i__ - 1;
 990         dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
 991 /* L230: */
 992     }
 993
 994 /*     IF METHOD = 0 THEN USE NEWTON STEP ONLY */
 995     if (*method == 0) {
 996
 997 /*       CHOLESKY DECOMPOSITION FOR H */
 998         choldr_(nr, n, &h__[h_offset], &wk2[1], &eps, &pivot[1], &e[1], &wk1[
 999                 1], &addmax);
1000
1001 /*       COMPUTE NEWTON STEP */
1002         i__1 = *n;
1003         for (i__ = 1; i__ <= i__1; ++i__) {
1004             wk1[i__] = -gpls[i__];
1005 /* L240: */
1006         }
1007         forslv_(nr, n, &h__[h_offset], &wk2[1], &wk1[1]);
1008         bakslv_(nr, n, &h__[h_offset], &d__[1], &wk2[1]);
1009
1010 /*       NO TENSOR STEP */
1011         nomin = TRUE_;
1012         goto L300;
1013
1014     }
1015
1016 /*     FORM TENSOR MODEL */
1017     mkmdl_(nr, n, &f, &fp, &gpls[1], &g[1], &s[1], &h__[h_offset], &alpha, &
1018             beta, &wk1[1], &d__[1]);
1019
1020 /*   SOLVE TENSOR MODEL AND COMPUTE NEWTON STEP */
1021 /*   ON INPUT : SH IS STORED IN WK1 */
1022 /*            A=(G-GPLS-SH-S*BETA/(6*STS)) IS STORED IN D */
1023 /*   ON OUTPUT: NEWTON STEP IS STORED IN DN */
1024 /*            TENSOR STEP IS STORED IN D */
1025     slvmdl_(nr, n, &h__[h_offset], &xpls[1], &wk2[1], &e[1], &g[1], &s[1], &
1026             gpls[1], &pivot[1], &d__[1], &wk1[1], &dn[1], &alpha, &beta, &
1027             nomin, &eps);
1028
1029 /*     IF TENSOR MODEL HAS NO MINIMIZER THEN USE NEWTON STEP */
1030     if (nomin) {
1031         dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
1032         goto L300;
1033     }
1034
1035 /*     IF TENSOR STEP IS NOT IN DESCENT DIRECTION THEN USE NEWTON STEP */
1036     gd = ddot_(n, &gpls[1], &c__1, &d__[1], &c__1);
1037     if (gd > 0.) {
1038         dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
1039         nomin = TRUE_;
1040     }
1041
1042 L300:
1043     ++(*itnno);
1044     dcopy_(n, &gpls[1], &c__1, &g[1], &c__1);
1045
1046 /*     APPLY LINESEARCH TO TENSOR (OR NEWTON) STEP */
1047     lnsrch_(nr, n, &x[1], &f, &g[1], &d__[1], &xpls[1], fpls, &mxtake, &
1048             iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
1049
1050     if (! nomin) {
1051 /*       TENSOR STEP IS FOUND AND IN DESCENT DIRECTION, */
1052 /*       APPLY LINESEARCH TO NEWTON STEP */
1053 /*       NEW NEWTON POINT IN WK2 */
1054         lnsrch_(nr, n, &x[1], &f, &g[1], &dn[1], &wk2[1], &fn, &mxtake, &
1055                 iretcd, stepmx, steptl, &typx[1], (S_fp)fcn, &wk1[1], &nfcnt);
1056
1057 /*       COMPARE TENSOR STEP TO NEWTON STEP */
1058 /*       IF NEWTON STEP IS BETTER, SET NEXT ITERATE TO NEW NEWTON POINT */
1059         if (fn < *fpls) {
1060             *fpls = fn;
1061             dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
1062             dcopy_(n, &wk2[1], &c__1, &xpls[1], &c__1);
1063         }
1064     }
1065     i__1 = *n;
1066     for (i__ = 1; i__ <= i__1; ++i__) {
1067         d__[i__] = xpls[i__] - x[i__];
1068 /* L320: */
1069     }
1070
1071 /*     UNSCALE XPLS, AND COMPUTE FPLS AND GPLS */
1072     i__1 = *n;
1073     for (i__ = 1; i__ <= i__1; ++i__) {
1074         wk1[i__] = xpls[i__] * typx[i__];
1075 /* L330: */
1076     }
1077     (*fcn)(n, &wk1[1], fpls);
1078     ++nfcnt;
1079     if (*grdflg >= 1) {
1080         (*grd)(n, &wk1[1], &gpls[1]);
1081     } else {
1082         fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
1083                  &rnf, &wk2[1], &c__1, &nfcnt);
1084     }
1085
1086 /*     CHECK STOPPING CONDITIONS */
1087     imsg = *msg;
1088     optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
1089             gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &rgx,
1090             &rsx);
1091
1092 /*     IF ITRMCD = 0 THEN NOT OVER YET */
1093     if (itrmcd == 0) {
1094         goto L500;
1095     }
1096
1097 /*     IF MSG >= 1 THEN PRINT OUT FINAL ITERATION */
1098 L350:
1099     if (imsg >= 1) {
1100
1101 /*       TRANSFORM X BACK TO ORIGINAL SPACE */
1102         i__1 = *n;
1103         for (i__ = 1; i__ <= i__1; ++i__) {
1104             xpls[i__] *= typx[i__];
1105 /* L360: */
1106         }
1107         io___93.ciunit = *ipr;
1108         s_wsfe(&io___93);
1109         i__1 = *n;
1110         for (i__ = 1; i__ <= i__1; ++i__) {
1111             do_fio(&c__1, (char *)&xpls[i__], (ftnlen)sizeof(doublereal));
1112         }
1113         e_wsfe();
1114         io___94.ciunit = *ipr;
1115         s_wsfe(&io___94);
1116         i__1 = *n;
1117         for (i__ = 1; i__ <= i__1; ++i__) {
1118             do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
1119         }
1120         e_wsfe();
1121         io___95.ciunit = *ipr;
1122         s_wsfe(&io___95);
1123         do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
1124         do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
1125         e_wsfe();
1126         if (imsg >= 3) {
1127             io___96.ciunit = *ipr;
1128             s_wsfe(&io___96);
1129             do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
1130             do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
1131             e_wsfe();
1132             io___97.ciunit = *ipr;
1133             s_wsfe(&io___97);
1134             do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
1135             e_wsfe();
1136         }
1137 /*       UPDATE THE HESSIAN */
1138         if (*hesflg == 0) {
1139             if (*grdflg == 1) {
1140                 fstofd_(nr, n, n, &xpls[1], (S_fp)grd, &gpls[1], &h__[
1141                         h_offset], &typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
1142             } else {
1143                 sndofd_(nr, n, &xpls[1], (S_fp)fcn, fpls, &h__[h_offset], &
1144                         typx[1], &rnf, &wk2[1], &d__[1], &nfcnt);
1145             }
1146         } else {
1147             (*hsn)(nr, n, &xpls[1], &h__[h_offset]);
1148         }
1149     }
1150     return 0;
1151
1152 /*     UPDATE INFORMATION AT THE CURRENT POINT */
1153 L500:
1154     dcopy_(n, &xpls[1], &c__1, &x[1], &c__1);
1155     i__1 = *n;
1156     for (i__ = 1; i__ <= i__1; ++i__) {
1157         s[i__] = -d__[i__];
1158 /* L550: */
1159     }
1160
1161 /*     IF TOO MANY ITERATIONS THEN RETURN */
1162     if (*itnno > *itnlim) {
1163         goto L350;
1164     }
1165
1166 /*     IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
1167     if (*msg >= 2) {
1168         io___98.ciunit = *ipr;
1169         s_wsfe(&io___98);
1170         do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
1171         e_wsfe();
1172         io___99.ciunit = *ipr;
1173         s_wsfe(&io___99);
1174         i__1 = *n;
1175         for (i__ = 1; i__ <= i__1; ++i__) {
1176             do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
1177         }
1178         e_wsfe();
1179         io___100.ciunit = *ipr;
1180         s_wsfe(&io___100);
1181         do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
1182         e_wsfe();
1183         io___101.ciunit = *ipr;
1184         s_wsfe(&io___101);
1185         i__1 = *n;
1186         for (i__ = 1; i__ <= i__1; ++i__) {
1187             do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
1188         }
1189         e_wsfe();
1190     }
1191     if (*msg >= 3) {
1192         io___102.ciunit = *ipr;
1193         s_wsfe(&io___102);
1194         do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
1195         do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
1196         e_wsfe();
1197         io___103.ciunit = *ipr;
1198         s_wsfe(&io___103);
1199         do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
1200         e_wsfe();
1201     }
1202 /*     UPDATE F */
1203     fp = f;
1204     f = *fpls;
1205
1206 /*     PERFORM NEXT ITERATION */
1207     goto L200;
1208
1209 /*     END OF ITERATION > 1 */
1210
1211 } /* opt_ */
1212
1213 /*  ---------------------- */
1214 /*  |  D F A U L T      | */
1215 /*  ---------------------- */
1216 /* Subroutine */ int dfault_(integer *n, doublereal *typx, doublereal *fscale,
1217          doublereal *gradtl, doublereal *steptl, integer *itnlim, doublereal *
1218         stepmx, integer *ipr, integer *method, integer *grdflg, integer *
1219         hesflg, integer *ndigit, integer *msg)
1220 {
1221     /* System generated locals */
1222     integer i__1;
1223
1224     /* Builtin functions */
1225     double pow_dd(doublereal *, doublereal *), d_lg10(doublereal *);
1226
1227     /* Local variables */
1228     integer i__;
1229     doublereal eps, temp;
1230     extern /* Subroutine */ int mcheps_(doublereal *);
1231
1232     /* Parameter adjustments */
1233     --typx;
1234
1235     /* Function Body */
1236     mcheps_(&eps);
1237     *method = 1;
1238     *fscale = 1.;
1239     *grdflg = 0;
1240     *hesflg = 0;
1241     i__1 = *n;
1242     for (i__ = 1; i__ <= i__1; ++i__) {
1243         typx[i__] = 1.;
1244 /* L1: */
1245     }
1246     temp = pow_dd(&eps, &c_b250);
1247     *gradtl = temp;
1248     *steptl = temp * temp;
1249     *ndigit = (integer) (-d_lg10(&eps));
1250 /*     SET ACTUAL DFAULT VALUE OF STEPMX IN OPTCHK */
1251     *stepmx = 0.;
1252     *itnlim = 100;
1253     *ipr = 6;
1254     *msg = 1;
1255     return 0;
1256 } /* dfault_ */
1257
1258 /*  ---------------------- */
1259 /*  |  O P T C H K       | */
1260 /*  ---------------------- */
1261 /* Subroutine */ int optchk_(integer *nr, integer *n, integer *msg,
1262         doublereal *x, doublereal *typx, doublereal *fscale, doublereal *
1263         gradtl, doublereal *steptl, integer *itnlim, integer *ndigit,
1264         doublereal *eps, integer *method, integer *grdflg, integer *hesflg,
1265         doublereal *stepmx, integer *ipr)
1266 {
1267     /* Format strings */
1268     static char fmt_901[] = "(\0020OPTCHK    ILLEGAL DIMENSION, N=\002,i5)";
1269     static char fmt_902[] = "(\002OPTCHK    ILLEGAL INPUT VALUES OF: NR=\002"
1270             ",i5,\002, N=\002,i5,\002, NR MUST BE <=  N.\002)";
1271
1272     /* System generated locals */
1273     integer i__1;
1274     doublereal d__1;
1275
1276     /* Builtin functions */
1277     double sqrt(doublereal), pow_dd(doublereal *, doublereal *), d_lg10(
1278             doublereal *), pow_di(doublereal *, integer *);
1279     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
1280
1281     /* Local variables */
1282     integer i__;
1283     doublereal temp, stpsiz;
1284
1285     /* Fortran I/O blocks */
1286     static cilist io___110 = { 0, 0, 0, fmt_901, 0 };
1287     static cilist io___111 = { 0, 0, 0, fmt_902, 0 };
1288
1289
1290
1291 /* PURPOSE */
1292 /* ------- */
1293 /* CHECK INPUT FOR REASONABLENESS */
1294
1295 /* PARAMETERS */
1296 /* ---------- */
1297 /* NR          <--> ROW DIMENSION OF H AND WRK */
1298 /* N            --> DIMENSION OF PROBLEM */
1299 /* X(N)         --> ON ENTRY, ESTIMATE TO ROOT OF FCN */
1300 /* TYPX(N)     <--> TYPICAL SIZE OF EACH COMPONENT OF X */
1301 /* FSCALE      <--> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
1302 /* GRADTL      <--> TOLERANCE AT WHICH GRADIENT CONSIDERED CLOSE */
1303 /*                  ENOUGH TO ZERO TO TERMINATE ALGORITHM */
1304 /* STEPTL      <--> TOLERANCE AT WHICH STEP LENGTH CONSIDERED CLOSE */
1305 /*                  ENOUGH TO ZERO TO TERMINATE ALGORITHM */
1306 /* ITNLIM      <--> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS */
1307 /* NDIGIT      <--> NUMBER OF GOOD DIGITS IN OPTIMIZATION FUNCTION FCN */
1308 /* EPS          --> MACHINE PRECISION */
1309 /* METHOD      <--> ALGORITHM INDICATOR */
1310 /* GRDFLG      <--> =1 IF ANALYTIC GRADIENT SUPPLIED */
1311 /* HESFLG      <--> =1 IF ANALYTIC HESSIAN SUPPLIED */
1312 /* STEPMX      <--> MAXIMUM STEP SIZE */
1313 /* MSG         <--> MESSAGE AND ERROR CODE */
1314 /* IPR          --> DEVICE TO WHICH TO SEND OUTPUT */
1315
1316
1317 /* CHECK DIMENSION OF PROBLEM */
1318
1319     /* Parameter adjustments */
1320     --typx;
1321     --x;
1322
1323     /* Function Body */
1324     if (*n <= 0) {
1325         goto L805;
1326     }
1327     if (*nr < *n) {
1328         goto L806;
1329     }
1330
1331 /* CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES. */
1332 /* IF NOT, SET THEM TO DEFAULT VALUES. */
1333     if (*method != 0) {
1334         *method = 1;
1335     }
1336     if (*grdflg != 2 && *grdflg != 1) {
1337         *grdflg = 0;
1338     }
1339     if (*hesflg != 2 && *hesflg != 1) {
1340         *hesflg = 0;
1341     }
1342     if (*msg > 3 || (real) (*msg) < 0.f) {
1343         *msg = 1;
1344     }
1345
1346 /* COMPUTE SCALE MATRIX */
1347
1348     i__1 = *n;
1349     for (i__ = 1; i__ <= i__1; ++i__) {
1350         if (typx[i__] == 0.f) {
1351             typx[i__] = 1.;
1352         }
1353         if (typx[i__] < 0.f) {
1354             typx[i__] = -typx[i__];
1355         }
1356 /* L10: */
1357     }
1358
1359 /* CHECK MAXIMUM STEP SIZE */
1360
1361     if (*stepmx > 0.) {
1362         goto L20;
1363     }
1364     stpsiz = 0.;
1365     i__1 = *n;
1366     for (i__ = 1; i__ <= i__1; ++i__) {
1367         stpsiz += x[i__] * x[i__] / typx[i__] * typx[i__];
1368 /* L15: */
1369     }
1370     stpsiz = sqrt(stpsiz);
1371 /* Computing MAX */
1372     d__1 = stpsiz * 1e3;
1373     *stepmx = max(d__1,1e3);
1374 L20:
1375 /* CHECK FUNCTION SCALE */
1376     if (*fscale == 0.f) {
1377         *fscale = 1.;
1378     }
1379     if (*fscale < 0.f) {
1380         *fscale = -(*fscale);
1381     }
1382 /* CHECK GRADIENT TOLERANCE */
1383     if (*gradtl < 0.f) {
1384         *gradtl = pow_dd(eps, &c_b250);
1385     }
1386 /* CHECK STEP TOLERANCE */
1387     if (*steptl < 0.f) {
1388         temp = pow_dd(eps, &c_b250);
1389         *steptl = temp * temp;
1390     }
1391
1392 /* CHECK ITERATION LIMIT */
1393     if (*itnlim <= 0) {
1394         *itnlim = 100;
1395     }
1396
1397 /* CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN */
1398     if (*ndigit <= 0) {
1399         *ndigit = (integer) (-d_lg10(eps));
1400     }
1401     i__1 = -(*ndigit);
1402     if (pow_di(&c_b134, &i__1) <= *eps) {
1403         *ndigit = (integer) (-d_lg10(eps));
1404     }
1405     return 0;
1406
1407 /* ERROR EXITS */
1408
1409 L805:
1410     io___110.ciunit = *ipr;
1411     s_wsfe(&io___110);
1412     do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
1413     e_wsfe();
1414     *msg = -20;
1415     return 0;
1416 L806:
1417     io___111.ciunit = *ipr;
1418     s_wsfe(&io___111);
1419     do_fio(&c__1, (char *)&(*nr), (ftnlen)sizeof(integer));
1420     do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
1421     e_wsfe();
1422     *msg = -20;
1423     return 0;
1424 } /* optchk_ */
1425
1426 /*  ---------------------- */
1427 /*  |  G R D C H K       | */
1428 /*  ---------------------- */
1429 /* Subroutine */ int grdchk_(integer *n, doublereal *x, S_fp fcn, doublereal *
1430         f, doublereal *g, doublereal *typsiz, doublereal *typx, doublereal *
1431         fscale, doublereal *rnf, doublereal *analtl, doublereal *wrk1,
1432         integer *msg, integer *ipr, integer *ifn)
1433 {
1434     /* Format strings */
1435     static char fmt_901[] = "(\0020GRDCHK    PROBABLE ERROR IN CODING OF ANA"
1436             "LYTIC\002,\002 GRADIENT FUNCTION.\002/\002 GRDCHK     COMP\002,1"
1437             "2x,\002ANALYTIC\002,12x,\002ESTIMATE\002)";
1438     static char fmt_902[] = "(\002 GRDCHK    \002,i5,3x,e20.13,3x,e20.13)";
1439
1440     /* System generated locals */
1441     integer i__1;
1442     doublereal d__1, d__2, d__3, d__4;
1443
1444     /* Builtin functions */
1445     integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
1446
1447     /* Local variables */
1448     integer i__;
1449     doublereal gs;
1450     integer ker;
1451     doublereal wrk;
1452     extern /* Subroutine */ int fstofd_(integer *, integer *, integer *,
1453             doublereal *, S_fp, doublereal *, doublereal *, doublereal *,
1454             doublereal *, doublereal *, integer *, integer *);
1455
1456     /* Fortran I/O blocks */
1457     static cilist io___116 = { 0, 0, 0, fmt_901, 0 };
1458     static cilist io___117 = { 0, 0, 0, fmt_902, 0 };
1459
1460
1461
1462 /* PURPOSE */
1463 /* ------- */
1464 /* CHECK ANALYTIC GRADIENT AGAINST ESTIMATED GRADIENT */
1465
1466 /* PARAMETERS */
1467 /* ---------- */
1468 /* N            --> DIMENSION OF PROBLEM */
1469 /* X(N)         --> ESTIMATE TO A ROOT OF FCN */
1470 /* FCN          --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION */
1471 /*                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
1472 /*                       FCN:  R(N) --> R(1) */
1473 /* F            --> FUNCTION VALUE:  FCN(X) */
1474 /* G(N)         --> GRADIENT:  G(X) */
1475 /* TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X */
1476 /* FSCALE       --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION FCN */
1477 /* RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN */
1478 /* ANALTL       --> TOLERANCE FOR COMPARISON OF ESTIMATED AND */
1479 /*                  ANALYTICAL GRADIENTS */
1480 /* WRK1(N)      --> WORKSPACE */
1481 /* MSG         <--  MESSAGE OR ERROR CODE */
1482 /*                    ON OUTPUT: =-21, PROBABLE CODING ERROR OF GRADIENT */
1483 /* IPR          --> DEVICE TO WHICH TO SEND OUTPUT */
1484 /* IFN         <--> NUMBER OF FUNCTION EVALUATIONS */
1485
1486 /* COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO */
1487 /* ANALYTIC GRADIENT. */
1488
1489     /* Parameter adjustments */
1490     --wrk1;
1491     --typx;
1492     --typsiz;
1493     --g;
1494     --x;
1495
1496     /* Function Body */
1497     fstofd_(&c__1, &c__1, n, &x[1], (S_fp)fcn, f, &wrk1[1], &typx[1], rnf, &
1498             wrk, &c__1, ifn);
1499     ker = 0;
1500     i__1 = *n;
1501     for (i__ = 1; i__ <= i__1; ++i__) {
1502 /* Computing MAX */
1503         d__2 = abs(*f);
1504 /* Computing MAX */
1505         d__3 = (d__1 = x[i__], abs(d__1)), d__4 = typsiz[i__];
1506         gs = max(d__2,*fscale) / max(d__3,d__4);
1507 /* Computing MAX */
1508         d__3 = (d__1 = g[i__], abs(d__1));
1509         if ((d__2 = g[i__] - wrk1[i__], abs(d__2)) > max(d__3,gs) * *analtl) {
1510             ker = 1;
1511         }
1512 /* L5: */
1513     }
1514     if (ker == 0) {
1515         goto L20;
1516     }
1517     io___116.ciunit = *ipr;
1518     s_wsfe(&io___116);
1519     e_wsfe();
1520     io___117.ciunit = *ipr;
1521     s_wsfe(&io___117);
1522     i__1 = *n;
1523     for (i__ = 1; i__ <= i__1; ++i__) {
1524         do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
1525         do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
1526         do_fio(&c__1, (char *)&wrk1[i__], (ftnlen)sizeof(doublereal));
1527     }
1528     e_wsfe();
1529     *msg = -21;
1530 L20:
1531     return 0;
1532 } /* grdchk_ */
1533
1534 /*  ---------------------- */
1535 /*  |  H E S C H K       | */
1536 /*  ---------------------- */
1537 /* Subroutine */ int heschk_(integer *nr, integer *n, doublereal *x, S_fp fcn,
1538          S_fp grd, S_fp hsn, doublereal *f, doublereal *g, doublereal *a,
1539         doublereal *typsiz, doublereal *typx, doublereal *rnf, doublereal *
1540         analtl, integer *iagflg, doublereal *udiag, doublereal *wrk1,
1541         doublereal *wrk2, integer *msg, integer *ipr, integer *ifn)
1542 {
1543     /* Format strings */
1544     static char fmt_901[] = "(\002 HESCHK    PROBABLE ERROR IN CODING OF ANA"
1545             "LYTIC\002,\002 HESSIAN FUNCTION.\002/\002 HESCHK      ROW  CO"
1546             "L\002,14x,\002ANALYTIC\002,14x,\002(ESTIMATE)\002)";
1547     static char fmt_902[] = "(\002 HESCHK    \002,2i5,2x,e20.13,2x,\002(\002"
1548             ",e20.13,\002)\002)";
1549
1550     /* System generated locals */
1551     integer a_dim1, a_offset, i__1, i__2;
1552     doublereal d__1, d__2, d__3, d__4, d__5;
1553
1554     /* Builtin functions */
1555     integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
1556
1557     /* Local variables */
1558     integer i__, j;
1559     doublereal hs;
1560     integer im1, jp1, ker;
1561     extern /* Subroutine */ int sndofd_(integer *, integer *, doublereal *,
1562             S_fp, doublereal *, doublereal *, doublereal *, doublereal *,
1563             doublereal *, doublereal *, integer *), fstofd_(integer *,
1564             integer *, integer *, doublereal *, S_fp, doublereal *,
1565             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
1566              integer *);
1567
1568     /* Fortran I/O blocks */
1569     static cilist io___123 = { 0, 0, 0, fmt_901, 0 };
1570     static cilist io___125 = { 0, 0, 0, fmt_902, 0 };
1571     static cilist io___126 = { 0, 0, 0, fmt_902, 0 };
1572
1573
1574
1575 /* PURPOSE */
1576 /* ------- */
1577 /* CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN */
1578 /*  (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN */
1579 /*   HSN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF A) */
1580
1581 /* PARAMETERS */
1582 /* ---------- */
1583 /* NR           --> ROW DIMENSION OF MATRIX */
1584 /* N            --> DIMENSION OF PROBLEM */
1585 /* X(N)         --> ESTIMATE TO A ROOT OF FCN */
1586 /* FCN          --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION */
1587 /*                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
1588 /*                       FCN:  R(N) --> R(1) */
1589 /* GRD          --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN. */
1590 /*                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
1591 /* HSN          --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN. */
1592 /*                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE */
1593 /*                  HSN SHOULD FILL ONLY THE LOWER TRIANGULAR PART */
1594 /*                  AND DIAGONAL OF A */
1595 /* F            --> FUNCTION VALUE:  FCN(X) */
1596 /* G(N)        <--  GRADIENT:  G(X) */
1597 /* A(N,N)      <--  ON EXIT:  HESSIAN IN LOWER TRIANGULAR PART AND DIAG */
1598 /* TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X */
1599 /* RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN */
1600 /* ANALTL       --> TOLERANCE FOR COMPARISON OF ESTIMATED AND */
1601 /*                  ANALYTICAL GRADIENTS */
1602 /* IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED */
1603 /* UDIAG(N)     --> WORKSPACE */
1604 /* WRK1(N)      --> WORKSPACE */
1605 /* WRK2(N)      --> WORKSPACE */
1606 /* MSG         <--> MESSAGE OR ERROR CODE */
1607 /*                    ON INPUT : IF =1XX DO NOT COMPARE ANAL + EST HESS */
1608 /*                    ON OUTPUT: =-22, PROBABLE CODING ERROR OF HESSIAN */
1609 /* IPR          --> DEVICE TO WHICH TO SEND OUTPUT */
1610 /* IFN         <--> NUMBER OF FUNCTION EVALUTATIONS */
1611
1612 /* COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN. */
1613
1614     /* Parameter adjustments */
1615     a_dim1 = *nr;
1616     a_offset = 1 + a_dim1;
1617     a -= a_offset;
1618     --wrk2;
1619     --wrk1;
1620     --udiag;
1621     --typx;
1622     --typsiz;
1623     --g;
1624     --x;
1625
1626     /* Function Body */
1627     if (*iagflg == 1) {
1628         fstofd_(nr, n, n, &x[1], (S_fp)grd, &g[1], &a[a_offset], &typx[1],
1629                 rnf, &wrk1[1], &c__3, ifn);
1630     }
1631     if (*iagflg != 1) {
1632         sndofd_(nr, n, &x[1], (S_fp)fcn, f, &a[a_offset], &typx[1], rnf, &
1633                 wrk1[1], &wrk2[1], ifn);
1634     }
1635
1636     ker = 0;
1637
1638 /* COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART */
1639 /* AND DIAGONAL OF "A" TO UDIAG */
1640
1641     i__1 = *n;
1642     for (j = 1; j <= i__1; ++j) {
1643         udiag[j] = a[j + j * a_dim1];
1644         if (j == *n) {
1645             goto L30;
1646         }
1647         jp1 = j + 1;
1648         i__2 = *n;
1649         for (i__ = jp1; i__ <= i__2; ++i__) {
1650             a[j + i__ * a_dim1] = a[i__ + j * a_dim1];
1651 /* L25: */
1652         }
1653 L30:
1654         ;
1655     }
1656
1657 /* COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE */
1658 /* APPROXIMATION. */
1659
1660     i__1 = *n;
1661     for (i__ = 1; i__ <= i__1; ++i__) {
1662         i__2 = *n;
1663         for (j = i__; j <= i__2; ++j) {
1664             a[j + i__ * a_dim1] = 0.;
1665 /* L32: */
1666         }
1667     }
1668
1669     (*hsn)(nr, n, &x[1], &a[a_offset]);
1670     i__2 = *n;
1671     for (j = 1; j <= i__2; ++j) {
1672 /* Computing MAX */
1673         d__3 = (d__1 = g[j], abs(d__1));
1674 /* Computing MAX */
1675         d__4 = (d__2 = x[j], abs(d__2)), d__5 = typsiz[j];
1676         hs = max(d__3,1.) / max(d__4,d__5);
1677 /* Computing MAX */
1678         d__3 = (d__1 = udiag[j], abs(d__1));
1679         if ((d__2 = a[j + j * a_dim1] - udiag[j], abs(d__2)) > max(d__3,hs) *
1680                 *analtl) {
1681             ker = 1;
1682         }
1683         if (j == *n) {
1684             goto L40;
1685         }
1686         jp1 = j + 1;
1687         i__1 = *n;
1688         for (i__ = jp1; i__ <= i__1; ++i__) {
1689 /* Computing MAX */
1690             d__3 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
1691             if ((d__2 = a[i__ + j * a_dim1] - a[j + i__ * a_dim1], abs(d__2))
1692                     > max(d__3,hs) * *analtl) {
1693                 ker = 1;
1694             }
1695 /* L35: */
1696         }
1697 L40:
1698         ;
1699     }
1700
1701     if (ker == 0) {
1702         goto L90;
1703     }
1704     io___123.ciunit = *ipr;
1705     s_wsfe(&io___123);
1706     e_wsfe();
1707     i__2 = *n;
1708     for (i__ = 1; i__ <= i__2; ++i__) {
1709         if (i__ == 1) {
1710             goto L45;
1711         }
1712         im1 = i__ - 1;
1713         i__1 = im1;
1714         for (j = 1; j <= i__1; ++j) {
1715             io___125.ciunit = *ipr;
1716             s_wsfe(&io___125);
1717             do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
1718             do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
1719             do_fio(&c__1, (char *)&a[i__ + j * a_dim1], (ftnlen)sizeof(
1720                     doublereal));
1721             do_fio(&c__1, (char *)&a[j + i__ * a_dim1], (ftnlen)sizeof(
1722                     doublereal));
1723             e_wsfe();
1724 /* L43: */
1725         }
1726 L45:
1727         io___126.ciunit = *ipr;
1728         s_wsfe(&io___126);
1729         do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
1730         do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(integer));
1731         do_fio(&c__1, (char *)&a[i__ + i__ * a_dim1], (ftnlen)sizeof(
1732                 doublereal));
1733         do_fio(&c__1, (char *)&udiag[i__], (ftnlen)sizeof(doublereal));
1734         e_wsfe();
1735 /* L50: */
1736     }
1737     *msg = -22;
1738 /*     ENDIF */
1739 L90:
1740     return 0;
1741 } /* heschk_ */
1742
1743 /*  ----------------------- */
1744 /*  |    M C H E P S    | */
1745 /*  ----------------------- */
1746 /* Subroutine */ int mcheps_(doublereal *eps)
1747 {
1748     doublereal temp, temp1;
1749
1750
1751 /* PURPOSE: */
1752 /*   COMPUTE MACHINE PRECISION */
1753
1754 /* ------------------------------------------------------------------------- */
1755
1756 /* PARAMETERS: */
1757
1758 /*   EPS <-- MACHINE PRECISION */
1759
1760     temp = 1.;
1761 L20:
1762     temp /= 2.;
1763     temp1 = temp + 1.;
1764     if (temp1 != 1.) {
1765         goto L20;
1766     }
1767     *eps = temp * 2.;
1768     return 0;
1769 } /* mcheps_ */
1770
1771
1772 doublereal twonrm_(integer *n, doublereal *v)
1773 {
1774     /* System generated locals */
1775     doublereal ret_val;
1776
1777     /* Builtin functions */
1778     double sqrt(doublereal);
1779
1780     /* Local variables */
1781     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
1782             integer *);
1783     doublereal temp;
1784
1785
1786 /* PURPOSE: */
1787 /*   COMPUTER L-2 NORM */
1788
1789 /* -------------------------------------------------------------------------- */
1790
1791 /* PARAMETER: */
1792
1793 /*   N       --> DIMENSION OF PROBLEM */
1794 /*   V(N)    --> VECTOR WHICH L-2 NORM IS EVALUATED */
1795
1796     /* Parameter adjustments */
1797     --v;
1798
1799     /* Function Body */
1800     temp = ddot_(n, &v[1], &c__1, &v[1], &c__1);
1801     ret_val = sqrt(temp);
1802     return ret_val;
1803 } /* twonrm_ */
1804
1805 /*  ---------------------- */
1806 /*  |  L N S R C H       | */
1807 /*  ---------------------- */
1808 /* Subroutine */ int lnsrch_(integer *nr, integer *n, doublereal *x,
1809         doublereal *f, doublereal *g, doublereal *p, doublereal *xpls,
1810         doublereal *fpls, logical *mxtake, integer *iretcd, doublereal *
1811         stepmx, doublereal *steptl, doublereal *typx, S_fp fcn, doublereal *
1812         w2, integer *nfcnt)
1813 {
1814     /* System generated locals */
1815     integer i__1;
1816     doublereal d__1, d__2;
1817
1818     /* Builtin functions */
1819     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
1820
1821     /* Local variables */
1822     doublereal a, b;
1823     integer i__, k;
1824     doublereal t1, t2, t3, scl, rln, sln, slp, tmp, disc;
1825     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
1826             integer *);
1827     doublereal temp, temp1, temp2, alpha;
1828     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
1829             integer *);
1830     doublereal pfpls, almbda, plmbda, tlmbda, almbmn;
1831
1832
1833 /*     THE ALPHA CONDITION ONLY LINE SEARCH */
1834
1835 /* PURPOSE */
1836 /* ------- */
1837 /* FIND A NEXT NEWTON ITERATE BY LINE SEARCH. */
1838
1839 /* PARAMETERS */
1840 /* ---------- */
1841 /* N            --> DIMENSION OF PROBLEM */
1842 /* X(N)         --> OLD ITERATE:   X[K-1] */
1843 /* F            --> FUNCTION VALUE AT OLD ITERATE, F(X) */
1844 /* G(N)         --> GRADIENT AT OLD ITERATE, G(X), OR APPROXIMATE */
1845 /* P(N)         --> NON-ZERO NEWTON STEP */
1846 /* XPLS(N)     <--  NEW ITERATE X[K] */
1847 /* FPLS        <--  FUNCTION VALUE AT NEW ITERATE, F(XPLS) */
1848 /* FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
1849 /* IRETCD      <--  RETURN CODE */
1850 /* MXTAKE      <--  BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED */
1851 /* STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE */
1852 /* STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES */
1853 /*                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM */
1854 /* TYPX(N)            --> DIAGONAL SCALING MATRIX FOR X (NOT IN UNCMIN) */
1855 /* IPR          --> DEVICE TO WHICH TO SEND OUTPUT */
1856 /* W2         --> WORKING SPACE */
1857 /* NFCNT      <--> NUMBER OF FUNCTION EVALUTIONS */
1858
1859 /* INTERNAL VARIABLES */
1860 /* ------------------ */
1861 /* SLN              NEWTON LENGTH */
1862 /* RLN              RELATIVE LENGTH OF NEWTON STEP */
1863
1864
1865     /* Parameter adjustments */
1866     --w2;
1867     --typx;
1868     --xpls;
1869     --p;
1870     --g;
1871     --x;
1872
1873     /* Function Body */
1874     *mxtake = FALSE_;
1875     *iretcd = 2;
1876     alpha = 1e-4;
1877 /* $    WRITE(IPR,954) */
1878 /* $    WRITE(IPR,955) (P(I),I=1,N) */
1879     tmp = 0.;
1880     i__1 = *n;
1881     for (i__ = 1; i__ <= i__1; ++i__) {
1882         tmp += p[i__] * p[i__];
1883 /* L5: */
1884     }
1885     sln = sqrt(tmp);
1886     if (sln <= *stepmx) {
1887         goto L10;
1888     }
1889
1890 /* NEWTON STEP LONGER THAN MAXIMUM ALLOWED */
1891     scl = *stepmx / sln;
1892     dscal_(n, &scl, &p[1], &c__1);
1893     sln = *stepmx;
1894 /* $     WRITE(IPR,954) */
1895 /* $     WRITE(IPR,955) (P(I),I=1,N) */
1896 L10:
1897     slp = ddot_(n, &g[1], &c__1, &p[1], &c__1);
1898     rln = 0.;
1899     i__1 = *n;
1900     for (i__ = 1; i__ <= i__1; ++i__) {
1901         temp = 1.;
1902         temp1 = (d__1 = x[i__], abs(d__1));
1903         temp2 = max(temp1,temp);
1904         temp1 = (d__1 = p[i__], abs(d__1));
1905 /* Computing MAX */
1906         d__1 = rln, d__2 = temp1 / temp2;
1907         rln = max(d__1,d__2);
1908 /* L15: */
1909     }
1910     almbmn = *steptl / rln;
1911     almbda = 1.;
1912 /* $    WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL */
1913
1914 /* LOOP */
1915 /* CHECK IF NEW ITERATE SATISFACTORY.  GENERATE NEW LAMBDA IF NECESSARY. */
1916
1917 L100:
1918     if (*iretcd < 2) {
1919         return 0;
1920     }
1921     i__1 = *n;
1922     for (i__ = 1; i__ <= i__1; ++i__) {
1923         xpls[i__] = x[i__] + almbda * p[i__];
1924 /* L105: */
1925     }
1926     i__1 = *n;
1927     for (k = 1; k <= i__1; ++k) {
1928         w2[k] = xpls[k] * typx[k];
1929 /* L101: */
1930     }
1931     (*fcn)(n, &w2[1], fpls);
1932     ++(*nfcnt);
1933 /* $    WRITE(IPR,956) ALMBDA */
1934 /* $    WRITE(IPR,951) */
1935 /* $    WRITE(IPR,955) (XPLS(I),I=1,N) */
1936 /* $    WRITE(IPR,953) FPLS */
1937     if (*fpls > *f + slp * alpha * almbda) {
1938         goto L130;
1939     }
1940 /*     IF(FPLS.LE. F+SLP*1.E-4*ALMBDA) */
1941 /*     THEN */
1942
1943 /* SOLUTION FOUND */
1944
1945     *iretcd = 0;
1946     if (almbda == 1. && sln > *stepmx * .99) {
1947         *mxtake = TRUE_;
1948     }
1949     goto L100;
1950
1951 /* SOLUTION NOT (YET) FOUND */
1952
1953 /*     ELSE */
1954 L130:
1955     if (almbda >= almbmn) {
1956         goto L140;
1957     }
1958 /*       IF(ALMBDA .LT. ALMBMN) */
1959 /*       THEN */
1960
1961 /* NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X */
1962
1963     *iretcd = 1;
1964     goto L100;
1965 /*       ELSE */
1966
1967 /* CALCULATE NEW LAMBDA */
1968
1969 L140:
1970     if (almbda != 1.) {
1971         goto L150;
1972     }
1973 /*         IF(ALMBDA.EQ.1.0) */
1974 /*         THEN */
1975
1976 /* FIRST BACKTRACK: QUADRATIC FIT */
1977
1978     tlmbda = -slp / ((*fpls - *f - slp) * 2.);
1979     goto L170;
1980 /*         ELSE */
1981
1982 /* ALL SUBSEQUENT BACKTRACKS: CUBIC FIT */
1983
1984 L150:
1985     t1 = *fpls - *f - almbda * slp;
1986     t2 = pfpls - *f - plmbda * slp;
1987     t3 = 1. / (almbda - plmbda);
1988     a = t3 * (t1 / (almbda * almbda) - t2 / (plmbda * plmbda));
1989     b = t3 * (t2 * almbda / (plmbda * plmbda) - t1 * plmbda / (almbda *
1990             almbda));
1991     disc = b * b - a * 3.f * slp;
1992     if (disc <= b * b) {
1993         goto L160;
1994     }
1995 /*           IF(DISC.GT. B*B) */
1996 /*           THEN */
1997
1998 /* ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM */
1999
2000     tlmbda = (-b + d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
2001     goto L165;
2002 /*           ELSE */
2003
2004 /* BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM */
2005
2006 L160:
2007     tlmbda = (-b - d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
2008 /*           ENDIF */
2009 L165:
2010     if (tlmbda > almbda * .5) {
2011         tlmbda = almbda * .5;
2012     }
2013 /*         ENDIF */
2014 L170:
2015     plmbda = almbda;
2016     pfpls = *fpls;
2017     if (tlmbda >= almbda * .1) {
2018         goto L180;
2019     }
2020 /*         IF(TLMBDA.LT.ALMBDA/10.) */
2021 /*         THEN */
2022     almbda *= .1f;
2023     goto L190;
2024 /*         ELSE */
2025 L180:
2026     almbda = tlmbda;
2027 /*         ENDIF */
2028 /*       ENDIF */
2029 /*     ENDIF */
2030 L190:
2031     goto L100;
2032 /* L956: */
2033 /* L951: */
2034 /* L952: */
2035 /* L953: */
2036 /* L954: */
2037 /* L955: */
2038 } /* lnsrch_ */
2039
2040 /*  ---------------------- */
2041 /*  |       Z H Z        | */
2042 /*  ---------------------- */
2043 /* Subroutine */ int zhz_(integer *nr, integer *n, doublereal *y, doublereal *
2044         h__, doublereal *u, doublereal *t)
2045 {
2046     /* System generated locals */
2047     integer h_dim1, h_offset, i__1, i__2;
2048     doublereal d__1;
2049
2050     /* Builtin functions */
2051     double sqrt(doublereal);
2052
2053     /* Local variables */
2054     doublereal d__;
2055     integer i__, j;
2056     doublereal s, sgn;
2057     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2058             integer *);
2059     doublereal temp1, temp2;
2060     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
2061             doublereal *, integer *);
2062     doublereal ynorm;
2063
2064
2065 /* PURPOSE: */
2066 /*   COMPUTE QTHQ(N,N) AND ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND */
2067 /*   FIRST N-1 COLUMNS OF QTHQ */
2068
2069 /* --------------------------------------------------------------------------- */
2070
2071 /* PARAMETERS: */
2072
2073 /*   NR            --> ROW DIMENSION OF MATRIX */
2074 /*   N       --> DIMENSION OF PROBLEM */
2075 /*   Y(N)    --> FIRST BASIS IN Q */
2076 /*   H(N,N)     <--> ON INPUT : HESSIAN */
2077 /*               ON OUTPUT: QTHQ (ZTHZ) */
2078 /*   U(N)       <--  VECTOR TO FORM Q AND Z */
2079 /*   T(N)        --> WORKSPACE */
2080
2081
2082 /*     U=Y+SGN(Y(N))||Y||E(N) */
2083     /* Parameter adjustments */
2084     --t;
2085     --u;
2086     h_dim1 = *nr;
2087     h_offset = 1 + h_dim1;
2088     h__ -= h_offset;
2089     --y;
2090
2091     /* Function Body */
2092     if (y[*n] != 0.) {
2093         sgn = y[*n] / (d__1 = y[*n], abs(d__1));
2094     } else {
2095         sgn = 1.;
2096     }
2097     ynorm = ddot_(n, &y[1], &c__1, &y[1], &c__1);
2098     ynorm = sqrt(ynorm);
2099     u[*n] = y[*n] + sgn * ynorm;
2100     i__1 = *n - 1;
2101     dcopy_(&i__1, &y[1], &c__1, &u[1], &c__1);
2102
2103 /*     D=UTU/2 */
2104     d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
2105     d__ /= 2.;
2106
2107 /*     T=2HU/UTU */
2108     i__1 = *n;
2109     for (i__ = 1; i__ <= i__1; ++i__) {
2110         t[i__] = 0.;
2111         i__2 = *n;
2112         for (j = 1; j <= i__2; ++j) {
2113             t[i__] = h__[i__ + j * h_dim1] * u[j] + t[i__];
2114 /* L30: */
2115         }
2116         t[i__] /= d__;
2117 /* L40: */
2118     }
2119
2120 /*     S=4UHU/(UTU)**2 */
2121     s = ddot_(n, &u[1], &c__1, &t[1], &c__1);
2122     s /= d__;
2123
2124 /*     COMPUTE QTHQ (ZTHZ) */
2125     i__1 = *n;
2126     for (i__ = 1; i__ <= i__1; ++i__) {
2127         temp1 = u[i__];
2128         i__2 = *n;
2129         for (j = 1; j <= i__2; ++j) {
2130             temp2 = u[j];
2131             h__[i__ + j * h_dim1] = h__[i__ + j * h_dim1] - u[i__] * t[j] - t[
2132                     i__] * u[j] + u[i__] * u[j] * s;
2133 /* L60: */
2134         }
2135 /* L70: */
2136     }
2137     return 0;
2138 } /* zhz_ */
2139
2140 /*  ---------------------- */
2141 /*  |    S O L V E W     | */
2142 /*  ---------------------- */
2143 /* Subroutine */ int solvew_(integer *nr, integer *n, doublereal *al,
2144         doublereal *u, doublereal *w, doublereal *b)
2145 {
2146     /* System generated locals */
2147     integer al_dim1, al_offset, i__1, i__2;
2148
2149     /* Local variables */
2150     doublereal d__;
2151     integer i__, j;
2152     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2153             integer *);
2154     extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
2155             doublereal *, doublereal *);
2156
2157
2158 /* PURPOSE: */
2159 /*   SOLVE L*W=ZT*V */
2160
2161 /* ---------------------------------------------------------------------- */
2162
2163 /* PARAMETERS: */
2164 /*   NR            --> ROW DIMENSION OF MATRIX */
2165 /*   N       --> DIMENSION OF PROBLEM */
2166 /*   AL(N-1,N-1)   --> LOWER TRIAGULAR MATRIX */
2167 /*   U(N)    --> VECTOR TO FORM Z */
2168 /*   W(N)    --> ON INPUT : VECTOR V IN SYSTEM OF LINEAR EQUATIONS */
2169 /*               ON OUTPUT: SOLUTION OF SYSTEM OF LINEAR EQUATIONS */
2170 /*   B(N)    --> WORKSPACE TO STORE ZT*V */
2171
2172 /*     FORM ZT*V (STORED IN B) */
2173     /* Parameter adjustments */
2174     al_dim1 = *nr;
2175     al_offset = 1 + al_dim1;
2176     al -= al_offset;
2177     --b;
2178     --w;
2179     --u;
2180
2181     /* Function Body */
2182     d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
2183     d__ /= 2.;
2184     i__1 = *n - 1;
2185     for (i__ = 1; i__ <= i__1; ++i__) {
2186         b[i__] = 0.;
2187         i__2 = *n;
2188         for (j = 1; j <= i__2; ++j) {
2189             b[i__] += u[j] * u[i__] * w[j] / d__;
2190 /* L15: */
2191         }
2192         b[i__] = w[i__] - b[i__];
2193 /* L20: */
2194     }
2195
2196 /*     SOLVE LW=ZT*V */
2197     i__1 = *n - 1;
2198     forslv_(nr, &i__1, &al[al_offset], &w[1], &b[1]);
2199     return 0;
2200 } /* solvew_ */
2201
2202 /*  ---------------------- */
2203 /*  |     D S T A R      | */
2204 /*  ---------------------- */
2205 /* Subroutine */ int dstar_(integer *nr, integer *n, doublereal *u,
2206         doublereal *s, doublereal *w1, doublereal *w2, doublereal *w3,
2207         doublereal *sigma, doublereal *al, doublereal *d__)
2208 {
2209     /* System generated locals */
2210     integer al_dim1, al_offset, i__1;
2211
2212     /* Local variables */
2213     integer i__;
2214     doublereal utt, utu;
2215     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2216             integer *);
2217     doublereal temp;
2218     extern /* Subroutine */ int bakslv_(integer *, integer *, doublereal *,
2219             doublereal *, doublereal *);
2220
2221
2222 /* PURPOSE: */
2223 /*   COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) */
2224
2225 /* ------------------------------------------------------------------------ */
2226
2227 /* PARAMETERS: */
2228 /*   NR            --> ROW DIMENSION OF MATRIX */
2229 /*   N       --> DIMENSION OF PROBLEM */
2230 /*   U(N)    --> VECTOR TO FORM Z */
2231 /*   S(N)    --> PREVIOUS STEP */
2232 /*   W1(N)   --> L**-1*ZT*A, WHERE A IS DESCRIBED IN SUBROUTINE SLVMDL */
2233 /*   W2(N)   --> L**-1*ZT*SH, WHERE H IS CURRENT HESSIAN */
2234 /*   W3(N)   --> L**-1*ZT*G, WHERE G IS CURRENT GRADIENT */
2235 /*   SIGMA   --> SOLUTION FOR REDUCED ONE VARIABLE MODEL */
2236 /*   AL(N-1,N-1) --> LOWER TRIANGULAR MATRIX L */
2237 /*   D(N)    --> TENSOR STEP */
2238
2239     /* Parameter adjustments */
2240     al_dim1 = *nr;
2241     al_offset = 1 + al_dim1;
2242     al -= al_offset;
2243     --d__;
2244     --w3;
2245     --w2;
2246     --w1;
2247     --s;
2248     --u;
2249
2250     /* Function Body */
2251     if (*n == 1) {
2252         d__[1] = *sigma * s[1];
2253     } else {
2254
2255 /*     COMPUTE T(SIGMA)=-(ZTHZ)*ZT*(G+SIGMA*SH+SIGMA**2*A/2) (STORED IN D) */
2256         i__1 = *n - 1;
2257         for (i__ = 1; i__ <= i__1; ++i__) {
2258             w2[i__] = w3[i__] + *sigma * w2[i__] + w1[i__] * .5 * *sigma * *
2259                     sigma;
2260 /* L10: */
2261         }
2262         i__1 = *n - 1;
2263         bakslv_(nr, &i__1, &al[al_offset], &d__[1], &w2[1]);
2264         d__[*n] = 0.;
2265
2266 /*     COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) */
2267         utu = ddot_(n, &u[1], &c__1, &u[1], &c__1);
2268         utt = ddot_(n, &u[1], &c__1, &d__[1], &c__1);
2269         temp = utt / utu;
2270         i__1 = *n;
2271         for (i__ = 1; i__ <= i__1; ++i__) {
2272             d__[i__] = *sigma * s[i__] - (d__[i__] - u[i__] * 2. * temp);
2273 /* L50: */
2274         }
2275     }
2276     return 0;
2277 } /* dstar_ */
2278
2279 /*  ---------------------- */
2280 /*  |     M K M D L      | */
2281 /*  ---------------------- */
2282 /* Subroutine */ int mkmdl_(integer *nr, integer *n, doublereal *f,
2283         doublereal *fp, doublereal *g, doublereal *gp, doublereal *s,
2284         doublereal *h__, doublereal *alpha, doublereal *beta, doublereal *sh,
2285         doublereal *a)
2286 {
2287     /* System generated locals */
2288     integer h_dim1, h_offset, i__1, i__2;
2289
2290     /* Local variables */
2291     integer i__, j;
2292     doublereal b1, b2, gs, gps, shs, sts;
2293     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2294             integer *);
2295
2296
2297 /* PURPOSE: */
2298 /*   FORM TENSOR MODEL */
2299
2300 /* ----------------------------------------------------------------------- */
2301
2302 /* PARAMETERS: */
2303 /*   NR            --> ROW DIMENSION OF MATRIX */
2304 /*   N       --> DIMENSION OF PROBLEM */
2305 /*   F       --> CURRENT FUNCTION VALUE */
2306 /*   FP            --> PREVIOUS FUNCTION VALUE */
2307 /*   G(N)    --> CURRENT GRADIENT */
2308 /*   GP(N)   --> PREVIOUS GRADIENT */
2309 /*   S(N)    --> STEP TO PREVIOUS POINT */
2310 /*   H(N,N)  --> HESSIAN */
2311 /*   ALPHA      <--  SCALAR TO FORM 3RD ORDER TERM OF TENSOR MODEL */
2312 /*   BETA       <--  SCALAR TO FORM 4TH ORDER TERM OF TENSOR MODEL */
2313 /*   SH(N)      <--  SH */
2314 /*   A(N)       <--  A=2*(GP-G-SH-S*BETA/(6*STS)) */
2315
2316
2317 /*     COMPUTE SH */
2318     /* Parameter adjustments */
2319     --a;
2320     --sh;
2321     h_dim1 = *nr;
2322     h_offset = 1 + h_dim1;
2323     h__ -= h_offset;
2324     --s;
2325     --gp;
2326     --g;
2327
2328     /* Function Body */
2329     i__1 = *n;
2330     for (i__ = 1; i__ <= i__1; ++i__) {
2331         sh[i__] = 0.;
2332         i__2 = *n;
2333         for (j = 1; j <= i__2; ++j) {
2334             sh[i__] += s[j] * h__[j + i__ * h_dim1];
2335 /* L10: */
2336         }
2337 /* L20: */
2338     }
2339     gs = ddot_(n, &g[1], &c__1, &s[1], &c__1);
2340     gps = ddot_(n, &gp[1], &c__1, &s[1], &c__1);
2341     shs = ddot_(n, &sh[1], &c__1, &s[1], &c__1);
2342     b1 = gps - gs - shs;
2343     b2 = *fp - *f - gs - shs * .5;
2344     *alpha = b2 * 24. - b1 * 6.;
2345     *beta = b1 * 24. - b2 * 72.;
2346
2347 /*     COMPUTE A */
2348     sts = ddot_(n, &s[1], &c__1, &s[1], &c__1);
2349     i__1 = *n;
2350     for (i__ = 1; i__ <= i__1; ++i__) {
2351         a[i__] = (gp[i__] - g[i__] - sh[i__] - s[i__] * *beta / (sts * 6.)) *
2352                 2.;
2353 /* L50: */
2354     }
2355     return 0;
2356 } /* mkmdl_ */
2357
2358 /*  ---------------------- */
2359 /*  |     S I G M A      | */
2360 /*  ---------------------- */
2361 /* Subroutine */ int sigma_(doublereal *sgstar, doublereal *a, doublereal *b,
2362         doublereal *c__, doublereal *d__)
2363 {
2364     doublereal s1, s2, s3;
2365     extern /* Subroutine */ int roots_(doublereal *, doublereal *, doublereal
2366             *, doublereal *, doublereal *, doublereal *, doublereal *),
2367             sortrt_(doublereal *, doublereal *, doublereal *);
2368
2369
2370 /* PURPOSE: */
2371 /*   COMPUTE DESIRABLE ROOT OF REDUCED ONE VARIABLE EQUATION */
2372
2373 /* ------------------------------------------------------------------------- */
2374
2375 /* PARAMETERS: */
2376 /*   SGSTAR     --> DESIRABLE ROOT */
2377 /*   A       --> COEFFICIENT OF 3RD ORDER TERM */
2378 /*   B       --> COEFFICIENT OF 2ND ORDER TERM */
2379 /*   C       --> COEFFICIENT OF 1ST ORDER TERM */
2380 /*   D       --> COEFFICIENT OF CONSTANT TERM */
2381
2382
2383 /*     COMPUTE ALL THREE ROOTS */
2384     roots_(&s1, &s2, &s3, a, b, c__, d__);
2385
2386 /*     SORT ROOTS */
2387     sortrt_(&s1, &s2, &s3);
2388
2389 /*     CHOOSE DESIRABLE ROOT */
2390     if (*a > 0.) {
2391         *sgstar = s3;
2392         if (s2 >= 0.) {
2393             *sgstar = s1;
2394         }
2395     } else {
2396 /*       IF  A  <  0  THEN */
2397         *sgstar = s2;
2398         if (s1 > 0. || s3 < 0.) {
2399             if (s1 > 0.) {
2400                 *sgstar = s1;
2401             } else {
2402                 *sgstar = s3;
2403             }
2404             *a = 0.;
2405         }
2406     }
2407     return 0;
2408 } /* sigma_ */
2409
2410 /*  ---------------------- */
2411 /*  |     R O O T S      | */
2412 /*  ---------------------- */
2413 /* Subroutine */ int roots_(doublereal *s1, doublereal *s2, doublereal *s3,
2414         doublereal *a, doublereal *b, doublereal *c__, doublereal *d__)
2415 {
2416     /* System generated locals */
2417     doublereal d__1;
2418
2419     /* Builtin functions */
2420     double sqrt(doublereal), pow_dd(doublereal *, doublereal *), acos(
2421             doublereal), cos(doublereal);
2422
2423     /* Local variables */
2424     doublereal q, r__, s, t, v, a1, a2, a3, pi, temp, theta;
2425
2426
2427 /* PURPOSE: */
2428 /*   COMPUTE ROOT(S) OF 3RD ORDER EQUATION */
2429
2430 /* --------------------------------------------------------------------------- */
2431
2432 /* PARAMETERS: */
2433 /*   S1             <--  ROOT   (IF THREE ROOTS ARE */
2434 /*   S2             <--  ROOT    EQUAL, THEN S1=S2=S3) */
2435 /*   S3             <--  ROOT */
2436 /*   A       --> COEFFICIENT OF 3RD ORDER TERM */
2437 /*   B       --> COEFFICIENT OF 2ND ORDER TERM */
2438 /*   C       --> COEFFICIENT OF 1ST ORDER TERM */
2439 /*   D       --> COEFFICIENT OF CONSTANT TERM */
2440
2441 /*     SET VALUE OF PI */
2442     pi = 3.141592653589793;
2443     a1 = *b / *a;
2444     a2 = *c__ / *a;
2445     a3 = *d__ / *a;
2446     q = (a2 * 3. - a1 * a1) / 9.;
2447     r__ = (a1 * 9. * a2 - a3 * 27. - a1 * 2. * a1 * a1) / 54.;
2448     v = q * q * q + r__ * r__;
2449     if (v > 0.) {
2450         s = r__ + sqrt(v);
2451         t = r__ - sqrt(v);
2452         if (t < 0.) {
2453             d__1 = -t;
2454             t = -pow_dd(&d__1, &c_b250);
2455         } else {
2456             t = pow_dd(&t, &c_b250);
2457         }
2458         if (s < 0.) {
2459             d__1 = -s;
2460             s = -pow_dd(&d__1, &c_b250);
2461         } else {
2462             s = pow_dd(&s, &c_b250);
2463         }
2464         *s1 = s + t - a1 / 3.;
2465         *s3 = *s1;
2466         *s2 = *s1;
2467     } else {
2468         temp = r__ / sqrt(-pow_dd(&q, &c_b384));
2469         theta = acos(temp);
2470         theta /= 3.;
2471         temp = sqrt(-q) * 2.;
2472         *s1 = temp * cos(theta) - a1 / 3.;
2473         *s2 = temp * cos(theta + pi * 2. / 3.) - a1 / 3.;
2474         *s3 = temp * cos(theta + pi * 4. / 3.) - a1 / 3.;
2475     }
2476     return 0;
2477 } /* roots_ */
2478
2479 /*  ---------------------- */
2480 /*  | S O R T R T         | */
2481 /*  ---------------------- */
2482 /* Subroutine */ int sortrt_(doublereal *s1, doublereal *s2, doublereal *s3)
2483 {
2484     doublereal t;
2485
2486
2487 /* PURPOSE: */
2488 /*   SORT ROOTS INTO ASCENDING ORDER */
2489
2490 /* ----------------------------------------------------------------------------- */
2491
2492 /* PARAMETERS: */
2493 /*   S1             <--> ROOT */
2494 /*   S2             <--> ROOT */
2495 /*   S3             <--> ROOT */
2496
2497     if (*s1 > *s2) {
2498         t = *s1;
2499         *s1 = *s2;
2500         *s2 = t;
2501     }
2502     if (*s2 > *s3) {
2503         t = *s2;
2504         *s2 = *s3;
2505         *s3 = t;
2506     }
2507     if (*s1 > *s2) {
2508         t = *s1;
2509         *s1 = *s2;
2510         *s2 = t;
2511     }
2512     return 0;
2513 } /* sortrt_ */
2514
2515 /*  ---------------------- */
2516 /*  |  F S T O F D       | */
2517 /*  ---------------------- */
2518 /* Subroutine */ int fstofd_(integer *nr, integer *m, integer *n, doublereal *
2519         xpls, S_fp fcn, doublereal *fpls, doublereal *a, doublereal *typx,
2520         doublereal *rnoise, doublereal *fhat, integer *icase, integer *nfcnt)
2521 {
2522     /* System generated locals */
2523     integer a_dim1, a_offset, i__1, i__2;
2524     doublereal d__1, d__2;
2525
2526     /* Builtin functions */
2527     double sqrt(doublereal);
2528
2529     /* Local variables */
2530     integer i__, j, jp1, nm1;
2531     doublereal xtmpj, stepsz;
2532
2533 /* PURPOSE */
2534 /* ------- */
2535 /* FIND FIRST ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" TO THE */
2536 /* FIRST DERIVATIVE OF THE FUNCTION DEFINED BY THE SUBPROGRAM "FNAME" */
2537 /* EVALUATED AT THE NEW ITERATE "XPLS". */
2538
2539
2540 /* FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE: */
2541 /* 1) THE FIRST DERIVATIVE (GRADIENT) OF THE OPTIMIZATION FUNCTION "FCN */
2542 /*    ANALYTIC USER ROUTINE HAS BEEN SUPPLIED; */
2543 /* 2) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION */
2544 /*    IF NO ANALYTIC USER ROUTINE HAS BEEN SUPPLIED FOR THE HESSIAN BUT */
2545 /*    ONE HAS BEEN SUPPLIED FOR THE GRADIENT ("FCN") AND IF THE */
2546 /*    OPTIMIZATION FUNCTION IS INEXPENSIVE TO EVALUATE */
2547
2548 /* NOTE */
2549 /* ---- */
2550 /* _M=1 (OPTIMIZATION) ALGORITHM ESTIMATES THE GRADIENT OF THE FUNCTION */
2551 /*      (FCN).   FCN(X) # F: R(N)-->R(1) */
2552 /* _M=N (SYSTEMS) ALGORITHM ESTIMATES THE JACOBIAN OF THE FUNCTION */
2553 /*      FCN(X) # F: R(N)-->R(N). */
2554 /* _M=N (OPTIMIZATION) ALGORITHM ESTIMATES THE HESSIAN OF THE OPTIMIZATIO */
2555 /*      FUNCTION, WHERE THE HESSIAN IS THE FIRST DERIVATIVE OF "FCN" */
2556
2557 /* PARAMETERS */
2558 /* ---------- */
2559 /* NR           --> ROW DIMENSION OF MATRIX */
2560 /* M            --> NUMBER OF ROWS IN A */
2561 /* N            --> NUMBER OF COLUMNS IN A; DIMENSION OF PROBLEM */
2562 /* XPLS(N)      --> NEW ITERATE:  X[K] */
2563 /* FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
2564 /* FPLS(M)      --> _M=1 (OPTIMIZATION) FUNCTION VALUE AT NEW ITERATE: */
2565 /*                       FCN(XPLS) */
2566 /*                  _M=N (OPTIMIZATION) VALUE OF FIRST DERIVATIVE */
2567 /*                       (GRADIENT) GIVEN BY USER FUNCTION FCN */
2568 /*                  _M=N (SYSTEMS)  FUNCTION VALUE OF ASSOCIATED */
2569 /*                       MINIMIZATION FUNCTION */
2570 /* A(NR,N)     <--  FINITE DIFFERENCE APPROXIMATION (SEE NOTE).  ONLY */
2571 /*                  LOWER TRIANGULAR MATRIX AND DIAGONAL ARE RETURNED */
2572 /* RNOISE       --> RELATIVE NOISE IN FCN [F(X)] */
2573 /* FHAT(M)      --> WORKSPACE */
2574 /* ICASE        --> =1 OPTIMIZATION (GRADIENT) */
2575 /*                  =2 SYSTEMS */
2576 /*                  =3 OPTIMIZATION (HESSIAN) */
2577 /* NFCNT       <--> NUMBER OF FUNCTION EVALUTIONS */
2578
2579 /* INTERNAL VARIABLES */
2580 /* ------------------ */
2581 /* STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION */
2582
2583
2584 /* FIND J-TH COLUMN OF A */
2585 /* EACH COLUMN IS DERIVATIVE OF F(FCN) WITH RESPECT TO XPLS(J) */
2586
2587     /* Parameter adjustments */
2588     --fhat;
2589     --fpls;
2590     --typx;
2591     a_dim1 = *nr;
2592     a_offset = 1 + a_dim1;
2593     a -= a_offset;
2594     --xpls;
2595
2596     /* Function Body */
2597     i__1 = *n;
2598     for (j = 1; j <= i__1; ++j) {
2599         xtmpj = xpls[j];
2600 /* Computing MAX */
2601         d__2 = (d__1 = xpls[j], abs(d__1));
2602         stepsz = sqrt(*rnoise) * max(d__2,1.);
2603         xpls[j] = xtmpj + stepsz;
2604         (*fcn)(n, &xpls[1], &fhat[1]);
2605         ++(*nfcnt);
2606         xpls[j] = xtmpj;
2607         i__2 = *m;
2608         for (i__ = 1; i__ <= i__2; ++i__) {
2609             a[i__ + j * a_dim1] = (fhat[i__] - fpls[i__]) / stepsz;
2610             a[i__ + j * a_dim1] *= typx[j];
2611 /* L20: */
2612         }
2613 /* L30: */
2614     }
2615     if (*icase != 3) {
2616         return 0;
2617     }
2618
2619 /* IF COMPUTING HESSIAN, A MUST BE SYMMETRIC */
2620
2621     if (*n == 1) {
2622         return 0;
2623     }
2624     nm1 = *n - 1;
2625     i__1 = nm1;
2626     for (j = 1; j <= i__1; ++j) {
2627         jp1 = j + 1;
2628         i__2 = *m;
2629         for (i__ = jp1; i__ <= i__2; ++i__) {
2630             a[i__ + j * a_dim1] = (a[i__ + j * a_dim1] + a[j + i__ * a_dim1])
2631                     / 2.f;
2632 /* L40: */
2633         }
2634 /* L50: */
2635     }
2636     return 0;
2637 } /* fstofd_ */
2638
2639 /*  ---------------------- */
2640 /*  |  S N D O F D       | */
2641 /*  ---------------------- */
2642 /* Subroutine */ int sndofd_(integer *nr, integer *n, doublereal *xpls, S_fp
2643         fcn, doublereal *fpls, doublereal *a, doublereal *typx, doublereal *
2644         rnoise, doublereal *stepsz, doublereal *anbr, integer *nfcnt)
2645 {
2646     /* System generated locals */
2647     integer a_dim1, a_offset, i__1, i__2;
2648     doublereal d__1, d__2;
2649
2650     /* Builtin functions */
2651     double pow_dd(doublereal *, doublereal *);
2652
2653     /* Local variables */
2654     integer i__, j, ip1;
2655     doublereal ov3, fhat, xtmpi, xtmpj;
2656
2657 /* PURPOSE */
2658 /* ------- */
2659 /* FIND SECOND ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A" */
2660 /* TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE SUBP */
2661 /* "FCN" EVALUATED AT THE NEW ITERATE "XPLS" */
2662
2663 /* FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE */
2664 /* 1) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION */
2665 /*    IF NO ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER */
2666 /*    THE GRADIENT OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION */
2667 /*    "FCN" IS INEXPENSIVE TO EVALUATE. */
2668
2669 /* PARAMETERS */
2670 /* ---------- */
2671 /* NR           --> ROW DIMENSION OF MATRIX */
2672 /* N            --> DIMENSION OF PROBLEM */
2673 /* XPLS(N)      --> NEW ITERATE:   X[K] */
2674 /* FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION */
2675 /* FPLS         --> FUNCTION VALUE AT NEW ITERATE, F(XPLS) */
2676 /* A(N,N)      <--  FINITE DIFFERENCE APPROXIMATION TO HESSIAN */
2677 /*                  ONLY LOWER TRIANGULAR MATRIX AND DIAGONAL */
2678 /*                  ARE RETURNED */
2679 /* RNOISE       --> RELATIVE NOISE IN FNAME [F(X)] */
2680 /* STEPSZ(N)    --> WORKSPACE (STEPSIZE IN I-TH COMPONENT DIRECTION) */
2681 /* ANBR(N)      --> WORKSPACE (NEIGHBOR IN I-TH DIRECTION) */
2682 /* NFCNT       <--> NUMBER OF FUNCTION EVALUATIONS */
2683
2684
2685
2686 /* FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION */
2687 /* OF I-TH UNIT VECTOR. */
2688
2689     /* Parameter adjustments */
2690     --anbr;
2691     --stepsz;
2692     --typx;
2693     a_dim1 = *nr;
2694     a_offset = 1 + a_dim1;
2695     a -= a_offset;
2696     --xpls;
2697
2698     /* Function Body */
2699     ov3 = .33333333333333331;
2700     i__1 = *n;
2701     for (i__ = 1; i__ <= i__1; ++i__) {
2702         xtmpi = xpls[i__];
2703 /* Computing MAX */
2704         d__2 = (d__1 = xpls[i__], abs(d__1));
2705         stepsz[i__] = pow_dd(rnoise, &ov3) * max(d__2,1.);
2706         xpls[i__] = xtmpi + stepsz[i__];
2707         (*fcn)(n, &xpls[1], &anbr[i__]);
2708         ++(*nfcnt);
2709         xpls[i__] = xtmpi;
2710 /* L10: */
2711     }
2712
2713 /* CALCULATE COLUMN I OF A */
2714
2715     i__1 = *n;
2716     for (i__ = 1; i__ <= i__1; ++i__) {
2717         xtmpi = xpls[i__];
2718         xpls[i__] = xtmpi + stepsz[i__] * 2.f;
2719         (*fcn)(n, &xpls[1], &fhat);
2720         ++(*nfcnt);
2721         a[i__ + i__ * a_dim1] = (*fpls - anbr[i__] + (fhat - anbr[i__])) / (
2722                 stepsz[i__] * stepsz[i__]);
2723         a[i__ + i__ * a_dim1] *= typx[i__] * typx[i__];
2724
2725 /* CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN */
2726         if (i__ == *n) {
2727             goto L25;
2728         }
2729         xpls[i__] = xtmpi + stepsz[i__];
2730         ip1 = i__ + 1;
2731         i__2 = *n;
2732         for (j = ip1; j <= i__2; ++j) {
2733             xtmpj = xpls[j];
2734             xpls[j] = xtmpj + stepsz[j];
2735             (*fcn)(n, &xpls[1], &fhat);
2736             ++(*nfcnt);
2737             a[j + i__ * a_dim1] = (*fpls - anbr[i__] + (fhat - anbr[j])) / (
2738                     stepsz[i__] * stepsz[j]);
2739             a[j + i__ * a_dim1] *= typx[i__] * typx[j];
2740             xpls[j] = xtmpj;
2741 /* L20: */
2742         }
2743 L25:
2744         xpls[i__] = xtmpi;
2745 /* L30: */
2746     }
2747     return 0;
2748 } /* sndofd_ */
2749
2750 /*  ---------------------- */
2751 /*  |  B A K S L V       | */
2752 /*  ---------------------- */
2753 /* Subroutine */ int bakslv_(integer *nr, integer *n, doublereal *a,
2754         doublereal *x, doublereal *b)
2755 {
2756     /* System generated locals */
2757     integer a_dim1, a_offset, i__1;
2758
2759     /* Local variables */
2760     integer i__, j, ip1;
2761     doublereal sum;
2762
2763
2764 /* PURPOSE */
2765 /* ------- */
2766 /* SOLVE  AX=B  WHERE A IS UPPER TRIANGULAR MATRIX. */
2767 /* NOTE THAT A IS INPUT AS A LOWER TRIANGULAR MATRIX AND */
2768 /* THAT THIS ROUTINE TAKES ITS TRANSPOSE IMPLICITLY. */
2769
2770 /* PARAMETERS */
2771 /* ---------- */
2772 /* NR           --> ROW DIMENSION OF MATRIX */
2773 /* N            --> DIMENSION OF PROBLEM */
2774 /* A(N,N)       --> LOWER TRIANGULAR MATRIX (PRESERVED) */
2775 /* X(N)        <--  SOLUTION VECTOR */
2776 /* B(N)         --> RIGHT-HAND SIDE VECTOR */
2777
2778 /* NOTE */
2779 /* ---- */
2780 /* IF B IS NO LONGER REQUIRED BY CALLING ROUTINE, */
2781 /* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE. */
2782
2783
2784 /* SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE) */
2785
2786     /* Parameter adjustments */
2787     --b;
2788     --x;
2789     a_dim1 = *nr;
2790     a_offset = 1 + a_dim1;
2791     a -= a_offset;
2792
2793     /* Function Body */
2794     i__ = *n;
2795     x[i__] = b[i__] / a[i__ + i__ * a_dim1];
2796     if (*n == 1) {
2797         return 0;
2798     }
2799 L30:
2800     ip1 = i__;
2801     --i__;
2802     sum = 0.f;
2803     i__1 = *n;
2804     for (j = ip1; j <= i__1; ++j) {
2805         sum += a[j + i__ * a_dim1] * x[j];
2806 /* L40: */
2807     }
2808     x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
2809     if (i__ > 1) {
2810         goto L30;
2811     }
2812     return 0;
2813 } /* bakslv_ */
2814
2815 /*  ---------------------- */
2816 /*  |  F O R S L V       | */
2817 /*  ---------------------- */
2818 /* Subroutine */ int forslv_(integer *nr, integer *n, doublereal *a,
2819         doublereal *x, doublereal *b)
2820 {
2821     /* System generated locals */
2822     integer a_dim1, a_offset, i__1, i__2;
2823
2824     /* Local variables */
2825     integer i__, j, im1;
2826     doublereal sum;
2827
2828
2829 /* PURPOSE */
2830 /* -------- */
2831 /* SOLVE  AX=B  WHERE A  IS LOWER TRIANGULAR  MATRIX */
2832
2833 /* PARAMETERS */
2834 /* --------- */
2835
2836 /* NR            -----> ROW DIMENSION OF MATRIX */
2837 /* N             -----> DIMENSION OF PROBLEM */
2838 /* A(N,N)        -----> LOWER TRIANGULAR MATRIX (PRESERVED) */
2839 /* X(N)          <----  SOLUTION VECTOR */
2840 /* B(N)           ----> RIGHT-HAND SIDE VECTOR */
2841
2842 /* NOTE */
2843 /* ----- */
2844 /* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE */
2845
2846
2847 /* SOLVE LX=B.  (FOREWARD  SOLVE) */
2848
2849     /* Parameter adjustments */
2850     --b;
2851     --x;
2852     a_dim1 = *nr;
2853     a_offset = 1 + a_dim1;
2854     a -= a_offset;
2855
2856     /* Function Body */
2857     x[1] = b[1] / a[a_dim1 + 1];
2858     if (*n == 1) {
2859         return 0;
2860     }
2861     i__1 = *n;
2862     for (i__ = 2; i__ <= i__1; ++i__) {
2863         sum = 0.f;
2864         im1 = i__ - 1;
2865         i__2 = im1;
2866         for (j = 1; j <= i__2; ++j) {
2867             sum += a[i__ + j * a_dim1] * x[j];
2868 /* L10: */
2869         }
2870         x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
2871 /* L20: */
2872     }
2873     return 0;
2874 } /* forslv_ */
2875
2876 /*  ---------------------- */
2877 /*  |  C H O L D R       | */
2878 /*  ---------------------- */
2879 /* Subroutine */ int choldr_(integer *nr, integer *n, doublereal *h__,
2880         doublereal *g, doublereal *eps, integer *pivot, doublereal *e,
2881         doublereal *diag, doublereal *addmax)
2882 {
2883     /* System generated locals */
2884     integer h_dim1, h_offset, i__1, i__2, i__3;
2885
2886     /* Builtin functions */
2887     double pow_dd(doublereal *, doublereal *), sqrt(doublereal);
2888
2889     /* Local variables */
2890     integer i__, j, k;
2891     doublereal tau1, tau2;
2892     logical redo;
2893     doublereal temp;
2894     extern /* Subroutine */ int modchl_(integer *, integer *, doublereal *,
2895             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
2896              doublereal *);
2897
2898
2899 /* PURPOSE: */
2900 /*   DRIVER FOR CHOLESKY DECOMPOSITION */
2901
2902 /* ---------------------------------------------------------------------- */
2903
2904 /* PARAMETERS: */
2905
2906 /*   NR            --> ROW DIMENSION */
2907 /*   N       --> DIMENSION OF PROBLEM */
2908 /*   H(N,N)  --> MATRIX */
2909 /*   G(N)    --> WORK SPACE */
2910 /*   EPS           --> MACHINE EPSILON */
2911 /*   PIVOT(N)      --> PIVOTING VECTOR */
2912 /*   E(N)    --> DIAGONAL MATRIX ADDED TO H FOR MAKING H P.D. */
2913 /*   DIAG(N) --> DIAGONAL OF H */
2914 /*   ADDMAX  --> ADDMAX * I  IS ADDED TO H */
2915     /* Parameter adjustments */
2916     --diag;
2917     --e;
2918     --pivot;
2919     --g;
2920     h_dim1 = *nr;
2921     h_offset = 1 + h_dim1;
2922     h__ -= h_offset;
2923
2924     /* Function Body */
2925     redo = FALSE_;
2926
2927 /*     SAVE DIAGONAL OF H */
2928     i__1 = *n;
2929     for (i__ = 1; i__ <= i__1; ++i__) {
2930         diag[i__] = h__[i__ + i__ * h_dim1];
2931 /* L10: */
2932     }
2933     tau1 = pow_dd(eps, &c_b250);
2934     tau2 = tau1;
2935     modchl_(nr, n, &h__[h_offset], &g[1], eps, &tau1, &tau2, &pivot[1], &e[1])
2936             ;
2937     *addmax = e[*n];
2938     i__1 = *n;
2939     for (i__ = 1; i__ <= i__1; ++i__) {
2940         if (pivot[i__] != i__) {
2941             redo = TRUE_;
2942         }
2943 /* L22: */
2944     }
2945     if (*addmax > 0. || redo) {
2946 /* ******************************** */
2947 /*                       * */
2948 /*       H IS NOT P.D.         * */
2949 /*                       * */
2950 /* ******************************** */
2951
2952 /*     H=H+UI */
2953         i__1 = *n;
2954         for (i__ = 2; i__ <= i__1; ++i__) {
2955             i__2 = i__ - 1;
2956             for (j = 1; j <= i__2; ++j) {
2957                 h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
2958 /* L32: */
2959             }
2960 /* L30: */
2961         }
2962         i__1 = *n;
2963         for (i__ = 1; i__ <= i__1; ++i__) {
2964             pivot[i__] = i__;
2965             h__[i__ + i__ * h_dim1] = diag[i__] + *addmax;
2966 /* L34: */
2967         }
2968 /* ******************************** */
2969 /*                       * */
2970 /*        COMPUTE L            * */
2971 /*                       * */
2972 /* ******************************** */
2973         i__1 = *n;
2974         for (j = 1; j <= i__1; ++j) {
2975
2976 /*     COMPUTE L(J,J) */
2977             temp = 0.;
2978             if (j > 1) {
2979                 i__2 = j - 1;
2980                 for (i__ = 1; i__ <= i__2; ++i__) {
2981                     temp += h__[j + i__ * h_dim1] * h__[j + i__ * h_dim1];
2982 /* L41: */
2983                 }
2984             }
2985             h__[j + j * h_dim1] -= temp;
2986             h__[j + j * h_dim1] = sqrt(h__[j + j * h_dim1]);
2987
2988 /*     COMPUTE L(I,J) */
2989             i__2 = *n;
2990             for (i__ = j + 1; i__ <= i__2; ++i__) {
2991                 temp = 0.;
2992                 if (j > 1) {
2993                     i__3 = j - 1;
2994                     for (k = 1; k <= i__3; ++k) {
2995                         temp += h__[i__ + k * h_dim1] * h__[j + k * h_dim1];
2996 /* L45: */
2997                     }
2998                 }
2999                 h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1] - temp;
3000                 h__[i__ + j * h_dim1] /= h__[j + j * h_dim1];
3001 /* L43: */
3002             }
3003 /* L40: */
3004         }
3005     }
3006     return 0;
3007 } /* choldr_ */
3008
3009 /*  ---------------------- */
3010 /*  |  M O D C H L       | */
3011 /*  ---------------------- */
3012 /* ********************************************************************* */
3013
3014 /*       SUBROUTINE NAME: MODCHL */
3015
3016 /*       AUTHORS :  ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
3017
3018 /*       DATE    : DECEMBER, 1988 */
3019
3020 /*       PURPOSE : PERFORM A MODIFIED CHOLESKY FACTORIZATION */
3021 /*                 OF THE FORM (PTRANSPOSE)AP  + E = L(LTRANSPOSE), */
3022 /*       WHERE L IS STORED IN THE LOWER TRIANGLE OF THE */
3023 /*       ORIGINAL MATRIX A. */
3024 /*       THE FACTORIZATION HAS 2 PHASES: */
3025 /*        PHASE 1: PIVOT ON THE MAXIMUM DIAGONAL ELEMENT. */
3026 /*            CHECK THAT THE NORMAL CHOLESKY UPDATE */
3027 /*            WOULD RESULT IN A POSITIVE DIAGONAL */
3028 /*            AT THE CURRENT ITERATION, AND */
3029 /*            IF SO, DO THE NORMAL CHOLESKY UPDATE, */
3030 /*            OTHERWISE SWITCH TO PHASE 2. */
3031 /*        PHASE 2: PIVOT ON THE MINIMUM OF THE NEGATIVES */
3032 /*            OF THE LOWER GERSCHGORIN BOUND */
3033 /*            ESTIMATES. */
3034 /*            COMPUTE THE AMOUNT TO ADD TO THE */
3035 /*            PIVOT ELEMENT AND ADD THIS */
3036 /*            TO THE PIVOT ELEMENT. */
3037 /*            DO THE CHOLESKY UPDATE. */
3038 /*            UPDATE THE ESTIMATES OF THE */
3039 /*            GERSCHGORIN BOUNDS. */
3040
3041 /*       INPUT   : NDIM    - LARGEST DIMENSION OF MATRIX THAT */
3042 /*                           WILL BE USED */
3043
3044 /*                 N       - DIMENSION OF MATRIX A */
3045
3046 /*                 A       - N*N SYMMETRIC MATRIX (ONLY LOWER TRIANGULAR */
3047 /*            PORTION OF A, INCLUDING THE MAIN DIAGONAL, IS USED) */
3048
3049 /*                 G       - N*1 WORK ARRAY */
3050
3051 /*                 MCHEPS - MACHINE PRECISION */
3052
3053 /*                TAU1    - TOLERANCE USED FOR DETERMINING WHEN TO SWITCH TO */
3054 /*                          PHASE 2 */
3055
3056 /*                TAU2    - TOLERANCE USED FOR DETERMINING THE MAXIMUM */
3057 /*                          CONDITION NUMBER OF THE FINAL 2X2 SUBMATRIX. */
3058
3059
3060 /*       OUTPUT  : L     - STORED IN THE MATRIX A (IN LOWER TRIANGULAR */
3061 /*                           PORTION OF A, INCLUDING THE MAIN DIAGONAL) */
3062
3063 /*                 P     - A RECORD OF HOW THE ROWS AND COLUMNS */
3064 /*                         OF THE MATRIX WERE PERMUTED WHILE */
3065 /*                         PERFORMING THE DECOMPOSITION */
3066
3067 /*                 E     - N*1 ARRAY, THE ITH ELEMENT IS THE */
3068 /*                         AMOUNT ADDED TO THE DIAGONAL OF A */
3069 /*                         AT THE ITH ITERATION */
3070
3071
3072 /* ************************************************************************ */
3073 /* Subroutine */ int modchl_(integer *ndim, integer *n, doublereal *a,
3074         doublereal *g, doublereal *mcheps, doublereal *tau1, doublereal *tau2,
3075          integer *p, doublereal *e)
3076 {
3077     /* System generated locals */
3078     integer a_dim1, a_offset, i__1, i__2, i__3;
3079     doublereal d__1;
3080
3081     /* Builtin functions */
3082     double sqrt(doublereal);
3083
3084     /* Local variables */
3085     integer i__, j, k, jp1;
3086     doublereal maxd, ming;
3087     extern /* Subroutine */ int init_(integer *, integer *, doublereal *,
3088             logical *, doublereal *, integer *, doublereal *, doublereal *,
3089             doublereal *, doublereal *, doublereal *, doublereal *);
3090     doublereal temp, gamma, delta, jdmin;
3091     integer imaxd, iming;
3092     extern /* Subroutine */ int fin2x2_(integer *, integer *, doublereal *,
3093             doublereal *, integer *, doublereal *, doublereal *, doublereal *)
3094             ;
3095     doublereal tdmin;
3096     integer itemp;
3097     doublereal normj, delta1;
3098     logical phase1;
3099     extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
3100             integer *, doublereal *);
3101     doublereal taugam, tempjj;
3102
3103
3104 /*  J        - CURRENT ITERATION NUMBER */
3105 /*  IMING    - INDEX OF THE ROW WITH THE MIN. OF THE */
3106 /*           NEG. LOWER GERSCH. BOUNDS */
3107 /*  IMAXD    - INDEX OF THE ROW WITH THE MAXIMUM DIAG. */
3108 /*           ELEMENT */
3109 /*  I,ITEMP,JPL,K  - TEMPORARY INTEGER VARIABLES */
3110 /*  DELTA    - AMOUNT TO ADD TO AJJ AT THE JTH ITERATION */
3111 /*  GAMMA    - THE MAXIMUM DIAGONAL ELEMENT OF THE ORIGINAL */
3112 /*           MATRIX A. */
3113 /*  NORMJ    - THE 1 NORM OF A(COLJ), ROWS J+1 --> N. */
3114 /*  MING     - THE MINIMUM OF THE NEG. LOWER GERSCH. BOUNDS */
3115 /*  MAXD     - THE MAXIMUM DIAGONAL ELEMENT */
3116 /*  TAUGAM - TAU1 * GAMMA */
3117 /*  PHASE1      - LOGICAL, TRUE IF IN PHASE1, OTHERWISE FALSE */
3118 /*  DELTA1,TEMP,JDMIN,TDMIN,TEMPJJ - TEMPORARY DOUBLE PRECISION VARS. */
3119
3120     /* Parameter adjustments */
3121     --e;
3122     --p;
3123     --g;
3124     a_dim1 = *ndim;
3125     a_offset = 1 + a_dim1;
3126     a -= a_offset;
3127
3128     /* Function Body */
3129     init_(n, ndim, &a[a_offset], &phase1, &delta, &p[1], &g[1], &e[1], &ming,
3130             tau1, &gamma, &taugam);
3131 /*     CHECK FOR N=1 */
3132     if (*n == 1) {
3133         delta = *tau2 * (d__1 = a[a_dim1 + 1], abs(d__1)) - a[a_dim1 + 1];
3134         if (delta > 0.) {
3135             e[1] = delta;
3136         }
3137         if (a[a_dim1 + 1] == 0.) {
3138             e[1] = *tau2;
3139         }
3140         a[a_dim1 + 1] = sqrt(a[a_dim1 + 1] + e[1]);
3141     }
3142
3143
3144     i__1 = *n - 1;
3145     for (j = 1; j <= i__1; ++j) {
3146
3147 /*        PHASE 1 */
3148
3149         if (phase1) {
3150
3151 /*           FIND INDEX OF MAXIMUM DIAGONAL ELEMENT A(I,I) WHERE I>=J */
3152
3153             maxd = a[j + j * a_dim1];
3154             imaxd = j;
3155             i__2 = *n;
3156             for (i__ = j + 1; i__ <= i__2; ++i__) {
3157                 if (maxd < a[i__ + i__ * a_dim1]) {
3158                     maxd = a[i__ + i__ * a_dim1];
3159                     imaxd = i__;
3160                 }
3161 /* L20: */
3162             }
3163
3164 /*           PIVOT TO THE TOP THE ROW AND COLUMN WITH THE MAX DIAG */
3165
3166             if (imaxd != j) {
3167
3168 /*              SWAP ROW J WITH ROW OF MAX DIAG */
3169
3170                 i__2 = j - 1;
3171                 for (i__ = 1; i__ <= i__2; ++i__) {
3172                     temp = a[j + i__ * a_dim1];
3173                     a[j + i__ * a_dim1] = a[imaxd + i__ * a_dim1];
3174                     a[imaxd + i__ * a_dim1] = temp;
3175 /* L30: */
3176                 }
3177
3178 /*              SWAP COLJ AND ROW MAXDIAG BETWEEN J AND MAXDIAG */
3179
3180                 i__2 = imaxd - 1;
3181                 for (i__ = j + 1; i__ <= i__2; ++i__) {
3182                     temp = a[i__ + j * a_dim1];
3183                     a[i__ + j * a_dim1] = a[imaxd + i__ * a_dim1];
3184                     a[imaxd + i__ * a_dim1] = temp;
3185 /* L35: */
3186                 }
3187
3188 /*              SWAP COLUMN J WITH COLUMN OF MAX DIAG */
3189
3190                 i__2 = *n;
3191                 for (i__ = imaxd + 1; i__ <= i__2; ++i__) {
3192                     temp = a[i__ + j * a_dim1];
3193                     a[i__ + j * a_dim1] = a[i__ + imaxd * a_dim1];
3194                     a[i__ + imaxd * a_dim1] = temp;
3195 /* L40: */
3196                 }
3197
3198 /*              SWAP DIAG ELEMENTS */
3199
3200                 temp = a[j + j * a_dim1];
3201                 a[j + j * a_dim1] = a[imaxd + imaxd * a_dim1];
3202                 a[imaxd + imaxd * a_dim1] = temp;
3203
3204 /*              SWAP ELEMENTS OF THE PERMUTATION VECTOR */
3205
3206                 itemp = p[j];
3207                 p[j] = p[imaxd];
3208                 p[imaxd] = itemp;
3209             }
3210 /*           CHECK TO SEE WHETHER THE NORMAL CHOLESKY UPDATE FOR THIS */
3211 /*           ITERATION WOULD RESULT IN A POSITIVE DIAGONAL, */
3212 /*           AND IF NOT THEN SWITCH TO PHASE 2. */
3213             jp1 = j + 1;
3214             tempjj = a[j + j * a_dim1];
3215             if (tempjj > 0.) {
3216                 jdmin = a[jp1 + jp1 * a_dim1];
3217                 i__2 = *n;
3218                 for (i__ = jp1; i__ <= i__2; ++i__) {
3219                     temp = a[i__ + j * a_dim1] * a[i__ + j * a_dim1] / tempjj;
3220                     tdmin = a[i__ + i__ * a_dim1] - temp;
3221                     jdmin = min(jdmin,tdmin);
3222 /* L60: */
3223                 }
3224                 if (jdmin < taugam) {
3225                     phase1 = FALSE_;
3226                 }
3227             } else {
3228                 phase1 = FALSE_;
3229             }
3230             if (phase1) {
3231
3232 /*              DO THE NORMAL CHOLESKY UPDATE IF STILL IN PHASE 1 */
3233
3234                 a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
3235                 tempjj = a[j + j * a_dim1];
3236                 i__2 = *n;
3237                 for (i__ = jp1; i__ <= i__2; ++i__) {
3238                     a[i__ + j * a_dim1] /= tempjj;
3239 /* L70: */
3240                 }
3241                 i__2 = *n;
3242                 for (i__ = jp1; i__ <= i__2; ++i__) {
3243                     temp = a[i__ + j * a_dim1];
3244                     i__3 = i__;
3245                     for (k = jp1; k <= i__3; ++k) {
3246                         a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
3247 /* L75: */
3248                     }
3249 /* L80: */
3250                 }
3251                 if (j == *n - 1) {
3252                     a[*n + *n * a_dim1] = sqrt(a[*n + *n * a_dim1]);
3253                 }
3254             } else {
3255
3256 /*              CALCULATE THE NEGATIVES OF THE LOWER GERSCHGORIN BOUNDS */
3257
3258                 gersch_(ndim, n, &a[a_offset], &j, &g[1]);
3259             }
3260         }
3261
3262 /*        PHASE 2 */
3263
3264         if (! phase1) {
3265             if (j != *n - 1) {
3266
3267 /*              FIND THE MINIMUM NEGATIVE GERSHGORIN BOUND */
3268
3269                 iming = j;
3270                 ming = g[j];
3271                 i__2 = *n;
3272                 for (i__ = j + 1; i__ <= i__2; ++i__) {
3273                     if (ming > g[i__]) {
3274                         ming = g[i__];
3275                         iming = i__;
3276                     }
3277 /* L90: */
3278                 }
3279
3280 /*               PIVOT TO THE TOP THE ROW AND COLUMN WITH THE */
3281 /*               MINIMUM NEGATIVE GERSCHGORIN BOUND */
3282
3283                 if (iming != j) {
3284
3285 /*                  SWAP ROW J WITH ROW OF MIN GERSCH BOUND */
3286
3287                     i__2 = j - 1;
3288                     for (i__ = 1; i__ <= i__2; ++i__) {
3289                         temp = a[j + i__ * a_dim1];
3290                         a[j + i__ * a_dim1] = a[iming + i__ * a_dim1];
3291                         a[iming + i__ * a_dim1] = temp;
3292 /* L100: */
3293                     }
3294
3295 /*                  SWAP COLJ WITH ROW IMING FROM J TO IMING */
3296
3297                     i__2 = iming - 1;
3298                     for (i__ = j + 1; i__ <= i__2; ++i__) {
3299                         temp = a[i__ + j * a_dim1];
3300                         a[i__ + j * a_dim1] = a[iming + i__ * a_dim1];
3301                         a[iming + i__ * a_dim1] = temp;
3302 /* L105: */
3303                     }
3304
3305 /*                 SWAP COLUMN J WITH COLUMN OF MIN GERSCH BOUND */
3306
3307                     i__2 = *n;
3308                     for (i__ = iming + 1; i__ <= i__2; ++i__) {
3309                         temp = a[i__ + j * a_dim1];
3310                         a[i__ + j * a_dim1] = a[i__ + iming * a_dim1];
3311                         a[i__ + iming * a_dim1] = temp;
3312 /* L110: */
3313                     }
3314
3315 /*                 SWAP DIAGONAL ELEMENTS */
3316
3317                     temp = a[j + j * a_dim1];
3318                     a[j + j * a_dim1] = a[iming + iming * a_dim1];
3319                     a[iming + iming * a_dim1] = temp;
3320
3321 /*                 SWAP ELEMENTS OF THE PERMUTATION VECTOR */
3322
3323                     itemp = p[j];
3324                     p[j] = p[iming];
3325                     p[iming] = itemp;
3326
3327 /*                 SWAP ELEMENTS OF THE NEGATIVE GERSCHGORIN BOUNDS VECTOR */
3328
3329                     temp = g[j];
3330                     g[j] = g[iming];
3331                     g[iming] = temp;
3332                 }
3333
3334 /*              CALCULATE DELTA AND ADD TO THE DIAGONAL. */
3335 /*              DELTA=MAX{0,-A(J,J) + MAX{NORMJ,TAUGAM},DELTA_PREVIOUS} */
3336 /*              WHERE NORMJ=SUM OF |A(I,J)|,FOR I=1,N, */
3337 /*              DELTA_PREVIOUS IS THE DELTA COMPUTED AT THE PREVIOUS ITERATION, */
3338 /*              AND TAUGAM IS TAU1*GAMMA. */
3339
3340                 normj = 0.f;
3341                 i__2 = *n;
3342                 for (i__ = j + 1; i__ <= i__2; ++i__) {
3343                     normj += (d__1 = a[i__ + j * a_dim1], abs(d__1));
3344 /* L140: */
3345                 }
3346                 temp = max(normj,taugam);
3347                 delta1 = temp - a[j + j * a_dim1];
3348                 temp = 0.f;
3349                 delta1 = max(temp,delta1);
3350                 delta = max(delta1,delta);
3351                 e[j] = delta;
3352                 a[j + j * a_dim1] += e[j];
3353
3354 /*              UPDATE THE GERSCHGORIN BOUND ESTIMATES */
3355 /*              (NOTE: G(I) IS THE NEGATIVE OF THE */
3356 /*               GERSCHGORIN LOWER BOUND.) */
3357
3358                 if (a[j + j * a_dim1] != normj) {
3359                     temp = normj / a[j + j * a_dim1] - 1.f;
3360                     i__2 = *n;
3361                     for (i__ = j + 1; i__ <= i__2; ++i__) {
3362                         g[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)) *
3363                                 temp;
3364 /* L150: */
3365                     }
3366                 }
3367
3368 /*              DO THE CHOLESKY UPDATE */
3369
3370                 a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
3371                 tempjj = a[j + j * a_dim1];
3372                 i__2 = *n;
3373                 for (i__ = j + 1; i__ <= i__2; ++i__) {
3374                     a[i__ + j * a_dim1] /= tempjj;
3375 /* L160: */
3376                 }
3377                 i__2 = *n;
3378                 for (i__ = j + 1; i__ <= i__2; ++i__) {
3379                     temp = a[i__ + j * a_dim1];
3380                     i__3 = i__;
3381                     for (k = j + 1; k <= i__3; ++k) {
3382                         a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
3383 /* L170: */
3384                     }
3385 /* L180: */
3386                 }
3387             } else {
3388                 fin2x2_(ndim, n, &a[a_offset], &e[1], &j, tau2, &delta, &
3389                         gamma);
3390             }
3391         }
3392 /* L200: */
3393     }
3394     return 0;
3395 } /* modchl_ */
3396
3397 /* ************************************************************************ */
3398 /*       SUBROUTINE NAME : INIT */
3399
3400 /*       PURPOSE : SET UP FOR START OF CHOLESKY FACTORIZATION */
3401
3402 /*       INPUT : N, NDIM, A, TAU1 */
3403
3404 /*       OUTPUT : PHASE1    - BOOLEAN VALUE SET TO TRUE IF IN PHASE ONE, */
3405 /*             OTHERWISE FALSE. */
3406 /*      DELTA     - AMOUNT TO ADD TO AJJ AT ITERATION J */
3407 /*      P,G,E - DESCRIBED ABOVE IN MODCHL */
3408 /*      MING      - THE MINIMUM NEGATIVE GERSCHGORIN BOUND */
3409 /*      GAMMA     - THE MAXIMUM DIAGONAL ELEMENT OF A */
3410 /*      TAUGAM  - TAU1 * GAMMA */
3411
3412 /* ************************************************************************ */
3413 /* Subroutine */ int init_(integer *n, integer *ndim, doublereal *a, logical *
3414         phase1, doublereal *delta, integer *p, doublereal *g, doublereal *e,
3415         doublereal *ming, doublereal *tau1, doublereal *gamma, doublereal *
3416         taugam)
3417 {
3418     /* System generated locals */
3419     integer a_dim1, a_offset, i__1;
3420     doublereal d__1, d__2, d__3;
3421
3422     /* Local variables */
3423     integer i__;
3424     extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
3425             integer *, doublereal *);
3426
3427     /* Parameter adjustments */
3428     --e;
3429     --g;
3430     --p;
3431     a_dim1 = *ndim;
3432     a_offset = 1 + a_dim1;
3433     a -= a_offset;
3434
3435     /* Function Body */
3436     *phase1 = TRUE_;
3437     *delta = 0.f;
3438     *ming = 0.f;
3439     i__1 = *n;
3440     for (i__ = 1; i__ <= i__1; ++i__) {
3441         p[i__] = i__;
3442         g[i__] = 0.f;
3443         e[i__] = 0.f;
3444 /* L10: */
3445     }
3446
3447 /*     FIND THE MAXIMUM MAGNITUDE OF THE DIAGONAL ELEMENTS. */
3448 /*     IF ANY DIAGONAL ELEMENT IS NEGATIVE, THEN PHASE1 IS FALSE. */
3449
3450     *gamma = 0.f;
3451     i__1 = *n;
3452     for (i__ = 1; i__ <= i__1; ++i__) {
3453 /* Computing MAX */
3454         d__2 = *gamma, d__3 = (d__1 = a[i__ + i__ * a_dim1], abs(d__1));
3455         *gamma = max(d__2,d__3);
3456         if (a[i__ + i__ * a_dim1] < 0.f) {
3457             *phase1 = FALSE_;
3458         }
3459 /* L20: */
3460     }
3461     *taugam = *tau1 * *gamma;
3462
3463 /*     IF NOT IN PHASE1, THEN CALCULATE THE INITIAL GERSCHGORIN BOUNDS */
3464 /*     NEEDED FOR THE START OF PHASE2. */
3465
3466     if (! (*phase1)) {
3467         gersch_(ndim, n, &a[a_offset], &c__1, &g[1]);
3468     }
3469     return 0;
3470 } /* init_ */
3471
3472 /* ************************************************************************ */
3473
3474 /*       SUBROUTINE NAME : GERSCH */
3475
3476 /*       PURPOSE : CALCULATE THE NEGATIVE OF THE GERSCHGORIN BOUNDS */
3477 /*                 CALLED ONCE AT THE START OF PHASE II. */
3478
3479 /*       INPUT   : NDIM, N, A, J */
3480
3481 /*       OUTPUT  : G - AN N VECTOR CONTAINING THE NEGATIVES OF THE */
3482 /*           GERSCHGORIN BOUNDS. */
3483
3484 /* ************************************************************************ */
3485 /* Subroutine */ int gersch_(integer *ndim, integer *n, doublereal *a,
3486         integer *j, doublereal *g)
3487 {
3488     /* System generated locals */
3489     integer a_dim1, a_offset, i__1, i__2;
3490     doublereal d__1;
3491
3492     /* Local variables */
3493     integer i__, k;
3494     doublereal offrow;
3495
3496     /* Parameter adjustments */
3497     --g;
3498     a_dim1 = *ndim;
3499     a_offset = 1 + a_dim1;
3500     a -= a_offset;
3501
3502     /* Function Body */
3503     i__1 = *n;
3504     for (i__ = *j; i__ <= i__1; ++i__) {
3505         offrow = 0.f;
3506         i__2 = i__ - 1;
3507         for (k = *j; k <= i__2; ++k) {
3508             offrow += (d__1 = a[i__ + k * a_dim1], abs(d__1));
3509 /* L10: */
3510         }
3511         i__2 = *n;
3512         for (k = i__ + 1; k <= i__2; ++k) {
3513             offrow += (d__1 = a[k + i__ * a_dim1], abs(d__1));
3514 /* L20: */
3515         }
3516         g[i__] = offrow - a[i__ + i__ * a_dim1];
3517 /* L30: */
3518     }
3519     return 0;
3520 } /* gersch_ */
3521
3522 /* ************************************************************************ */
3523
3524 /*  SUBROUTINE NAME : FIN2X2 */
3525
3526 /*  PURPOSE : HANDLES FINAL 2X2 SUBMATRIX IN PHASE II. */
3527 /*            FINDS EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX, */
3528 /*            CALCULATES THE AMOUNT TO ADD TO THE DIAGONAL, */
3529 /*            ADDS TO THE FINAL 2 DIAGONAL ELEMENTS, */
3530 /*            AND DOES THE FINAL UPDATE. */
3531
3532 /*  INPUT : NDIM, N, A, E, J, TAU2, */
3533 /*          DELTA - AMOUNT ADDED TO THE DIAGONAL IN THE */
3534 /*                  PREVIOUS ITERATION */
3535
3536 /*  OUTPUT : A - MATRIX WITH COMPLETE L FACTOR IN THE LOWER TRIANGLE, */
3537 /*           E - N*1 VECTOR CONTAINING THE AMOUNT ADDED TO THE DIAGONAL */
3538 /*               AT EACH ITERATION, */
3539 /*           DELTA - AMOUNT ADDED TO DIAGONAL ELEMENTS N-1 AND N. */
3540
3541 /* ************************************************************************ */
3542 /* Subroutine */ int fin2x2_(integer *ndim, integer *n, doublereal *a,
3543         doublereal *e, integer *j, doublereal *tau2, doublereal *delta,
3544         doublereal *gamma)
3545 {
3546     /* System generated locals */
3547     integer a_dim1, a_offset;
3548     doublereal d__1;
3549
3550     /* Builtin functions */
3551     double sqrt(doublereal);
3552
3553     /* Local variables */
3554     doublereal t1, t2, t3, t1a, t2a, temp, lmbd1, lmbd2, delta1, lmbdhi,
3555             lmbdlo;
3556
3557
3558 /*     FIND EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX */
3559
3560     /* Parameter adjustments */
3561     --e;
3562     a_dim1 = *ndim;
3563     a_offset = 1 + a_dim1;
3564     a -= a_offset;
3565
3566     /* Function Body */
3567     t1 = a[*n - 1 + (*n - 1) * a_dim1] + a[*n + *n * a_dim1];
3568     t2 = a[*n - 1 + (*n - 1) * a_dim1] - a[*n + *n * a_dim1];
3569     t1a = abs(t2);
3570     t2a = (d__1 = a[*n + (*n - 1) * a_dim1], abs(d__1)) * 2.;
3571     if (t1a >= t2a) {
3572         if (t1a > 0.) {
3573             t2a /= t1a;
3574         }
3575 /* Computing 2nd power */
3576         d__1 = t2a;
3577         t3 = t1a * sqrt(d__1 * d__1 + 1.);
3578     } else {
3579         t1a /= t2a;
3580 /* Computing 2nd power */
3581         d__1 = t1a;
3582         t3 = t2a * sqrt(d__1 * d__1 + 1.);
3583     }
3584     lmbd1 = (t1 - t3) / 2.f;
3585     lmbd2 = (t1 + t3) / 2.f;
3586     lmbdhi = max(lmbd1,lmbd2);
3587     lmbdlo = min(lmbd1,lmbd2);
3588
3589 /*     FIND DELTA SUCH THAT: */
3590 /*     1.  THE L2 CONDITION NUMBER OF THE FINAL */
3591 /*     2X2 SUBMATRIX + DELTA*I <= TAU2 */
3592 /*     2. DELTA >= PREVIOUS DELTA, */
3593 /*     3. LMBDLO + DELTA >= TAU2 * GAMMA, */
3594 /*     WHERE LMBDLO IS THE SMALLEST EIGENVALUE OF THE FINAL */
3595 /*     2X2 SUBMATRIX */
3596
3597     delta1 = (lmbdhi - lmbdlo) / (1.f - *tau2);
3598     delta1 = max(delta1,*gamma);
3599     delta1 = *tau2 * delta1 - lmbdlo;
3600     temp = 0.f;
3601     *delta = max(*delta,temp);
3602     *delta = max(*delta,delta1);
3603     if (*delta > 0.f) {
3604         a[*n - 1 + (*n - 1) * a_dim1] += *delta;
3605         a[*n + *n * a_dim1] += *delta;
3606         e[*n - 1] = *delta;
3607         e[*n] = *delta;
3608     }
3609
3610 /*     FINAL UPDATE */
3611
3612     a[*n - 1 + (*n - 1) * a_dim1] = sqrt(a[*n - 1 + (*n - 1) * a_dim1]);
3613     a[*n + (*n - 1) * a_dim1] /= a[*n - 1 + (*n - 1) * a_dim1];
3614 /* Computing 2nd power */
3615     d__1 = a[*n + (*n - 1) * a_dim1];
3616     a[*n + *n * a_dim1] -= d__1 * d__1;
3617     a[*n + *n * a_dim1] = sqrt(a[*n + *n * a_dim1]);
3618     return 0;
3619 } /* fin2x2_ */
3620
3621 /*  ---------------------- */
3622 /*  |    S L V M D L     | */
3623 /*  ---------------------- */
3624 /* Subroutine */ int slvmdl_(integer *nr, integer *n, doublereal *h__,
3625         doublereal *u, doublereal *t, doublereal *e, doublereal *diag,
3626         doublereal *s, doublereal *g, integer *pivot, doublereal *w1,
3627         doublereal *w2, doublereal *w3, doublereal *alpha, doublereal *beta,
3628         logical *nomin, doublereal *eps)
3629 {
3630     /* System generated locals */
3631     integer h_dim1, h_offset, i__1, i__2;
3632
3633     /* Builtin functions */
3634     double sqrt(doublereal);
3635
3636     /* Local variables */
3637     integer i__, j;
3638     doublereal r__, r1, r2, ca, cb, cc, cd, w11, sg, w22, w12, w33, w13, w23,
3639             ss, uu, shs;
3640     extern /* Subroutine */ int zhz_(integer *, integer *, doublereal *,
3641             doublereal *, doublereal *, doublereal *);
3642     doublereal temp;
3643     extern /* Subroutine */ int sigma_(doublereal *, doublereal *, doublereal
3644             *, doublereal *, doublereal *), dstar_(integer *, integer *,
3645             doublereal *, doublereal *, doublereal *, doublereal *,
3646             doublereal *, doublereal *, doublereal *, doublereal *);
3647     doublereal addmax;
3648     extern /* Subroutine */ int choldr_(integer *, integer *, doublereal *,
3649             doublereal *, doublereal *, integer *, doublereal *, doublereal *,
3650              doublereal *), bakslv_(integer *, integer *, doublereal *,
3651             doublereal *, doublereal *);
3652     doublereal sgstar;
3653     extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
3654             doublereal *, doublereal *), solvew_(integer *, integer *,
3655             doublereal *, doublereal *, doublereal *, doublereal *);
3656
3657
3658 /* PURPOSE: */
3659 /*   COMPUTE TENSOR AND NEWTON STEPS */
3660
3661 /* ---------------------------------------------------------------------------- */
3662
3663 /* PARAMETERS: */
3664
3665 /*   NR            --> ROW DIMENSION OF MATRIX */
3666 /*   N       --> DIMENSION OF PROBLEM */
3667 /*   H(N,N)  --> HESSIAN */
3668 /*   U(N)    --> VECTOR TO FORM Q IN QR */
3669 /*   T(N)    --> WORKSPACE */
3670 /*   E(N)    --> DIAGONAL ADDED TO HESSIAN IN CHOLESKY DECOMPOSITION */
3671 /*   DIAG(N) --> DIAGONAL OF HESSIAN */
3672 /*   S(N)    --> STEP TO PREVIOUS POINT (FOR TENSOR MODEL) */
3673 /*   G(N)    --> CURRENT GRADIENT */
3674 /*   PIVOT(N)      --> PIVOT VECTOR FOR CHOLESKY DECOMPOSITION */
3675 /*   W1(N)      <--> ON INPUT: A=2*(GP-G-HS-S*BETA/(6*STS)) */
3676 /*               ON OUTPUT: TENSOR STEP */
3677 /*   W2(N)   --> SH */
3678 /*   W3(N)      <--  NEWTON STEP */
3679 /*   ALPHA   --> SCALAR FOR 3RD ORDER TERM OF TENSOR MODEL */
3680 /*   BETA    --> SCALAR FOR 4TH ORDER TERM OF TENSOR MODEL */
3681 /*   NOMIN      <--  =.TRUE. IF TENSOR MODEL HAS NO MINIMIZER */
3682 /*   EPS           --> MACHINE EPSILON */
3683
3684
3685 /*     S O L V E    M O D E L */
3686
3687     /* Parameter adjustments */
3688     --w3;
3689     --w2;
3690     --w1;
3691     --pivot;
3692     --g;
3693     --s;
3694     --diag;
3695     --e;
3696     --t;
3697     --u;
3698     h_dim1 = *nr;
3699     h_offset = 1 + h_dim1;
3700     h__ -= h_offset;
3701
3702     /* Function Body */
3703     *nomin = FALSE_;
3704
3705 /*     COMPUTE QTHQ(N,N), ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND N-1 */
3706 /*     COLUMNS OF QTHQ */
3707     if (*n > 1) {
3708         zhz_(nr, n, &s[1], &h__[h_offset], &u[1], &t[1]);
3709
3710 /*     IN CHOLESKY DECOMPOSITION WILL STORE H(1,1) ... H(N-1,N-1) */
3711 /*     IN DIAG(1) ... DIAG(N-1), STORE H(N,N) IN DIAG(N) FIRST */
3712         diag[*n] = h__[*n + *n * h_dim1];
3713
3714 /*     COLESKY DECOMPOSITION FOR FIRST N-1 ROWS AND N-1 COLUMNS OF ZTHZ */
3715 /*     ZTHZ(N-1,N-1)=LLT */
3716         i__1 = *n - 1;
3717         choldr_(nr, &i__1, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
3718                 diag[1], &addmax);
3719     }
3720
3721 /*     ON INPUT: SH IS STORED IN W2 */
3722     shs = 0.;
3723     i__1 = *n;
3724     for (i__ = 1; i__ <= i__1; ++i__) {
3725         shs += w2[i__] * s[i__];
3726         w3[i__] = g[i__];
3727 /* L100: */
3728     }
3729
3730 /*   COMPUTE W1,W2,W3 */
3731 /*   W1=L**-1*ZT*A */
3732 /*   W2=L**-1*ZT*SH */
3733 /*   W3=L**-1*ZT*G */
3734
3735     if (*n > 1) {
3736         solvew_(nr, n, &h__[h_offset], &u[1], &w1[1], &t[1]);
3737         solvew_(nr, n, &h__[h_offset], &u[1], &w2[1], &t[1]);
3738         solvew_(nr, n, &h__[h_offset], &u[1], &w3[1], &t[1]);
3739     }
3740
3741 /*     COMPUTE COEFFICIENTS CA,CB,CC AND CD FOR REDUCED ONE VARIABLE */
3742 /*     3RD ORDER EQUATION */
3743     w11 = 0.;
3744     i__1 = *n - 1;
3745     for (i__ = 1; i__ <= i__1; ++i__) {
3746         w11 += w1[i__] * w1[i__];
3747 /* L110: */
3748     }
3749     ca = *beta / 6. - w11 / 2.;
3750     w12 = 0.;
3751     i__1 = *n - 1;
3752     for (i__ = 1; i__ <= i__1; ++i__) {
3753         w12 += w1[i__] * w2[i__];
3754 /* L120: */
3755     }
3756     cb = *alpha / 2. - w12 * 3. / 2.;
3757     w13 = 0.;
3758     i__1 = *n - 1;
3759     for (i__ = 1; i__ <= i__1; ++i__) {
3760         w13 += w1[i__] * w3[i__];
3761 /* L130: */
3762     }
3763     w22 = 0.;
3764     i__1 = *n - 1;
3765     for (i__ = 1; i__ <= i__1; ++i__) {
3766         w22 += w2[i__] * w2[i__];
3767 /* L133: */
3768     }
3769     cc = shs - w22 - w13;
3770     sg = 0.;
3771     i__1 = *n;
3772     for (i__ = 1; i__ <= i__1; ++i__) {
3773         sg += s[i__] * g[i__];
3774 /* L140: */
3775     }
3776     w23 = 0.;
3777     i__1 = *n - 1;
3778     for (i__ = 1; i__ <= i__1; ++i__) {
3779         w23 += w2[i__] * w3[i__];
3780 /* L145: */
3781     }
3782     cd = sg - w23;
3783     w33 = 0.;
3784     i__1 = *n - 1;
3785     for (i__ = 1; i__ <= i__1; ++i__) {
3786         w33 += w3[i__] * w3[i__];
3787 /* L147: */
3788     }
3789
3790 /*     COMPUTE DESIRABLE ROOT, SGSTAR, OF 3RD ORDER EQUATION */
3791     if (ca != 0.) {
3792         sigma_(&sgstar, &ca, &cb, &cc, &cd);
3793         if (ca == 0.) {
3794             *nomin = TRUE_;
3795             goto L200;
3796         }
3797     } else {
3798 /*       2ND ORDER ( CA=0 ) */
3799         if (cb != 0.) {
3800             r__ = cc * cc - cb * 4. * cd;
3801             if (r__ < 0.) {
3802                 *nomin = TRUE_;
3803                 goto L200;
3804             } else {
3805                 r1 = (-cc + sqrt(r__)) / (cb * 2.);
3806                 r2 = (-cc - sqrt(r__)) / (cb * 2.);
3807                 if (r2 < r1) {
3808                     temp = r1;
3809                     r1 = r2;
3810                     r2 = temp;
3811                 }
3812                 if (cb > 0.) {
3813                     sgstar = r2;
3814                 } else {
3815                     sgstar = r1;
3816                 }
3817                 if (r1 > 0. && sgstar == r2 || r2 < 0. && sgstar == r1) {
3818                     *nomin = TRUE_;
3819                     goto L200;
3820                 }
3821             }
3822         } else {
3823 /*         1ST ORDER (CA=0,CB=0) */
3824             if (cc > 0.) {
3825                 sgstar = -cd / cc;
3826             } else {
3827                 *nomin = TRUE_;
3828                 goto L200;
3829             }
3830         }
3831     }
3832
3833 /*     FIND TENSOR STEP, W1 (FUNCTION OF SGSTAR) */
3834     dstar_(nr, n, &u[1], &s[1], &w1[1], &w2[1], &w3[1], &sgstar, &h__[
3835             h_offset], &w1[1]);
3836
3837 /*     COMPUTE DN */
3838 L200:
3839     uu = 0.;
3840     ss = 0.;
3841     i__1 = *n;
3842     for (i__ = 1; i__ <= i__1; ++i__) {
3843         uu += u[i__] * u[i__];
3844         ss += s[i__] * s[i__];
3845 /* L202: */
3846     }
3847     uu /= 2.;
3848     ss = sqrt(ss);
3849     if (*n == 1) {
3850         choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &diag[1],
3851                  &addmax);
3852     } else {
3853
3854 /* COMPUTE LAST ROW OF L(N,N) */
3855         i__1 = *n - 1;
3856         for (i__ = 1; i__ <= i__1; ++i__) {
3857             temp = 0.;
3858             if (i__ > 1) {
3859                 i__2 = i__ - 1;
3860                 for (j = 1; j <= i__2; ++j) {
3861                     temp += h__[*n + j * h_dim1] * h__[i__ + j * h_dim1];
3862 /* L210: */
3863                 }
3864             }
3865             h__[*n + i__ * h_dim1] = (h__[i__ + *n * h_dim1] - temp) / h__[
3866                     i__ + i__ * h_dim1];
3867 /* L220: */
3868         }
3869         temp = 0.;
3870         i__1 = *n - 1;
3871         for (i__ = 1; i__ <= i__1; ++i__) {
3872             temp += h__[*n + i__ * h_dim1] * h__[*n + i__ * h_dim1];
3873 /* L224: */
3874         }
3875         h__[*n + *n * h_dim1] = diag[*n] - temp + addmax;
3876         if (h__[*n + *n * h_dim1] > 0.) {
3877             h__[*n + *n * h_dim1] = sqrt(h__[*n + *n * h_dim1]);
3878         } else {
3879 /*     AFTER ADDING THE LAST COLUMN AND ROW */
3880 /*     QTHQ IS NOT P.D., NEED TO REDO CHOLESKY DECOMPOSITION */
3881             i__1 = *n;
3882             for (i__ = 2; i__ <= i__1; ++i__) {
3883                 i__2 = i__ - 1;
3884                 for (j = 1; j <= i__2; ++j) {
3885                     h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
3886 /* L230: */
3887                 }
3888                 h__[i__ + i__ * h_dim1] = diag[i__];
3889 /* L232: */
3890             }
3891             h__[h_dim1 + 1] = diag[1];
3892             choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
3893                     diag[1], &addmax);
3894         }
3895     }
3896
3897 /*   SOLVE QTHQ*QT*W3=-QT*G, WHERE W3 IS NEWTON STEP */
3898 /*   W2=-QT*G */
3899     i__1 = *n;
3900     for (i__ = 1; i__ <= i__1; ++i__) {
3901         w2[i__] = 0.;
3902         i__2 = *n;
3903         for (j = 1; j <= i__2; ++j) {
3904             w2[i__] += u[j] * u[i__] * g[j] / uu;
3905 /* L300: */
3906         }
3907         w2[i__] -= g[i__];
3908 /* L302: */
3909     }
3910     forslv_(nr, n, &h__[h_offset], &w3[1], &w2[1]);
3911     bakslv_(nr, n, &h__[h_offset], &w2[1], &w3[1]);
3912 /*     W2=QT*W3 => W3=Q*W2 --- NEWTON STEP */
3913     i__1 = *n;
3914     for (i__ = 1; i__ <= i__1; ++i__) {
3915         w3[i__] = 0.;
3916         i__2 = *n;
3917         for (j = 1; j <= i__2; ++j) {
3918             w3[i__] += u[i__] * u[j] * w2[j] / uu;
3919 /* L310: */
3920         }
3921         w3[i__] = w2[i__] - w3[i__];
3922 /* L312: */
3923     }
3924     return 0;
3925 } /* slvmdl_ */
3926
3927 /*  ---------------------- */
3928 /*  |  O P T S T P       | */
3929 /*  ---------------------- */
3930 /* Subroutine */ int optstp_(integer *n, doublereal *xpls, doublereal *fpls,
3931         doublereal *gpls, doublereal *x, integer *itncnt, integer *icscmx,
3932         integer *itrmcd, doublereal *gradtl, doublereal *steptl, doublereal *
3933         fscale, integer *itnlim, integer *iretcd, logical *mxtake, integer *
3934         ipr, integer *msg, doublereal *rgx, doublereal *rsx)
3935 {
3936     /* Format strings */
3937     static char fmt_900[] = "(\002 OPTSTP    STEP OF MAXIMUM LENGTH (STEPMX)"
3938             " TAKEN\002)";
3939     static char fmt_901[] = "(\002 OPTSTP    RELATIVE GRADIENT CLOSE TO ZE"
3940             "RO.\002/\002 OPTSTP    CURRENT ITERATE IS PROBABLY SOLUTION.\002)"
3941             ;
3942     static char fmt_902[] = "(\002 OPTSTP    SUCCESSIVE ITERATES WITHIN TOLE"
3943             "RANCE.\002/\002 OPTSTP    CURRENT ITERATE IS PROBABLY SOLUTION"
3944             ".\002)";
3945     static char fmt_903[] = "(\002 OPTSTP    LAST GLOBAL STEP FAILED TO LOCA"
3946             "TE A POINT\002,\002 LOWER THAN X.\002/\002 OPTSTP    EITHER X IS"
3947             " AN APPROXIMATE LOCAL MINIMUM\002,\002 OF THE FUNCTION,\002/\002"
3948             " OPTSTP    THE FUNCTION IS TOO NON-LINEAR FOR THIS\002,\002 ALGO"
3949             "RITHM,\002/\002 OPTSTP    OR STEPTL IS TOO LARGE.\002)";
3950     static char fmt_904[] = "(\002 OPTSTP    ITERATION LIMIT EXCEEDED.\002"
3951             "/\002 OPTSTP    ALGORITHM FAILED.\002)";
3952     static char fmt_905[] = "(\002 OPTSTP    MAXIMUM STEP SIZE EXCEEDED 5"
3953             "\002,\002 CONSECUTIVE TIMES.\002/\002 OPTSTP    EITHER THE FUNCT"
3954             "ION IS UNBOUNDED BELOW,\002/\002 OPTSTP    BECOMES ASYMPTOTIC TO"
3955             " A FINITE VALUE\002,\002 FROM ABOVE IN SOME DIRECTION,\002/\002 "
3956             "OPTSTP    OR STEPMX IS TOO SMALL\002)";
3957
3958     /* System generated locals */
3959     integer i__1;
3960     doublereal d__1, d__2, d__3;
3961
3962     /* Builtin functions */
3963     integer s_wsfe(cilist *), e_wsfe(void);
3964
3965     /* Local variables */
3966     doublereal d__;
3967     integer i__, imsg;
3968     doublereal relgrd;
3969     integer jtrmcd;
3970     doublereal relstp;
3971
3972     /* Fortran I/O blocks */
3973     static cilist io___280 = { 0, 0, 0, fmt_900, 0 };
3974     static cilist io___281 = { 0, 0, 0, fmt_901, 0 };
3975     static cilist io___282 = { 0, 0, 0, fmt_902, 0 };
3976     static cilist io___283 = { 0, 0, 0, fmt_903, 0 };
3977     static cilist io___284 = { 0, 0, 0, fmt_904, 0 };
3978     static cilist io___285 = { 0, 0, 0, fmt_905, 0 };
3979
3980
3981
3982 /* UNCONSTRAINED MINIMIZATION STOPPING CRITERIA */
3983 /* -------------------------------------------- */
3984 /* FIND WHETHER THE ALGORITHM SHOULD TERMINATE, DUE TO ANY */
3985 /* OF THE FOLLOWING: */
3986 /* 1) PROBLEM SOLVED WITHIN USER TOLERANCE */
3987 /* 2) CONVERGENCE WITHIN USER TOLERANCE */
3988 /* 3) ITERATION LIMIT REACHED */
3989 /* 4) DIVERGENCE OR TOO RESTRICTIVE MAXIMUM STEP (STEPMX) SUSPECTED */
3990
3991
3992 /* PARAMETERS */
3993 /* ---------- */
3994 /* N            --> DIMENSION OF PROBLEM */
3995 /* XPLS(N)      --> NEW ITERATE X[K] */
3996 /* FPLS         --> FUNCTION VALUE AT NEW ITERATE F(XPLS) */
3997 /* GPLS(N)      --> GRADIENT AT NEW ITERATE, G(XPLS), OR APPROXIMATE */
3998 /* X(N)         --> OLD ITERATE X[K-1] */
3999 /* ITNCNT       --> CURRENT ITERATION K */
4000 /* ICSCMX      <--> NUMBER CONSECUTIVE STEPS .GE. STEPMX */
4001 /*                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS] */
4002 /* ITRMCD      <--  TERMINATION CODE */
4003 /* GRADTL       --> TOLERANCE AT WHICH RELATIVE GRADIENT CONSIDERED CLOSE */
4004 /*                  ENOUGH TO ZERO TO TERMINATE ALGORITHM */
4005 /* STEPTL       --> RELATIVE STEP SIZE AT WHICH SUCCESSIVE ITERATES */
4006 /*                  CONSIDERED CLOSE ENOUGH TO TERMINATE ALGORITHM */
4007 /* FSCALE       --> ESTIMATE OF SCALE OF OBJECTIVE FUNCTION */
4008 /* ITNLIM       --> MAXIMUM NUMBER OF ALLOWABLE ITERATIONS */
4009 /* IRETCD       --> RETURN CODE */
4010 /* MXTAKE       --> BOOLEAN FLAG INDICATING STEP OF MAXIMUM LENGTH USED */
4011 /* IPR          --> DEVICE TO WHICH TO SEND OUTPUT */
4012 /* MSG         <--> CONTROL OUTPUT ON INPUT AND CONTAIN STOPPING */
4013 /*                  CONDITION ON OUTPUT */
4014
4015
4016
4017     /* Parameter adjustments */
4018     --x;
4019     --gpls;
4020     --xpls;
4021
4022     /* Function Body */
4023     *itrmcd = 0;
4024     imsg = *msg;
4025     *rgx = 0.;
4026     *rsx = 0.;
4027
4028 /* LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X */
4029     if (*iretcd != 1) {
4030         goto L50;
4031     }
4032 /*     IF(IRETCD.EQ.1) */
4033 /*     THEN */
4034     jtrmcd = 3;
4035     goto L600;
4036 /*     ENDIF */
4037 L50:
4038
4039 /* FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM. */
4040 /* CHECK WHETHER WITHIN TOLERANCE */
4041
4042 /* Computing MAX */
4043     d__1 = abs(*fpls);
4044     d__ = max(d__1,*fscale);
4045 /*     D=1D0 */
4046     *rgx = 0.f;
4047     i__1 = *n;
4048     for (i__ = 1; i__ <= i__1; ++i__) {
4049 /* Computing MAX */
4050         d__3 = (d__2 = xpls[i__], abs(d__2));
4051         relgrd = (d__1 = gpls[i__], abs(d__1)) * max(d__3,1.) / d__;
4052         *rgx = max(*rgx,relgrd);
4053 /* L100: */
4054     }
4055     jtrmcd = 1;
4056     if (*rgx <= *gradtl) {
4057         goto L600;
4058     }
4059
4060     if (*itncnt == 0) {
4061         return 0;
4062     }
4063
4064 /* FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM */
4065 /* CHECK WHETHER WITHIN TOLERANCE. */
4066
4067     *rsx = 0.f;
4068     i__1 = *n;
4069     for (i__ = 1; i__ <= i__1; ++i__) {
4070 /* Computing MAX */
4071         d__3 = (d__2 = xpls[i__], abs(d__2));
4072         relstp = (d__1 = xpls[i__] - x[i__], abs(d__1)) / max(d__3,1.);
4073         *rsx = max(*rsx,relstp);
4074 /* L120: */
4075     }
4076     jtrmcd = 2;
4077     if (*rsx <= *steptl) {
4078         goto L600;
4079     }
4080
4081 /* CHECK ITERATION LIMIT */
4082
4083     jtrmcd = 4;
4084     if (*itncnt >= *itnlim) {
4085         goto L600;
4086     }
4087
4088 /* CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX */
4089
4090     if (*mxtake) {
4091         goto L140;
4092     }
4093 /*     IF(.NOT.MXTAKE) */
4094 /*     THEN */
4095     *icscmx = 0;
4096     return 0;
4097 /*     ELSE */
4098 L140:
4099 /*       IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900) */
4100     if (imsg >= 1) {
4101         io___280.ciunit = *ipr;
4102         s_wsfe(&io___280);
4103         e_wsfe();
4104     }
4105     ++(*icscmx);
4106     if (*icscmx < 5) {
4107         return 0;
4108     }
4109     jtrmcd = 5;
4110 /*     ENDIF */
4111
4112
4113 /* PRINT TERMINATION CODE */
4114
4115 L600:
4116     *itrmcd = jtrmcd;
4117 /*     IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD */
4118     if (imsg >= 1) {
4119         switch (*itrmcd) {
4120             case 1:  goto L601;
4121             case 2:  goto L602;
4122             case 3:  goto L603;
4123             case 4:  goto L604;
4124             case 5:  goto L605;
4125         }
4126     }
4127     goto L700;
4128 L601:
4129     io___281.ciunit = *ipr;
4130     s_wsfe(&io___281);
4131     e_wsfe();
4132     goto L700;
4133 L602:
4134     io___282.ciunit = *ipr;
4135     s_wsfe(&io___282);
4136     e_wsfe();
4137     goto L700;
4138 L603:
4139     io___283.ciunit = *ipr;
4140     s_wsfe(&io___283);
4141     e_wsfe();
4142     goto L700;
4143 L604:
4144     io___284.ciunit = *ipr;
4145     s_wsfe(&io___284);
4146     e_wsfe();
4147     goto L700;
4148 L605:
4149     io___285.ciunit = *ipr;
4150     s_wsfe(&io___285);
4151     e_wsfe();
4152
4153 L700:
4154     *msg = -(*itrmcd);
4155     return 0;
4156
4157 } /* optstp_ */
4158
4159
4160 /* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,
4161         doublereal *dy, integer *incy)
4162 {
4163     /* System generated locals */
4164     integer i__1;
4165
4166     /* Local variables */
4167     integer i__, m, ix, iy, mp1;
4168
4169
4170 /*     COPIES A VECTOR, X, TO A VECTOR, Y. */
4171 /*     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
4172 /*     JACK DONGARRA, LINPACK, 3/11/78. */
4173
4174
4175     /* Parameter adjustments */
4176     --dy;
4177     --dx;
4178
4179     /* Function Body */
4180     if (*n <= 0) {
4181         return 0;
4182     }
4183     if (*incx == 1 && *incy == 1) {
4184         goto L20;
4185     }
4186
4187 /*        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
4188 /*          NOT EQUAL TO 1 */
4189
4190     ix = 1;
4191     iy = 1;
4192     if (*incx < 0) {
4193         ix = (-(*n) + 1) * *incx + 1;
4194     }
4195     if (*incy < 0) {
4196         iy = (-(*n) + 1) * *incy + 1;
4197     }
4198     i__1 = *n;
4199     for (i__ = 1; i__ <= i__1; ++i__) {
4200         dy[iy] = dx[ix];
4201         ix += *incx;
4202         iy += *incy;
4203 /* L10: */
4204     }
4205     return 0;
4206
4207 /*        CODE FOR BOTH INCREMENTS EQUAL TO 1 */
4208
4209
4210 /*        CLEAN-UP LOOP */
4211
4212 L20:
4213     m = *n % 7;
4214     if (m == 0) {
4215         goto L40;
4216     }
4217     i__1 = m;
4218     for (i__ = 1; i__ <= i__1; ++i__) {
4219         dy[i__] = dx[i__];
4220 /* L30: */
4221     }
4222     if (*n < 7) {
4223         return 0;
4224     }
4225 L40:
4226     mp1 = m + 1;
4227     i__1 = *n;
4228     for (i__ = mp1; i__ <= i__1; i__ += 7) {
4229         dy[i__] = dx[i__];
4230         dy[i__ + 1] = dx[i__ + 1];
4231         dy[i__ + 2] = dx[i__ + 2];
4232         dy[i__ + 3] = dx[i__ + 3];
4233         dy[i__ + 4] = dx[i__ + 4];
4234         dy[i__ + 5] = dx[i__ + 5];
4235         dy[i__ + 6] = dx[i__ + 6];
4236 /* L50: */
4237     }
4238     return 0;
4239 } /* dcopy_ */
4240
4241
4242 doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
4243         integer *incy)
4244 {
4245     /* System generated locals */
4246     integer i__1;
4247     doublereal ret_val;
4248
4249     /* Local variables */
4250     integer i__, m, ix, iy, mp1;
4251     doublereal dtemp;
4252
4253
4254 /*     FORMS THE DOT PRODUCT OF TWO VECTORS. */
4255 /*     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
4256 /*     JACK DONGARRA, LINPACK, 3/11/78. */
4257
4258
4259     /* Parameter adjustments */
4260     --dy;
4261     --dx;
4262
4263     /* Function Body */
4264     ret_val = 0.;
4265     dtemp = 0.;
4266     if (*n <= 0) {
4267         return ret_val;
4268     }
4269     if (*incx == 1 && *incy == 1) {
4270         goto L20;
4271     }
4272
4273 /*        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
4274 /*          NOT EQUAL TO 1 */
4275
4276     ix = 1;
4277     iy = 1;
4278     if (*incx < 0) {
4279         ix = (-(*n) + 1) * *incx + 1;
4280     }
4281     if (*incy < 0) {
4282         iy = (-(*n) + 1) * *incy + 1;
4283     }
4284     i__1 = *n;
4285     for (i__ = 1; i__ <= i__1; ++i__) {
4286         dtemp += dx[ix] * dy[iy];
4287         ix += *incx;
4288         iy += *incy;
4289 /* L10: */
4290     }
4291     ret_val = dtemp;
4292     return ret_val;
4293
4294 /*        CODE FOR BOTH INCREMENTS EQUAL TO 1 */
4295
4296
4297 /*        CLEAN-UP LOOP */
4298
4299 L20:
4300     m = *n % 5;
4301     if (m == 0) {
4302         goto L40;
4303     }
4304     i__1 = m;
4305     for (i__ = 1; i__ <= i__1; ++i__) {
4306         dtemp += dx[i__] * dy[i__];
4307 /* L30: */
4308     }
4309     if (*n < 5) {
4310         goto L60;
4311     }
4312 L40:
4313     mp1 = m + 1;
4314     i__1 = *n;
4315     for (i__ = mp1; i__ <= i__1; i__ += 5) {
4316         dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
4317                 i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
4318                 4] * dy[i__ + 4];
4319 /* L50: */
4320     }
4321 L60:
4322     ret_val = dtemp;
4323     return ret_val;
4324 } /* ddot_ */
4325
4326
4327 /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,
4328         integer *incx)
4329 {
4330     /* System generated locals */
4331     integer i__1, i__2;
4332
4333     /* Local variables */
4334     integer i__, m, mp1, nincx;
4335
4336
4337 /*     SCALES A VECTOR BY A CONSTANT. */
4338 /*     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
4339 /*     JACK DONGARRA, LINPACK, 3/11/78. */
4340
4341
4342     /* Parameter adjustments */
4343     --dx;
4344
4345     /* Function Body */
4346     if (*n <= 0) {
4347         return 0;
4348     }
4349     if (*incx == 1) {
4350         goto L20;
4351     }
4352
4353 /*        CODE FOR INCREMENT NOT EQUAL TO 1 */
4354
4355     nincx = *n * *incx;
4356     i__1 = nincx;
4357     i__2 = *incx;
4358     for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
4359         dx[i__] = *da * dx[i__];
4360 /* L10: */
4361     }
4362     return 0;
4363
4364 /*        CODE FOR INCREMENT EQUAL TO 1 */
4365
4366
4367 /*        CLEAN-UP LOOP */
4368
4369 L20:
4370     m = *n % 5;
4371     if (m == 0) {
4372         goto L40;
4373     }
4374     i__2 = m;
4375     for (i__ = 1; i__ <= i__2; ++i__) {
4376         dx[i__] = *da * dx[i__];
4377 /* L30: */
4378     }
4379     if (*n < 5) {
4380         return 0;
4381     }
4382 L40:
4383     mp1 = m + 1;
4384     i__2 = *n;
4385     for (i__ = mp1; i__ <= i__2; i__ += 5) {
4386         dx[i__] = *da * dx[i__];
4387         dx[i__ + 1] = *da * dx[i__ + 1];
4388         dx[i__ + 2] = *da * dx[i__ + 2];
4389         dx[i__ + 3] = *da * dx[i__ + 3];
4390         dx[i__ + 4] = *da * dx[i__ + 4];
4391 /* L50: */
4392     }
4393     return 0;
4394 } /* dscal_ */
4395
4396 /* Main program alias */ int driver_ () { MAIN__ (); return 0; }