/* * ME -- * * $Id: lbfgsb.c,v 1.1.1.1 2004/02/16 23:45:44 taku-ku Exp $; * * Copyright (C) 2001-2002 Taku Kudo * All rights reserved. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later verjsion. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * */ /* NEOS, November 1994. (Latest revision April 1997.) Optimization Technology Center. Argonne National Laboratory and Northwestern University. Written by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */ #include #include typedef long int integer; typedef unsigned long int uinteger; typedef char *address; typedef short int shortint; typedef float real; typedef double doublereal; typedef long int logical; typedef short int shortlogical; typedef char logical1; typedef char integer1; typedef short flag; typedef short ftnlen; typedef short ftnint; #define TRUE_ (1) #define FALSE_ (0) #define abs(x) ((x) >= 0 ? (x) : -(x)) #define min(a,b) ((a) <= (b) ? (a) : (b)) #define max(a,b) ((a) >= (b) ? (a) : (b)) #define s_cmp(a,b,c,d) (strncmp ((a),(b), (c) >= (d) ? (d) : (c))) #define s_copy(a,b,c,d) (strcpy((a),(b))) static doublereal c_b9 = 0.; static integer c__1 = 1; static integer c__11 = 11; static doublereal c_b275 = .001; static doublereal c_b276 = .9; static doublereal c_b277 = .1; static /* Subroutine */ int mainlb_(); static /* Subroutine */ doublereal ddot_(); static /* Subroutine */ int dscal_(); static /* Subroutine */ int freev_(); static /* Subroutine */ int dcopy_(); static /* Subroutine */ int formk_(); static /* Subroutine */ int formt_(); static /* Subroutine */ int subsm_(); static /* Subroutine */ int prn2lb_(); static /* Subroutine */ int errclb_(); static /* Subroutine */ int active_(); static /* Subroutine */ int cauchy_(); static /* Subroutine */ doublereal dpmeps_(); static /* Subroutine */ int cmprlb_(); static /* Subroutine */ int lnsrlb_(); static /* Subroutine */ int matupd_(); static /* Subroutine */ int projgr_(); static /* Subroutine */ int dcopy_(); static /* Subroutine */ int daxpy_(); static /* Subroutine */ int hpsolb_(); static /* Subroutine */ int bmv_(); static /* Subroutine */ int dpofa_(); static /* Subroutine */ int dtrsl_(); static /* Subroutine */ int dcsrch_(); static /* Subroutine */ int dcopy_(); static /* Subroutine */ int dcstep_(); /* ================ L-BFGS-B (version 2.3) ========================== */ /* Subroutine */ int lbfgsb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave, task_len, csave_len) integer *n, *m; doublereal *x, *l, *u; integer *nbd; doublereal *f, *g, *factr, *pgtol, *wa; integer *iwa; char *task; integer *iprint; char *csave; logical *lsave; integer *isave; doublereal *dsave; ftnlen task_len; ftnlen csave_len; { /* System generated locals */ integer i__1; /* Local variables */ static integer lsnd, l1, l2, l3, ld, lr, lt; static integer lz, lwa, lwn, lss, lws, lwt, lsy, lwy; /* Parameter adjustments */ --iwa; --g; --nbd; --u; --l; --x; --wa; --lsave; --isave; --dsave; /* Function Body */ if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) { isave[1] = *m * *n; /* Computing 2nd power */ i__1 = *m; isave[2] = i__1 * i__1; /* Computing 2nd power */ i__1 = *m; isave[3] = i__1 * i__1 << 2; isave[4] = 1; isave[5] = isave[4] + isave[1]; isave[6] = isave[5] + isave[1]; isave[7] = isave[6] + isave[2]; isave[8] = isave[7] + isave[2]; isave[9] = isave[8]; isave[10] = isave[9] + isave[2]; isave[11] = isave[10] + isave[3]; isave[12] = isave[11] + isave[3]; isave[13] = isave[12] + *n; isave[14] = isave[13] + *n; isave[15] = isave[14] + *n; isave[16] = isave[15] + *n; } l1 = isave[1]; l2 = isave[2]; l3 = isave[3]; lws = isave[4]; lwy = isave[5]; lsy = isave[6]; lss = isave[7]; lwt = isave[9]; lwn = isave[10]; lsnd = isave[11]; lz = isave[12]; lr = isave[13]; ld = isave[14]; lt = isave[15]; lwa = isave[16]; mainlb_(n, m, &x[1], &l[1], &u[1], &nbd[1], f, &g[1], factr, pgtol, &wa[ lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lwt], &wa[lwn], &wa[lsnd], &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lwa], &iwa[1], &iwa[*n + 1], &iwa[(*n << 1) + 1], task, iprint, csave, &lsave[1], &isave[ 22], &dsave[1], (ftnlen)60, (ftnlen)60); return 0; } /* setulb_ */ /* ======================= The end of setulb ============================= */ /* Subroutine */ int mainlb_(n, m, x, l, u, nbd, f, g, factr, pgtol, ws, wy, sy, ss, wt, wn, snd, z__, r__, d__, t, wa, index, iwhere, indx2, task, iprint, csave, lsave, isave, dsave, task_len, csave_len) integer *n, *m; doublereal *x, *l, *u; integer *nbd; doublereal *f, *g, *factr, *pgtol, *ws, *wy, *sy, *ss, *wt, *wn, *snd, *z__, * r__, *d__, *t, *wa; integer *index, *iwhere, *indx2; char *task; integer *iprint; char *csave; logical *lsave; integer *isave; doublereal *dsave; ftnlen task_len; ftnlen csave_len; { /* System generated locals */ integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, wt_dim1, wt_offset, wn_dim1, wn_offset, snd_dim1, snd_offset, i__1; doublereal d__1, d__2; /* Local variables */ static integer head; static doublereal fold; static integer nact; static doublereal ddum; static integer info; static integer nfgv, ifun, iter, nint; static char word[3]; static integer i__, iback, k; static doublereal gdold; static integer nfree; static logical boxed; static integer itail; static doublereal theta; static doublereal dnorm; static integer nskip, iword; static doublereal xstep, stpmx; static doublereal gd, dr, rr; static integer ileave; static integer itfile; static doublereal cachyt, epsmch; static logical updatd; static logical prjctd; static integer iupdat; static logical cnstnd; static doublereal sbgnrm; static integer nenter; static doublereal lnscht; static integer nintol; static doublereal dtd; static integer col; static doublereal tol; static logical wrk; static doublereal stp, cpu1, cpu2; /* Parameter adjustments */ --indx2; --iwhere; --index; --t; --d__; --r__; --z__; --g; --nbd; --u; --l; --x; --wa; snd_dim1 = 2 * *m; snd_offset = 1 + snd_dim1 * 1; snd -= snd_offset; wn_dim1 = 2 * *m; wn_offset = 1 + wn_dim1 * 1; wn -= wn_offset; wt_dim1 = *m; wt_offset = 1 + wt_dim1 * 1; wt -= wt_offset; ss_dim1 = *m; ss_offset = 1 + ss_dim1 * 1; ss -= ss_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; --lsave; --isave; --dsave; /* Function Body */ if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) { /* Generate the current machine precision. */ epsmch = dpmeps_(); fold = 0.; dnorm = 0.; cpu1 = 0.; gd = 0.; sbgnrm = 0.; stp = 0.; stpmx = 0.; gdold = 0.; dtd = 0.; /* Initialize counters and scalars when task='START'. */ /* for the limited memory BFGS matrices: */ col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; iback = 0; itail = 0; ifun = 0; iword = 0; nact = 0; ileave = 0; nenter = 0; /* for operation counts: */ iter = 0; nfgv = 0; nint = 0; nintol = 0; nskip = 0; nfree = *n; /* for stopping tolerance: */ tol = *factr * epsmch; cachyt = 0.; lnscht = 0.; /* 'word' records the status of subspace solutions. */ s_copy(word, "---", (ftnlen)3, (ftnlen)3); /* 'info' records the termination information. */ info = 0; itfile = 0; /* Check the input arguments for errors. */ errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (ftnlen) 60); if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) { return 0; } /* Initialize iwhere & project x onto the feasible set. */ active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd, &cnstnd, &boxed); /* The end of the initialization. */ } else { /* restore local variables. */ prjctd = lsave[1]; cnstnd = lsave[2]; boxed = lsave[3]; updatd = lsave[4]; nintol = isave[1]; itfile = isave[3]; iback = isave[4]; nskip = isave[5]; head = isave[6]; col = isave[7]; itail = isave[8]; iter = isave[9]; iupdat = isave[10]; nint = isave[12]; nfgv = isave[13]; info = isave[14]; ifun = isave[15]; iword = isave[16]; nfree = isave[17]; nact = isave[18]; ileave = isave[19]; nenter = isave[20]; theta = dsave[1]; fold = dsave[2]; tol = dsave[3]; dnorm = dsave[4]; epsmch = dsave[5]; cpu1 = dsave[6]; cachyt = dsave[7]; lnscht = dsave[9]; gd = dsave[11]; stpmx = dsave[12]; sbgnrm = dsave[13]; stp = dsave[14]; gdold = dsave[15]; dtd = dsave[16]; /* After returning from the driver go to the point where execution */ /* is to resume. */ if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) { goto L666; } if (s_cmp(task, "NEW_X", (ftnlen)5, (ftnlen)5) == 0) { goto L777; } if (s_cmp(task, "FG_ST", (ftnlen)5, (ftnlen)5) == 0) { goto L111; } if (s_cmp(task, "STOP", (ftnlen)4, (ftnlen)4) == 0) { if (s_cmp(task + 6, "CPU", (ftnlen)3, (ftnlen)3) == 0) { /* restore the previous iterate. */ dcopy_(n, &t[1], &c__1, &x[1], &c__1); dcopy_(n, &r__[1], &c__1, &g[1], &c__1); *f = fold; } goto L999; } } /* Compute f0 and g0. */ s_copy(task, "FG_START", (ftnlen)60, (ftnlen)8); /* return to the driver to calculate f and g; reenter at 111. */ goto L1000; L111: nfgv = 1; /* Compute the infinity norm of the (-) projected gradient. */ projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm); if (sbgnrm <= *pgtol) { /* terminate the algorithm. */ s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", ( ftnlen)60, (ftnlen)48); goto L999; } /* ----------------- the beginning of the loop -------------------------- */ L222: iword = -1; if (! cnstnd && col > 0) { /* skip the search for GCP. */ dcopy_(n, &x[1], &c__1, &z__[1], &c__1); wrk = updatd; nint = 0; goto L333; } /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Compute the Generalized Cauchy Point (GCP). */ /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[ 1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[ sy_offset], &wt[wt_offset], &theta, &col, &head, &wa[1], &wa[(*m << 1) + 1], &wa[(*m << 2) + 1], &wa[*m * 6 + 1], &nint, iprint, & sbgnrm, &info, &epsmch); if (info != 0) { /* singular triangular system detected; refresh the lbfgs memory. */ info = 0; col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; cachyt = cachyt + cpu2 - cpu1; goto L222; } cachyt = cachyt + cpu2 - cpu1; nintol += nint; /* Count the entering and leaving variables for iter > 0; */ /* find the index set of free and active variables at the GCP. */ freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], & wrk, &updatd, &cnstnd, iprint, &iter); nact = *n - nfree; L333: /* If there are no free variables or B=theta*I, then */ /* skip the subspace minimization. */ if (nfree == 0 || col == 0) { goto L555; } /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Subspace minimization. */ /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Form the LEL^T factorization of the indefinite */ /* matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */ /* [L_a -R_z theta*S'AA'S ] */ /* where E = [-I 0] */ /* [ 0 I] */ if (wrk) { formk_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iupdat, & updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], & wy[wy_offset], &sy[sy_offset], &theta, &col, &head, &info); } if (info != 0) { /* nonpositive definiteness in Cholesky factorization; */ /* refresh the lbfgs memory and restart the iteration. */ info = 0; col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; goto L222; } /* compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */ /* from 'cauchy'). */ cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset] , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], &theta, & col, &head, &nfree, &cnstnd, &info); if (info != 0) { goto L444; } /* call the direct method. */ subsm_(n, m, &nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], & ws[ws_offset], &wy[wy_offset], &theta, &col, &head, &iword, &wa[1] , &wn[wn_offset], iprint, &info); L444: if (info != 0) { /* singular triangular system detected; */ /* refresh the lbfgs memory and restart the iteration. */ info = 0; col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; goto L222; } L555: /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Line search and optimality tests. */ /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Generate the search direction d:=z-x. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { d__[i__] = z__[i__] - x[i__]; /* L40: */ } L666: lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], & d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, & stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd, csave, &isave[22], &dsave[17], (ftnlen)60, (ftnlen)60); if (info != 0 || iback >= 20) { /* restore the previous iterate. */ dcopy_(n, &t[1], &c__1, &x[1], &c__1); dcopy_(n, &r__[1], &c__1, &g[1], &c__1); *f = fold; if (col == 0) { /* abnormal termination. */ if (info == 0) { info = -9; /* restore the actual number of f and g evaluations etc. */ --nfgv; --ifun; --iback; } s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (ftnlen)60, ( ftnlen)30); ++iter; goto L999; } else { /* refresh the lbfgs memory and restart the iteration. */ if (info == 0) { --nfgv; } info = 0; col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; s_copy(task, "RESTART_FROM_LNSRCH", (ftnlen)60, (ftnlen)19); lnscht = lnscht + cpu2 - cpu1; goto L222; } } else if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) { /* return to the driver for calculating f and g; reenter at 666. */ goto L1000; } else { /* calculate and print out the quantities related to the new X. */ lnscht = lnscht + cpu2 - cpu1; ++iter; /* Compute the infinity norm of the projected (-)gradient. */ projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm); /* Print iteration information. */ prn2lb_(n, &x[1], f, &g[1], iprint, &itfile, &iter, &nfgv, &nact, & sbgnrm, &nint, word, &iword, &iback, &stp, &xstep, (ftnlen)3); goto L1000; } L777: /* Test for termination. */ if (sbgnrm <= *pgtol) { /* terminate the algorithm. */ s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", ( ftnlen)60, (ftnlen)48); goto L999; } /* Computing MAX */ d__1 = abs(fold), d__2 = abs(*f), d__1 = max(d__1,d__2); ddum = max(d__1,1.); if (fold - *f <= tol * ddum) { /* terminate the algorithm. */ s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", ( ftnlen)60, (ftnlen)47); if (iback >= 10) { info = -5; } /* i.e., to issue a warning if iback>10 in the line search. */ goto L999; } /* Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { r__[i__] = g[i__] - r__[i__]; /* L42: */ } rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1); if (stp == 1.) { dr = gd - gdold; ddum = -gdold; } else { dr = (gd - gdold) * stp; dscal_(n, &stp, &d__[1], &c__1); ddum = -gdold * stp; } if (dr <= epsmch * ddum) { /* skip the L-BFGS update. */ ++nskip; updatd = FALSE_; goto L888; } /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ /* Update the L-BFGS matrix. */ /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */ updatd = TRUE_; ++iupdat; /* Update matrices WS and WY and form the middle matrix in B. */ matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[ ss_offset], &d__[1], &r__[1], &itail, &iupdat, &col, &head, & theta, &rr, &dr, &stp, &dtd); /* Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */ /* Store T in the upper triangular of the array wt; */ /* Cholesky factorize T to J*J' with */ /* J' stored in the upper triangular of wt. */ formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], &col, &theta, & info); if (info != 0) { /* nonpositive definiteness in Cholesky factorization; */ /* refresh the lbfgs memory and restart the iteration. */ info = 0; col = 0; head = 1; theta = 1.; iupdat = 0; updatd = FALSE_; goto L222; } /* Now the inverse of the middle matrix in B is */ /* [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] */ /* [ -L*D^(-1/2) J ] [ 0 J' ] */ L888: /* -------------------- the end of the loop ----------------------------- */ goto L222; L999: L1000: /* Save local variables. */ lsave[1] = prjctd; lsave[2] = cnstnd; lsave[3] = boxed; lsave[4] = updatd; isave[1] = nintol; isave[3] = itfile; isave[4] = iback; isave[5] = nskip; isave[6] = head; isave[7] = col; isave[8] = itail; isave[9] = iter; isave[10] = iupdat; isave[12] = nint; isave[13] = nfgv; isave[14] = info; isave[15] = ifun; isave[16] = iword; isave[17] = nfree; isave[18] = nact; isave[19] = ileave; isave[20] = nenter; dsave[1] = theta; dsave[2] = fold; dsave[3] = tol; dsave[4] = dnorm; dsave[5] = epsmch; dsave[6] = cpu1; dsave[7] = cachyt; dsave[9] = lnscht; dsave[11] = gd; dsave[12] = stpmx; dsave[13] = sbgnrm; dsave[14] = stp; dsave[15] = gdold; dsave[16] = dtd; return 0; } /* mainlb_ */ /* ======================= The end of mainlb ============================= */ /* Subroutine */ int active_(n, l, u, nbd, x, iwhere, iprint, prjctd, cnstnd, boxed) integer *n; doublereal *l, *u; integer *nbd; doublereal *x; integer *iwhere, *iprint; logical *prjctd, *cnstnd, *boxed; { /* System generated locals */ integer i__1; /* Local variables */ static integer nbdd, i__; /* Initialize nbdd, prjctd, cnstnd and boxed. */ /* Parameter adjustments */ --iwhere; --x; --nbd; --u; --l; /* Function Body */ nbdd = 0; *prjctd = FALSE_; *cnstnd = FALSE_; *boxed = TRUE_; /* Project the initial x to the easible set if necessary. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (nbd[i__] > 0) { if (nbd[i__] <= 2 && x[i__] <= l[i__]) { if (x[i__] < l[i__]) { *prjctd = TRUE_; x[i__] = l[i__]; } ++nbdd; } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) { if (x[i__] > u[i__]) { *prjctd = TRUE_; x[i__] = u[i__]; } ++nbdd; } } /* L10: */ } /* Initialize iwhere and assign values to cnstnd and boxed. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (nbd[i__] != 2) { *boxed = FALSE_; } if (nbd[i__] == 0) { /* this variable is always free */ iwhere[i__] = -1; /* otherwise set x(i)=mid(x(i), u(i), l(i)). */ } else { *cnstnd = TRUE_; if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) { /* this variable is always fixed */ iwhere[i__] = 3; } else { iwhere[i__] = 0; } } /* L20: */ } return 0; } /* active_ */ /* ======================= The end of active ============================= */ /* Subroutine */ int bmv_(m, sy, wt, col, v, p, info) integer *m; doublereal *sy, *wt; integer *col; doublereal *v, *p; integer *info; { /* System generated locals */ integer sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2; /* Builtin functions */ double sqrt(); /* Local variables */ static integer i__, k; static integer i2; static doublereal sum; /* Parameter adjustments */ wt_dim1 = *m; wt_offset = 1 + wt_dim1 * 1; wt -= wt_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; --p; --v; /* Function Body */ if (*col == 0) { return 0; } /* PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] */ /* [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. */ /* solve Jp2=v2+LD^(-1)v1. */ p[*col + 1] = v[*col + 1]; i__1 = *col; for (i__ = 2; i__ <= i__1; ++i__) { i2 = *col + i__; sum = 0.; i__2 = i__ - 1; for (k = 1; k <= i__2; ++k) { sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1]; /* L10: */ } p[i2] = v[i2] + sum; /* L20: */ } /* Solve the triangular system */ dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info); if (*info != 0) { return 0; } /* solve D^(1/2)p1=v1. */ i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]); /* L30: */ } /* PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] */ /* [ 0 J' ] [ p2 ] [ p2 ]. */ /* solve J^Tp2=p2. */ dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info); if (*info != 0) { return 0; } /* compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */ /* =-D^(-1/2)p1+D^(-1)L'p2. */ i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { p[i__] = -p[i__] / sqrt(sy[i__ + i__ * sy_dim1]); /* L40: */ } i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; i__2 = *col; for (k = i__ + 1; k <= i__2; ++k) { sum += sy[k + i__ * sy_dim1] * p[*col + k] / sy[i__ + i__ * sy_dim1]; /* L50: */ } p[i__] += sum; /* L60: */ } return 0; } /* bmv_ */ /* ======================== The end of bmv =============================== */ /* Subroutine */ int cauchy_(n, x, l, u, nbd, g, iorder, iwhere, t, d__, xcp, m, wy, ws, sy, wt, theta, col, head, p, c__, wbp, v, nint, iprint, sbgnrm, info, epsmch) integer *n; doublereal *x, *l, *u; integer *nbd; doublereal *g; integer *iorder, *iwhere; doublereal *t, *d__, *xcp; integer *m; doublereal *wy, *ws, *sy, *wt, *theta; integer *col, *head; doublereal *p, *c__, *wbp, *v; integer *nint, *iprint; doublereal *sbgnrm; integer *info; doublereal *epsmch; { /* System generated locals */ integer wy_dim1, wy_offset, ws_dim1, ws_offset, sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2; doublereal d__1; /* Local variables */ static doublereal dibp; static integer iter; static doublereal zibp, tsum, dibp2; static integer i__, j; static logical bnded; static doublereal neggi; static integer nfree; static doublereal bkmin; static integer nleft; static doublereal f1, f2, f2_org__, dt, tj, tl; static integer nbreak, ibkmin; static doublereal tu; static integer pointr; static doublereal tj0; static logical xlower, xupper; static integer ibp; static doublereal dtm; static doublereal wmc, wmp, wmw; static integer col2; /* Parameter adjustments */ --xcp; --d__; --t; --iwhere; --iorder; --g; --nbd; --u; --l; --x; --v; --wbp; --c__; --p; wt_dim1 = *m; wt_offset = 1 + wt_dim1 * 1; wt -= wt_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; /* Function Body */ if (*sbgnrm <= 0.) { dcopy_(n, &x[1], &c__1, &xcp[1], &c__1); return 0; } bnded = TRUE_; nfree = *n + 1; nbreak = 0; ibkmin = 0; bkmin = 0.; col2 = *col << 1; f1 = 0.; /* We set p to zero and build it up as we determine d. */ i__1 = col2; for (i__ = 1; i__ <= i__1; ++i__) { p[i__] = 0.; /* L20: */ } /* In the following loop we determine for each variable its bound */ /* status and its breakpoint, and update p accordingly. */ /* Smallest breakpoint is identified. */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { neggi = -g[i__]; if (iwhere[i__] != 3 && iwhere[i__] != -1) { /* if x(i) is not a constant and has bounds, */ /* compute the difference between x(i) and its bounds. */ if (nbd[i__] <= 2) { tl = x[i__] - l[i__]; } if (nbd[i__] >= 2) { tu = u[i__] - x[i__]; } /* If a variable is close enough to a bound */ /* we treat it as at bound. */ xlower = nbd[i__] <= 2 && tl <= 0.; xupper = nbd[i__] >= 2 && tu <= 0.; /* reset iwhere(i). */ iwhere[i__] = 0; if (xlower) { if (neggi <= 0.) { iwhere[i__] = 1; } } else if (xupper) { if (neggi >= 0.) { iwhere[i__] = 2; } } else { if (abs(neggi) <= 0.) { iwhere[i__] = -3; } } } pointr = *head; if (iwhere[i__] != 0 && iwhere[i__] != -1) { d__[i__] = 0.; } else { d__[i__] = neggi; f1 -= neggi * neggi; /* calculate p := p - W'e_i* (g_i). */ i__2 = *col; for (j = 1; j <= i__2; ++j) { p[j] += wy[i__ + pointr * wy_dim1] * neggi; p[*col + j] += ws[i__ + pointr * ws_dim1] * neggi; pointr = pointr % *m + 1; /* L40: */ } if (nbd[i__] <= 2 && nbd[i__] != 0 && neggi < 0.) { /* x(i) + d(i) is bounded; compute t(i). */ ++nbreak; iorder[nbreak] = i__; t[nbreak] = tl / (-neggi); if (nbreak == 1 || t[nbreak] < bkmin) { bkmin = t[nbreak]; ibkmin = nbreak; } } else if (nbd[i__] >= 2 && neggi > 0.) { /* x(i) + d(i) is bounded; compute t(i). */ ++nbreak; iorder[nbreak] = i__; t[nbreak] = tu / neggi; if (nbreak == 1 || t[nbreak] < bkmin) { bkmin = t[nbreak]; ibkmin = nbreak; } } else { /* x(i) + d(i) is not bounded. */ --nfree; iorder[nfree] = i__; if (abs(neggi) > 0.) { bnded = FALSE_; } } } /* L50: */ } /* The indices of the nonzero components of d are now stored */ /* in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */ /* The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */ if (*theta != 1.) { /* complete the initialization of p for theta not= one. */ dscal_(col, theta, &p[*col + 1], &c__1); } /* Initialize GCP xcp = x. */ dcopy_(n, &x[1], &c__1, &xcp[1], &c__1); if (nbreak == 0 && nfree == *n + 1) { /* is a zero vector, return with the initial xcp as GCP. */ return 0; } /* Initialize c = W'(xcp - x) = 0. */ i__1 = col2; for (j = 1; j <= i__1; ++j) { c__[j] = 0.; /* L60: */ } /* Initialize derivative f2. */ f2 = -(*theta) * f1; f2_org__ = f2; if (*col > 0) { bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info); if (*info != 0) { return 0; } f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1); } dtm = -f1 / f2; tsum = 0.; *nint = 1; /* If there are no breakpoints, locate the GCP and return. */ if (nbreak == 0) { goto L888; } nleft = nbreak; iter = 1; tj = 0.; /* ------------------- the beginning of the loop ------------------------- */ L777: /* Find the next smallest breakpoint; */ /* compute dt = t(nleft) - t(nleft + 1). */ tj0 = tj; if (iter == 1) { /* Since we already have the smallest breakpoint we need not do */ /* heapsort yet. Often only one breakpoint is used and the */ /* cost of heapsort is avoided. */ tj = bkmin; ibp = iorder[ibkmin]; } else { if (iter == 2) { /* Replace the already used smallest breakpoint with the */ /* breakpoint numbered nbreak > nlast, before heapsort call. */ if (ibkmin != nbreak) { t[ibkmin] = t[nbreak]; iorder[ibkmin] = iorder[nbreak]; } /* Update heap structure of breakpoints */ /* (if iter=2, initialize heap). */ } i__1 = iter - 2; hpsolb_(&nleft, &t[1], &iorder[1], &i__1); tj = t[nleft]; ibp = iorder[nleft]; } dt = tj - tj0; /* If a minimizer is within this interval, */ /* locate the GCP and return. */ if (dtm < dt) { goto L888; } /* Otherwise fix one variable and */ /* reset the corresponding component of d to zero. */ tsum += dt; --nleft; ++iter; dibp = d__[ibp]; d__[ibp] = 0.; if (dibp > 0.) { zibp = u[ibp] - x[ibp]; xcp[ibp] = u[ibp]; iwhere[ibp] = 2; } else { zibp = l[ibp] - x[ibp]; xcp[ibp] = l[ibp]; iwhere[ibp] = 1; } if (nleft == 0 && nbreak == *n) { /* all n variables are fixed, */ /* return with xcp as GCP. */ dtm = dt; goto L999; } /* Update the derivative information. */ ++(*nint); /* Computing 2nd power */ d__1 = dibp; dibp2 = d__1 * d__1; /* Update f1 and f2. */ /* temporarily set f1 and f2 for col=0. */ f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp; f2 -= *theta * dibp2; if (*col > 0) { /* update c = c + dt*p. */ daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1); /* choose wbp, */ /* the row of W corresponding to the breakpoint encountered. */ pointr = *head; i__1 = *col; for (j = 1; j <= i__1; ++j) { wbp[j] = wy[ibp + pointr * wy_dim1]; wbp[*col + j] = *theta * ws[ibp + pointr * ws_dim1]; pointr = pointr % *m + 1; /* L70: */ } /* compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */ bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info); if (*info != 0) { return 0; } wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1); wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1); wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1); /* update p = p - dibp*wbp. */ d__1 = -dibp; daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1); /* complete updating f1 and f2 while col > 0. */ f1 += dibp * wmc; f2 = f2 + dibp * 2. * wmp - dibp2 * wmw; } /* Computing MAX */ d__1 = *epsmch * f2_org__; f2 = max(d__1,f2); if (nleft > 0) { dtm = -f1 / f2; goto L777; /* to repeat the loop for unsearched intervals. */ } else if (bnded) { f1 = 0.; f2 = 0.; dtm = 0.; } else { dtm = -f1 / f2; } /* ------------------- the end of the loop ------------------------------- */ L888: if (dtm <= 0.) { dtm = 0.; } tsum += dtm; /* Move free variables (i.e., the ones w/o breakpoints) and */ /* the variables whose breakpoints haven't been reached. */ daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1); L999: /* Update c = c + dtm*p = W'(x^c - x) */ /* which will be used in computing r = Z'(B(x^c - x) + g). */ if (*col > 0) { daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1); } return 0; } /* cauchy_ */ /* ====================== The end of cauchy ============================== */ /* Subroutine */ int cmprlb_(n, m, x, g, ws, wy, sy, wt, z__, r__, wa, index, theta, col, head, nfree, cnstnd, info) integer *n, *m; doublereal *x, *g, *ws, *wy, *sy, *wt, *z__, *r__, *wa; integer *index; doublereal *theta; integer *col, *head, *nfree; logical *cnstnd; integer *info; { /* System generated locals */ integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2; /* Local variables */ static integer i__, j, k; static doublereal a1, a2; static integer pointr; /* Parameter adjustments */ --index; --r__; --z__; --g; --x; --wa; wt_dim1 = *m; wt_offset = 1 + wt_dim1 * 1; wt -= wt_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; /* Function Body */ if (! (*cnstnd) && *col > 0) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { r__[i__] = -g[i__]; /* L26: */ } } else { i__1 = *nfree; for (i__ = 1; i__ <= i__1; ++i__) { k = index[i__]; r__[i__] = -(*theta) * (z__[k] - x[k]) - g[k]; /* L30: */ } bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wa[(*m << 1) + 1], &wa[ 1], info); if (*info != 0) { *info = -8; return 0; } pointr = *head; i__1 = *col; for (j = 1; j <= i__1; ++j) { a1 = wa[j]; a2 = *theta * wa[*col + j]; i__2 = *nfree; for (i__ = 1; i__ <= i__2; ++i__) { k = index[i__]; r__[i__] = r__[i__] + wy[k + pointr * wy_dim1] * a1 + ws[k + pointr * ws_dim1] * a2; /* L32: */ } pointr = pointr % *m + 1; /* L34: */ } } return 0; } /* cmprlb_ */ /* ======================= The end of cmprlb ============================= */ /* Subroutine */ int errclb_(n, m, factr, l, u, nbd, task, info, k, task_len) integer *n, *m; doublereal *factr, *l, *u; integer *nbd; char *task; integer *info, *k; ftnlen task_len; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__; /* Parameter adjustments */ --nbd; --u; --l; /* Function Body */ if (*n <= 0) { s_copy(task, "ERROR: N .LE. 0", (ftnlen)60, (ftnlen)15); } if (*m <= 0) { s_copy(task, "ERROR: M .LE. 0", (ftnlen)60, (ftnlen)15); } if (*factr < 0.) { s_copy(task, "ERROR: FACTR .LT. 0", (ftnlen)60, (ftnlen)19); } /* Check the validity of the arrays nbd(i), u(i), and l(i). */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (nbd[i__] < 0 || nbd[i__] > 3) { /* return */ s_copy(task, "ERROR: INVALID NBD", (ftnlen)60, (ftnlen)18); *info = -6; *k = i__; } if (nbd[i__] == 2) { if (l[i__] > u[i__]) { /* return */ s_copy(task, "ERROR: NO FEASIBLE SOLUTION", (ftnlen)60, ( ftnlen)27); *info = -7; *k = i__; } } /* L10: */ } return 0; } /* errclb_ */ /* ======================= The end of errclb ============================= */ /* Subroutine */ int formk_(n, nsub, ind, nenter, ileave, indx2, iupdat, updatd, wn, wn1, m, ws, wy, sy, theta, col, head, info) integer *n, *nsub, *ind, *nenter, *ileave, *indx2, *iupdat; logical *updatd; doublereal *wn, *wn1; integer *m; doublereal *ws, *wy, *sy, *theta; integer *col, *head, *info; { /* System generated locals */ integer wn_dim1, wn_offset, wn1_dim1, wn1_offset, ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, i__1, i__2, i__3; /* Local variables */ static integer dend, pend; static integer upcl; static doublereal temp1, temp2, temp3, temp4; static integer i__, k; static integer ipntr, jpntr, k1, m2, dbegin, is, js, iy, jy, pbegin, is1, js1, col2; /* Parameter adjustments */ --indx2; --ind; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; wn1_dim1 = 2 * *m; wn1_offset = 1 + wn1_dim1 * 1; wn1 -= wn1_offset; wn_dim1 = 2 * *m; wn_offset = 1 + wn_dim1 * 1; wn -= wn_offset; /* Function Body */ if (*updatd) { if (*iupdat > *m) { /* shift old part of WN1. */ i__1 = *m - 1; for (jy = 1; jy <= i__1; ++jy) { js = *m + jy; i__2 = *m - jy; dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[ jy + jy * wn1_dim1], &c__1); i__2 = *m - jy; dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[ js + js * wn1_dim1], &c__1); i__2 = *m - 1; dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[ *m + 1 + jy * wn1_dim1], &c__1); /* L10: */ } } /* put new rows in blocks (1,1), (2,1) and (2,2). */ pbegin = 1; pend = *nsub; dbegin = *nsub + 1; dend = *n; iy = *col; is = *m + *col; ipntr = *head + *col - 1; if (ipntr > *m) { ipntr -= *m; } jpntr = *head; i__1 = *col; for (jy = 1; jy <= i__1; ++jy) { js = *m + jy; temp1 = 0.; temp2 = 0.; temp3 = 0.; /* compute element jy of row 'col' of Y'ZZ'Y */ i__2 = pend; for (k = pbegin; k <= i__2; ++k) { k1 = ind[k]; temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1]; /* L15: */ } /* compute elements jy of row 'col' of L_a and S'AA'S */ i__2 = dend; for (k = dbegin; k <= i__2; ++k) { k1 = ind[k]; temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1]; temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1]; /* L16: */ } wn1[iy + jy * wn1_dim1] = temp1; wn1[is + js * wn1_dim1] = temp2; wn1[is + jy * wn1_dim1] = temp3; jpntr = jpntr % *m + 1; /* L20: */ } /* put new column in block (2,1). */ jy = *col; jpntr = *head + *col - 1; if (jpntr > *m) { jpntr -= *m; } ipntr = *head; i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { is = *m + i__; temp3 = 0.; /* compute element i of column 'col' of R_z */ i__2 = pend; for (k = pbegin; k <= i__2; ++k) { k1 = ind[k]; temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1]; /* L25: */ } ipntr = ipntr % *m + 1; wn1[is + jy * wn1_dim1] = temp3; /* L30: */ } upcl = *col - 1; } else { upcl = *col; } /* modify the old parts in blocks (1,1) and (2,2) due to changes */ /* in the set of free variables. */ ipntr = *head; i__1 = upcl; for (iy = 1; iy <= i__1; ++iy) { is = *m + iy; jpntr = *head; i__2 = iy; for (jy = 1; jy <= i__2; ++jy) { js = *m + jy; temp1 = 0.; temp2 = 0.; temp3 = 0.; temp4 = 0.; i__3 = *nenter; for (k = 1; k <= i__3; ++k) { k1 = indx2[k]; temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1]; temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1]; /* L35: */ } i__3 = *n; for (k = *ileave; k <= i__3; ++k) { k1 = indx2[k]; temp3 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1]; temp4 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1]; /* L36: */ } wn1[iy + jy * wn1_dim1] = wn1[iy + jy * wn1_dim1] + temp1 - temp3; wn1[is + js * wn1_dim1] = wn1[is + js * wn1_dim1] - temp2 + temp4; jpntr = jpntr % *m + 1; /* L40: */ } ipntr = ipntr % *m + 1; /* L45: */ } /* modify the old parts in block (2,1). */ ipntr = *head; i__1 = *m + upcl; for (is = *m + 1; is <= i__1; ++is) { jpntr = *head; i__2 = upcl; for (jy = 1; jy <= i__2; ++jy) { temp1 = 0.; temp3 = 0.; i__3 = *nenter; for (k = 1; k <= i__3; ++k) { k1 = indx2[k]; temp1 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1]; /* L50: */ } i__3 = *n; for (k = *ileave; k <= i__3; ++k) { k1 = indx2[k]; temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1]; /* L51: */ } if (is <= jy + *m) { wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] + temp1 - temp3; } else { wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] - temp1 + temp3; } jpntr = jpntr % *m + 1; /* L55: */ } ipntr = ipntr % *m + 1; /* L60: */ } /* Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] */ /* [-L_a +R_z S'AA'S*theta] */ m2 = *m << 1; i__1 = *col; for (iy = 1; iy <= i__1; ++iy) { is = *col + iy; is1 = *m + iy; i__2 = iy; for (jy = 1; jy <= i__2; ++jy) { js = *col + jy; js1 = *m + jy; wn[jy + iy * wn_dim1] = wn1[iy + jy * wn1_dim1] / *theta; wn[js + is * wn_dim1] = wn1[is1 + js1 * wn1_dim1] * *theta; /* L65: */ } i__2 = iy - 1; for (jy = 1; jy <= i__2; ++jy) { wn[jy + is * wn_dim1] = -wn1[is1 + jy * wn1_dim1]; /* L66: */ } i__2 = *col; for (jy = iy; jy <= i__2; ++jy) { wn[jy + is * wn_dim1] = wn1[is1 + jy * wn1_dim1]; /* L67: */ } wn[iy + iy * wn_dim1] += sy[iy + iy * sy_dim1]; /* L70: */ } /* Form the upper triangle of */ /* WN= [ LL' L^-1(-L_a'+R_z')] */ /* [(-L_a +R_z)L'^-1 S'AA'S*theta ] */ /* first Cholesky factor (1,1) block of wn to get LL' */ /* with L' stored in the upper triangle of wn. */ dpofa_(&wn[wn_offset], &m2, col, info); if (*info != 0) { *info = -1; return 0; } /* then form L^-1(-L_a'+R_z') in the (1,2) block. */ col2 = *col << 1; i__1 = col2; for (js = *col + 1; js <= i__1; ++js) { dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info); /* L71: */ } /* Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the */ /* upper triangle of (2,2) block of wn. */ i__1 = col2; for (is = *col + 1; is <= i__1; ++is) { i__2 = col2; for (js = is; js <= i__2; ++js) { wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], &c__1, &wn[js * wn_dim1 + 1], &c__1); /* L74: */ } /* L72: */ } /* Cholesky factorization of (2,2) block of wn. */ dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info); if (*info != 0) { *info = -2; return 0; } return 0; } /* formk_ */ /* ======================= The end of formk ============================== */ /* Subroutine */ int formt_(m, wt, sy, ss, col, theta, info) integer *m; doublereal *wt, *sy, *ss; integer *col; doublereal *theta; integer *info; { /* System generated locals */ integer wt_dim1, wt_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, i__1, i__2, i__3; /* Local variables */ static doublereal ddum; static integer i__, j, k; static integer k1; /* Parameter adjustments */ ss_dim1 = *m; ss_offset = 1 + ss_dim1 * 1; ss -= ss_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; wt_dim1 = *m; wt_offset = 1 + wt_dim1 * 1; wt -= wt_offset; /* Function Body */ i__1 = *col; for (j = 1; j <= i__1; ++j) { wt[j * wt_dim1 + 1] = *theta * ss[j * ss_dim1 + 1]; /* L52: */ } i__1 = *col; for (i__ = 2; i__ <= i__1; ++i__) { i__2 = *col; for (j = i__; j <= i__2; ++j) { k1 = min(i__,j) - 1; ddum = 0.; i__3 = k1; for (k = 1; k <= i__3; ++k) { ddum += sy[i__ + k * sy_dim1] * sy[j + k * sy_dim1] / sy[k + k * sy_dim1]; /* L53: */ } wt[i__ + j * wt_dim1] = ddum + *theta * ss[i__ + j * ss_dim1]; /* L54: */ } /* L55: */ } /* Cholesky factorize T to J*J' with */ /* J' stored in the upper triangle of wt. */ dpofa_(&wt[wt_offset], m, col, info); if (*info != 0) { *info = -3; } return 0; } /* formt_ */ /* ======================= The end of formt ============================== */ /* Subroutine */ int freev_(n, nfree, index, nenter, ileave, indx2, iwhere, wrk, updatd, cnstnd, iprint, iter) integer *n, *nfree, *index, *nenter, *ileave, *indx2, *iwhere; logical *wrk, *updatd, *cnstnd; integer *iprint, *iter; { /* System generated locals */ integer i__1; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(); /* Local variables */ static integer iact, i__, k; /* Parameter adjustments */ --iwhere; --indx2; --index; /* Function Body */ *nenter = 0; *ileave = *n + 1; if (*iter > 0 && *cnstnd) { /* count the entering and leaving variables. */ i__1 = *nfree; for (i__ = 1; i__ <= i__1; ++i__) { k = index[i__]; if (iwhere[k] > 0) { --(*ileave); indx2[*ileave] = k; } /* L20: */ } i__1 = *n; for (i__ = *nfree + 1; i__ <= i__1; ++i__) { k = index[i__]; if (iwhere[k] <= 0) { ++(*nenter); indx2[*nenter] = k; } /* L22: */ } } *wrk = *ileave < *n + 1 || *nenter > 0 || *updatd; /* Find the index set of free and active variables at the GCP. */ *nfree = 0; iact = *n + 1; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (iwhere[i__] <= 0) { ++(*nfree); index[*nfree] = i__; } else { --iact; index[iact] = i__; } /* L24: */ } return 0; } /* freev_ */ /* ======================= The end of freev ============================== */ /* Subroutine */ int hpsolb_(n, t, iorder, iheap) integer *n; doublereal *t; integer *iorder, *iheap; { /* System generated locals */ integer i__1; /* Local variables */ static doublereal ddum; static integer i__, j, k, indxin, indxou; static doublereal out; /* ************ */ /* Parameter adjustments */ --iorder; --t; /* Function Body */ if (*iheap == 0) { /* Rearrange the elements t(1) to t(n) to form a heap. */ i__1 = *n; for (k = 2; k <= i__1; ++k) { ddum = t[k]; indxin = iorder[k]; /* Add ddum to the heap. */ i__ = k; L10: if (i__ > 1) { j = i__ / 2; if (ddum < t[j]) { t[i__] = t[j]; iorder[i__] = iorder[j]; i__ = j; goto L10; } } t[i__] = ddum; iorder[i__] = indxin; /* L20: */ } } /* Assign to 'out' the value of t(1), the least member of the heap, */ /* and rearrange the remaining members to form a heap as */ /* elements 1 to n-1 of t. */ if (*n > 1) { i__ = 1; out = t[1]; indxou = iorder[1]; ddum = t[*n]; indxin = iorder[*n]; /* Restore the heap */ L30: j = i__ + i__; if (j <= *n - 1) { if (t[j + 1] < t[j]) { ++j; } if (t[j] < ddum) { t[i__] = t[j]; iorder[i__] = iorder[j]; i__ = j; goto L30; } } t[i__] = ddum; iorder[i__] = indxin; /* Put the least member in t(n). */ t[*n] = out; iorder[*n] = indxou; } return 0; } /* hpsolb_ */ /* ====================== The end of hpsolb ============================== */ /* Subroutine */ int lnsrlb_(n, l, u, nbd, x, f, fold, gd, gdold, g, d__, r__, t, z__, stp, dnorm, dtd, xstep, stpmx, iter, ifun, iback, nfgv, info, task, boxed, cnstnd, csave, isave, dsave, task_len, csave_len) integer *n; doublereal *l, *u; integer *nbd; doublereal *x, *f, *fold, *gd, *gdold, *g, *d__, *r__, *t, *z__, *stp, *dnorm, *dtd, *xstep, *stpmx; integer *iter, *ifun, *iback, *nfgv, *info; char *task; logical *boxed, *cnstnd; char *csave; integer *isave; doublereal *dsave; ftnlen task_len; ftnlen csave_len; { /* System generated locals */ integer i__1; doublereal d__1; /* Builtin functions */ double sqrt(); /* Local variables */ static integer i__; static doublereal a1, a2; /* ********** */ /* Parameter adjustments */ --z__; --t; --r__; --d__; --g; --x; --nbd; --u; --l; --isave; --dsave; /* Function Body */ if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) { goto L556; } *dtd = ddot_(n, &d__[1], &c__1, &d__[1], &c__1); *dnorm = sqrt(*dtd); /* Determine the maximum step length. */ *stpmx = 1e10; if (*cnstnd) { if (*iter == 0) { *stpmx = 1.; } else { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { a1 = d__[i__]; if (nbd[i__] != 0) { if (a1 < 0. && nbd[i__] <= 2) { a2 = l[i__] - x[i__]; if (a2 >= 0.) { *stpmx = 0.; } else if (a1 * *stpmx < a2) { *stpmx = a2 / a1; } } else if (a1 > 0. && nbd[i__] >= 2) { a2 = u[i__] - x[i__]; if (a2 <= 0.) { *stpmx = 0.; } else if (a1 * *stpmx > a2) { *stpmx = a2 / a1; } } } /* L43: */ } } } if (*iter == 0 && ! (*boxed)) { /* Computing MIN */ d__1 = 1. / *dnorm; *stp = min(d__1,*stpmx); } else { *stp = 1.; } dcopy_(n, &x[1], &c__1, &t[1], &c__1); dcopy_(n, &g[1], &c__1, &r__[1], &c__1); *fold = *f; *ifun = 0; *iback = 0; s_copy(csave, "START", (ftnlen)60, (ftnlen)5); L556: *gd = ddot_(n, &g[1], &c__1, &d__[1], &c__1); if (*ifun == 0) { *gdold = *gd; if (*gd >= 0.) { /* the directional derivative >=0. */ /* Line search is impossible. */ *info = -4; return 0; } } dcsrch_(f, gd, stp, &c_b275, &c_b276, &c_b277, &c_b9, stpmx, csave, & isave[1], &dsave[1], (ftnlen)60); *xstep = *stp * *dnorm; if (s_cmp(csave, "CONV", (ftnlen)4, (ftnlen)4) != 0 && s_cmp(csave, "WARN" , (ftnlen)4, (ftnlen)4) != 0) { s_copy(task, "FG_LNSRCH", (ftnlen)60, (ftnlen)9); ++(*ifun); ++(*nfgv); *iback = *ifun - 1; if (*stp == 1.) { dcopy_(n, &z__[1], &c__1, &x[1], &c__1); } else { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { x[i__] = *stp * d__[i__] + t[i__]; /* L41: */ } } } else { s_copy(task, "NEW_X", (ftnlen)60, (ftnlen)5); } return 0; } /* lnsrlb_ */ /* ======================= The end of lnsrlb ============================= */ /* Subroutine */ int matupd_(n, m, ws, wy, sy, ss, d__, r__, itail, iupdat, col, head, theta, rr, dr, stp, dtd) integer *n, *m; doublereal *ws, *wy, *sy, *ss, *d__, *r__; integer *itail, *iupdat, *col, *head; doublereal *theta, *rr, *dr, *stp, *dtd; { /* System generated locals */ integer ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, i__1, i__2; /* Local variables */ static integer j; static integer pointr; /* Parameter adjustments */ --r__; --d__; ss_dim1 = *m; ss_offset = 1 + ss_dim1 * 1; ss -= ss_offset; sy_dim1 = *m; sy_offset = 1 + sy_dim1 * 1; sy -= sy_offset; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; /* Function Body */ if (*iupdat <= *m) { *col = *iupdat; *itail = (*head + *iupdat - 2) % *m + 1; } else { *itail = *itail % *m + 1; *head = *head % *m + 1; } /* Update matrices WS and WY. */ dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1); dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1); /* Set theta=yy/ys. */ *theta = *rr / *dr; /* Form the middle matrix in B. */ /* update the upper triangle of SS, */ /* and the lower triangle of SY: */ if (*iupdat > *m) { /* move old information */ i__1 = *col - 1; for (j = 1; j <= i__1; ++j) { dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1] , &c__1); i__2 = *col - j; dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j * sy_dim1], &c__1); /* L50: */ } } /* add new information: the last row of SY */ /* and the last column of SS: */ pointr = *head; i__1 = *col - 1; for (j = 1; j <= i__1; ++j) { sy[*col + j * sy_dim1] = ddot_(n, &d__[1], &c__1, &wy[pointr * wy_dim1 + 1], &c__1); ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, & d__[1], &c__1); pointr = pointr % *m + 1; /* L51: */ } if (*stp == 1.) { ss[*col + *col * ss_dim1] = *dtd; } else { ss[*col + *col * ss_dim1] = *stp * *stp * *dtd; } sy[*col + *col * sy_dim1] = *dr; return 0; } /* matupd_ */ /* Subroutine */ int prn2lb_(n, x, f, g, iprint, itfile, iter, nfgv, nact, sbgnrm, nint, word, iword, iback, stp, xstep, word_len) integer *n; doublereal *x, *f, *g; integer *iprint, *itfile, *iter, *nfgv, *nact; doublereal *sbgnrm; integer *nint; char *word; integer *iword, *iback; doublereal *stp, *xstep; ftnlen word_len; { /* Parameter adjustments */ --g; --x; /* Function Body */ if (*iword == 0) { /* the subspace minimization converged. */ s_copy(word, "con", (ftnlen)3, (ftnlen)3); } else if (*iword == 1) { /* the subspace minimization stopped at a bound. */ s_copy(word, "bnd", (ftnlen)3, (ftnlen)3); } else if (*iword == 5) { /* the truncated Newton step has been used. */ s_copy(word, "TNT", (ftnlen)3, (ftnlen)3); } else { s_copy(word, "---", (ftnlen)3, (ftnlen)3); } return 0; } /* prn2lb_ */ /* ======================= The end of prn2lb ============================= */ /* Subroutine */ int projgr_(n, l, u, nbd, x, g, sbgnrm) integer *n; doublereal *l, *u; integer *nbd; doublereal *x, *g, *sbgnrm; { /* System generated locals */ integer i__1; doublereal d__1, d__2; /* Local variables */ static integer i__; static doublereal gi; /* Parameter adjustments */ --g; --x; --nbd; --u; --l; /* Function Body */ *sbgnrm = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { gi = g[i__]; if (nbd[i__] != 0) { if (gi < 0.) { if (nbd[i__] >= 2) { /* Computing MAX */ d__1 = x[i__] - u[i__]; gi = max(d__1,gi); } } else { if (nbd[i__] <= 2) { /* Computing MIN */ d__1 = x[i__] - l[i__]; gi = min(d__1,gi); } } } /* Computing MAX */ d__1 = *sbgnrm, d__2 = abs(gi); *sbgnrm = max(d__1,d__2); /* L15: */ } return 0; } /* projgr_ */ /* ======================= The end of projgr ============================= */ /* Subroutine */ int subsm_(n, m, nsub, ind, l, u, nbd, x, d__, ws, wy, theta, col, head, iword, wv, wn, iprint, info) integer *n, *m, *nsub, *ind; doublereal *l, *u; integer *nbd; doublereal *x, *d__, *ws, *wy, *theta; integer *col, *head, *iword; doublereal *wv, *wn; integer *iprint, *info; { /* System generated locals */ integer ws_dim1, ws_offset, wy_dim1, wy_offset, wn_dim1, wn_offset, i__1, i__2; /* Local variables */ static doublereal temp1, temp2; static integer i__, j, k; static doublereal alpha; static integer m2; static doublereal dk; static integer js, jy, pointr, ibd, col2; /* Parameter adjustments */ --d__; --x; --nbd; --u; --l; wn_dim1 = 2 * *m; wn_offset = 1 + wn_dim1 * 1; wn -= wn_offset; --wv; wy_dim1 = *n; wy_offset = 1 + wy_dim1 * 1; wy -= wy_offset; ws_dim1 = *n; ws_offset = 1 + ws_dim1 * 1; ws -= ws_offset; --ind; /* Function Body */ if (*nsub <= 0) { return 0; } /* Compute wv = W'Zd. */ pointr = *head; i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { temp1 = 0.; temp2 = 0.; i__2 = *nsub; for (j = 1; j <= i__2; ++j) { k = ind[j]; temp1 += wy[k + pointr * wy_dim1] * d__[j]; temp2 += ws[k + pointr * ws_dim1] * d__[j]; /* L10: */ } wv[i__] = temp1; wv[*col + i__] = *theta * temp2; pointr = pointr % *m + 1; /* L20: */ } /* Compute wv:=K^(-1)wv. */ m2 = *m << 1; col2 = *col << 1; dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info); if (*info != 0) { return 0; } i__1 = *col; for (i__ = 1; i__ <= i__1; ++i__) { wv[i__] = -wv[i__]; /* L25: */ } dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info); if (*info != 0) { return 0; } /* Compute d = (1/theta)d + (1/theta**2)Z'W wv. */ pointr = *head; i__1 = *col; for (jy = 1; jy <= i__1; ++jy) { js = *col + jy; i__2 = *nsub; for (i__ = 1; i__ <= i__2; ++i__) { k = ind[i__]; d__[i__] = d__[i__] + wy[k + pointr * wy_dim1] * wv[jy] / *theta + ws[k + pointr * ws_dim1] * wv[js]; /* L30: */ } pointr = pointr % *m + 1; /* L40: */ } i__1 = *nsub; for (i__ = 1; i__ <= i__1; ++i__) { d__[i__] /= *theta; /* L50: */ } /* Backtrack to the feasible region. */ alpha = 1.; temp1 = alpha; i__1 = *nsub; for (i__ = 1; i__ <= i__1; ++i__) { k = ind[i__]; dk = d__[i__]; if (nbd[k] != 0) { if (dk < 0. && nbd[k] <= 2) { temp2 = l[k] - x[k]; if (temp2 >= 0.) { temp1 = 0.; } else if (dk * alpha < temp2) { temp1 = temp2 / dk; } } else if (dk > 0. && nbd[k] >= 2) { temp2 = u[k] - x[k]; if (temp2 <= 0.) { temp1 = 0.; } else if (dk * alpha > temp2) { temp1 = temp2 / dk; } } if (temp1 < alpha) { alpha = temp1; ibd = i__; } } /* L60: */ } if (alpha < 1.) { dk = d__[ibd]; k = ind[ibd]; if (dk > 0.) { x[k] = u[k]; d__[ibd] = 0.; } else if (dk < 0.) { x[k] = l[k]; d__[ibd] = 0.; } } i__1 = *nsub; for (i__ = 1; i__ <= i__1; ++i__) { k = ind[i__]; x[k] += alpha * d__[i__]; /* L70: */ } if (alpha < 1.) { *iword = 1; } else { *iword = 0; } return 0; } /* subsm_ */ /* ====================== The end of subsm =============================== */ /* Subroutine */ int dcsrch_(f, g, stp, ftol, gtol, xtol, stpmin, stpmax, task, isave, dsave, task_len) doublereal *f, *g, *stp, *ftol, *gtol, *xtol, *stpmin, *stpmax; char *task; integer *isave; doublereal *dsave; ftnlen task_len; { /* System generated locals */ doublereal d__1; /* Local variables */ static integer stage; static doublereal finit, ginit, width, ftest, gtest, stmin, stmax, width1, fm, gm, fx, fy, gx, gy; static logical brackt; static doublereal fxm, fym, gxm, gym, stx, sty; /* Parameter adjustments */ --dsave; --isave; /* Function Body */ if (s_cmp(task, "START", (ftnlen)5, (ftnlen)5) == 0) { /* Check the input arguments for errors. */ if (*stp < *stpmin) { s_copy(task, "ERROR: STP .LT. STPMIN", task_len, (ftnlen)22); } if (*stp > *stpmax) { s_copy(task, "ERROR: STP .GT. STPMAX", task_len, (ftnlen)22); } if (*g >= 0.) { s_copy(task, "ERROR: INITIAL G .GE. ZERO", task_len, (ftnlen)26); } if (*ftol < 0.) { s_copy(task, "ERROR: FTOL .LT. ZERO", task_len, (ftnlen)21); } if (*gtol < 0.) { s_copy(task, "ERROR: GTOL .LT. ZERO", task_len, (ftnlen)21); } if (*xtol < 0.) { s_copy(task, "ERROR: XTOL .LT. ZERO", task_len, (ftnlen)21); } if (*stpmin < 0.) { s_copy(task, "ERROR: STPMIN .LT. ZERO", task_len, (ftnlen)23); } if (*stpmax < *stpmin) { s_copy(task, "ERROR: STPMAX .LT. STPMIN", task_len, (ftnlen)25); } /* Exit if there are errors on input. */ if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) { return 0; } /* Initialize local variables. */ brackt = FALSE_; stage = 1; finit = *f; ginit = *g; gtest = *ftol * ginit; width = *stpmax - *stpmin; width1 = width / .5; /* The variables stx, fx, gx contain the values of the step, */ /* function, and derivative at the best step. */ /* The variables sty, fy, gy contain the value of the step, */ /* function, and derivative at sty. */ /* The variables stp, f, g contain the values of the step, */ /* function, and derivative at stp. */ stx = 0.; fx = finit; gx = ginit; sty = 0.; fy = finit; gy = ginit; stmin = 0.; stmax = *stp + *stp * 4.; s_copy(task, "FG", task_len, (ftnlen)2); goto L1000; } else { /* Restore local variables. */ if (isave[1] == 1) { brackt = TRUE_; } else { brackt = FALSE_; } stage = isave[2]; ginit = dsave[1]; gtest = dsave[2]; gx = dsave[3]; gy = dsave[4]; finit = dsave[5]; fx = dsave[6]; fy = dsave[7]; stx = dsave[8]; sty = dsave[9]; stmin = dsave[10]; stmax = dsave[11]; width = dsave[12]; width1 = dsave[13]; } /* If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */ /* algorithm enters the second stage. */ ftest = finit + *stp * gtest; if (stage == 1 && *f <= ftest && *g >= 0.) { stage = 2; } /* Test for warnings. */ if (brackt && (*stp <= stmin || *stp >= stmax)) { s_copy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS", task_len, ( ftnlen)41); } if (brackt && stmax - stmin <= *xtol * stmax) { s_copy(task, "WARNING: XTOL TEST SATISFIED", task_len, (ftnlen)28); } if (*stp == *stpmax && *f <= ftest && *g <= gtest) { s_copy(task, "WARNING: STP = STPMAX", task_len, (ftnlen)21); } if (*stp == *stpmin && (*f > ftest || *g >= gtest)) { s_copy(task, "WARNING: STP = STPMIN", task_len, (ftnlen)21); } /* Test for convergence. */ if (*f <= ftest && abs(*g) <= *gtol * (-ginit)) { s_copy(task, "CONVERGENCE", task_len, (ftnlen)11); } /* Test for termination. */ if (s_cmp(task, "WARN", (ftnlen)4, (ftnlen)4) == 0 || s_cmp(task, "CONV", (ftnlen)4, (ftnlen)4) == 0) { goto L1000; } /* A modified function is used to predict the step during the */ /* first stage if a lower function value has been obtained but */ /* the decrease is not sufficient. */ if (stage == 1 && *f <= fx && *f > ftest) { /* Define the modified function and derivative values. */ fm = *f - *stp * gtest; fxm = fx - stx * gtest; fym = fy - sty * gtest; gm = *g - gtest; gxm = gx - gtest; gym = gy - gtest; /* Call dcstep to update stx, sty, and to compute the new step. */ dcstep_(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, &fm, &gm, &brackt, & stmin, &stmax); /* Reset the function and derivative values for f. */ fx = fxm + stx * gtest; fy = fym + sty * gtest; gx = gxm + gtest; gy = gym + gtest; } else { /* Call dcstep to update stx, sty, and to compute the new step. */ dcstep_(&stx, &fx, &gx, &sty, &fy, &gy, stp, f, g, &brackt, &stmin, & stmax); } /* Decide if a bisection step is needed. */ if (brackt) { if ((d__1 = sty - stx, abs(d__1)) >= width1 * .66) { *stp = stx + (sty - stx) * .5; } width1 = width; width = (d__1 = sty - stx, abs(d__1)); } /* Set the minimum and maximum steps allowed for stp. */ if (brackt) { stmin = min(stx,sty); stmax = max(stx,sty); } else { stmin = *stp + (*stp - stx) * 1.1; stmax = *stp + (*stp - stx) * 4.; } /* Force the step to be within the bounds stpmax and stpmin. */ *stp = max(*stp,*stpmin); *stp = min(*stp,*stpmax); /* If further progress is not possible, let stp be the best */ /* point obtained during the search. */ if ((brackt && (*stp <= stmin || *stp >= stmax)) || ((brackt) && (stmax - stmin <= *xtol * stmax))) { *stp = stx; } /* Obtain another function and derivative. */ s_copy(task, "FG", task_len, (ftnlen)2); L1000: /* Save local variables. */ if (brackt) { isave[1] = 1; } else { isave[1] = 0; } isave[2] = stage; dsave[1] = ginit; dsave[2] = gtest; dsave[3] = gx; dsave[4] = gy; dsave[5] = finit; dsave[6] = fx; dsave[7] = fy; dsave[8] = stx; dsave[9] = sty; dsave[10] = stmin; dsave[11] = stmax; dsave[12] = width; dsave[13] = width1; return 0; } /* dcsrch_ */ /* ====================== The end of dcsrch ============================== */ /* Subroutine */ int dcstep_(stx, fx, dx, sty, fy, dy, stp, fp, dp, brackt, stpmin, stpmax) doublereal *stx, *fx, *dx, *sty, *fy, *dy, *stp, *fp, *dp; logical *brackt; doublereal *stpmin, *stpmax; { /* System generated locals */ doublereal d__1, d__2, d__3; /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal sgnd, stpc, stpf, stpq, p, q, gamma, r__, s, theta; sgnd = *dp * (*dx / abs(*dx)); /* First case: A higher function value. The minimum is bracketed. */ /* If the cubic step is closer to stx than the quadratic step, the */ /* cubic step is taken, otherwise the average of the cubic and */ /* quadratic steps is taken. */ if (*fp > *fx) { theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp; /* Computing MAX */ d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs( *dp); s = max(d__1,d__2); /* Computing 2nd power */ d__1 = theta / s; gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s)); if (*stp < *stx) { gamma = -gamma; } p = gamma - *dx + theta; q = gamma - *dx + gamma + *dp; r__ = p / q; stpc = *stx + r__ * (*stp - *stx); stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2. * (*stp - *stx); if ((d__1 = stpc - *stx, abs(d__1)) < (d__2 = stpq - *stx, abs(d__2))) { stpf = stpc; } else { stpf = stpc + (stpq - stpc) / 2.; } *brackt = TRUE_; /* Second case: A lower function value and derivatives of opposite */ /* sign. The minimum is bracketed. If the cubic step is farther from */ /* stp than the secant step, the cubic step is taken, otherwise the */ /* secant step is taken. */ } else if (sgnd < 0.) { theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp; /* Computing MAX */ d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs( *dp); s = max(d__1,d__2); /* Computing 2nd power */ d__1 = theta / s; gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s)); if (*stp > *stx) { gamma = -gamma; } p = gamma - *dp + theta; q = gamma - *dp + gamma + *dx; r__ = p / q; stpc = *stp + r__ * (*stx - *stp); stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp); if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs(d__2))) { stpf = stpc; } else { stpf = stpq; } *brackt = TRUE_; /* Third case: A lower function value, derivatives of the same sign, */ /* and the magnitude of the derivative decreases. */ } else if (abs(*dp) < abs(*dx)) { /* The cubic step is computed only if the cubic tends to infinity */ /* in the direction of the step or if the minimum of the cubic */ /* is beyond stp. Otherwise the cubic step is defined to be the */ /* secant step. */ theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp; /* Computing MAX */ d__1 = abs(theta), d__2 = abs(*dx), d__1 = max(d__1,d__2), d__2 = abs( *dp); s = max(d__1,d__2); /* The case gamma = 0 only arises if the cubic does not tend */ /* to infinity in the direction of the step. */ /* Computing MAX */ /* Computing 2nd power */ d__3 = theta / s; d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s); gamma = s * sqrt((max(d__1,d__2))); if (*stp > *stx) { gamma = -gamma; } p = gamma - *dp + theta; q = gamma + (*dx - *dp) + gamma; r__ = p / q; if (r__ < 0. && gamma != 0.) { stpc = *stp + r__ * (*stx - *stp); } else if (*stp > *stx) { stpc = *stpmax; } else { stpc = *stpmin; } stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp); if (*brackt) { /* A minimizer has been bracketed. If the cubic step is */ /* closer to stp than the secant step, the cubic step is */ /* taken, otherwise the secant step is taken. */ if ((d__1 = stpc - *stp, abs(d__1)) < (d__2 = stpq - *stp, abs( d__2))) { stpf = stpc; } else { stpf = stpq; } if (*stp > *stx) { /* Computing MIN */ d__1 = *stp + (*sty - *stp) * .66; stpf = min(d__1,stpf); } else { /* Computing MAX */ d__1 = *stp + (*sty - *stp) * .66; stpf = max(d__1,stpf); } } else { /* A minimizer has not been bracketed. If the cubic step is */ /* farther from stp than the secant step, the cubic step is */ /* taken, otherwise the secant step is taken. */ if ((d__1 = stpc - *stp, abs(d__1)) > (d__2 = stpq - *stp, abs( d__2))) { stpf = stpc; } else { stpf = stpq; } stpf = min(*stpmax,stpf); stpf = max(*stpmin,stpf); } /* Fourth case: A lower function value, derivatives of the */ /* same sign, and the magnitude of the derivative does not */ /* decrease. If the minimum is not bracketed, the step is either */ /* stpmin or stpmax, otherwise the cubic step is taken. */ } else { if (*brackt) { theta = (*fp - *fy) * 3. / (*sty - *stp) + *dy + *dp; /* Computing MAX */ d__1 = abs(theta), d__2 = abs(*dy), d__1 = max(d__1,d__2), d__2 = abs(*dp); s = max(d__1,d__2); /* Computing 2nd power */ d__1 = theta / s; gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s)); if (*stp > *sty) { gamma = -gamma; } p = gamma - *dp + theta; q = gamma - *dp + gamma + *dy; r__ = p / q; stpc = *stp + r__ * (*sty - *stp); stpf = stpc; } else if (*stp > *stx) { stpf = *stpmax; } else { stpf = *stpmin; } } /* Update the interval which contains a minimizer. */ if (*fp > *fx) { *sty = *stp; *fy = *fp; *dy = *dp; } else { if (sgnd < 0.) { *sty = *stx; *fy = *fx; *dy = *dx; } *stx = *stp; *fx = *fp; *dx = *dp; } /* Compute the new step. */ *stp = stpf; return 0; } /* dcstep_ */ doublereal dnrm2_(n, x, incx) integer *n; doublereal *x; integer *incx; { /* System generated locals */ integer i__1, i__2; doublereal ret_val, d__1, d__2, d__3; /* Builtin functions */ double sqrt(); /* Local variables */ static integer i__; static doublereal scale; /* Parameter adjustments */ --x; /* Function Body */ ret_val = 0.; scale = 0.; i__1 = *n; i__2 = *incx; for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { /* Computing MAX */ d__2 = scale, d__3 = (d__1 = x[i__], abs(d__1)); scale = max(d__2,d__3); /* L10: */ } if (scale == 0.) { return ret_val; } i__2 = *n; i__1 = *incx; for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) { /* Computing 2nd power */ d__1 = x[i__] / scale; ret_val += d__1 * d__1; /* L20: */ } ret_val = scale * sqrt(ret_val); return ret_val; } /* dnrm2_ */ /* ====================== The end of dnrm2 =============================== */ doublereal dpmeps_() { /* Initialized data */ static doublereal zero = 0.; static doublereal one = 1.; static doublereal two = 2.; /* System generated locals */ integer i__1; doublereal ret_val; /* Local variables */ static doublereal beta; static integer irnd; static doublereal temp, temp1, a, b; static integer i__; static doublereal betah; static integer ibeta, negep; static doublereal tempa; static integer itemp, it; static doublereal betain; /* determine ibeta, beta ala malcolm. */ a = one; b = one; L10: a += a; temp = a + one; temp1 = temp - a; if (temp1 - one == zero) { goto L10; } L20: b += b; temp = a + b; itemp = (integer) (temp - a); if (itemp == 0) { goto L20; } ibeta = itemp; beta = (doublereal) ibeta; /* determine it, irnd. */ it = 0; b = one; L30: ++it; b *= beta; temp = b + one; temp1 = temp - b; if (temp1 - one == zero) { goto L30; } irnd = 0; betah = beta / two; temp = a + betah; if (temp - a != zero) { irnd = 1; } tempa = a + beta; temp = tempa + betah; if (irnd == 0 && temp - tempa != zero) { irnd = 2; } /* determine dpmeps. */ negep = it + 3; betain = one / beta; a = one; i__1 = negep; for (i__ = 1; i__ <= i__1; ++i__) { a *= betain; /* L40: */ } L50: temp = one + a; if (temp - one != zero) { goto L60; } a *= beta; goto L50; L60: ret_val = a; if (ibeta == 2 || irnd == 0) { goto L70; } a = a * (one + a) / two; temp = one + a; if (temp - one != zero) { ret_val = a; } L70: return ret_val; } /* dpmeps_ */ /* ====================== The end of dpmeps ============================== */ /* Subroutine */ int daxpy_(n, da, dx, incx, dy, incy) integer *n; doublereal *da, *dx; integer *incx; doublereal *dy; integer *incy; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, m, ix, iy, mp1; /* Parameter adjustments */ --dy; --dx; /* Function Body */ if (*n <= 0) { return 0; } if (*da == 0.) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { dy[iy] += *da * dx[ix]; ix += *incx; iy += *incy; /* L10: */ } return 0; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 4; if (m == 0) { goto L40; } i__1 = m; for (i__ = 1; i__ <= i__1; ++i__) { dy[i__] += *da * dx[i__]; /* L30: */ } if (*n < 4) { return 0; } L40: mp1 = m + 1; i__1 = *n; for (i__ = mp1; i__ <= i__1; i__ += 4) { dy[i__] += *da * dx[i__]; dy[i__ + 1] += *da * dx[i__ + 1]; dy[i__ + 2] += *da * dx[i__ + 2]; dy[i__ + 3] += *da * dx[i__ + 3]; /* L50: */ } return 0; } /* daxpy_ */ /* ====================== The end of daxpy =============================== */ /* Subroutine */ int dcopy_(n, dx, incx, dy, incy) integer *n; doublereal *dx; integer *incx; doublereal *dy; integer *incy; { /* System generated locals */ integer i__1; /* Local variables */ static integer i__, m, ix, iy, mp1; /* copies a vector, x, to a vector, y. */ /* uses unrolled loops for increments equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* Parameter adjustments */ --dy; --dx; /* Function Body */ if (*n <= 0) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { dy[iy] = dx[ix]; ix += *incx; iy += *incy; /* L10: */ } return 0; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 7; if (m == 0) { goto L40; } i__1 = m; for (i__ = 1; i__ <= i__1; ++i__) { dy[i__] = dx[i__]; /* L30: */ } if (*n < 7) { return 0; } L40: mp1 = m + 1; i__1 = *n; for (i__ = mp1; i__ <= i__1; i__ += 7) { dy[i__] = dx[i__]; dy[i__ + 1] = dx[i__ + 1]; dy[i__ + 2] = dx[i__ + 2]; dy[i__ + 3] = dx[i__ + 3]; dy[i__ + 4] = dx[i__ + 4]; dy[i__ + 5] = dx[i__ + 5]; dy[i__ + 6] = dx[i__ + 6]; /* L50: */ } return 0; } /* dcopy_ */ /* ====================== The end of dcopy =============================== */ doublereal ddot_(n, dx, incx, dy, incy) integer *n; doublereal *dx; integer *incx; doublereal *dy; integer *incy; { /* System generated locals */ integer i__1; doublereal ret_val; /* Local variables */ static integer i__, m; static doublereal dtemp; static integer ix, iy, mp1; /* forms the dot product of two vectors. */ /* uses unrolled loops for increments equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* Parameter adjustments */ --dy; --dx; /* Function Body */ ret_val = 0.; dtemp = 0.; if (*n <= 0) { return ret_val; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { dtemp += dx[ix] * dy[iy]; ix += *incx; iy += *incy; /* L10: */ } ret_val = dtemp; return ret_val; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 5; if (m == 0) { goto L40; } i__1 = m; for (i__ = 1; i__ <= i__1; ++i__) { dtemp += dx[i__] * dy[i__]; /* L30: */ } if (*n < 5) { goto L60; } L40: mp1 = m + 1; i__1 = *n; for (i__ = mp1; i__ <= i__1; i__ += 5) { dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[ i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 4] * dy[i__ + 4]; /* L50: */ } L60: ret_val = dtemp; return ret_val; } /* ddot_ */ /* ====================== The end of ddot ================================ */ /* Subroutine */ int dpofa_(a, lda, n, info) doublereal *a; integer *lda, *n, *info; { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; /* Builtin functions */ double sqrt(); /* Local variables */ static integer j, k; static doublereal s, t; static integer jm1; /* Parameter adjustments */ a_dim1 = *lda; a_offset = 1 + a_dim1 * 1; a -= a_offset; /* Function Body */ i__1 = *n; for (j = 1; j <= i__1; ++j) { *info = j; s = 0.; jm1 = j - 1; if (jm1 < 1) { goto L20; } i__2 = jm1; for (k = 1; k <= i__2; ++k) { i__3 = k - 1; t = a[k + j * a_dim1] - ddot_(&i__3, &a[k * a_dim1 + 1], &c__1, & a[j * a_dim1 + 1], &c__1); t /= a[k + k * a_dim1]; a[k + j * a_dim1] = t; s += t * t; /* L10: */ } L20: s = a[j + j * a_dim1] - s; /* ......exit */ if (s <= 0.) { goto L40; } a[j + j * a_dim1] = sqrt(s); /* L30: */ } *info = 0; L40: return 0; } /* dpofa_ */ /* ====================== The end of dpofa =============================== */ /* Subroutine */ int dscal_(n, da, dx, incx) integer *n; doublereal *da, *dx; integer *incx; { /* System generated locals */ integer i__1, i__2; /* Local variables */ static integer i__, m, nincx, mp1; /* scales a vector by a constant. */ /* uses unrolled loops for increment equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* modified 3/93 to return if incx .le. 0. */ /* Parameter adjustments */ --dx; /* Function Body */ if (*n <= 0 || *incx <= 0) { return 0; } if (*incx == 1) { goto L20; } /* code for increment not equal to 1 */ nincx = *n * *incx; i__1 = nincx; i__2 = *incx; for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { dx[i__] = *da * dx[i__]; /* L10: */ } return 0; /* code for increment equal to 1 */ /* clean-up loop */ L20: m = *n % 5; if (m == 0) { goto L40; } i__2 = m; for (i__ = 1; i__ <= i__2; ++i__) { dx[i__] = *da * dx[i__]; /* L30: */ } if (*n < 5) { return 0; } L40: mp1 = m + 1; i__2 = *n; for (i__ = mp1; i__ <= i__2; i__ += 5) { dx[i__] = *da * dx[i__]; dx[i__ + 1] = *da * dx[i__ + 1]; dx[i__ + 2] = *da * dx[i__ + 2]; dx[i__ + 3] = *da * dx[i__ + 3]; dx[i__ + 4] = *da * dx[i__ + 4]; /* L50: */ } return 0; } /* dscal_ */ /* ====================== The end of dscal =============================== */ /* Subroutine */ int dtrsl_(t, ldt, n, b, job, info) doublereal *t; integer *ldt, *n; doublereal *b; integer *job, *info; { /* System generated locals */ integer t_dim1, t_offset, i__1, i__2; /* Local variables */ static integer case__; static doublereal temp; static integer j; static integer jj; /* Parameter adjustments */ t_dim1 = *ldt; t_offset = 1 + t_dim1 * 1; t -= t_offset; --b; /* Function Body */ i__1 = *n; for (*info = 1; *info <= i__1; ++(*info)) { /* ......exit */ if (t[*info + *info * t_dim1] == 0.) { goto L150; } /* L10: */ } *info = 0; /* determine the task and go to it. */ case__ = 1; if (*job % 10 != 0) { case__ = 2; } if (*job % 100 / 10 != 0) { case__ += 2; } switch ((int)case__) { case 1: goto L20; case 2: goto L50; case 3: goto L80; case 4: goto L110; } /* solve t*x=b for t lower triangular */ L20: b[1] /= t[t_dim1 + 1]; if (*n < 2) { goto L40; } i__1 = *n; for (j = 2; j <= i__1; ++j) { temp = -b[j - 1]; i__2 = *n - j + 1; daxpy_(&i__2, &temp, &t[j + (j - 1) * t_dim1], &c__1, &b[j], &c__1); b[j] /= t[j + j * t_dim1]; /* L30: */ } L40: goto L140; /* solve t*x=b for t upper triangular. */ L50: b[*n] /= t[*n + *n * t_dim1]; if (*n < 2) { goto L70; } i__1 = *n; for (jj = 2; jj <= i__1; ++jj) { j = *n - jj + 1; temp = -b[j + 1]; daxpy_(&j, &temp, &t[(j + 1) * t_dim1 + 1], &c__1, &b[1], &c__1); b[j] /= t[j + j * t_dim1]; /* L60: */ } L70: goto L140; /* solve trans(t)*x=b for t lower triangular. */ L80: b[*n] /= t[*n + *n * t_dim1]; if (*n < 2) { goto L100; } i__1 = *n; for (jj = 2; jj <= i__1; ++jj) { j = *n - jj + 1; i__2 = jj - 1; b[j] -= ddot_(&i__2, &t[j + 1 + j * t_dim1], &c__1, &b[j + 1], &c__1); b[j] /= t[j + j * t_dim1]; /* L90: */ } L100: goto L140; /* solve trans(t)*x=b for t upper triangular. */ L110: b[1] /= t[t_dim1 + 1]; if (*n < 2) { goto L130; } i__1 = *n; for (j = 2; j <= i__1; ++j) { i__2 = j - 1; b[j] -= ddot_(&i__2, &t[j * t_dim1 + 1], &c__1, &b[1], &c__1); b[j] /= t[j + j * t_dim1]; /* L120: */ } L130: L140: L150: return 0; } /* dtrsl_ */