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
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
10 http://www.netlib.org/f2c/libf2c.zip
15 /* Table of constant values */
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.;
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. */
31 /* Main program */ int MAIN__(void)
34 static char fmt_211[] = "(\002 G=\002,999e20.13)";
36 /* System generated locals */
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);
50 doublereal h__[200] /* was [50][4] */;
54 extern /* Subroutine */ int fcn_(), grd_();
56 extern /* Subroutine */ int hsn_();
58 doublereal wrk[400] /* was [50][8] */, fpls, gpls[50];
60 doublereal xpls[50], typx[50];
63 integer grdflg, hesflg;
66 extern /* Subroutine */ int dfault_(integer *, doublereal *, doublereal *,
67 doublereal *, doublereal *, integer *, doublereal *, integer *,
68 integer *, integer *, integer *, integer *, integer *);
70 extern /* Subroutine */ int prtfcn_(integer *);
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;
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 };
114 o__1.ofnm = "wood.inp";
124 for (i__ = 1; i__ <= i__1; ++i__) {
125 do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
129 do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
134 for (i__ = 1; i__ <= i__1; ++i__) {
135 do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
140 do_lio(&c__9, &c__1, " ", (ftnlen)1);
143 do_lio(&c__9, &c__1, "CALLING TENSRD - ALL INPUT PARAMETERS ARE SET", (
147 do_lio(&c__9, &c__1, " TO THEIR DEFAULT VALUES.", (ftnlen)
151 do_lio(&c__9, &c__1, " OUTPUT WILL BE WRITTEN TO THE STANDARD OUTPUT.", (
155 do_lio(&c__9, &c__1, " ", (ftnlen)1);
157 tensrd_(&nr, &n, x, (U_fp)fcn_, &msg, xpls, &fpls, gpls, h__, &itnno, wrk,
160 do_lio(&c__9, &c__1, "RESULTS OF TENSRD:", (ftnlen)18);
164 for (i__ = 1; i__ <= i__1; ++i__) {
165 do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
169 do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
171 for (i__ = 1; i__ <= i__1; ++i__) {
172 do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
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));
191 for (i__ = 1; i__ <= i__1; ++i__) {
192 do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
196 do_lio(&c__9, &c__1, "INITIAL APPROXIMATION TO THE SOLUTION X*:", (ftnlen)
201 for (i__ = 1; i__ <= i__1; ++i__) {
202 do_lio(&c__5, &c__1, (char *)&x[i__ - 1], (ftnlen)sizeof(doublereal));
206 do_lio(&c__9, &c__1, " ", (ftnlen)1);
209 do_lio(&c__9, &c__1, "CALLING TENSOR AFTER RESETTING INPUT PARAMETERS", (
213 do_lio(&c__9, &c__1, "IPR, MSG AND METHOD.", (ftnlen)20);
216 do_lio(&c__9, &c__1, "THE STANDARD METHOD IS USED IN THIS EXAMPLE.", (
220 do_lio(&c__9, &c__1, "ADDITIONAL OUTPUT WILL BE WRITTEN TO FILE 'FORT8'.",
224 do_lio(&c__9, &c__1, " ", (ftnlen)1);
227 dfault_(&n, typx, &fscale, &gradtl, &steptl, &itnlim, &stepmx, &ipr, &
228 method, &grdflg, &hesflg, &ndigit, &msg);
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);
246 do_lio(&c__9, &c__1, "RESULTS OF TENSOR:", (ftnlen)18);
250 for (i__ = 1; i__ <= i__1; ++i__) {
251 do_fio(&c__1, (char *)&gpls[i__ - 1], (ftnlen)sizeof(doublereal));
255 do_lio(&c__9, &c__1, "XPLS=", (ftnlen)5);
257 for (i__ = 1; i__ <= i__1; ++i__) {
258 do_lio(&c__5, &c__1, (char *)&xpls[i__ - 1], (ftnlen)sizeof(
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));
281 s_stop("", (ftnlen)0);
285 /* Subroutine */ int fcn_(integer *n, doublereal *x, doublereal *f)
287 /* System generated locals */
288 doublereal d__1, d__2, d__3, d__4, d__5, d__6;
290 /* Builtin functions */
291 double pow_dd(doublereal *, doublereal *);
293 /* Parameter adjustments */
300 d__1 = x[1] * x[1] - x[2];
302 d__3 = x[3] * x[3] - 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]) *
314 /* Subroutine */ int prtfcn_(integer *n)
316 /* Builtin functions */
317 integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen),
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 };
328 do_lio(&c__9, &c__1, "__________________________________________", (
332 do_lio(&c__9, &c__1, "f(x)= Wood function", (ftnlen)19);
335 do_lio(&c__9, &c__1, "__________________________________________", (
341 /* ---------------------- */
343 /* ---------------------- */
344 /* Subroutine */ int grd_(integer *n, real *x, real *g)
350 /* ---------------------- */
352 /* ---------------------- */
353 /* Subroutine */ int hsn_(integer *nr, integer *n, real *x, real *h__)
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,
371 /* System generated locals */
372 integer h_dim1, h_offset, wrk_dim1, wrk_offset;
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 *);
384 /* AUTHORS: TA-TUNG CHOW, ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
386 /* DATE: MARCH 29, 1991 */
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 */
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 */
401 /* ----------------------------------------------------------------------- */
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 */
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 */
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) */
445 /* Parameter adjustments */
447 wrk_offset = 1 + wrk_dim1;
451 h_offset = 1 + h_dim1;
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]);
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)
475 /* System generated locals */
476 integer h_dim1, h_offset, wrk_dim1, wrk_offset;
478 /* Local variables */
479 extern /* Subroutine */ int grd_(integer *, real *, real *), hsn_(integer
480 *, integer *, real *, real *);
483 integer grdflg, hesflg;
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 *,
496 doublereal steptl, stepmx;
500 /* SHORT CALLING SEQUENCE SUBROUTINE */
502 /* ----------------------------------------------------------------------- */
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 */
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) */
529 /* Parameter adjustments */
531 wrk_offset = 1 + wrk_dim1;
535 h_offset = 1 + h_dim1;
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]);
551 /* ---------------------- */
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)
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"
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)";
590 /* System generated locals */
591 integer h_dim1, h_offset, i__1, i__2;
592 doublereal d__1, d__2;
594 /* Builtin functions */
595 double pow_di(doublereal *, integer *), sqrt(doublereal);
596 integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
598 /* Local variables */
601 doublereal gd, fn, fp, rnf, eps, rgx, rsx, beta;
602 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
605 doublereal temp, alpha;
606 extern /* Subroutine */ int mkmdl_(integer *, integer *, doublereal *,
607 doublereal *, doublereal *, doublereal *, doublereal *,
608 doublereal *, doublereal *, doublereal *, doublereal *,
612 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
613 doublereal *, integer *);
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 *,
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 *);
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 *);
637 extern /* Subroutine */ int optchk_(integer *, integer *, integer *,
638 doublereal *, doublereal *, doublereal *, doublereal *,
639 doublereal *, integer *, integer *, doublereal *, integer *,
640 integer *, integer *, doublereal *, integer *);
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 *);
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 };
681 /* DERIVATIVE TENSOR METHODS FOR UNCONSTRAINED OPTIMIZATION */
683 /* ----------------------------------------------------------------------- */
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 */
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 */
725 /* -------------------------------- */
726 /* INITIALIZATION | */
727 /* -------------------------------- */
729 /* Parameter adjustments */
739 h_offset = 1 + h_dim1;
751 optchk_(nr, n, msg, &x[1], &typx[1], fscale, gradtl, steptl, itnlim,
752 ndigit, &eps, method, grdflg, hesflg, stepmx, ipr);
759 for (i__ = 1; i__ <= i__1; ++i__) {
764 /* -------------------------------- */
765 /* INITIAL ITERATION | */
766 /* -------------------------------- */
768 /* COMPUTE TYPICAL SIZE OF X */
770 for (i__ = 1; i__ <= i__1; ++i__) {
771 dn[i__] = 1. / typx[i__];
776 d__1 = pow_di(&c_b134, &i__1);
779 d__1 = .01, d__2 = sqrt(rnf);
780 analtl = max(d__1,d__2);
781 /* UNSCALE X AND COMPUTE F AND G */
783 for (i__ = 1; i__ <= i__1; ++i__) {
784 wk1[i__] = x[i__] * typx[i__];
787 (*fcn)(n, &wk1[1], &f);
790 (*grd)(n, &wk1[1], &g[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);
797 fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, &f, &g[1], &typx[1], &
798 rnf, &wk2[1], &c__1, &nfcnt);
804 gnorm = twonrm_(n, &g[1]);
806 /* PRINT OUT 1ST ITERATION? */
808 io___70.ciunit = *ipr;
810 do_fio(&c__1, (char *)&f, (ftnlen)sizeof(doublereal));
812 io___71.ciunit = *ipr;
815 for (i__ = 1; i__ <= i__1; ++i__) {
816 do_fio(&c__1, (char *)&g[i__], (ftnlen)sizeof(doublereal));
821 /* TEST WHETHER THE INITIAL GUESS SATISFIES THE STOPPING CRITERIA */
822 if (gnorm <= *gradtl) {
825 for (i__ = 1; i__ <= i__1; ++i__) {
830 optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
831 gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &
837 /* ------------------------ */
839 /* ------------------------ */
843 /* COMPUTE HESSIAN */
846 fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
847 typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
849 sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
850 rnf, &wk2[1], &d__[1], &nfcnt);
854 (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
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);
868 for (i__ = 2; i__ <= i__1; ++i__) {
870 dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
874 /* CHOLESKY DECOMPOSITION FOR H (H=LLT) */
875 choldr_(nr, n, &h__[h_offset], &d__[1], &eps, &pivot[1], &e[1], &wk1[1], &
878 /* SOLVE FOR NEWTON STEP D */
880 for (i__ = 1; i__ <= i__1; ++i__) {
884 forslv_(nr, n, &h__[h_offset], &wk1[1], &wk2[1]);
885 bakslv_(nr, n, &h__[h_offset], &d__[1], &wk1[1]);
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);
892 /* CALL DCOPY(N,GPLS(1),1,GP(1),1) */
894 /* UNSCALE XPLS AND COMPUTE GPLS */
896 for (i__ = 1; i__ <= i__1; ++i__) {
897 wk1[i__] = xpls[i__] * typx[i__];
901 (*grd)(n, &wk1[1], &gpls[1]);
903 fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
904 &rnf, &wk2[1], &c__1, &nfcnt);
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,
912 /* IF ITRMCD > 0 THEN STOPPING CONDITIONS SATISFIED */
917 /* UPDATE X,F AND S FOR TENSOR MODEL */
921 for (i__ = 1; i__ <= i__1; ++i__) {
923 s[i__] = x[i__] - temp;
928 /* IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
930 io___81.ciunit = *ipr;
932 do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
934 io___82.ciunit = *ipr;
937 for (i__ = 1; i__ <= i__1; ++i__) {
938 do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
941 io___83.ciunit = *ipr;
943 do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
945 io___84.ciunit = *ipr;
948 for (i__ = 1; i__ <= i__1; ++i__) {
949 do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
954 io___85.ciunit = *ipr;
956 do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
957 do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
959 io___86.ciunit = *ipr;
961 do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
965 /* ------------------------ */
966 /* ITERATION > 1 | */
967 /* ------------------------ */
969 /* UNSCALE X AND COMPUTE H */
972 for (i__ = 1; i__ <= i__1; ++i__) {
973 wk1[i__] = x[i__] * typx[i__];
978 fstofd_(nr, n, n, &wk1[1], (S_fp)grd, &g[1], &h__[h_offset], &
979 typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
981 sndofd_(nr, n, &wk1[1], (S_fp)fcn, &f, &h__[h_offset], &typx[1], &
982 rnf, &wk2[1], &d__[1], &nfcnt);
985 (*hsn)(nr, n, &wk1[1], &h__[h_offset]);
988 for (i__ = 2; i__ <= i__1; ++i__) {
990 dcopy_(&i__2, &h__[i__ + h_dim1], nr, &h__[i__ * h_dim1 + 1], &c__1);
994 /* IF METHOD = 0 THEN USE NEWTON STEP ONLY */
997 /* CHOLESKY DECOMPOSITION FOR H */
998 choldr_(nr, n, &h__[h_offset], &wk2[1], &eps, &pivot[1], &e[1], &wk1[
1001 /* COMPUTE NEWTON STEP */
1003 for (i__ = 1; i__ <= i__1; ++i__) {
1004 wk1[i__] = -gpls[i__];
1007 forslv_(nr, n, &h__[h_offset], &wk2[1], &wk1[1]);
1008 bakslv_(nr, n, &h__[h_offset], &d__[1], &wk2[1]);
1010 /* NO TENSOR STEP */
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]);
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, &
1029 /* IF TENSOR MODEL HAS NO MINIMIZER THEN USE NEWTON STEP */
1031 dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
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);
1038 dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
1044 dcopy_(n, &gpls[1], &c__1, &g[1], &c__1);
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);
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);
1057 /* COMPARE TENSOR STEP TO NEWTON STEP */
1058 /* IF NEWTON STEP IS BETTER, SET NEXT ITERATE TO NEW NEWTON POINT */
1061 dcopy_(n, &dn[1], &c__1, &d__[1], &c__1);
1062 dcopy_(n, &wk2[1], &c__1, &xpls[1], &c__1);
1066 for (i__ = 1; i__ <= i__1; ++i__) {
1067 d__[i__] = xpls[i__] - x[i__];
1071 /* UNSCALE XPLS, AND COMPUTE FPLS AND GPLS */
1073 for (i__ = 1; i__ <= i__1; ++i__) {
1074 wk1[i__] = xpls[i__] * typx[i__];
1077 (*fcn)(n, &wk1[1], fpls);
1080 (*grd)(n, &wk1[1], &gpls[1]);
1082 fstofd_(&c__1, &c__1, n, &wk1[1], (S_fp)fcn, fpls, &gpls[1], &typx[1],
1083 &rnf, &wk2[1], &c__1, &nfcnt);
1086 /* CHECK STOPPING CONDITIONS */
1088 optstp_(n, &xpls[1], fpls, &gpls[1], &x[1], itnno, &icscmx, &itrmcd,
1089 gradtl, steptl, fscale, itnlim, &iretcd, &mxtake, ipr, msg, &rgx,
1092 /* IF ITRMCD = 0 THEN NOT OVER YET */
1097 /* IF MSG >= 1 THEN PRINT OUT FINAL ITERATION */
1101 /* TRANSFORM X BACK TO ORIGINAL SPACE */
1103 for (i__ = 1; i__ <= i__1; ++i__) {
1104 xpls[i__] *= typx[i__];
1107 io___93.ciunit = *ipr;
1110 for (i__ = 1; i__ <= i__1; ++i__) {
1111 do_fio(&c__1, (char *)&xpls[i__], (ftnlen)sizeof(doublereal));
1114 io___94.ciunit = *ipr;
1117 for (i__ = 1; i__ <= i__1; ++i__) {
1118 do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
1121 io___95.ciunit = *ipr;
1123 do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
1124 do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
1127 io___96.ciunit = *ipr;
1129 do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
1130 do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
1132 io___97.ciunit = *ipr;
1134 do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
1137 /* UPDATE THE HESSIAN */
1140 fstofd_(nr, n, n, &xpls[1], (S_fp)grd, &gpls[1], &h__[
1141 h_offset], &typx[1], &rnf, &wk2[1], &c__3, &nfcnt);
1143 sndofd_(nr, n, &xpls[1], (S_fp)fcn, fpls, &h__[h_offset], &
1144 typx[1], &rnf, &wk2[1], &d__[1], &nfcnt);
1147 (*hsn)(nr, n, &xpls[1], &h__[h_offset]);
1152 /* UPDATE INFORMATION AT THE CURRENT POINT */
1154 dcopy_(n, &xpls[1], &c__1, &x[1], &c__1);
1156 for (i__ = 1; i__ <= i__1; ++i__) {
1161 /* IF TOO MANY ITERATIONS THEN RETURN */
1162 if (*itnno > *itnlim) {
1166 /* IF MSG >= 2 THEN PRINT OUT EACH ITERATION */
1168 io___98.ciunit = *ipr;
1170 do_fio(&c__1, (char *)&(*itnno), (ftnlen)sizeof(integer));
1172 io___99.ciunit = *ipr;
1175 for (i__ = 1; i__ <= i__1; ++i__) {
1176 do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(doublereal));
1179 io___100.ciunit = *ipr;
1181 do_fio(&c__1, (char *)&(*fpls), (ftnlen)sizeof(doublereal));
1183 io___101.ciunit = *ipr;
1186 for (i__ = 1; i__ <= i__1; ++i__) {
1187 do_fio(&c__1, (char *)&gpls[i__], (ftnlen)sizeof(doublereal));
1192 io___102.ciunit = *ipr;
1194 do_fio(&c__1, (char *)&nfcnt, (ftnlen)sizeof(integer));
1195 do_fio(&c__1, (char *)&rgx, (ftnlen)sizeof(doublereal));
1197 io___103.ciunit = *ipr;
1199 do_fio(&c__1, (char *)&rsx, (ftnlen)sizeof(doublereal));
1206 /* PERFORM NEXT ITERATION */
1209 /* END OF ITERATION > 1 */
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)
1221 /* System generated locals */
1224 /* Builtin functions */
1225 double pow_dd(doublereal *, doublereal *), d_lg10(doublereal *);
1227 /* Local variables */
1229 doublereal eps, temp;
1230 extern /* Subroutine */ int mcheps_(doublereal *);
1232 /* Parameter adjustments */
1242 for (i__ = 1; i__ <= i__1; ++i__) {
1246 temp = pow_dd(&eps, &c_b250);
1248 *steptl = temp * temp;
1249 *ndigit = (integer) (-d_lg10(&eps));
1250 /* SET ACTUAL DFAULT VALUE OF STEPMX IN OPTCHK */
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)
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)";
1272 /* System generated locals */
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);
1281 /* Local variables */
1283 doublereal temp, stpsiz;
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 };
1293 /* CHECK INPUT FOR REASONABLENESS */
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 */
1317 /* CHECK DIMENSION OF PROBLEM */
1319 /* Parameter adjustments */
1331 /* CHECK THAT PARAMETERS ONLY TAKE ON ACCEPTABLE VALUES. */
1332 /* IF NOT, SET THEM TO DEFAULT VALUES. */
1336 if (*grdflg != 2 && *grdflg != 1) {
1339 if (*hesflg != 2 && *hesflg != 1) {
1342 if (*msg > 3 || (real) (*msg) < 0.f) {
1346 /* COMPUTE SCALE MATRIX */
1349 for (i__ = 1; i__ <= i__1; ++i__) {
1350 if (typx[i__] == 0.f) {
1353 if (typx[i__] < 0.f) {
1354 typx[i__] = -typx[i__];
1359 /* CHECK MAXIMUM STEP SIZE */
1366 for (i__ = 1; i__ <= i__1; ++i__) {
1367 stpsiz += x[i__] * x[i__] / typx[i__] * typx[i__];
1370 stpsiz = sqrt(stpsiz);
1372 d__1 = stpsiz * 1e3;
1373 *stepmx = max(d__1,1e3);
1375 /* CHECK FUNCTION SCALE */
1376 if (*fscale == 0.f) {
1379 if (*fscale < 0.f) {
1380 *fscale = -(*fscale);
1382 /* CHECK GRADIENT TOLERANCE */
1383 if (*gradtl < 0.f) {
1384 *gradtl = pow_dd(eps, &c_b250);
1386 /* CHECK STEP TOLERANCE */
1387 if (*steptl < 0.f) {
1388 temp = pow_dd(eps, &c_b250);
1389 *steptl = temp * temp;
1392 /* CHECK ITERATION LIMIT */
1397 /* CHECK NUMBER OF DIGITS OF ACCURACY IN FUNCTION FCN */
1399 *ndigit = (integer) (-d_lg10(eps));
1402 if (pow_di(&c_b134, &i__1) <= *eps) {
1403 *ndigit = (integer) (-d_lg10(eps));
1410 io___110.ciunit = *ipr;
1412 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
1417 io___111.ciunit = *ipr;
1419 do_fio(&c__1, (char *)&(*nr), (ftnlen)sizeof(integer));
1420 do_fio(&c__1, (char *)&(*n), (ftnlen)sizeof(integer));
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)
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)";
1440 /* System generated locals */
1442 doublereal d__1, d__2, d__3, d__4;
1444 /* Builtin functions */
1445 integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
1447 /* Local variables */
1452 extern /* Subroutine */ int fstofd_(integer *, integer *, integer *,
1453 doublereal *, S_fp, doublereal *, doublereal *, doublereal *,
1454 doublereal *, doublereal *, integer *, integer *);
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 };
1464 /* CHECK ANALYTIC GRADIENT AGAINST ESTIMATED GRADIENT */
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 */
1486 /* COMPUTE FIRST ORDER FINITE DIFFERENCE GRADIENT AND COMPARE TO */
1487 /* ANALYTIC GRADIENT. */
1489 /* Parameter adjustments */
1497 fstofd_(&c__1, &c__1, n, &x[1], (S_fp)fcn, f, &wrk1[1], &typx[1], rnf, &
1501 for (i__ = 1; i__ <= i__1; ++i__) {
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);
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) {
1517 io___116.ciunit = *ipr;
1520 io___117.ciunit = *ipr;
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));
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)
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)";
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;
1554 /* Builtin functions */
1555 integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen);
1557 /* Local variables */
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 *,
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 };
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) */
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 */
1612 /* COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN. */
1614 /* Parameter adjustments */
1616 a_offset = 1 + a_dim1;
1628 fstofd_(nr, n, n, &x[1], (S_fp)grd, &g[1], &a[a_offset], &typx[1],
1629 rnf, &wrk1[1], &c__3, ifn);
1632 sndofd_(nr, n, &x[1], (S_fp)fcn, f, &a[a_offset], &typx[1], rnf, &
1633 wrk1[1], &wrk2[1], ifn);
1638 /* COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART */
1639 /* AND DIAGONAL OF "A" TO UDIAG */
1642 for (j = 1; j <= i__1; ++j) {
1643 udiag[j] = a[j + j * a_dim1];
1649 for (i__ = jp1; i__ <= i__2; ++i__) {
1650 a[j + i__ * a_dim1] = a[i__ + j * a_dim1];
1657 /* COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE */
1658 /* APPROXIMATION. */
1661 for (i__ = 1; i__ <= i__1; ++i__) {
1663 for (j = i__; j <= i__2; ++j) {
1664 a[j + i__ * a_dim1] = 0.;
1669 (*hsn)(nr, n, &x[1], &a[a_offset]);
1671 for (j = 1; j <= i__2; ++j) {
1673 d__3 = (d__1 = g[j], abs(d__1));
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);
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) *
1688 for (i__ = jp1; i__ <= i__1; ++i__) {
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) {
1704 io___123.ciunit = *ipr;
1708 for (i__ = 1; i__ <= i__2; ++i__) {
1714 for (j = 1; j <= i__1; ++j) {
1715 io___125.ciunit = *ipr;
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(
1721 do_fio(&c__1, (char *)&a[j + i__ * a_dim1], (ftnlen)sizeof(
1727 io___126.ciunit = *ipr;
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(
1733 do_fio(&c__1, (char *)&udiag[i__], (ftnlen)sizeof(doublereal));
1743 /* ----------------------- */
1744 /* | M C H E P S | */
1745 /* ----------------------- */
1746 /* Subroutine */ int mcheps_(doublereal *eps)
1748 doublereal temp, temp1;
1752 /* COMPUTE MACHINE PRECISION */
1754 /* ------------------------------------------------------------------------- */
1758 /* EPS <-- MACHINE PRECISION */
1772 doublereal twonrm_(integer *n, doublereal *v)
1774 /* System generated locals */
1777 /* Builtin functions */
1778 double sqrt(doublereal);
1780 /* Local variables */
1781 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
1787 /* COMPUTER L-2 NORM */
1789 /* -------------------------------------------------------------------------- */
1793 /* N --> DIMENSION OF PROBLEM */
1794 /* V(N) --> VECTOR WHICH L-2 NORM IS EVALUATED */
1796 /* Parameter adjustments */
1800 temp = ddot_(n, &v[1], &c__1, &v[1], &c__1);
1801 ret_val = sqrt(temp);
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 *
1814 /* System generated locals */
1816 doublereal d__1, d__2;
1818 /* Builtin functions */
1819 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
1821 /* Local variables */
1824 doublereal t1, t2, t3, scl, rln, sln, slp, tmp, disc;
1825 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
1827 doublereal temp, temp1, temp2, alpha;
1828 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
1830 doublereal pfpls, almbda, plmbda, tlmbda, almbmn;
1833 /* THE ALPHA CONDITION ONLY LINE SEARCH */
1837 /* FIND A NEXT NEWTON ITERATE BY LINE SEARCH. */
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 */
1859 /* INTERNAL VARIABLES */
1860 /* ------------------ */
1861 /* SLN NEWTON LENGTH */
1862 /* RLN RELATIVE LENGTH OF NEWTON STEP */
1865 /* Parameter adjustments */
1877 /* $ WRITE(IPR,954) */
1878 /* $ WRITE(IPR,955) (P(I),I=1,N) */
1881 for (i__ = 1; i__ <= i__1; ++i__) {
1882 tmp += p[i__] * p[i__];
1886 if (sln <= *stepmx) {
1890 /* NEWTON STEP LONGER THAN MAXIMUM ALLOWED */
1891 scl = *stepmx / sln;
1892 dscal_(n, &scl, &p[1], &c__1);
1894 /* $ WRITE(IPR,954) */
1895 /* $ WRITE(IPR,955) (P(I),I=1,N) */
1897 slp = ddot_(n, &g[1], &c__1, &p[1], &c__1);
1900 for (i__ = 1; i__ <= i__1; ++i__) {
1902 temp1 = (d__1 = x[i__], abs(d__1));
1903 temp2 = max(temp1,temp);
1904 temp1 = (d__1 = p[i__], abs(d__1));
1906 d__1 = rln, d__2 = temp1 / temp2;
1907 rln = max(d__1,d__2);
1910 almbmn = *steptl / rln;
1912 /* $ WRITE(IPR,952) SLN,SLP,RMNLMB,STEPMX,STEPTL */
1915 /* CHECK IF NEW ITERATE SATISFACTORY. GENERATE NEW LAMBDA IF NECESSARY. */
1922 for (i__ = 1; i__ <= i__1; ++i__) {
1923 xpls[i__] = x[i__] + almbda * p[i__];
1927 for (k = 1; k <= i__1; ++k) {
1928 w2[k] = xpls[k] * typx[k];
1931 (*fcn)(n, &w2[1], fpls);
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) {
1940 /* IF(FPLS.LE. F+SLP*1.E-4*ALMBDA) */
1943 /* SOLUTION FOUND */
1946 if (almbda == 1. && sln > *stepmx * .99) {
1951 /* SOLUTION NOT (YET) FOUND */
1955 if (almbda >= almbmn) {
1958 /* IF(ALMBDA .LT. ALMBMN) */
1961 /* NO SATISFACTORY XPLS FOUND SUFFICIENTLY DISTINCT FROM X */
1967 /* CALCULATE NEW LAMBDA */
1973 /* IF(ALMBDA.EQ.1.0) */
1976 /* FIRST BACKTRACK: QUADRATIC FIT */
1978 tlmbda = -slp / ((*fpls - *f - slp) * 2.);
1982 /* ALL SUBSEQUENT BACKTRACKS: CUBIC FIT */
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 *
1991 disc = b * b - a * 3.f * slp;
1992 if (disc <= b * b) {
1995 /* IF(DISC.GT. B*B) */
1998 /* ONLY ONE POSITIVE CRITICAL POINT, MUST BE MINIMUM */
2000 tlmbda = (-b + d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
2004 /* BOTH CRITICAL POINTS POSITIVE, FIRST IS MINIMUM */
2007 tlmbda = (-b - d_sign(&c_b324, &a) * sqrt(disc)) / (a * 3.f);
2010 if (tlmbda > almbda * .5) {
2011 tlmbda = almbda * .5;
2017 if (tlmbda >= almbda * .1) {
2020 /* IF(TLMBDA.LT.ALMBDA/10.) */
2040 /* ---------------------- */
2042 /* ---------------------- */
2043 /* Subroutine */ int zhz_(integer *nr, integer *n, doublereal *y, doublereal *
2044 h__, doublereal *u, doublereal *t)
2046 /* System generated locals */
2047 integer h_dim1, h_offset, i__1, i__2;
2050 /* Builtin functions */
2051 double sqrt(doublereal);
2053 /* Local variables */
2057 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2059 doublereal temp1, temp2;
2060 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
2061 doublereal *, integer *);
2066 /* COMPUTE QTHQ(N,N) AND ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND */
2067 /* FIRST N-1 COLUMNS OF QTHQ */
2069 /* --------------------------------------------------------------------------- */
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 */
2082 /* U=Y+SGN(Y(N))||Y||E(N) */
2083 /* Parameter adjustments */
2087 h_offset = 1 + h_dim1;
2093 sgn = y[*n] / (d__1 = y[*n], abs(d__1));
2097 ynorm = ddot_(n, &y[1], &c__1, &y[1], &c__1);
2098 ynorm = sqrt(ynorm);
2099 u[*n] = y[*n] + sgn * ynorm;
2101 dcopy_(&i__1, &y[1], &c__1, &u[1], &c__1);
2104 d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
2109 for (i__ = 1; i__ <= i__1; ++i__) {
2112 for (j = 1; j <= i__2; ++j) {
2113 t[i__] = h__[i__ + j * h_dim1] * u[j] + t[i__];
2120 /* S=4UHU/(UTU)**2 */
2121 s = ddot_(n, &u[1], &c__1, &t[1], &c__1);
2124 /* COMPUTE QTHQ (ZTHZ) */
2126 for (i__ = 1; i__ <= i__1; ++i__) {
2129 for (j = 1; j <= i__2; ++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;
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)
2146 /* System generated locals */
2147 integer al_dim1, al_offset, i__1, i__2;
2149 /* Local variables */
2152 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2154 extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
2155 doublereal *, doublereal *);
2159 /* SOLVE L*W=ZT*V */
2161 /* ---------------------------------------------------------------------- */
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 */
2172 /* FORM ZT*V (STORED IN B) */
2173 /* Parameter adjustments */
2175 al_offset = 1 + al_dim1;
2182 d__ = ddot_(n, &u[1], &c__1, &u[1], &c__1);
2185 for (i__ = 1; i__ <= i__1; ++i__) {
2188 for (j = 1; j <= i__2; ++j) {
2189 b[i__] += u[j] * u[i__] * w[j] / d__;
2192 b[i__] = w[i__] - b[i__];
2198 forslv_(nr, &i__1, &al[al_offset], &w[1], &b[1]);
2202 /* ---------------------- */
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__)
2209 /* System generated locals */
2210 integer al_dim1, al_offset, i__1;
2212 /* Local variables */
2214 doublereal utt, utu;
2215 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2218 extern /* Subroutine */ int bakslv_(integer *, integer *, doublereal *,
2219 doublereal *, doublereal *);
2223 /* COMPUTE TENSOR STEP D=SIGMA*S+ZT*T(SIGMA) */
2225 /* ------------------------------------------------------------------------ */
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 */
2239 /* Parameter adjustments */
2241 al_offset = 1 + al_dim1;
2252 d__[1] = *sigma * s[1];
2255 /* COMPUTE T(SIGMA)=-(ZTHZ)*ZT*(G+SIGMA*SH+SIGMA**2*A/2) (STORED IN D) */
2257 for (i__ = 1; i__ <= i__1; ++i__) {
2258 w2[i__] = w3[i__] + *sigma * w2[i__] + w1[i__] * .5 * *sigma * *
2263 bakslv_(nr, &i__1, &al[al_offset], &d__[1], &w2[1]);
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);
2271 for (i__ = 1; i__ <= i__1; ++i__) {
2272 d__[i__] = *sigma * s[i__] - (d__[i__] - u[i__] * 2. * temp);
2279 /* ---------------------- */
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,
2287 /* System generated locals */
2288 integer h_dim1, h_offset, i__1, i__2;
2290 /* Local variables */
2292 doublereal b1, b2, gs, gps, shs, sts;
2293 extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
2298 /* FORM TENSOR MODEL */
2300 /* ----------------------------------------------------------------------- */
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 */
2314 /* A(N) <-- A=2*(GP-G-SH-S*BETA/(6*STS)) */
2318 /* Parameter adjustments */
2322 h_offset = 1 + h_dim1;
2330 for (i__ = 1; i__ <= i__1; ++i__) {
2333 for (j = 1; j <= i__2; ++j) {
2334 sh[i__] += s[j] * h__[j + i__ * h_dim1];
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.;
2348 sts = ddot_(n, &s[1], &c__1, &s[1], &c__1);
2350 for (i__ = 1; i__ <= i__1; ++i__) {
2351 a[i__] = (gp[i__] - g[i__] - sh[i__] - s[i__] * *beta / (sts * 6.)) *
2358 /* ---------------------- */
2360 /* ---------------------- */
2361 /* Subroutine */ int sigma_(doublereal *sgstar, doublereal *a, doublereal *b,
2362 doublereal *c__, doublereal *d__)
2364 doublereal s1, s2, s3;
2365 extern /* Subroutine */ int roots_(doublereal *, doublereal *, doublereal
2366 *, doublereal *, doublereal *, doublereal *, doublereal *),
2367 sortrt_(doublereal *, doublereal *, doublereal *);
2371 /* COMPUTE DESIRABLE ROOT OF REDUCED ONE VARIABLE EQUATION */
2373 /* ------------------------------------------------------------------------- */
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 */
2383 /* COMPUTE ALL THREE ROOTS */
2384 roots_(&s1, &s2, &s3, a, b, c__, d__);
2387 sortrt_(&s1, &s2, &s3);
2389 /* CHOOSE DESIRABLE ROOT */
2398 if (s1 > 0. || s3 < 0.) {
2410 /* ---------------------- */
2412 /* ---------------------- */
2413 /* Subroutine */ int roots_(doublereal *s1, doublereal *s2, doublereal *s3,
2414 doublereal *a, doublereal *b, doublereal *c__, doublereal *d__)
2416 /* System generated locals */
2419 /* Builtin functions */
2420 double sqrt(doublereal), pow_dd(doublereal *, doublereal *), acos(
2421 doublereal), cos(doublereal);
2423 /* Local variables */
2424 doublereal q, r__, s, t, v, a1, a2, a3, pi, temp, theta;
2428 /* COMPUTE ROOT(S) OF 3RD ORDER EQUATION */
2430 /* --------------------------------------------------------------------------- */
2433 /* S1 <-- ROOT (IF THREE ROOTS ARE */
2434 /* S2 <-- ROOT EQUAL, THEN S1=S2=S3) */
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 */
2441 /* SET VALUE OF PI */
2442 pi = 3.141592653589793;
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__;
2454 t = -pow_dd(&d__1, &c_b250);
2456 t = pow_dd(&t, &c_b250);
2460 s = -pow_dd(&d__1, &c_b250);
2462 s = pow_dd(&s, &c_b250);
2464 *s1 = s + t - a1 / 3.;
2468 temp = r__ / sqrt(-pow_dd(&q, &c_b384));
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.;
2479 /* ---------------------- */
2480 /* | S O R T R T | */
2481 /* ---------------------- */
2482 /* Subroutine */ int sortrt_(doublereal *s1, doublereal *s2, doublereal *s3)
2488 /* SORT ROOTS INTO ASCENDING ORDER */
2490 /* ----------------------------------------------------------------------------- */
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)
2522 /* System generated locals */
2523 integer a_dim1, a_offset, i__1, i__2;
2524 doublereal d__1, d__2;
2526 /* Builtin functions */
2527 double sqrt(doublereal);
2529 /* Local variables */
2530 integer i__, j, jp1, nm1;
2531 doublereal xtmpj, stepsz;
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". */
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 */
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" */
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: */
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) */
2576 /* =3 OPTIMIZATION (HESSIAN) */
2577 /* NFCNT <--> NUMBER OF FUNCTION EVALUTIONS */
2579 /* INTERNAL VARIABLES */
2580 /* ------------------ */
2581 /* STEPSZ - STEPSIZE IN THE J-TH VARIABLE DIRECTION */
2584 /* FIND J-TH COLUMN OF A */
2585 /* EACH COLUMN IS DERIVATIVE OF F(FCN) WITH RESPECT TO XPLS(J) */
2587 /* Parameter adjustments */
2592 a_offset = 1 + a_dim1;
2598 for (j = 1; j <= i__1; ++j) {
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]);
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];
2619 /* IF COMPUTING HESSIAN, A MUST BE SYMMETRIC */
2626 for (j = 1; j <= i__1; ++j) {
2629 for (i__ = jp1; i__ <= i__2; ++i__) {
2630 a[i__ + j * a_dim1] = (a[i__ + j * a_dim1] + a[j + i__ * a_dim1])
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)
2646 /* System generated locals */
2647 integer a_dim1, a_offset, i__1, i__2;
2648 doublereal d__1, d__2;
2650 /* Builtin functions */
2651 double pow_dd(doublereal *, doublereal *);
2653 /* Local variables */
2654 integer i__, j, ip1;
2655 doublereal ov3, fhat, xtmpi, xtmpj;
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" */
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. */
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 */
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 */
2686 /* FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION */
2687 /* OF I-TH UNIT VECTOR. */
2689 /* Parameter adjustments */
2694 a_offset = 1 + a_dim1;
2699 ov3 = .33333333333333331;
2701 for (i__ = 1; i__ <= i__1; ++i__) {
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__]);
2713 /* CALCULATE COLUMN I OF A */
2716 for (i__ = 1; i__ <= i__1; ++i__) {
2718 xpls[i__] = xtmpi + stepsz[i__] * 2.f;
2719 (*fcn)(n, &xpls[1], &fhat);
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__];
2725 /* CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN */
2729 xpls[i__] = xtmpi + stepsz[i__];
2732 for (j = ip1; j <= i__2; ++j) {
2734 xpls[j] = xtmpj + stepsz[j];
2735 (*fcn)(n, &xpls[1], &fhat);
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];
2750 /* ---------------------- */
2751 /* | B A K S L V | */
2752 /* ---------------------- */
2753 /* Subroutine */ int bakslv_(integer *nr, integer *n, doublereal *a,
2754 doublereal *x, doublereal *b)
2756 /* System generated locals */
2757 integer a_dim1, a_offset, i__1;
2759 /* Local variables */
2760 integer i__, j, ip1;
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. */
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 */
2780 /* IF B IS NO LONGER REQUIRED BY CALLING ROUTINE, */
2781 /* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE. */
2784 /* SOLVE (L-TRANSPOSE)X=B. (BACK SOLVE) */
2786 /* Parameter adjustments */
2790 a_offset = 1 + a_dim1;
2795 x[i__] = b[i__] / a[i__ + i__ * a_dim1];
2804 for (j = ip1; j <= i__1; ++j) {
2805 sum += a[j + i__ * a_dim1] * x[j];
2808 x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
2815 /* ---------------------- */
2816 /* | F O R S L V | */
2817 /* ---------------------- */
2818 /* Subroutine */ int forslv_(integer *nr, integer *n, doublereal *a,
2819 doublereal *x, doublereal *b)
2821 /* System generated locals */
2822 integer a_dim1, a_offset, i__1, i__2;
2824 /* Local variables */
2825 integer i__, j, im1;
2831 /* SOLVE AX=B WHERE A IS LOWER TRIANGULAR MATRIX */
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 */
2844 /* THEN VECTORS B AND X MAY SHARE THE SAME STORAGE */
2847 /* SOLVE LX=B. (FOREWARD SOLVE) */
2849 /* Parameter adjustments */
2853 a_offset = 1 + a_dim1;
2857 x[1] = b[1] / a[a_dim1 + 1];
2862 for (i__ = 2; i__ <= i__1; ++i__) {
2866 for (j = 1; j <= i__2; ++j) {
2867 sum += a[i__ + j * a_dim1] * x[j];
2870 x[i__] = (b[i__] - sum) / a[i__ + i__ * a_dim1];
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)
2883 /* System generated locals */
2884 integer h_dim1, h_offset, i__1, i__2, i__3;
2886 /* Builtin functions */
2887 double pow_dd(doublereal *, doublereal *), sqrt(doublereal);
2889 /* Local variables */
2891 doublereal tau1, tau2;
2894 extern /* Subroutine */ int modchl_(integer *, integer *, doublereal *,
2895 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
2900 /* DRIVER FOR CHOLESKY DECOMPOSITION */
2902 /* ---------------------------------------------------------------------- */
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 */
2921 h_offset = 1 + h_dim1;
2927 /* SAVE DIAGONAL OF H */
2929 for (i__ = 1; i__ <= i__1; ++i__) {
2930 diag[i__] = h__[i__ + i__ * h_dim1];
2933 tau1 = pow_dd(eps, &c_b250);
2935 modchl_(nr, n, &h__[h_offset], &g[1], eps, &tau1, &tau2, &pivot[1], &e[1])
2939 for (i__ = 1; i__ <= i__1; ++i__) {
2940 if (pivot[i__] != i__) {
2945 if (*addmax > 0. || redo) {
2946 /* ******************************** */
2948 /* H IS NOT P.D. * */
2950 /* ******************************** */
2954 for (i__ = 2; i__ <= i__1; ++i__) {
2956 for (j = 1; j <= i__2; ++j) {
2957 h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
2963 for (i__ = 1; i__ <= i__1; ++i__) {
2965 h__[i__ + i__ * h_dim1] = diag[i__] + *addmax;
2968 /* ******************************** */
2972 /* ******************************** */
2974 for (j = 1; j <= i__1; ++j) {
2976 /* COMPUTE L(J,J) */
2980 for (i__ = 1; i__ <= i__2; ++i__) {
2981 temp += h__[j + i__ * h_dim1] * h__[j + i__ * h_dim1];
2985 h__[j + j * h_dim1] -= temp;
2986 h__[j + j * h_dim1] = sqrt(h__[j + j * h_dim1]);
2988 /* COMPUTE L(I,J) */
2990 for (i__ = j + 1; i__ <= i__2; ++i__) {
2994 for (k = 1; k <= i__3; ++k) {
2995 temp += h__[i__ + k * h_dim1] * h__[j + k * h_dim1];
2999 h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1] - temp;
3000 h__[i__ + j * h_dim1] /= h__[j + j * h_dim1];
3009 /* ---------------------- */
3010 /* | M O D C H L | */
3011 /* ---------------------- */
3012 /* ********************************************************************* */
3014 /* SUBROUTINE NAME: MODCHL */
3016 /* AUTHORS : ELIZABETH ESKOW AND ROBERT B. SCHNABEL */
3018 /* DATE : DECEMBER, 1988 */
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 */
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. */
3041 /* INPUT : NDIM - LARGEST DIMENSION OF MATRIX THAT */
3044 /* N - DIMENSION OF MATRIX A */
3046 /* A - N*N SYMMETRIC MATRIX (ONLY LOWER TRIANGULAR */
3047 /* PORTION OF A, INCLUDING THE MAIN DIAGONAL, IS USED) */
3049 /* G - N*1 WORK ARRAY */
3051 /* MCHEPS - MACHINE PRECISION */
3053 /* TAU1 - TOLERANCE USED FOR DETERMINING WHEN TO SWITCH TO */
3056 /* TAU2 - TOLERANCE USED FOR DETERMINING THE MAXIMUM */
3057 /* CONDITION NUMBER OF THE FINAL 2X2 SUBMATRIX. */
3060 /* OUTPUT : L - STORED IN THE MATRIX A (IN LOWER TRIANGULAR */
3061 /* PORTION OF A, INCLUDING THE MAIN DIAGONAL) */
3063 /* P - A RECORD OF HOW THE ROWS AND COLUMNS */
3064 /* OF THE MATRIX WERE PERMUTED WHILE */
3065 /* PERFORMING THE DECOMPOSITION */
3067 /* E - N*1 ARRAY, THE ITH ELEMENT IS THE */
3068 /* AMOUNT ADDED TO THE DIAGONAL OF A */
3069 /* AT THE ITH ITERATION */
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)
3077 /* System generated locals */
3078 integer a_dim1, a_offset, i__1, i__2, i__3;
3081 /* Builtin functions */
3082 double sqrt(doublereal);
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 *)
3097 doublereal normj, delta1;
3099 extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
3100 integer *, doublereal *);
3101 doublereal taugam, tempjj;
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. */
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 */
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. */
3120 /* Parameter adjustments */
3125 a_offset = 1 + a_dim1;
3129 init_(n, ndim, &a[a_offset], &phase1, &delta, &p[1], &g[1], &e[1], &ming,
3130 tau1, &gamma, &taugam);
3133 delta = *tau2 * (d__1 = a[a_dim1 + 1], abs(d__1)) - a[a_dim1 + 1];
3137 if (a[a_dim1 + 1] == 0.) {
3140 a[a_dim1 + 1] = sqrt(a[a_dim1 + 1] + e[1]);
3145 for (j = 1; j <= i__1; ++j) {
3151 /* FIND INDEX OF MAXIMUM DIAGONAL ELEMENT A(I,I) WHERE I>=J */
3153 maxd = a[j + j * a_dim1];
3156 for (i__ = j + 1; i__ <= i__2; ++i__) {
3157 if (maxd < a[i__ + i__ * a_dim1]) {
3158 maxd = a[i__ + i__ * a_dim1];
3164 /* PIVOT TO THE TOP THE ROW AND COLUMN WITH THE MAX DIAG */
3168 /* SWAP ROW J WITH ROW OF MAX DIAG */
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;
3178 /* SWAP COLJ AND ROW MAXDIAG BETWEEN J AND MAXDIAG */
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;
3188 /* SWAP COLUMN J WITH COLUMN OF MAX DIAG */
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;
3198 /* SWAP DIAG ELEMENTS */
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;
3204 /* SWAP ELEMENTS OF THE PERMUTATION VECTOR */
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. */
3214 tempjj = a[j + j * a_dim1];
3216 jdmin = a[jp1 + jp1 * a_dim1];
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);
3224 if (jdmin < taugam) {
3232 /* DO THE NORMAL CHOLESKY UPDATE IF STILL IN PHASE 1 */
3234 a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
3235 tempjj = a[j + j * a_dim1];
3237 for (i__ = jp1; i__ <= i__2; ++i__) {
3238 a[i__ + j * a_dim1] /= tempjj;
3242 for (i__ = jp1; i__ <= i__2; ++i__) {
3243 temp = a[i__ + j * a_dim1];
3245 for (k = jp1; k <= i__3; ++k) {
3246 a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
3252 a[*n + *n * a_dim1] = sqrt(a[*n + *n * a_dim1]);
3256 /* CALCULATE THE NEGATIVES OF THE LOWER GERSCHGORIN BOUNDS */
3258 gersch_(ndim, n, &a[a_offset], &j, &g[1]);
3267 /* FIND THE MINIMUM NEGATIVE GERSHGORIN BOUND */
3272 for (i__ = j + 1; i__ <= i__2; ++i__) {
3273 if (ming > g[i__]) {
3280 /* PIVOT TO THE TOP THE ROW AND COLUMN WITH THE */
3281 /* MINIMUM NEGATIVE GERSCHGORIN BOUND */
3285 /* SWAP ROW J WITH ROW OF MIN GERSCH BOUND */
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;
3295 /* SWAP COLJ WITH ROW IMING FROM J TO IMING */
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;
3305 /* SWAP COLUMN J WITH COLUMN OF MIN GERSCH BOUND */
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;
3315 /* SWAP DIAGONAL ELEMENTS */
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;
3321 /* SWAP ELEMENTS OF THE PERMUTATION VECTOR */
3327 /* SWAP ELEMENTS OF THE NEGATIVE GERSCHGORIN BOUNDS VECTOR */
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. */
3342 for (i__ = j + 1; i__ <= i__2; ++i__) {
3343 normj += (d__1 = a[i__ + j * a_dim1], abs(d__1));
3346 temp = max(normj,taugam);
3347 delta1 = temp - a[j + j * a_dim1];
3349 delta1 = max(temp,delta1);
3350 delta = max(delta1,delta);
3352 a[j + j * a_dim1] += e[j];
3354 /* UPDATE THE GERSCHGORIN BOUND ESTIMATES */
3355 /* (NOTE: G(I) IS THE NEGATIVE OF THE */
3356 /* GERSCHGORIN LOWER BOUND.) */
3358 if (a[j + j * a_dim1] != normj) {
3359 temp = normj / a[j + j * a_dim1] - 1.f;
3361 for (i__ = j + 1; i__ <= i__2; ++i__) {
3362 g[i__] += (d__1 = a[i__ + j * a_dim1], abs(d__1)) *
3368 /* DO THE CHOLESKY UPDATE */
3370 a[j + j * a_dim1] = sqrt(a[j + j * a_dim1]);
3371 tempjj = a[j + j * a_dim1];
3373 for (i__ = j + 1; i__ <= i__2; ++i__) {
3374 a[i__ + j * a_dim1] /= tempjj;
3378 for (i__ = j + 1; i__ <= i__2; ++i__) {
3379 temp = a[i__ + j * a_dim1];
3381 for (k = j + 1; k <= i__3; ++k) {
3382 a[i__ + k * a_dim1] -= temp * a[k + j * a_dim1];
3388 fin2x2_(ndim, n, &a[a_offset], &e[1], &j, tau2, &delta, &
3397 /* ************************************************************************ */
3398 /* SUBROUTINE NAME : INIT */
3400 /* PURPOSE : SET UP FOR START OF CHOLESKY FACTORIZATION */
3402 /* INPUT : N, NDIM, A, TAU1 */
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 */
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 *
3418 /* System generated locals */
3419 integer a_dim1, a_offset, i__1;
3420 doublereal d__1, d__2, d__3;
3422 /* Local variables */
3424 extern /* Subroutine */ int gersch_(integer *, integer *, doublereal *,
3425 integer *, doublereal *);
3427 /* Parameter adjustments */
3432 a_offset = 1 + a_dim1;
3440 for (i__ = 1; i__ <= i__1; ++i__) {
3447 /* FIND THE MAXIMUM MAGNITUDE OF THE DIAGONAL ELEMENTS. */
3448 /* IF ANY DIAGONAL ELEMENT IS NEGATIVE, THEN PHASE1 IS FALSE. */
3452 for (i__ = 1; i__ <= i__1; ++i__) {
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) {
3461 *taugam = *tau1 * *gamma;
3463 /* IF NOT IN PHASE1, THEN CALCULATE THE INITIAL GERSCHGORIN BOUNDS */
3464 /* NEEDED FOR THE START OF PHASE2. */
3467 gersch_(ndim, n, &a[a_offset], &c__1, &g[1]);
3472 /* ************************************************************************ */
3474 /* SUBROUTINE NAME : GERSCH */
3476 /* PURPOSE : CALCULATE THE NEGATIVE OF THE GERSCHGORIN BOUNDS */
3477 /* CALLED ONCE AT THE START OF PHASE II. */
3479 /* INPUT : NDIM, N, A, J */
3481 /* OUTPUT : G - AN N VECTOR CONTAINING THE NEGATIVES OF THE */
3482 /* GERSCHGORIN BOUNDS. */
3484 /* ************************************************************************ */
3485 /* Subroutine */ int gersch_(integer *ndim, integer *n, doublereal *a,
3486 integer *j, doublereal *g)
3488 /* System generated locals */
3489 integer a_dim1, a_offset, i__1, i__2;
3492 /* Local variables */
3496 /* Parameter adjustments */
3499 a_offset = 1 + a_dim1;
3504 for (i__ = *j; i__ <= i__1; ++i__) {
3507 for (k = *j; k <= i__2; ++k) {
3508 offrow += (d__1 = a[i__ + k * a_dim1], abs(d__1));
3512 for (k = i__ + 1; k <= i__2; ++k) {
3513 offrow += (d__1 = a[k + i__ * a_dim1], abs(d__1));
3516 g[i__] = offrow - a[i__ + i__ * a_dim1];
3522 /* ************************************************************************ */
3524 /* SUBROUTINE NAME : FIN2X2 */
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. */
3532 /* INPUT : NDIM, N, A, E, J, TAU2, */
3533 /* DELTA - AMOUNT ADDED TO THE DIAGONAL IN THE */
3534 /* PREVIOUS ITERATION */
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. */
3541 /* ************************************************************************ */
3542 /* Subroutine */ int fin2x2_(integer *ndim, integer *n, doublereal *a,
3543 doublereal *e, integer *j, doublereal *tau2, doublereal *delta,
3546 /* System generated locals */
3547 integer a_dim1, a_offset;
3550 /* Builtin functions */
3551 double sqrt(doublereal);
3553 /* Local variables */
3554 doublereal t1, t2, t3, t1a, t2a, temp, lmbd1, lmbd2, delta1, lmbdhi,
3558 /* FIND EIGENVALUES OF FINAL 2 BY 2 SUBMATRIX */
3560 /* Parameter adjustments */
3563 a_offset = 1 + a_dim1;
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];
3570 t2a = (d__1 = a[*n + (*n - 1) * a_dim1], abs(d__1)) * 2.;
3575 /* Computing 2nd power */
3577 t3 = t1a * sqrt(d__1 * d__1 + 1.);
3580 /* Computing 2nd power */
3582 t3 = t2a * sqrt(d__1 * d__1 + 1.);
3584 lmbd1 = (t1 - t3) / 2.f;
3585 lmbd2 = (t1 + t3) / 2.f;
3586 lmbdhi = max(lmbd1,lmbd2);
3587 lmbdlo = min(lmbd1,lmbd2);
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 */
3597 delta1 = (lmbdhi - lmbdlo) / (1.f - *tau2);
3598 delta1 = max(delta1,*gamma);
3599 delta1 = *tau2 * delta1 - lmbdlo;
3601 *delta = max(*delta,temp);
3602 *delta = max(*delta,delta1);
3604 a[*n - 1 + (*n - 1) * a_dim1] += *delta;
3605 a[*n + *n * a_dim1] += *delta;
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]);
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)
3630 /* System generated locals */
3631 integer h_dim1, h_offset, i__1, i__2;
3633 /* Builtin functions */
3634 double sqrt(doublereal);
3636 /* Local variables */
3638 doublereal r__, r1, r2, ca, cb, cc, cd, w11, sg, w22, w12, w33, w13, w23,
3640 extern /* Subroutine */ int zhz_(integer *, integer *, doublereal *,
3641 doublereal *, doublereal *, doublereal *);
3643 extern /* Subroutine */ int sigma_(doublereal *, doublereal *, doublereal
3644 *, doublereal *, doublereal *), dstar_(integer *, integer *,
3645 doublereal *, doublereal *, doublereal *, doublereal *,
3646 doublereal *, doublereal *, doublereal *, doublereal *);
3648 extern /* Subroutine */ int choldr_(integer *, integer *, doublereal *,
3649 doublereal *, doublereal *, integer *, doublereal *, doublereal *,
3650 doublereal *), bakslv_(integer *, integer *, doublereal *,
3651 doublereal *, doublereal *);
3653 extern /* Subroutine */ int forslv_(integer *, integer *, doublereal *,
3654 doublereal *, doublereal *), solvew_(integer *, integer *,
3655 doublereal *, doublereal *, doublereal *, doublereal *);
3659 /* COMPUTE TENSOR AND NEWTON STEPS */
3661 /* ---------------------------------------------------------------------------- */
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 */
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 */
3685 /* S O L V E M O D E L */
3687 /* Parameter adjustments */
3699 h_offset = 1 + h_dim1;
3705 /* COMPUTE QTHQ(N,N), ZTHZ(N-1,N-1) = FIRST N-1 ROWS AND N-1 */
3706 /* COLUMNS OF QTHQ */
3708 zhz_(nr, n, &s[1], &h__[h_offset], &u[1], &t[1]);
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];
3714 /* COLESKY DECOMPOSITION FOR FIRST N-1 ROWS AND N-1 COLUMNS OF ZTHZ */
3715 /* ZTHZ(N-1,N-1)=LLT */
3717 choldr_(nr, &i__1, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
3721 /* ON INPUT: SH IS STORED IN W2 */
3724 for (i__ = 1; i__ <= i__1; ++i__) {
3725 shs += w2[i__] * s[i__];
3730 /* COMPUTE W1,W2,W3 */
3732 /* W2=L**-1*ZT*SH */
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]);
3741 /* COMPUTE COEFFICIENTS CA,CB,CC AND CD FOR REDUCED ONE VARIABLE */
3742 /* 3RD ORDER EQUATION */
3745 for (i__ = 1; i__ <= i__1; ++i__) {
3746 w11 += w1[i__] * w1[i__];
3749 ca = *beta / 6. - w11 / 2.;
3752 for (i__ = 1; i__ <= i__1; ++i__) {
3753 w12 += w1[i__] * w2[i__];
3756 cb = *alpha / 2. - w12 * 3. / 2.;
3759 for (i__ = 1; i__ <= i__1; ++i__) {
3760 w13 += w1[i__] * w3[i__];
3765 for (i__ = 1; i__ <= i__1; ++i__) {
3766 w22 += w2[i__] * w2[i__];
3769 cc = shs - w22 - w13;
3772 for (i__ = 1; i__ <= i__1; ++i__) {
3773 sg += s[i__] * g[i__];
3778 for (i__ = 1; i__ <= i__1; ++i__) {
3779 w23 += w2[i__] * w3[i__];
3785 for (i__ = 1; i__ <= i__1; ++i__) {
3786 w33 += w3[i__] * w3[i__];
3790 /* COMPUTE DESIRABLE ROOT, SGSTAR, OF 3RD ORDER EQUATION */
3792 sigma_(&sgstar, &ca, &cb, &cc, &cd);
3798 /* 2ND ORDER ( CA=0 ) */
3800 r__ = cc * cc - cb * 4. * cd;
3805 r1 = (-cc + sqrt(r__)) / (cb * 2.);
3806 r2 = (-cc - sqrt(r__)) / (cb * 2.);
3817 if (r1 > 0. && sgstar == r2 || r2 < 0. && sgstar == r1) {
3823 /* 1ST ORDER (CA=0,CB=0) */
3833 /* FIND TENSOR STEP, W1 (FUNCTION OF SGSTAR) */
3834 dstar_(nr, n, &u[1], &s[1], &w1[1], &w2[1], &w3[1], &sgstar, &h__[
3842 for (i__ = 1; i__ <= i__1; ++i__) {
3843 uu += u[i__] * u[i__];
3844 ss += s[i__] * s[i__];
3850 choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &diag[1],
3854 /* COMPUTE LAST ROW OF L(N,N) */
3856 for (i__ = 1; i__ <= i__1; ++i__) {
3860 for (j = 1; j <= i__2; ++j) {
3861 temp += h__[*n + j * h_dim1] * h__[i__ + j * h_dim1];
3865 h__[*n + i__ * h_dim1] = (h__[i__ + *n * h_dim1] - temp) / h__[
3866 i__ + i__ * h_dim1];
3871 for (i__ = 1; i__ <= i__1; ++i__) {
3872 temp += h__[*n + i__ * h_dim1] * h__[*n + i__ * h_dim1];
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]);
3879 /* AFTER ADDING THE LAST COLUMN AND ROW */
3880 /* QTHQ IS NOT P.D., NEED TO REDO CHOLESKY DECOMPOSITION */
3882 for (i__ = 2; i__ <= i__1; ++i__) {
3884 for (j = 1; j <= i__2; ++j) {
3885 h__[i__ + j * h_dim1] = h__[j + i__ * h_dim1];
3888 h__[i__ + i__ * h_dim1] = diag[i__];
3891 h__[h_dim1 + 1] = diag[1];
3892 choldr_(nr, n, &h__[h_offset], &t[1], eps, &pivot[1], &e[1], &
3897 /* SOLVE QTHQ*QT*W3=-QT*G, WHERE W3 IS NEWTON STEP */
3900 for (i__ = 1; i__ <= i__1; ++i__) {
3903 for (j = 1; j <= i__2; ++j) {
3904 w2[i__] += u[j] * u[i__] * g[j] / uu;
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 */
3914 for (i__ = 1; i__ <= i__1; ++i__) {
3917 for (j = 1; j <= i__2; ++j) {
3918 w3[i__] += u[i__] * u[j] * w2[j] / uu;
3921 w3[i__] = w2[i__] - w3[i__];
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)
3936 /* Format strings */
3937 static char fmt_900[] = "(\002 OPTSTP STEP OF MAXIMUM LENGTH (STEPMX)"
3939 static char fmt_901[] = "(\002 OPTSTP RELATIVE GRADIENT CLOSE TO ZE"
3940 "RO.\002/\002 OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION.\002)"
3942 static char fmt_902[] = "(\002 OPTSTP SUCCESSIVE ITERATES WITHIN TOLE"
3943 "RANCE.\002/\002 OPTSTP CURRENT ITERATE IS PROBABLY SOLUTION"
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)";
3958 /* System generated locals */
3960 doublereal d__1, d__2, d__3;
3962 /* Builtin functions */
3963 integer s_wsfe(cilist *), e_wsfe(void);
3965 /* Local variables */
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 };
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 */
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 */
4017 /* Parameter adjustments */
4028 /* LAST GLOBAL STEP FAILED TO LOCATE A POINT LOWER THAN X */
4032 /* IF(IRETCD.EQ.1) */
4039 /* FIND DIRECTION IN WHICH RELATIVE GRADIENT MAXIMUM. */
4040 /* CHECK WHETHER WITHIN TOLERANCE */
4044 d__ = max(d__1,*fscale);
4048 for (i__ = 1; i__ <= i__1; ++i__) {
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);
4056 if (*rgx <= *gradtl) {
4064 /* FIND DIRECTION IN WHICH RELATIVE STEPSIZE MAXIMUM */
4065 /* CHECK WHETHER WITHIN TOLERANCE. */
4069 for (i__ = 1; i__ <= i__1; ++i__) {
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);
4077 if (*rsx <= *steptl) {
4081 /* CHECK ITERATION LIMIT */
4084 if (*itncnt >= *itnlim) {
4088 /* CHECK NUMBER OF CONSECUTIVE STEPS \ STEPMX */
4093 /* IF(.NOT.MXTAKE) */
4099 /* IF (MOD(MSG/8,2) .EQ. 0) WRITE(IPR,900) */
4101 io___280.ciunit = *ipr;
4113 /* PRINT TERMINATION CODE */
4117 /* IF (MOD(MSG/8,2) .EQ. 0) GO TO(601,602,603,604,605), ITRMCD */
4129 io___281.ciunit = *ipr;
4134 io___282.ciunit = *ipr;
4139 io___283.ciunit = *ipr;
4144 io___284.ciunit = *ipr;
4149 io___285.ciunit = *ipr;
4160 /* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx,
4161 doublereal *dy, integer *incy)
4163 /* System generated locals */
4166 /* Local variables */
4167 integer i__, m, ix, iy, mp1;
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. */
4175 /* Parameter adjustments */
4183 if (*incx == 1 && *incy == 1) {
4187 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
4188 /* NOT EQUAL TO 1 */
4193 ix = (-(*n) + 1) * *incx + 1;
4196 iy = (-(*n) + 1) * *incy + 1;
4199 for (i__ = 1; i__ <= i__1; ++i__) {
4207 /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
4218 for (i__ = 1; i__ <= i__1; ++i__) {
4228 for (i__ = mp1; i__ <= i__1; i__ += 7) {
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];
4242 doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy,
4245 /* System generated locals */
4249 /* Local variables */
4250 integer i__, m, ix, iy, mp1;
4254 /* FORMS THE DOT PRODUCT OF TWO VECTORS. */
4255 /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
4256 /* JACK DONGARRA, LINPACK, 3/11/78. */
4259 /* Parameter adjustments */
4269 if (*incx == 1 && *incy == 1) {
4273 /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
4274 /* NOT EQUAL TO 1 */
4279 ix = (-(*n) + 1) * *incx + 1;
4282 iy = (-(*n) + 1) * *incy + 1;
4285 for (i__ = 1; i__ <= i__1; ++i__) {
4286 dtemp += dx[ix] * dy[iy];
4294 /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
4305 for (i__ = 1; i__ <= i__1; ++i__) {
4306 dtemp += dx[i__] * dy[i__];
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__ +
4327 /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx,
4330 /* System generated locals */
4333 /* Local variables */
4334 integer i__, m, mp1, nincx;
4337 /* SCALES A VECTOR BY A CONSTANT. */
4338 /* USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. */
4339 /* JACK DONGARRA, LINPACK, 3/11/78. */
4342 /* Parameter adjustments */
4353 /* CODE FOR INCREMENT NOT EQUAL TO 1 */
4358 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
4359 dx[i__] = *da * dx[i__];
4364 /* CODE FOR INCREMENT EQUAL TO 1 */
4375 for (i__ = 1; i__ <= i__2; ++i__) {
4376 dx[i__] = *da * dx[i__];
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];
4396 /* Main program alias */ int driver_ () { MAIN__ (); return 0; }