/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:42 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dmlc.h" #include #include #include /* PARAMETER translations */ #define IFVAL 1 #define IG 3 #define IIACT 5 #define IINFO 1 #define IITERC 2 #define INACT 4 #define INFVAL 3 #define KTNORM 2 /* end of PARAMETER translations */ void /*FUNCTION*/ dmlc01( void (*dmlcfg)(long,double[],double*,double[],LOGICAL32*), long n, long m, long meq, double *a, long lda, double b[], double xl[], double xu[], double x[], double acc, long iprint, long nfmax, long iw[], long liw, double w[], long lw) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, ibres, id, igm, igs, info, ipar, ireskt, iu, iw20, ixbig, ixs, iz, iztg; double temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; long *const Iw = &iw[0] - 1; double *const W = &w[0] - 1; double *const X = &x[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. * File: DMLC.for Minimization with linear constraints. *>> 2001-10-05 DMLC Krogh Minor format changes -- 1PE => 1P,E *>> 2001-06-18 DMLC Krogh Changed ". LT." to " .LT." *>> 1996-07-08 DMLC Krogh Multiple changes for C conversion. *>> 1996-03-30 DMLC Krogh Removed MAX from type stmt., added external. *>> 1995-11-15 DMLC Krogh Moved formats up for C conversion. *>> 1994-11-11 DMLC Krogh Declared all vars. *>> 1994-11-02 DMLC Krogh Changes to use M77CON *>> 1992-04-27 DMLC CLL *>> 1992-04-17 CLL *>> 1992-03-27 CLL *>> 1992-02-05 CLL Declared all floating point variables. *>> 1992-01-16 CLL *>> 1991-06-10 CLL & FTK Editing comments. *>> 1991-04-30 CLL & FTK *>> 1991-04-19 CLL *>> 1991-04-15 FTK & CLL Made FIRST an argument of DMLC20 *>> 1990-09-12 C. L. Lawson and F. T. Krogh, JPL *>> 1990-07-25 C. L. Lawson, JPL *>> 1990-07-12 C. L. Lawson, JPL *>> 1990-04-03 C. L. Lawson, JPL *--D replaces "?": ?MLC,?MLC01,?MLC02,?MLCFG,?MLC20,?MLC21,?MLC04,?MLC05 *--& ?MLC13,?MLC15,?MLC03,?MLC12,?MLC06,?MLC11,?MLC09,?MLC16,?MLC19 *--& ?MLC07,?MLC18,?MLC08,?MLC14,?MLC10,?MLC17,?ERV1 * Also calls ERMOR, ERMSG, IERM1, IERV1 * ------------------------------------------------------------------ * This algorithm and the original code is due to * M. J. D. Powell, Cambridge Univ., 1989. * 1990-04-03 This version adapted for usage at JPL by * C. L. Lawson. Instrumented with "c--" lines to automate type * conversions. * 1990-07-11 CLL Added new argument NFMAX to DMLC01. Added new * subroutine DMLC20 to compute gradient by finite differences. * Added argument HAVEG in the user-provided subr DMLCFG. * Added args XL and XU to subr _MLC09 in order to be able to pass * them on to DMLC20. * Added reference to [D/R]1MACH in _MLC15. * 1990-07-12 CLL Transposed indices in the array A(,). * The constraints are now Ax .le. b rather than (A**t)x .le. b. * 1990-07-19 CLL Reorganized the arg list of DMLC01. * Prev args INFO, NACT, IACT() are embedded in new arg IW(). * New item, FVAL, and previous args W() and PAR() are combined * into new arg W(). The new item FVAL is the final value of the * objective function. * 1990-07-24 CLL Setting INFO = 0 initially. It keeps this value * till it is changed for some reason. Previously it was set to * 4 initially. * Changed specification of IPRINT. It now contains 5 print * parameters. This gives more flexibility, such as requesting * printing only at termination. * Moved the main block of intermediate printing into a new subr, * DMLC21. This allows cleaner logic for testing when to print. * DMLC21 is called both from DMLC02 and from _MLC05. In revising * the print control we removed variables NFVALK and ITERP, and * added variables NFVPR and TOLPRT. * Printing of error messages is now done using JPL MATH77 error * printing subroutines. Requires linking of * DERV1, ERFIN, ERMOR, ERMSG, IERM1, and IERV1. The calls are in * DMLC01, DMLC02, and _MLC15. ERFIN is not called directly. * Setting INFO = 9 in DMLC02 when terminate due to solution being * determined just by the constraints. Also compute and return * the func value in this case. * 1990-09-12 CLL and FTK Added arg. W20 for _MLC20. * 1991-04-19 CLL Change to have NFVALS count every call to _MLCFG. * Move test of IPFREQ to be sure last iteration gets printed if * it is a multiple of IPFREQ. Changed usage of NFVPR. * 1992 March/April CLL Additions to the code in MINFUN (C05) and * LSRCH (C09). The code will now return INFO = 2 in some cases * in which it previously returned INFO = 3. It will also finish in * fewer iterations in some cases. Added more to error msg on quit- * ting with INFO = 3. Added HAVEG into arg list of _MLC05. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *-- Begin Mask code changes. * Names of subroutines have been changed as follows: * Original names New DP names New SP names * GETMIN DMLC01 SMLC01 * MINFLC DMLC02 SMLC02 * EQCONS DMLC03 SMLC03 * GETFES DMLC04 SMLC04 * MINFUN DMLC05 SMLC05 * CONRES DMLC06 SMLC06 * GETD DMLC07 SMLC07 * SDEGEN DMLC08 SMLC08 * LSRCH DMLC09 SMLC09 * NEWCON DMLC10 SMLC10 * SATACT DMLC11 SMLC11 * ADDCON DMLC12 SMLC12 * ADJTOL DMLC13 SMLC13 * DELCON DMLC14 SMLC14 * INITZU DMLC15 SMLC15 * KTVEC DMLC16 SMLC16 * SDIRN DMLC17 SMLC17 * STEPBD DMLC18 SMLC18 * ZBFGS DMLC19 SMLC19 * FGCALC DMLCFG SMLCFG * New code by FTK --> DMLC20 SMLC20 * New code by CLL --> DMLC21 SMLC21 */ /*-- End Mask code changes. */ /* This is main entry point to a package of subroutines that calculate * the least value of a differentiable function of several variables * subject to linear constraints on the values of the variables. * ------------------------------------------------------------------ * Subroutine Arguments * * DMLCFG [in] Name of subroutine having an interface of the form * SUBROUTINE DMLCFG (N, X, F, G, HAVEG) * INTEGER N * [DOUBLE PRECISION or REAL] X(N), F, G(N) * LOGICAL HAVEG * provided by the user to compute values of the objective function, * F, and optionally its gradient vector, G(), at the argument value * given by X(). The user may choose either to write code to compute * G() or not. If DMLCFG computes G() it must also set HAVEG = .true * If DMLCFG does not compute G() it must set HAVEG = .false. * In this latter case the package will make additional calls to * DMLCFG to compute a finite difference approximation to the * gradient vector. * The vectors X() and G() will be of length N. The subroutine must * not change the values of N or X(). * DMLCFG may use COMMON block(s) if neccessary to communicate * data between the user's main program and DMLCFG. * N is the number of variables and must be set by the user. * M is the number of linear constraints (excluding simple bounds) and * must be set by the user. * MEQ is the number of constraints that are equalities and must be set * by the user. * A(.,.) is a 2-dimensional array whose rows are the gradients of * the M constraint functions. Its entries must be set by the user * and its dimensions must be at least M by N. * LDA is the actual first dimension of the array A that is supplied by * the user, so its value may not be less than M. * B(.) is a vector of constraint right hand sides that must also be set * by the user. Specifically the constraints on the variables X(I) * I=1(1)N are * A(K,1)*X(1)+...+A(K,N)*X(N) .EQ. B(K) K=1,...,MEQ * A(K,1)*X(1)+...+A(K,N)*X(N) .LE. B(K) K=MEQ+1,...,M . * Note that the data that define the equality constraints come * before the data of the inequalities. * XL(.) and XU(.) are vectors whose components must be set to lower and * upper bounds on the variables. Choose very large negative and * positive entries if a component should be unconstrained, or set * XL(I)=XU(I) to freeze the I-th variable. Specifically these * simple bounds are * XL(I) .LE. X(I) and X(I) .LE. XU(I) I=1,...,N . * X(.) is the vector of variables of the optimization calculation. Its * initial elements must be set by the user to an estimate of the * required solution. The subroutines can usually cope with poor * estimates, and there is no need for X(.) to be feasible initially. * These variables are adjusted automatically and the values that * give the least feasible calculated value of the objective function * are available in X(.) on the return from DMLC01. * ACC is a tolerance on the first order conditions at the calculated * solution of the optimization problem. These first order * conditions state that, if X(.) is a solution, then there is a set * of active constraints with indices IACT(K) K=1(1)NACT, say, such * that X(.) is on the boundaries of these constraints, and the * gradient of the objective function can be expressed in the form * GRAD(F)=PAR(1)*GRAD(C(IACT(1)))+... * ...+PAR(NACT)*GRAD(C(IACT(NACT))) . * Here PAR(K) K=1(1)NACT are Lagrange multipliers that are * nonpositive * for inequality constraints, and GRAD(C(IACT(K))) is the gradient * of the IACT(K)-th constraint function, so it is A(IACT(K),.) if * IACT(K) .LE. M, and it is minus or plus the J-th coordinate vector * if the constraint is the lower or upper bound on X(J) * respectively. The normal return from the calculation occurs when * X(.) is feasible and the sum of squares of components of the * vector RESKT(.) is at most ACC**2, where RESKT(.) is the * N-component vector of residuals of the first order condition that * is displayed above. Sometimes the package cannot satisfy this * condition, because noise in the function values can prevent a * change to the variables, no line search being allowed to increase * the objective function. * IPRINT [in, integer] Contains print control parameters, packed as * follows: * IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST * Printing only begins after a first feasible point is found. * The items to be printed are selected by IPMORE. The other * parameters determine when to print. * IPFREQ is zero or positive. If positive, printing will be done on * iterations 0, IPFREQ, 2*IPFREQ, etc. * IPTOL = 0 or 1. 1 means to print the new TOL value * and the standard items each time TOL is changed. * IPFRST = 0 or 1. 1 means to print on the first iteration, i.e. on * iteration No. 0. * IPMORE = 0 or 1. 0 means the items to be printed are ITERC, * NFVALS, F, X(1:N), and G(1:N). 1 means to also print * IACT(1:NAXT), PAR(1:NACT) and RESKT(1:N). * IPLAST = 0 or 1. 1 means to print the final results, and the * reason for termination. * NFMAX [in, integer] If positive this sets an upper limit on the * number of function evaluations. When gradients are estimated by * differences, the count includes the function evaluations for this * purpose also. If zero there is no upper limit. * IW(.) [work, out, integer] Integer work space. Also used to return */ /* status information. Its length must be at least LIW. * On return: * IW(1) contains INFO. Indicates reason for termination. * See below for explanation of values. * IW(2) contains ITERC. No. of iterations. * IW(3) contains NFVALS. No. of function evaluations, not * counting extra func evals done to estimate the gradient * when the gradient is not computed explicitly. In this * latter case the actual No. of func evals will be * (K+1)*NFVALS, where K is the number of solution * components whose lower and upper bounds are not equal. * IW(4) contains NACT. The number of active constraints at the * solution point. Will be in the interval [0, N]. * IW(5:4+NACT) contains IACT(1:NACT). IACT() is used as work * space of length M + 2*N. On return the first NACT * locations contain the indices of the constraints that * are active at the final point. Indices [1:M] refer to * rows of the system Ax .le. b. Indices [M+1:M+N] * refer to component lower bounds. Indices [M+N+1:M+2*N] * refer to component upper bounds. * * LIW [in, integer] Dimension of the IW() array. Require * LIW .ge. 4 + M + 2*N. * * W(.) [work, out, float] Working space array of float variables. * Also used to return results. * Its length must be at least LW. * On return: * W(1) contains FVAL, the final value of the objective function. * W(2:1+N) contains the final gradient vector, GRAD(1:N). * W(2+N:1+2*N) contains the Kuhn-Tucker residual vector, RESKT(1:N). * W(2+2*N:1+2*N+NACT) contains the Lagrange multipliers, PAR(1:NACT) * where NACT has a value in [0,N]. GRAD(), RESKT(), and PAR() * are mentioned above in the description of ACC. * LW [in, integer] Dimension of the W() array. Require * LW .ge. 3 + M + 16*N + N**2. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Return values of INFO * * INFO indicates the reason for termination. * * INFO=1 X(.) is feasible and the condition that depends on * ACC is satisfied. * INFO=2 X(.) is feasible and rounding errors are preventing * further progress. * INFO=3 X(.) is feasible but the objective function fails to * decrease although a decrease is predicted by the current gradient * vector. If this return occurs and RESKT(.) has large components * then the user's calculation of the gradient of the objective * function may be incorrect. One should also question the coding of * the gradient when the final rate of convergence is slow. * INFO=4 In this case the calculation cannot begin because IA * is less than N or because the lower bound on a variable is greater * than the upper bound. * INFO=5 This value indicates that the equality constraints * are inconsistent. These constraints include any components of * X(.) that are frozen by setting XL(I)=XU(I). * INFO=6 In this case there is an error return because the * equality constraints and the bounds on the variables are found to * be inconsistent. * INFO=7 This value indicates that there is no vector of * variables that satisfies all of the constraints. Specifically, * when this return or an INFO=6 return occurs, the current active * constraints (whose indices are IACT(K) K=1(1)NACT) prevent any * change in X(.) that reduces the sum of constraint violations, * where only bounds are included in this sum if INFO=6. * INFO=8 In this case the limit on the number of calls of * subroutine DMLCFG (see below) has been reached, and there would * have been further calculation otherwise. * INFO = 9 The solution is determined just by the constraints. * In this case PAR() and RESKT() are not defined on return. * FVAL will be defined on return. G() will only be defined on * return if DMLCFG returns HAVEG = .true. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Note 1. The variables N, M, MEQ, LDA, ACC, IPRINT, NFMAX and the * elements of the arrays A(.,.), B(.), XL(.) and XU(.) are not * altered by the optimization procedure. Their values, and the * initial components of X(.) must be set on entry to DMLC01. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Note 2. A paper on the method of calculation and a report on the * main features of the computer code are available from the author * M.J.D.Powell (D.A.M.T.P., University of Cambridge, Silver Street, * Cambridge CB3 9EW, England). * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Check sizes of LIW and LW */ info = 0; if (liw < 4 + m + 2*n) { info = 4; ierm1( "DMLC01", info, 0, "Dimension LIW is too small.", "LIW" , liw, ',' ); ierv1( "Need", 4 + m + 2*n, '.' ); } if (lw < 3 + m + n*(16 + n)) { info = 4; ierm1( "DMLC01", info, 0, "Dimension LW is too small.", "LW" , lw, ',' ); ierv1( "Need", 3 + m + n*(16 + n), '.' ); } if (info != 0) { Iw[1] = info; return; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Partition the array W() * parameter(IFVAL = 1, KTNORM = 2, IG = 3) */ ireskt = IG + n; ipar = ireskt + n; iz = ipar + n; iu = iz + n*n; ixbig = iu + n; ibres = ixbig + n; id = ibres + m + n + n; iztg = id + n; igm = iztg + n; ixs = igm + n; igs = ixs + n; /* Need 4*N+1 locations beginning at IW20. */ iw20 = igs + n; /* Partition the array IW() * parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5) * * Zero the RESKT vector so it will be safe to compute its norm on * return from DMLC02 even in cases when DMLC02 does not compute * RESKT. * */ for (i = ireskt; i <= (ireskt + n - 1); i++) { W[i] = 0.0e0; } /* Call the optimization package. * */ dmlc02( dmlcfg, n, m, meq, a, lda, b, xl, xu, x, acc, &Iw[IIACT], &Iw[INACT], &W[ipar], iprint, nfmax, &Iw[IINFO], &Iw[IITERC], &Iw[INFVAL], &W[IG], &W[iz], &W[iu], &W[ixbig], &W[ireskt], &W[KTNORM], &W[ibres], &W[id], &W[iztg], &W[igm], &W[ixs], &W[igs], &W[IFVAL], &W[iw20] ); temp = 0.0e0; for (i = ireskt; i <= (ireskt + n - 1); i++) { temp += SQ(W[i]); } W[KTNORM] = sqrt( temp ); return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc02( void (*dmlcfg)(long,double[],double*,double[],LOGICAL32*), long n, long m, long meq, double *a, long lda, double b[], double xl[], double xu[], double x[], double acc, long iact[], long *nact, double par[], long iprint, long nfmax, long *info, long *iterc, long *nfvals, double g[], double z[], double u[], double xbig[], double reskt[], double *enrmkt, double bres[], double d[], double ztg[], double gm[], double xs[], double gs[], double *fval, double w20[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) LOGICAL32 haveg; long int gmode, i, ipfreq, ipfrst, iplast, ipmore, iptol, k, meql, mp, msat, mtot, nfvpr; double relacc, ssqkt, tol, zznorm; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; double *const Bres = &bres[0] - 1; double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gm = &gm[0] - 1; double *const Gs = &gs[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Reskt = &reskt[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xs = &xs[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* Main subroutine to control the Powell constrained minimization * package. * 1990-07-19 CLL Added NFMAX, ITERC, NFVALS, and FVAL to arg list. * Also added FVAL to arg list of DMLC05 that is called * from this subr. * NFMAX provides limit on no. of func evaluations, but zero means * no limit. This was previously overloaded onto INFO. * ITERC & NFVALS give output of counts of iterations and * function evaluations. * FVAL provides output of the final objective function value. * ------------------------------------------------------------------ * Arguments * * DMLCFG, N,M,MEQ,A(),LDA,B(),XL(),XU(),X(), * ACC,IACT(),NACT,PAR(), * IPRINT [in, integer] Contains print parameters packed as follows: * IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST * NFMAX [in, integer] Limit on No. of function evaluations. But zero * means no limit. * INFO [out, integer] Termination status flag. * ITERC [out, integer] Count of No. of iterations. Initialized to 0 * in this subr. Incremented in DMLC05. * NFVALS [out, integer] Count of number of function evaluations. * G(),Z(),U(),XBIG(), * RESKT() * ENRMKT Euclidean norm of RESKT(). * BRES(),D(), * ZTG(),GM(),XS(),GS(), * FVAL [out, float] Final (best) value of F. * W20 [out, work, float] Working storage for DMLC20. * ------------------------------------------------------------------ * Descriptions of some of the internal variables. * * GMODE [integer] Initialized to 0 here. Sent to DMLC05 for eventual * use by DMLC20. * HAVEG [logical] Arg in call to DMLCFG, not used in this subroutine. * NFVPR [integer] Used to avoid printing results from the same * function evaluation more than once. Initialized to -1 in this * subr. Modified in DMLC05. Used in DMLC02 and DNLC05. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Initialize ZZNORM, ITERC, NFVPR, NFVALS, INFO * */ zznorm = -1.0e0; w20[0] = -1.0e0; w20[1] = -1.0e0; gmode = 0; *iterc = 0; nfvpr = -1; *nfvals = 0; *info = 1; /* Decompose IPRINT into 5 integers: * IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST * */ ipfreq = iprint/100; i = iprint - ipfreq*100; iptol = i/8; i += -iptol*8; ipfrst = i/4; i += -ipfrst*4; ipmore = i/2; iplast = i - ipmore*2; /* Check the bounds on N, M and MEQ. * */ if (((n < 1 || m < 0) || meq < 0) || m < meq) { *info = 4; ierm1( "DMLC02", *info, 0, "Bad value for M, N, or MEQ.", "M", m, ',' ); ierv1( "N", n, ',' ); ierv1( "MEQ", meq, '.' ); goto L_40; } /* Initialize RELACC, Z, U and TOL. * */ dmlc15( n, m, xl, xu, x, iact, &meql, info, z, u, xbig, &relacc ); if (*info == 4) goto L_40; tol = fmax( 0.01e0, 10.0e0*relacc ); if (((ipfreq != 0 || ipfrst != 0) || iplast != 0) || iptol != 0) { printf("\n Beginning subroutine DMLC02, called from DMLC01:\n N =%4ld, M =%4ld, MEQ =%4ld, ACC =%11.3g\n RELACC =%11.3g, TOL = %11.3g\n", n, m, meq, acc, relacc, tol); } /* Add any equality constraints to the active set. * */ if (meq > 0) { dmlc03( n, m, meq, a, lda, b, xu, iact, &meql, info, z, u, relacc, xs, gs ); if (*info == 5) { ermsg( "DMLC02", *info, 0, "Equality constraints are inconsistent." , '.' ); goto L_40; } } *nact = meql; msat = meql; /* Add the bounds to the list of constraints. * */ mtot = *nact; for (i = 1; i <= n; i++) { if (Xl[i] < Xu[i]) { mtot += 2; Iact[mtot - 1] = m + i; Iact[mtot] = m + n + i; } } /* Try to satisfy the bound constraints. * */ dmlc04( n, m, a, lda, b, xl, xu, x, iact, nact, par, info, g, z, u, xbig, relacc, &tol, meql, &msat, mtot, bres, d, ztg, gm, reskt, xs, gs, w20 ); if (msat < mtot) { *info = 6; ermsg( "DMLC02", *info, 0, "Equalities and bounds are inconsistent." , '.' ); goto L_40; } /* Add the ordinary inequalities to the list of constraints. * */ if (m > meq) { mp = meq + 1; for (k = mp; k <= m; k++) { mtot += 1; Iact[mtot] = k; } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ L_30: ; /* Begin main loop * * Correct any constraint violations. * */ dmlc04( n, m, a, lda, b, xl, xu, x, iact, nact, par, info, g, z, u, xbig, relacc, &tol, meql, &msat, mtot, bres, d, ztg, gm, reskt, xs, gs, w20 ); if (msat < mtot) { *info = 7; ermsg( "DMLC02", *info, 0, "Constraints are inconsistent." , '.' ); goto L_40; } else if (meql == n) { /* Here the soln is determined just by the constraints. * PAR() and RESKT() are not defined on return in this case. * FVAL has not yet been evaluated. We evaluate FVAL here in * order to return its value. We leave G() not computed if HAVEG * is false, however it would be easy to compute G by calling * DMLC20 if that seemed preferable. * */ (*dmlcfg)( n, x, fval, g, &haveg ); *nfvals += 1; /* call DMLCFG(N, X, FVAL, G, HAVEG) */ *enrmkt = 0.0e0; *info = 9; if (((ipfreq != 0 || ipfrst != 0) || iplast != 0) || iptol != 0) { printf("\n [1] DMLC02..\n The solution is determined by the equality constraints.\n"); dmlc21( ipmore, FALSE, n, *iterc, *nfvals, *fval, tol, x, g, *nact, iact, par, reskt, *enrmkt ); } goto L_40; } /* Minimize the objective function in the case when constraints are * treated as degenerate if their residuals are less than TOL. * * DMLC05 is MINFUN */ dmlc05( dmlcfg, n, m, a, lda, b, xl, xu, x, acc, iact, nact, par, ipfreq, iptol, ipfrst, ipmore, info, g, z, u, xbig, relacc, &zznorm, tol, meql, mtot, iterc, &nfvpr, nfvals, nfmax, reskt, &ssqkt, bres, d, ztg, gm, xs, gs, fval, &gmode, &haveg, w20 ); /* Reduce TOL if necessary. * */ if (tol > relacc && *nact > 0) { if (nfmax <= 0 || *nfvals < nfmax) { /* DMLC13 is ADJTOL */ dmlc13( n, m, a, lda, b, xl, xu, x, iact, *nact, xbig, relacc, &tol, meql ); goto L_30; } else { *info = 8; } } /* Here when INFO = 1, 2, 3, or 8. * We treat 1 and 2 as normal conditions and 3 and 8 as errors. * */ if (*info == 1) { if (iplast != 0) { printf("[1] DMLC01 quitting because requested accuracy has been attained.\n"); } } else if (*info == 2) { if (iplast != 0) { printf("[2] DMLC01 quitting due to limitation of computational precision.\n"); } } else if (*info == 3) { ermsg( "DMLC02", *info, 0, "F not decreasing, although gradient indicates it should." , ',' ); if (haveg) { ermor( "Could be due to limitation of precision", ',' ); ermor( "or to incorrect code for the gradient.", '.' ); } else { ermor( "Could be due to limitation of precision.", '.' ); } } else if (*info == 8) { ierm1( "DMLC02", *info, 0, "Reached specified max number of function evaluations." , "Count", *nfvals, ',' ); ierv1( "NFMAX", nfmax, '.' ); } /* Test to print results before quitting. * */ *enrmkt = sqrt( ssqkt ); if (iplast != 0 && *nfvals != nfvpr) { /* print*,'DMLC02.. Debug.. Call C21 on leaving C02.' */ dmlc21( ipmore, TRUE, n, *iterc, *nfvals, *fval, tol, x, g, *nact, iact, par, reskt, *enrmkt ); } L_40: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc03( long n, long m, long meq, double *a, long lda, double b[], double xu[], long iact[], long *meql, long *info, double z[], double u[], double relacc, double am[], double cgrad[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, iz, j, jm, k, keq, np; double rhs, sum, sumabs, vmult; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Am = &am[0] - 1; double *const B = &b[0] - 1; double *const Cgrad = &cgrad[0] - 1; long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* Try to add the next equality constraint to the active set. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ for (keq = 1; keq <= meq; keq++) { if (*meql < n) { np = *meql + 1; Iact[np] = keq; dmlc12( n, m, a, lda, iact, meql, z, u, relacc, np, am, cgrad ); if (*meql == np) goto L_50; } /* If linear dependence occurs then find the multipliers of the * dependence relation and apply them to the right hand sides. * */ sum = B[keq]; sumabs = fabs( B[keq] ); if (*meql > 0) { for (i = 1; i <= n; i++) { Am[i] = A(i - 1,keq - 1); } k = *meql; L_20: vmult = 0.0e0; iz = k; for (i = 1; i <= n; i++) { vmult += Z[iz]*Am[i]; iz += n; } vmult *= U[k]; j = Iact[k]; if (j <= m) { for (i = 1; i <= n; i++) { Am[i] += -vmult*A(i - 1,j - 1); } rhs = B[j]; } else { jm = j - m - n; Am[jm] -= vmult; rhs = Xu[jm]; } sum += -rhs*vmult; sumabs += fabs( rhs*vmult ); k -= 1; if (k >= 1) goto L_20; } /* Error return if the constraints are inconsistent. * */ if (fabs( sum ) > relacc*sumabs) { *info = 5; goto L_60; } L_50: ; } L_60: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc04( long n, long m, double *a, long lda, double b[], double xl[], double xu[], double x[], long iact[], long *nact, double par[], long *info, double g[], double z[], double u[], double xbig[], double relacc, double *tol, long meql, long *msat, long mtot, double bres[], double d[], double ztg[], double gm[], double gmnew[], double parnew[], double cgrad[], double w20[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, indxbd, itest, msatk; double stepcb, sumres, sumrsk; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; double *const Bres = &bres[0] - 1; double *const Cgrad = &cgrad[0] - 1; double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gm = &gm[0] - 1; double *const Gmnew = &gmnew[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Parnew = &parnew[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* Make the correction to X for the active constraints. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ *info = 0; L_10: dmlc11( n, m, a, lda, b, xl, xu, x, iact, nact, info, z, u, xbig, relacc, *tol, meql ); if (*info > 0) *msat = *nact; if (*msat == mtot) goto L_60; /* Try to correct the infeasibility. * */ L_20: msatk = *msat; sumrsk = 0.0e0; L_30: dmlc06( n, m, a, lda, b, xl, xu, x, iact, nact, par, g, z, u, xbig, bres, d, ztg, relacc, *tol, &stepcb, &sumres, meql, msat, mtot, &indxbd, gm, gmnew, parnew, cgrad, w20 ); /* Include the new constraint in the active set. * */ if (stepcb > 0.0e0) { for (i = 1; i <= n; i++) { X[i] += stepcb*D[i]; Xbig[i] = fmax( Xbig[i], fabs( X[i] ) ); } dmlc12( n, m, a, lda, iact, nact, z, u, relacc, indxbd, gmnew, cgrad ); } /* Test whether to continue the search for feasibility. * */ if (*msat < mtot) { if (stepcb == 0.0e0) goto L_50; if (msatk < *msat) goto L_20; if (sumrsk == 0.0e0 || sumres < sumrsk) { sumrsk = sumres; itest = 0; } itest += 1; if (itest <= 2) goto L_30; /* Reduce TOL if it may be too large to allow feasibility. * */ L_50: if (*tol > relacc) { dmlc13( n, m, a, lda, b, xl, xu, x, iact, *nact, xbig, relacc, tol, meql ); goto L_10; } } L_60: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc05( void (*dmlcfg)(long,double[],double*,double[],LOGICAL32*), long n, long m, double *a, long lda, double b[], double xl[], double xu[], double x[], double acc, long iact[], long *nact, double par[], long ipfreq, long iptol, long ipfrst, long ipmore, long *info, double g[], double z[], double u[], double xbig[], double relacc, double *zznorm, double tol, long meql, long mtot, long *iterc, long *nfvpr, long *nfvals, long nfmax, double reskt[], double *ssqkt, double bres[], double d[], double ztg[], double gm[], double xs[], double gs[], double *fval, long *gmode, LOGICAL32 *haveg, double w20[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) LOGICAL32 doprnt; long int _l0, i, indxbd, iterk, k, msat; double ddotg, diff, fprev, relaxf, step, stepcb, sum; static double f; static double eps = 0.0e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; double *const Bres = &bres[0] - 1; double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gm = &gm[0] - 1; double *const Gs = &gs[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Reskt = &reskt[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xs = &xs[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* MINFUN: Solve the minimization problem using the current value of * TOL. * 1990-07-19 CLL Added FVAL to argument list to return final value * of F. * 1990-07-23 CLL Removed IPRINT from arg list and added IPFREQ, * IPTOL, IPFRST, IPMORE, and NFVPR. NFVPR keeps a record of * last function evaluation which printing has been done. * 1990-07-23 CLL Added logical variable TOLPRT. * ------------------------------------------------------------------ * Arguments * * DMLCFG, N,M,A(),LDA,B(),XL(),XU(),X(),ACC,IACT(),NACT, * PAR() * IPFREQ [in, integer] Zero or positive. If positive, printing of * intermediate results using DMLC21 will be done on iterations * 0, IPFREQ, 2*IPFREQ, etc. * IPTOL [in, integer] = 0 or 1. 1 means to print the new TOL value * and items printed by DMLC21 each time this subr is entered. * IPFRST [in, integer] = 0 or 1. 1 means to print using DMLC21 on * the first iteration, i.e. iteration No. 0. * IPMORE [in, integer] = 0 or 1. Will be passed to DMLC21 to cause * printing to include more items when = 1 than when = 0. * INFO [inout, integer] Status flag. * G(),Z(),U(),XBIG(),RELACC,ZZNORM,TOL,MEQL,MTOT, * NFVALS [inout, integer] Count of calls for function evaluation. * Only counts calls to DMLCFG when the algorithm directly needs a * function value, not calls to DMLCFG that occur in DMLC20 for * estimation of the gradient. * NFMAX [in, integer] Limit on No. of function evaluations, but zero * means no limit. * RESKT() * SSQKT Sum of Squares of elts of RESKT(). * BRES(),D(),ZTG(),GM(),XS(),GS(), * FVAL [out, float] On return set to current (best) value of F. * GMODE [inout, integer] For use by DMLC20. * HAVEG [out, logical] Primarily an internal variable in this subr, * but returned so calling subr can use it to select proper error msg * when INFO = 3. * Set true by DMLCFG if DMLCFG provides the gradient * vector and false otherwise. If false, this subr calls DMLC20 to * compute a finite difference approximation for the gradient. * W20 [out, work, float] Working storage for DMLC20. * ------------------------------------------------------------------ * Description of some of the internal variables. * * DOPRNT [logical] When IPTOL = 1 this subr prints TOL on entry and * sets DOPRNT = true. Also if IPFRST = 1 and ITERC = 0 this subr * sets DOPRNT = true. * Having DOPRNT = true causes subsequent printing of intermediate * results for the current iteration, after PAR() and RESKT() have * been computed. * HAVEG [logical] Set true by DMLCFG if DMLCFG provides the gradient * vector and false otherwise. If false, this subr calls DMLC20 to * compute a finite difference approximation for the gradient. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Initialize the minimization calculation. * * print'(/''DMLC05.. Debug.. On entry, INFO = '',i6,'', ITERC ='', * i6)', INFO, ITERC */ if (eps == 0.0e0) { eps = DBL_EPSILON; } msat = mtot; iterk = *iterc; if (*nfvals == 0 || *info == 1) { /*++ CODE for ~.C. is inactive * CALL DMLCFG (N,X,F,G, HAVEG) *++ CODE for .C. is active */ (*dmlcfg)( n, x, &f, g, haveg ); *haveg = *haveg; /*++ END * print'('' DMLC05 798.. X()='',2g23.15/(19x,2g23.15))', * * (X(I),I=1,N) * print'('' .. F ='',g23.15)', F * print'('' ............ G()='',2g23.15/(19x,2g23.15))', * * (G(I),I=1,N) */ *nfvals += 1; if (!*haveg) dmlc20( dmlcfg, n, x, f, g, xl, xu, gmode, nfvals, w20 ); } fprev = fabs( f + f + 1.0e0 ); /* Setting DIFF here is not needed for the * algorithm, but is convenient when doing * debug printing. CLL 3/26/92 */ diff = 1.0e0; if (iptol != 0) { printf("\n New value of TOL = %13.5g\n", tol); doprnt = TRUE; } else if (ipfrst != 0 && *iterc == 0) { doprnt = TRUE; } else { doprnt = FALSE; } /* Calculate the next search direction. * DMLC06 is CONRES */ L_10: dmlc06( n, m, a, lda, b, xl, xu, x, iact, nact, par, g, z, u, xbig, bres, d, ztg, relacc, tol, &stepcb, &ddotg, meql, &msat, mtot, &indxbd, gm, reskt, xs, gs, w20 ); /* Calculate the Kuhn Tucker residual vector. * DMLC16 is KTVEC */ dmlc16( n, m, a, lda, iact, *nact, par, g, reskt, z, u, bres, &relaxf, meql, ssqkt, xs, gs ); /* Test for printing results every IPFREQ iterations. * */ if (ipfreq != 0) { if ((*iterc%ipfreq) == 0) doprnt = TRUE; } if (doprnt) { dmlc21( ipmore, TRUE, n, *iterc, *nfvals, f, tol, x, g, *nact, iact, par, reskt, sqrt( *ssqkt ) ); *nfvpr = *nfvals; doprnt = FALSE; } /* Test for convergence. * * print'(/'' DMLC05 849..''/'' DMLC05 849..'',3g13.5/ * * '' F, FPREV, FPREV-F='',3g13.5/ * * '' DIFF='',g13.5/ * * '' TOL, RELACC, NACT='',2g13.5,i5)', * * ACC**2, SSQKT, DDOTG, F, FPREV, FPREV-F, DIFF, * * TOL, RELACC, NACT */ if (*ssqkt <= acc*acc) { *info = 1; goto L_70; } if (ddotg >= 0.0e0) { *info = 2; goto L_70; } /* Test for termination due to no decrease in F. * */ if (f >= fprev) { if (tol == relacc || *nact == 0) { if (diff > 0.0e0) goto L_20; } /* print'(/'' DMLC05 869.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F */ if (fabs( ddotg ) <= eps*fabs( f )) { *info = 2; } else { *info = 3; } goto L_70; } L_20: ; diff = fprev - f; fprev = f; /* Test that more calls of DMLCFG are allowed. * */ if (nfmax > 0 && *nfvals >= nfmax) { *info = 8; goto L_70; } /* Test whether to reduce TOL. * */ if ((tol > relacc && *iterc > iterk) && 0.1e0*relaxf >= fmax( diff, -0.5e0*ddotg )) { /* print'('' DMLC05.. Debug.. return to reduce TOL.''/ * * '' ........ TOL='',g13.5,'', RELACC='',g13.5,'', INFO='', * * i5)', TOL, RELACC, INFO */ goto L_70; } /* Calculate the step along the search direction. * */ *iterc += 1; /* print'(/'' DMLC05 895.. ITERC = '',i6)', ITERC * DMLC09 is LSRCH */ dmlc09( dmlcfg, n, x, g, d, xs, gs, relacc, stepcb, ddotg, &f, &step, nfvals, nfmax, bres, xl, xu, gmode, w20 ); if (step == 0.0e0) { /* print'(/'' DMLC05 918.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F */ *info = 3; if (fabs( ddotg ) <= eps*fabs( f )) { *info = 2; } else { sum = 0.0e0; for (i = 1; i <= n; i++) { sum += fabs( D[i]*Gs[i] ); } /* print'(/'' DMLC05 890.. DDOTG, SUM='',2g13.5/ * * '' RELACC*SUM, DDOTG+RELACC*SUM='',2g13.5)', * * DDOTG, SUM, RELACC*SUM, DDOTG+RELACC*SUM */ if (ddotg + relacc*sum >= 0.0e0) *info = 2; } goto L_70; } /* Revise XBIG. * */ for (i = 1; i <= n; i++) { Xbig[i] = fmax( Xbig[i], fabs( X[i] ) ); } /* Revise the second derivative approximation. * */ dmlc19( n, x, *nact, g, z, ztg, xs, gs, zznorm ); /* Add a constraint to the active set if it restricts the step. * */ if (step == stepcb) { k = Iact[indxbd]; if (k > m) { k -= m; if (k <= n) { X[k] = Xl[k]; } else { X[k - n] = Xu[k - n]; } } dmlc12( n, m, a, lda, iact, nact, z, u, relacc, indxbd, xs, gs ); } goto L_10; L_70: ; *fval = f; return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc06( long n, long m, double *a, long lda, double b[], double xl[], double xu[], double x[], long iact[], long *nact, double par[], double g[], double z[], double u[], double xbig[], double bres[], double d[], double ztg[], double relacc, double tol, double *stepcb, double *sumres, long meql, long *msat, long mtot, long *indxbd, double gm[], double gmnew[], double parnew[], double cgrad[], double w20[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, idiff, j, jm, k, kl, mdeg, msatk; double ddotg, res, resabs, sum, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; double *const Bres = &bres[0] - 1; double *const Cgrad = &cgrad[0] - 1; double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gm = &gm[0] - 1; double *const Gmnew = &gmnew[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Parnew = &parnew[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* CONRES * Calculate and partition the residuals of the inactive constraints, * and set the gradient vector when seeking feasibility. * W20 has length 4*N+1. Only locations 1 to N are used here. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ idiff = mtot - *msat; if (idiff > 0) { for (i = 1; i <= n; i++) { G[i] = 0.0e0; } *sumres = 0.0e0; } msatk = *msat; mdeg = *nact; *msat = *nact; kl = meql + 1; for (k = kl; k <= mtot; k++) { j = Iact[k]; /* Calculate the residual of the current constraint. * */ if (j <= m) { res = B[j]; resabs = fabs( B[j] ); for (i = 1; i <= n; i++) { res += -X[i]*A(i - 1,j - 1); resabs += fabs( Xbig[i]*A(i - 1,j - 1) ); } } else { jm = j - m; if (jm <= n) { res = X[jm] - Xl[jm]; resabs = fabs( Xbig[jm] ) + fabs( Xl[jm] ); } else { jm -= n; res = Xu[jm] - X[jm]; resabs = fabs( Xbig[jm] ) + fabs( Xu[jm] ); } } Bres[j] = res; /* Set TEMP to the relative residual. * */ temp = 0.0e0; if (resabs != 0.0e0) temp = res/resabs; if (k > msatk && temp < 0.0e0) { if (temp + relacc >= 0.0e0) { if (j <= m) { sum = fabs( B[j] ); for (i = 1; i <= n; i++) { sum += fabs( X[i]*A(i - 1,j - 1) ); } } else { jm = j - m; if (jm <= n) { sum = fabs( X[jm] ) + fabs( Xl[jm] ); } else { sum = fabs( X[jm - n] ) + fabs( Xu[jm - n] ); } } if (fabs( res ) <= sum*relacc) temp = 0.0e0; } } /* Place the residual in the appropriate position. * */ if (k <= *nact) goto L_50; if (k <= msatk || temp >= 0.0e0) { *msat += 1; if (*msat < k) { Iact[k] = Iact[*msat]; } if (temp > tol) { Iact[*msat] = j; } else { mdeg += 1; Iact[*msat] = Iact[mdeg]; Iact[mdeg] = j; } /* Update the gradient and SUMRES if the constraint is violated when * seeking feasibility. * */ } else { if (j <= m) { for (i = 1; i <= n; i++) { G[i] += A(i - 1,j - 1); } } else { j -= m; if (j <= n) { G[j] -= 1.0e0; } else { G[j - n] += 1.0e0; } } *sumres += fabs( res ); } L_50: ; } /* Seek the next search direction unless DMLC06 was called from * DMLC04 [GETFES] and feasibility has been achieved. * */ *stepcb = 0.0e0; if (idiff > 0 && *msat == mtot) goto L_60; /* GETD */ dmlc07( n, m, a, lda, iact, nact, par, g, z, u, d, ztg, relacc, &ddotg, meql, mdeg, gm, gmnew, parnew, cgrad, w20 ); /* Calculate the (bound on the) step-length due to the constraints. * */ if (ddotg < 0.0e0) { dmlc18( n, m, a, lda, iact, bres, d, stepcb, &ddotg, mdeg, msat, mtot, indxbd ); } if (idiff == 0) *sumres = ddotg; L_60: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc07( long n, long m, double *a, long lda, long iact[], long *nact, double par[], double g[], double z[], double u[], double d[], double ztg[], double relacc, double *ddotg, long meql, long mdeg, double gm[], double gmnew[], double parnew[], double cgrad[], double w20[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, iz, j, jm, k; double ddotgm, size, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cgrad = &cgrad[0] - 1; double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gm = &gm[0] - 1; double *const Gmnew = &gmnew[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Parnew = &parnew[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* GETD * Initialize GM and cycle backwards through the active set. * W20 has length 4*N+1. Only locations 0 to N are used here. * If W20(1) is nonnegative it means [D/S]MLC20 has been called and * has stored into W20(1:N) estimates of the errors in the gradient * values stored in G(1:N). * If W20(1) is negative the user is computing the gradients and * W20(1:N) is not being otherwise used. * This subr stores a value into W20(0). This will only be used if * [D/S]MLC20 is called. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ L_10: for (i = 1; i <= n; i++) { Gm[i] = G[i]; } k = *nact; L_30: if (k > 0) { /* Set TEMP to the next multiplier, but reduce the active set if * TEMP has an unacceptable sign. * */ temp = 0.0e0; iz = k; for (i = 1; i <= n; i++) { temp += Z[iz]*Gm[i]; iz += n; } temp *= U[k]; if (k > meql && temp > 0.0e0) { dmlc14( n, m, a, lda, iact, nact, z, u, relacc, k ); goto L_10; } /* Update GM using the multiplier that has just been calculated. * */ j = Iact[k]; if (j <= m) { for (i = 1; i <= n; i++) { Gm[i] += -temp*A(i - 1,j - 1); } } else { jm = j - m; if (jm <= n) { Gm[jm] += temp; } else { Gm[jm - n] -= temp; } } Par[k] = temp; k -= 1; goto L_30; } /* Calculate the search direction and DDOTG. * */ *ddotg = 0.0e0; if (*nact < n) { /* SDEGEN */ dmlc08( n, m, a, lda, iact, nact, par, z, u, d, ztg, gm, relacc, &ddotgm, meql, mdeg, gmnew, parnew, cgrad ); if (ddotgm < 0.0e0) { size = 0.0e0; for (i = 1; i <= n; i++) { temp = D[i]*G[i]; *ddotg += temp; if (w20[1] >= 0.0e0) { size += fabs( D[i]*w20[i] ); } else { size += fabs( temp ); } } if (w20[1] < 0.0e0) size *= relacc; if (*ddotg >= -size) { *ddotg = 0.0e0; w20[0] = 0.0e0; } else { w20[0] = size/fabs( *ddotg ); } } else { w20[0] = 0.0e0; } } return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc08( long n, long m, double *a, long lda, long iact[], long *nact, double par[], double z[], double u[], double d[], double ztg[], double gm[], double relacc, double *ddotgm, long meql, long mdeg, double gmnew[], double parnew[], double cgrad[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, idrop, itest, iz, j, jm, k, ku, mp, np; double dtest, ratio, sum, temp, theta; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cgrad = &cgrad[0] - 1; double *const D = &d[0] - 1; double *const Gm = &gm[0] - 1; double *const Gmnew = &gmnew[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Parnew = &parnew[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* SDEGEN * Calculate the search direction and branch if it is not downhill. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ mp = meql + 1; dtest = 0.0e0; L_10: dmlc17( n, *nact, z, d, ztg, gm, relacc, ddotgm ); if (*ddotgm == 0.0e0) goto L_120; /* Branch if there is no need to consider any degenerate constraints. * The test gives termination if two consecutive additions to the * active set fail to increase the predicted new value of F. * */ if (*nact == mdeg) goto L_120; np = *nact + 1; sum = 0.0e0; for (j = np; j <= n; j++) { sum += SQ(Ztg[j]); } if (dtest > 0.0e0 && sum >= dtest) { if (itest == 1) goto L_120; itest = 1; } else { dtest = sum; itest = 0; } /* Add a constraint to the active set if there are any significant * violations of degenerate constraints. * */ k = *nact; dmlc10( n, m, a, lda, iact, nact, z, u, d, relacc, mdeg, gmnew, parnew, cgrad ); if (*nact == k) goto L_120; Par[*nact] = 0.0e0; /* Calculate the new reduced gradient and Lagrange parameters. * */ L_30: for (i = 1; i <= n; i++) { Gmnew[i] = Gm[i]; } k = *nact; L_50: temp = 0.0e0; iz = k; for (i = 1; i <= n; i++) { temp += Z[iz]*Gmnew[i]; iz += n; } temp *= U[k]; Parnew[k] = Par[k] + temp; if (k == *nact) Parnew[k] = fmin( Parnew[k], 0.0e0 ); j = Iact[k]; if (j <= m) { for (i = 1; i <= n; i++) { Gmnew[i] += -temp*A(i - 1,j - 1); } } else { jm = j - m; if (jm <= n) { Gmnew[jm] += temp; } else { Gmnew[jm - n] -= temp; } } k -= 1; if (k > meql) goto L_50; /* Set RATIO for linear interpolation between PAR and PARNEW. * */ ratio = 0.0e0; if (mp < *nact) { ku = *nact - 1; for (k = mp; k <= ku; k++) { if (Parnew[k] > 0.0e0) { ratio = Parnew[k]/(Parnew[k] - Par[k]); idrop = k; } } } /* Apply the linear interpolation. * */ theta = 1.0e0 - ratio; for (k = mp; k <= *nact; k++) { Par[k] = fmin( theta*Parnew[k] + ratio*Par[k], 0.0e0 ); } for (i = 1; i <= n; i++) { Gm[i] = theta*Gmnew[i] + ratio*Gm[i]; } /* Drop a constraint if RATIO is positive. * */ if (ratio > 0.0e0) { dmlc14( n, m, a, lda, iact, nact, z, u, relacc, idrop ); for (k = idrop; k <= *nact; k++) { Par[k] = Par[k + 1]; } goto L_30; } /* Return if there is no freedom for a new search direction. * */ if (*nact < n) goto L_10; *ddotgm = 0.0e0; L_120: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc09( void (*dmlcfg)(long,double[],double*,double[],LOGICAL32*), long n, double x[], double g[], double d[], double xs[], double gs[], double relacc, double stepcb, double ddotg, double *f, double *step, long *nfvals, long nfmax, double gopt[], double xl[], double xu[], long *gmode, double w20[]) { LOGICAL32 haveg, keep, noback; long int i, icount; double ddotgb, dghgh, dgknot, dglow, dgmid, dgopt, fbase, fhgh, flow, fopt, ratio, relint, sbase, size, stphgh, stplow, stpmin, stpopt, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const G = &g[0] - 1; double *const Gopt = &gopt[0] - 1; double *const Gs = &gs[0] - 1; double *const X = &x[0] - 1; double *const Xl = &xl[0] - 1; double *const Xs = &xs[0] - 1; double *const Xu = &xu[0] - 1; /* end of OFFSET VECTORS */ /* LSRCH Line search * 4/18/91 CLL changed handling of ICOUNT and NFVALS. Added GMODE. * ------------------------------------------------------------------ * Subroutine Arguments * * Input: DMLCFG, N, D(), RELACC, STEPCB, DDOTG, NFMAX, XL(), XU() * Inout: X(), G(), F, NFVALS, GMODE, W20() * Out: GS(), STEP * Work space: XS(), GOPT() * * DMLCFG, N,X(),G(),D() * XS() [out] Copy of the initial X(). * GS() [out] Copy of the initial G(). * RELACC,STEPCB, DDOTG, * F,STEP, NFVALS,NFMAX,GOPT(), XL(), XU(), * GMODE [integer,inout] For use by DMLC20. * W20() [float,inout] For use by DMLC20. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * * Initialization. * * print'(/'' DMLC09 1290.. STEPCB='',g13.5/ * * '' .. D()='',3g13.5/(22x,3g13.5))', * * STEPCB, (D(I),I=1,N) */ relint = 0.9e0; icount = 0; ratio = -1.0e0; for (i = 1; i <= n; i++) { Xs[i] = X[i]; Gs[i] = G[i]; Gopt[i] = G[i]; if (D[i] != 0.0e0) { temp = fabs( X[i]/D[i] ); if (ratio < 0.0e0 || temp < ratio) ratio = temp; } } *step = fmin( 1.0e0, stepcb ); stpmin = fmax( relacc*ratio, 1.0e-12**step ); *step = fmax( stpmin, *step ); sbase = 0.0e0; fbase = *f; ddotgb = ddotg; stplow = 0.0e0; flow = *f; dglow = ddotg; stphgh = 0.0e0; stpopt = 0.0e0; fopt = *f; dgopt = fabs( ddotg ); noback = FALSE; /* print'('' DMLC09 1318.. RATIO,STPMIN='',2g13.5/ * * '' DDOTG ='',g13.5)', * * RATIO, STPMIN, DDOTG * * Calculate another function and gradient value. * */ L_20: ; if (noback && *step < stpopt) { /* print*,' DMLC09 1355.. Ending DMLC09 on NOBACK.' */ goto L_70; } noback = FALSE; for (i = 1; i <= n; i++) { X[i] = Xs[i] + *step*D[i]; } (*dmlcfg)( n, x, f, g, &haveg ); *nfvals += 1; /* call DMLCFG(N, X, F, G, HAVEG) */ if (!haveg) dmlc20( dmlcfg, n, x, *f, g, xl, xu, gmode, nfvals, w20 ); icount += 1; /* print'(/'' DMLC09 1321.. ICOUNT='',i3/ * * '' .. STEP,STPOPT ='',2g23.15)', ICOUNT, STEP, STPOPT * print'('' .. X()='',2g23.15/(19x,2g23.15))', * * (X(I),I=1,N) * print'('' .. F ='',g23.15)',F * print'('' .. G()='',2g23.15/(19x,2g23.15))', * * (G(I),I=1,N) */ dgmid = 0.0e0; size = 0.0e0; for (i = 1; i <= n; i++) { temp = D[i]*G[i]; dgmid += temp; size += fabs( temp ); } size *= relacc; /* print'(/'' DMLC09 1335.. FOPT - F ='',g13.5/ * * '' .. DGMID ='',g13.5/ * * '' .. SIZE ='',g13.5/ * * '' .. DGOPT - abs(DGMID) ='',g13.5)', * * FOPT - F, DGMID, SIZE, DGOPT - abs(DGMID) */ if (*f <= fopt) { if (icount == 1 && dgmid <= 0.0e0) { /* print '('' DMLC09 1381.. Setting NOBACK on new test.'')' */ noback = TRUE; keep = TRUE; } else { keep = *f < fopt || fabs( dgmid ) < dgopt; } if (keep) { stpopt = *step; /* print'(/'' DMLC09 1384.. New STPOPT='',g13.5)', STPOPT */ fopt = *f; for (i = 1; i <= n; i++) { Gopt[i] = G[i]; } dgopt = fabs( dgmid ); } } if (nfmax > 0 && *nfvals >= nfmax) goto L_70; /* Modify the bounds on the steplength or convergence. * */ if (*f >= fbase + 0.1e0*(*step - sbase)*ddotgb) { if ((stphgh > 0.0e0 || *f > fbase) || dgmid > 0.5e0*ddotg) { stphgh = *step; fhgh = *f; dghgh = dgmid; goto L_60; } sbase = *step; fbase = *f; ddotgb = dgmid; } if (dgmid >= 0.7e0*ddotgb) goto L_70; stplow = *step; flow = *f; dglow = dgmid; L_60: if (stphgh > 0.0e0 && stplow >= relint*stphgh) goto L_70; /* Calculate the next step length or end the iterations. * * print'(/'' DMLC09 1366.. ICOUNT, STEP='',i23,g23.15/ * * '' STPLOW, STPHGH='',2g23.15)', * * ICOUNT, STEP, STPLOW, STPHGH */ if (stphgh == 0.0e0) { if (*step == stepcb) goto L_70; temp = 10.0e0; if (dgmid > 0.9e0*ddotg) temp = ddotg/(ddotg - dgmid); *step = fmin( temp**step, stepcb ); /* print'('' DMLC09 1374.. To 20 with STEP='',g23.15)', STEP */ goto L_20; } else if (icount == 1 || stplow > 0.0e0) { dgknot = 2.0e0*(fhgh - flow)/(stphgh - stplow) - 0.5e0*(dglow + dghgh); if (dgknot >= 0.0e0) { ratio = fmax( 0.1e0, 0.5e0*dglow/(dglow - dgknot) ); } else { ratio = (0.5e0*dghgh - dgknot)/(dghgh - dgknot); } *step = stplow + ratio*(stphgh - stplow); /* print'('' DMLC09 1384.. To 20 with STEP='',g23.15)', STEP */ goto L_20; } else { *step *= 0.1e0; if (*step >= stpmin) { /* print'('' DMLC09 1389.. To 20 with STEP='',g23.15)', STEP */ goto L_20; } } /* Return from subroutine. * */ L_70: ; if (*step != stpopt) { *step = stpopt; *f = fopt; for (i = 1; i <= n; i++) { X[i] = Xs[i] + *step*D[i]; G[i] = Gopt[i]; } } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc10( long n, long m, double *a, long lda, long iact[], long *nact, double z[], double u[], double d[], double relacc, long mdeg, double zzdiag[], double gmnew[], double cgrad[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, iadd, iz, j, jm, jmv, k, khigh, np; double cviol, cvmax, savabs, savsum, sum, sumabs, sumd, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cgrad = &cgrad[0] - 1; double *const D = &d[0] - 1; double *const Gmnew = &gmnew[0] - 1; long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; double *const Zzdiag = &zzdiag[0] - 1; /* end of OFFSET VECTORS */ /* NEWCON * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * * Initialization. * */ np = *nact + 1; khigh = mdeg; iz = 0; for (i = 1; i <= n; i++) { Zzdiag[i] = 0.0e0; for (j = np; j <= n; j++) { Zzdiag[i] += SQ(Z[iz + j]); } iz += n; } /* Calculate the scalar products of D with its constraints. * */ L_30: cvmax = 0.0e0; for (k = np; k <= khigh; k++) { j = Iact[k]; if (j <= m) { sum = 0.0e0; sumabs = 0.0e0; sumd = 0.0e0; for (i = 1; i <= n; i++) { temp = D[i]*A(i - 1,j - 1); sum += temp; sumabs += fabs( temp ); sumd += Zzdiag[i]*SQ(A(i - 1,j - 1)); } } else { jm = j - m; if (jm <= n) { sum = -D[jm]; } else { jm -= n; sum = D[jm]; } sumabs = fabs( sum ); sumd = Zzdiag[jm]; } /* Pick out the most violated constraint, or return if the * violation is negligible. * */ if (sum > relacc*sumabs) { cviol = sum*sum/sumd; if (cviol > cvmax) { cvmax = cviol; iadd = k; savsum = sum; savabs = sumabs; } } } if (cvmax <= 0.0e0) goto L_140; if (*nact == 0) goto L_120; /* Set GMNEW to the gradient of the most violated constraint. * */ j = Iact[iadd]; if (j <= m) { jmv = 0; for (i = 1; i <= n; i++) { Gmnew[i] = A(i - 1,j - 1); } } else { jmv = j - m; for (i = 1; i <= n; i++) { Gmnew[i] = 0.0e0; } if (jmv <= n) { Gmnew[jmv] = -1.0e0; } else { jmv -= n; Gmnew[jmv] = 1.0e0; } } /* Modify GMNEW for the next active constraint. * */ k = *nact; L_80: temp = 0.0e0; iz = k; for (i = 1; i <= n; i++) { temp += Z[iz]*Gmnew[i]; iz += n; } temp *= U[k]; j = Iact[k]; if (j <= m) { for (i = 1; i <= n; i++) { Gmnew[i] += -temp*A(i - 1,j - 1); } } else { jm = j - m; if (jm <= n) { Gmnew[jm] += temp; } else { Gmnew[jm - n] -= temp; } } /* Revise the values of SAVSUM and SAVABS. * */ sum = 0.0e0; sumabs = 0.0e0; for (i = 1; i <= n; i++) { temp = D[i]*Gmnew[i]; sum += temp; sumabs += fabs( temp ); } savsum = fmin( savsum, sum ); savabs = fmax( savabs, sumabs ); k -= 1; if (k >= 1) goto L_80; /* Add the new constraint to the active set if the constraint * violation is still significant. * */ if (jmv > 0) D[jmv] = 0.0e0; if (savsum <= relacc*savabs) goto L_130; L_120: k = *nact; dmlc12( n, m, a, lda, iact, nact, z, u, relacc, iadd, gmnew, cgrad ); if (*nact > k) goto L_140; /* Seek another constraint violation. * */ iadd = np; L_130: if (np < khigh) { k = Iact[khigh]; Iact[khigh] = Iact[iadd]; Iact[iadd] = k; khigh -= 1; goto L_30; } L_140: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc11( long n, long m, double *a, long lda, double b[], double xl[], double xu[], double x[], long iact[], long *nact, long *info, double z[], double u[], double xbig[], double relacc, double tol, long meql) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, idrop, iz, j, jx, k; double res, resabs, resbig, savex, scale, temp, tempa; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* SATACT * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (*nact == 0) goto L_50; for (k = 1; k <= *nact; k++) { /* Calculate the next constraint residual. * */ j = Iact[k]; if (j <= m) { res = B[j]; resabs = fabs( B[j] ); resbig = resabs; for (i = 1; i <= n; i++) { tempa = A(i - 1,j - 1); temp = tempa*X[i]; res -= temp; resabs += fabs( temp ); resbig += fabs( tempa )*Xbig[i]; } } else { jx = j - m; if (jx <= n) { res = X[jx] - Xl[jx]; resabs = fabs( X[jx] ) + fabs( Xl[jx] ); resbig = Xbig[jx] + fabs( Xl[jx] ); savex = Xl[jx]; } else { jx -= n; res = Xu[jx] - X[jx]; resabs = fabs( X[jx] ) + fabs( Xu[jx] ); resbig = Xbig[jx] + fabs( Xu[jx] ); savex = Xu[jx]; } } /* Shift X if necessary. * */ if (res != 0.0e0) { temp = res/resabs; if (k <= meql) temp = -fabs( temp ); if (tol == relacc || temp + relacc < 0.0e0) { *info = 1; scale = res*U[k]; iz = k; for (i = 1; i <= n; i++) { X[i] += scale*Z[iz]; iz += n; Xbig[i] = fmax( Xbig[i], fabs( X[i] ) ); } if (j > m) X[jx] = savex; /* Else flag a constraint deletion if necessary. * */ } else if (res/resbig > tol) { Iact[k] = -Iact[k]; } } } /* Delete any flagged constraints and then return. * */ idrop = *nact; L_40: if (Iact[idrop] < 0) { Iact[idrop] = -Iact[idrop]; dmlc14( n, m, a, lda, iact, nact, z, u, relacc, idrop ); } idrop -= 1; if (idrop > meql) goto L_40; L_50: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc12( long n, long m, double *a, long lda, long iact[], long *nact, double z[], double u[], double relacc, long indxbd, double ztc[], double cgrad[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, icon, inewbd, ipiv, iz, iznbd, j, jp, np; double sum, sumabs, temp, tempa, tempb, wcos, wpiv, wsin; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cgrad = &cgrad[0] - 1; long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; double *const Ztc = &ztc[0] - 1; /* end of OFFSET VECTORS */ /* ADDCON * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ np = *nact + 1; icon = Iact[indxbd]; Iact[indxbd] = Iact[np]; Iact[np] = icon; /* Form ZTC when the new constraint is a bound. * */ if (icon > m) { inewbd = icon - m; if (inewbd <= n) { temp = -1.0e0; } else { inewbd -= n; temp = 1.0e0; } iznbd = inewbd*n - n; for (j = 1; j <= n; j++) { Ztc[j] = temp*Z[iznbd + j]; } /* Else form ZTC for an ordinary constraint. * */ } else { for (i = 1; i <= n; i++) { Cgrad[i] = A(i - 1,icon - 1); } for (j = 1; j <= n; j++) { Ztc[j] = 0.0e0; iz = j; for (i = 1; i <= n; i++) { Ztc[j] += Z[iz]*Cgrad[i]; iz += n; } } } /* Find any Givens rotations to apply to the last columns of Z. * */ j = n; L_40: jp = j; j -= 1; if (j > *nact) { if (Ztc[jp] == 0.0e0) goto L_40; if (fabs( Ztc[jp] ) <= relacc*fabs( Ztc[j] )) { temp = fabs( Ztc[j] ); } else if (fabs( Ztc[j] ) <= relacc*fabs( Ztc[jp] )) { temp = fabs( Ztc[jp] ); } else { temp = fabs( Ztc[jp] )*sqrt( 1.0e0 + powi(Ztc[j]/Ztc[jp],2) ); } wcos = Ztc[j]/temp; wsin = Ztc[jp]/temp; Ztc[j] = temp; /* Apply the rotation when the new constraint is a bound. * */ iz = j; if (icon > m) { for (i = 1; i <= n; i++) { temp = wcos*Z[iz + 1] - wsin*Z[iz]; Z[iz] = wcos*Z[iz] + wsin*Z[iz + 1]; Z[iz + 1] = temp; iz += n; } Z[iznbd + jp] = 0.0e0; /* Else apply the rotation for an ordinary constraint. * */ } else { wpiv = 0.0e0; for (i = 1; i <= n; i++) { tempa = wcos*Z[iz + 1]; tempb = wsin*Z[iz]; temp = fabs( Cgrad[i] )*(fabs( tempa ) + fabs( tempb )); if (temp > wpiv) { wpiv = temp; ipiv = i; } Z[iz] = wcos*Z[iz] + wsin*Z[iz + 1]; Z[iz + 1] = tempa - tempb; iz += n; } /* Ensure orthogonality of Z(.,JP) to CGRAD. * */ sum = 0.0e0; iz = jp; for (i = 1; i <= n; i++) { sum += Z[iz]*Cgrad[i]; iz += n; } if (sum != 0.0e0) { iz = ipiv*n - n + jp; Z[iz] += -sum/Cgrad[ipiv]; } } goto L_40; } /* Test for linear independence in the proposed new active set. * */ if (Ztc[np] == 0.0e0) goto L_90; if (icon <= m) { sum = 0.0e0; sumabs = 0.0e0; iz = np; for (i = 1; i <= n; i++) { temp = Z[iz]*Cgrad[i]; sum += temp; sumabs += fabs( temp ); iz += n; } if (fabs( sum ) <= relacc*sumabs) goto L_90; } /* Set the new diagonal element of U() and return. * */ U[np] = 1.0e0/Ztc[np]; *nact = np; L_90: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc13( long n, long m, double *a, long lda, double b[], double xl[], double xu[], double x[], long iact[], long nact, double xbig[], double relacc, double *tol, long meql) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, j, jm, k, kl; double res, resabs, viol; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const B = &b[0] - 1; long *const Iact = &iact[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; /* end of OFFSET VECTORS */ /* ADJTOL: Change the tolerance TOL to a smaller value, if possible. * May also recompute XBIG(). * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Set VIOL to the greatest relative constraint residual of the first * NACT constraints. */ viol = 0.0e0; if (nact > meql) { kl = meql + 1; for (k = kl; k <= nact; k++) { j = Iact[k]; if (j <= m) { res = B[j]; resabs = fabs( B[j] ); for (i = 1; i <= n; i++) { res += -A(i - 1,j - 1)*X[i]; resabs += fabs( A(i - 1,j - 1)*Xbig[i] ); } } else { jm = j - m; if (jm <= n) { res = X[jm] - Xl[jm]; resabs = Xbig[jm] + fabs( Xl[jm] ); } else { jm -= n; res = Xu[jm] - X[jm]; resabs = Xbig[jm] + fabs( Xu[jm] ); } } if (res > 0.0e0) viol = fmax( viol, res/resabs ); } } /* Adjust TOL. * */ *tol = 0.1e0*fmin( *tol, viol ); if (*tol <= relacc + relacc) { *tol = relacc; for (i = 1; i <= n; i++) { Xbig[i] = fabs( X[i] ); } } return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc14( long n, long m, double *a, long lda, long iact[], long *nact, double z[], double u[], double relacc, long idrop) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, ibd, icon, ipiv, isave, iz, izbd, j, jp, nm; double denom, rjjp, sum, temp, tempa, tempb, ujp, wcos, wpiv, wsin; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* DELCON: Cycle through the constraint exchanges that are needed. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ nm = *nact - 1; if (idrop == *nact) goto L_60; isave = Iact[idrop]; /* Cycle through the constraint exchanges that are needed. * */ for (j = idrop; j <= nm; j++) { jp = j + 1; icon = Iact[jp]; Iact[j] = icon; /* Calculate the (J,JP) element of R. * */ if (icon <= m) { rjjp = 0.0e0; iz = j; for (i = 1; i <= n; i++) { rjjp += Z[iz]*A(i - 1,icon - 1); iz += n; } } else { ibd = icon - m; if (ibd <= n) { izbd = ibd*n - n; rjjp = -Z[izbd + j]; } else { ibd -= n; izbd = ibd*n - n; rjjp = Z[izbd + j]; } } /* Calculate the parameters of the next rotation. * */ ujp = U[jp]; temp = rjjp*ujp; denom = fabs( temp ); if (denom*relacc < 1.0e0) denom = sqrt( 1.0e0 + denom*denom ); wcos = temp/denom; wsin = 1.0e0/denom; /* Rotate Z when a bound constraint is promoted. * */ iz = j; if (icon > m) { for (i = 1; i <= n; i++) { temp = wcos*Z[iz + 1] - wsin*Z[iz]; Z[iz] = wcos*Z[iz] + wsin*Z[iz + 1]; Z[iz + 1] = temp; iz += n; } Z[izbd + jp] = 0.0e0; /* Rotate Z when an ordinary constraint is promoted. * */ } else { wpiv = 0.0e0; for (i = 1; i <= n; i++) { tempa = wcos*Z[iz + 1]; tempb = wsin*Z[iz]; temp = fabs( A(i - 1,icon - 1) )*(fabs( tempa ) + fabs( tempb )); if (temp > wpiv) { wpiv = temp; ipiv = i; } Z[iz] = wcos*Z[iz] + wsin*Z[iz + 1]; Z[iz + 1] = tempa - tempb; iz += n; } /* Ensure orthogonality to promoted constraint. * */ sum = 0.0e0; iz = jp; for (i = 1; i <= n; i++) { sum += Z[iz]*A(i - 1,icon - 1); iz += n; } if (sum != 0.0e0) { iz = ipiv*n - n + jp; Z[iz] += -sum/A(ipiv - 1,icon - 1); } } /* Set the new diagonal elements of U. * */ U[jp] = -denom*U[j]; U[j] = ujp/denom; } /* Return. * */ Iact[*nact] = isave; L_60: *nact = nm; return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc15( long n, long m, double xl[], double xu[], double x[], long iact[], long *meql, long *info, double z[], double u[], double xbig[], double *relacc) { long int _l0, i, iz, j, jact, nn; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iact = &iact[0] - 1; double *const U = &u[0] - 1; double *const X = &x[0] - 1; double *const Xbig = &xbig[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* INITZU: Initialize RELACC, Z(), U(). Adjust X(). Set XBIG(). * 1990-07-20 CLL Changed error handling and setting of INFO. * Expect INFO to be set to 1 before entering this subr. This subr * sets INFO = 4 if it finds XL(j) > XU(j) for some j. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * * Set RELACC. * July 1990 CLL Commented out following 6 lines and replaced with * reference to D1MACH(4) which is the smallest x such that the * stored value of 1+x is greater than 1. Efforts to determine * such an x with portable code such as the following 6 lines have * generally eventually failed on some new computer. Also using * D1MACH allows adjustments to be made for known deficiencies * in particular computers, for example for the Cray X/MP & Y/MP. * * ZTPAR=100.0D0 * RELACC=1.0D0 * 10 RELACC=0.5D0*RELACC * TEMPA=ZTPAR+0.5D0*RELACC * TEMPB=ZTPAR+RELACC * IF (ZTPAR .LT. TEMPA .AND. TEMPA .LT. TEMPB) GOTO 10 */ *relacc = 100.0e0*DBL_EPSILON; /* Seek bound inconsistencies and bound equality constraints. * */ *meql = 0; for (j = 1; j <= n; j++) { if (Xl[j] > Xu[j]) { *info = 4; ierm1( "DMLC15", *info, 0, "Bad bounds: XL(j) > XU(j)" , "j", j, ',' ); derv1( "XL(j)", Xl[j], ',' ); derv1( "XU(j)", Xu[j], '.' ); return; } if (Xl[j] == Xu[j]) *meql += 1; } /* Initialize U, Z and XBIG. * */ jact = 0; nn = n*n; for (i = 1; i <= nn; i++) { Z[i] = 0.0e0; } iz = 0; for (i = 1; i <= n; i++) { if (Xl[i] == Xu[i]) { X[i] = Xu[i]; jact += 1; U[jact] = 1.0e0; Iact[jact] = i + m + n; j = jact; } else { j = i + *meql - jact; } Z[iz + j] = 1.0e0; iz += n; Xbig[i] = fabs( X[i] ); } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc16( long n, long m, double *a, long lda, long iact[], long nact, double par[], double g[], double reskt[], double z[], double u[], double bres[], double *relaxf, long meql, double *ssqkt, double parw[], double resktw[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, icase, iz, j, jm, k, kk, kl; double ssqktw, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Bres = &bres[0] - 1; double *const G = &g[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Parw = &parw[0] - 1; double *const Reskt = &reskt[0] - 1; double *const Resktw = &resktw[0] - 1; double *const U = &u[0] - 1; double *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ /* Calculate the Lagrange parameters and the residual vector. * * Input: N, M, A(), LDA, IACT, NACT, G(), Z(), U(), BRES(), MEQL * Output: PAR(), RESKT, RELAXF, SSQKT * Work space: PARW(), RESKTW() * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ for (i = 1; i <= n; i++) { Reskt[i] = G[i]; } /* print'(/'' DMLC16 1993.. G()='',2g23.15/(19x,2g23.15))', * * (G(I),I=1,N) * print'('' .. U()='', 3g23.15/(8x,3g23.15))', (U(I),I=1,3) * print'('' .. Z()='', 2g23.15/(8x,2g23.15))', (Z(I),I=1,16) * print'('' .. A(,)='',2g23.15/(9x,2g23.15))',((A(I,J),J=1,4),I=1,3) */ if (nact > 0) { icase = 0; L_20: for (kk = 1; kk <= nact; kk++) { k = nact + 1 - kk; j = Iact[k]; temp = 0.0e0; iz = k; for (i = 1; i <= n; i++) { temp += Z[iz]*Reskt[i]; /* print'('' DMLC16 1991.. Z(),RESKT(),TEMP='',3g13.5)', * * Z(IZ),RESKT(I),TEMP */ iz += n; } temp *= U[k]; /* print'('' DMLC16 1995.. U(K),TEMP='', 2g13.5)', U(K), TEMP */ if (icase == 0) Par[k] = 0.0e0; if (k <= meql || Par[k] + temp < 0.0e0) { Par[k] += temp; /* print'('' DMLC16 1999.. TEMP,PAR(K)='',2g13.5)', * * TEMP,PAR(K) */ } else { temp = -Par[k]; Par[k] = 0.0e0; } if (temp != 0.0e0) { if (j <= m) { for (i = 1; i <= n; i++) { Reskt[i] += -temp*A(i - 1,j - 1); /* print'('' DMLC16 2009.. TEMP,A='',2g13.5/ * * '' RESKT(I) ='',g13.5)', * * TEMP, A(J,I), RESKT(I) */ } } else { jm = j - m; if (jm <= n) { Reskt[jm] += temp; /* print'('' DMLC16 2018.. TEMP,RESKT(JM)='',2g13.5)', * * TEMP,RESKT(JM) */ } else { Reskt[jm - n] -= temp; /* print'('' DMLC16 2022.. TEMP,RESKT(JM-N)='',2g13.5)', * * TEMP,RESKT(JM-N) */ } } } } /* Calculate the sum of squares of the KT residual vector. * */ *ssqkt = 0.0e0; if (nact == n) goto L_130; /* print'(/'' DMLC16 2014.. RESKT()='',4g12.4/(23x,4g12.4))', * * (RESKT(I),I=1,N) */ for (i = 1; i <= n; i++) { *ssqkt += SQ(Reskt[i]); } /* print'('' DMLC16 2018.. SSQKT='',g13.5)', SSQKT * * Apply iterative refinement to the residual vector. * */ if (icase == 0) { icase = 1; for (k = 1; k <= nact; k++) { Parw[k] = Par[k]; } for (i = 1; i <= n; i++) { Resktw[i] = Reskt[i]; } ssqktw = *ssqkt; goto L_20; } /* Undo the iterative refinement if it does not reduce SSQKT. * */ if (ssqktw < *ssqkt) { for (k = 1; k <= nact; k++) { Par[k] = Parw[k]; } for (i = 1; i <= n; i++) { Reskt[i] = Resktw[i]; } *ssqkt = ssqktw; } /* Calculate SSQKT when there are no active constraints. * */ } else { *ssqkt = 0.0e0; for (i = 1; i <= n; i++) { *ssqkt += SQ(G[i]); } } /* Predict the reduction in F if one corrects any positive residuals * of active inequality constraints. * */ *relaxf = 0.0e0; if (meql < nact) { kl = meql + 1; for (k = kl; k <= nact; k++) { j = Iact[k]; if (Bres[j] > 0.0e0) { *relaxf += -Par[k]*Bres[j]; } } } L_130: return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc17( long n, long nact, double z[], double d[], double ztg[], double gm[], double relacc, double *ddotgm) { long int i, iz, j, np; double sum, sumabs, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const D = &d[0] - 1; double *const Gm = &gm[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* SDIRN: Compute a search direction. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ *ddotgm = 0.0e0; if (nact >= n) goto L_60; /* Premultiply GM by the transpose of Z. * */ np = nact + 1; for (j = np; j <= n; j++) { sum = 0.0e0; sumabs = 0.0e0; iz = j; for (i = 1; i <= n; i++) { temp = Z[iz]*Gm[i]; sum += temp; sumabs += fabs( temp ); iz += n; } if (fabs( sum ) <= relacc*sumabs) sum = 0.0e0; Ztg[j] = sum; } /* Form D by premultiplying ZTG by -Z. * */ iz = 0; for (i = 1; i <= n; i++) { sum = 0.0e0; sumabs = 0.0e0; for (j = np; j <= n; j++) { temp = Z[iz + j]*Ztg[j]; sum -= temp; sumabs += fabs( temp ); } if (fabs( sum ) <= relacc*sumabs) sum = 0.0e0; D[i] = sum; iz += n; } /* Test that the search direction is downhill. * */ sumabs = 0.0e0; for (i = 1; i <= n; i++) { temp = D[i]*Gm[i]; *ddotgm += temp; sumabs += fabs( temp ); } if (*ddotgm + relacc*sumabs >= 0.0e0) *ddotgm = 0.0e0; L_60: return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc18( long n, long m, double *a, long lda, long iact[], double bres[], double d[], double *stepcb, double *ddotg, long mdeg, long *msat, long mtot, long *indxbd) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int i, iflag, j, jm, k, kl; double sp, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Bres = &bres[0] - 1; double *const D = &d[0] - 1; long *const Iact = &iact[0] - 1; /* end of OFFSET VECTORS */ /* STEPBD: * Set steps to constraint boundaries and find the least positive * one. * *>> 1990-04-02 CLL Changes to avoid jumping into scope of Block If. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ iflag = 0; *stepcb = 0.0e0; *indxbd = 0; k = mdeg; L_10: k += 1; if (k > mtot) goto L_30; /* Use remote code to compute scalar product. */ goto L_80; L_20: ; /* Set BRES(J) to indicate the status of the j-th constraint. * */ if (sp*Bres[j] <= 0.0e0) { Bres[j] = 0.0e0; } else { Bres[j] /= sp; if (*stepcb == 0.0e0 || Bres[j] < *stepcb) { *stepcb = Bres[j]; *indxbd = k; } } goto L_10; L_30: ; /* Try to pass through the boundary of a violated constraint. * */ L_40: ; if (*indxbd <= *msat) goto L_70; iflag = 1; k = *indxbd; /* Use remote code to compute scalar product. */ goto L_80; L_50: ; *msat += 1; Iact[*indxbd] = Iact[*msat]; Iact[*msat] = j; Bres[j] = 0.0e0; *indxbd = *msat; *ddotg -= sp; if (*ddotg < 0.0e0 && *msat < mtot) { /* Seek the next constraint boundary along the search direction. * */ temp = 0.0e0; kl = mdeg + 1; for (k = kl; k <= mtot; k++) { j = Iact[k]; if (Bres[j] > 0.0e0) { if (temp == 0.0e0 || Bres[j] < temp) { temp = Bres[j]; *indxbd = k; } } } if (temp > 0.0e0) { *stepcb = temp; goto L_40; } } L_70: ; return; /* ------------------------------------------------------------------ * This is a remote block of code to compute the * scalar product of D with the current constraint normal. * */ L_80: j = Iact[k]; if (j <= m) { sp = 0.0e0; for (i = 1; i <= n; i++) { sp += D[i]*A(i - 1,j - 1); } } else { jm = j - m; if (jm <= n) { sp = -D[jm]; } else { sp = D[jm - n]; } } /* Return from remote block is selected by IFLAG. */ if (iflag == 0) { goto L_20; } else { goto L_50; } return; #undef A } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc19( long n, double x[], long nact, double g[], double z[], double ztg[], double xs[], double gs[], double *zznorm) { long int i, iz, k, km, kp, np; double dd, dg, sum, temp, wcos, wsin; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const G = &g[0] - 1; double *const Gs = &gs[0] - 1; double *const X = &x[0] - 1; double *const Xs = &xs[0] - 1; double *const Z = &z[0] - 1; double *const Ztg = &ztg[0] - 1; /* end of OFFSET VECTORS */ /* Test if there is sufficient convexity for the update. * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ dd = 0.0e0; dg = 0.0e0; temp = 0.0e0; for (i = 1; i <= n; i++) { Xs[i] = X[i] - Xs[i]; dd += SQ(Xs[i]); temp += Gs[i]*Xs[i]; Gs[i] = G[i] - Gs[i]; dg += Gs[i]*Xs[i]; } if (dg < 0.1e0*fabs( temp )) goto L_90; /* Transform the Z matrix. * */ k = n; L_20: kp = k; k -= 1; if (k > nact) { if (Ztg[kp] == 0.0e0) goto L_20; temp = fabs( Ztg[kp] )*sqrt( 1.0e0 + powi(Ztg[k]/Ztg[kp],2) ); wcos = Ztg[k]/temp; wsin = Ztg[kp]/temp; Ztg[k] = temp; iz = k; for (i = 1; i <= n; i++) { temp = wcos*Z[iz + 1] - wsin*Z[iz]; Z[iz] = wcos*Z[iz] + wsin*Z[iz + 1]; Z[iz + 1] = temp; iz += n; } goto L_20; } /* Update the value of ZZNORM. * */ if (*zznorm < 0.0e0) { *zznorm = dd/dg; } else { temp = sqrt( *zznorm*dd/dg ); *zznorm = fmin( *zznorm, temp ); *zznorm = fmax( *zznorm, 0.1e0*temp ); } /* Complete the updating of Z. * */ np = nact + 1; temp = sqrt( dg ); iz = np; for (i = 1; i <= n; i++) { Z[iz] = Xs[i]/temp; iz += n; } if (np < n) { km = np + 1; for (k = km; k <= n; k++) { temp = 0.0e0; iz = k; for (i = 1; i <= n; i++) { temp += Gs[i]*Z[iz]; iz += n; } temp /= dg; sum = 0.0e0; iz = k; for (i = 1; i <= n; i++) { Z[iz] += -temp*Xs[i]; sum += SQ(Z[iz]); iz += n; } if (sum < *zznorm) { temp = sqrt( *zznorm/sum ); iz = k; for (i = 1; i <= n; i++) { Z[iz] *= temp; iz += n; } } } } L_90: return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc20( void (*dmlcfg)(long,double[],double*,double[],LOGICAL32*), long n, double x[], double f, double g[], double xl[], double xu[], long *gmode, long *nfvals, double w[]) { LOGICAL32 mode1, nopp; static LOGICAL32 lnopp; long int _l0, i, lgl, lpp, lxl; double d, dxn, fnew, h, hb, tp, xi, xnew; static double epsr, facx; static long kdebug = 0; static double fac = 0.0e0; static long kount = 0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const G = &g[0] - 1; double *const X = &x[0] - 1; double *const Xl = &xl[0] - 1; double *const Xu = &xu[0] - 1; /* end of OFFSET VECTORS */ /*>> 1991-04-15 FTK & CLL Made FIRST an argument of DMLC20 *>> 1990-09-11 Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA. * * This subr computes a finite difference approximation to a gradient * vector. Values of F are not computed outside the bounds given by XL() * and XU(). For any I for which XL(I) = XU(I), G(I) will be set to zero * and computation to estimate G(I) will be skipped. * On entry F must have already been evaluated at X(). * The algorithm used is based on the following observations. * 1. Since the gradient vector is being used in an iteration, there is * no point to computing it more accurately than is useful in the * iteration. An extra N function values to get better error * estimates is not likely to pay for itself. * 2. If the delta x used for the difference is in the direction of the * next move, then discretization error in approximating the gradient * is likely to help rather than hurt. * 3. The sign of the change in X(I) for the next move is likely to be * the same as it was on the last move. * 4. Because of 3 and 4, the increment used is larger than would be * suggested from a consideration of round off and discretization * errors. * 5. When the distance between values on successive iterations gets very * small it is useful to have second derivative information in order * to balance the effects of round off and discretization errors. * * *************************** Variables ******************************** * * D The difference between the current X(I) and the value on the * previous iteration. * DXN L2 norm of the move just taken. * DMLCFG [in, subroutine name] Name of user-provided subroutine for * computing function values. When called from this subr, DMLCFG must * compute F, set HAVEG = .false., and not store into G(). * EPSR Save variable = relative machine precision. * F [in, float] On entry F will contain the function value * evaluated at the given X(). * FAC Save variable = sqrt(relative machine precision) * FACX Save variable, 4. * FAC. * G() [out, float] On return G() will contain an approximation of * the N-dimensional gradient of f with respect to x at X() computed * using one-sided finite differences. * GMODE [inout, integer] Has value 0, or 1. * 0 means initial entry. Compute 1-sided differences, assuming no * saved info, and set GMODE to 1. * 1 Compute differences, using saved info. * H Increment in X(I) used for computing gradient. * HB Lower bound desired for value of abs(H). * KOUNT Counts interations. If multiple of 10, two sided differences * are computed. * LGL W(LGL+I) contains the value of G(I) from previous iteration. * LNOPP Value of NOPP for start of next iteration. * LPP W(LPP+I) contains estimate of the second derivative (d/dv)G(I), * where v is the variable consisting of the linear combinations of * X's that gave the last move. * LXL W(LXL+I) contains X() from the previous iteration. * MODE1 .true. if have 1-sided difference, = .false. if have 2-sided. * N [in, integer] Number of components in X(), G(), XL(), & XU(). * NFVALS [inout, integer] Counter of number of calls to DMLCFG. * NOPP Logical variable set .true. if second derivative information is * not available. * TP Temporary storage. * X() [in, work, float] On entry contains the current N-dimensional * parameter vector. This subr will alter components in X() but will * reset X() to its original contents on return. * XI Current base value for X(I) when computing F for a gradient. * XL() [in, float] Lower bounds for components of X(). * XNEW Current value for X(I) in computing F for a gradient. Also * used for temporary storage of 1-sided gradient. * XU() [in, float] Upper bounds for components of X(). * W() [inout, float] Work array of dimension at least 4N+1. * W(0) and W(1) are initialized to -1.0 in [D/S]MLC02. If this * subr is never called, i.e., if the user provides computed * gradients, then W(1) will remain at this initial value. * If this subr is called, it will (generallly) set W(1) positive. * [D/S]MLC07 will (generally) set W(0) positive for subsequent * testing when and if this subr is called. * W(1:N) contains an estimate of the error in G(I) for use in * [D/S]MLC07. * W(N+1 : 4*N) is used as temporary work space in this subr. See the * usage of the indices LGL, LPP, LXL. * ------------------------------------------------------------------ * * ********************** Variable Declarations ************************* * */ LOGICAL32 haveg; /* logical HAVEG */ /* ********************** Start of Executable Code ********************** * */ if (fac == 0.0e0) { /* Get machine constants */ epsr = DBL_EPSILON; fac = sqrt( epsr ); facx = 4.e0*fac; } /* Set up base locations in W() */ lxl = n; lgl = lxl + lxl; lpp = lgl + lxl; dxn = 0.e0; if (*gmode == 0) { *gmode = 1; kount = 1; nopp = TRUE; } else { kount += 1; if ((kount%10) == 1) { if (w[0] > 1.e-5) kount = 0; } else { if (w[0] > 0.01e0) kount = 0; } for (i = 1; i <= lxl; i++) { dxn += powi(X[i] - w[lxl + i],2); } dxn = sqrt( dxn ); nopp = lnopp; } lnopp = FALSE; if (kdebug != 0) { printf("0Gradient #=%3ld F=%24.17e DXN=%15.8e\n I X DX G ERR H dG/dV\n", kdebug, f, dxn); kdebug += 1; } for (i = 1; i <= lxl; i++) { xi = X[i]; d = xi - w[lxl + i]; if (nopp) { hb = fac*(1.e0 + fabs( xi )); } else { hb = fac*(sqrt( fabs( f/w[lpp + i] ) ) + facx*(fabs( X[i] ) + facx)); } L_20: if ((w[lpp + i] < 0.e0) || ((kount%10) == 0)) { if ((xi + hb <= Xu[i]) && (xi - hb >= Xl[i])) { mode1 = FALSE; h = sign( fmin( 100.e0*hb, fmin( Xu[i] - xi, xi - Xl[i] ) ), d ); goto L_30; } } mode1 = TRUE; h = sign( fmax( hb, fabs( 1.e-4*d ) ), d ); L_30: xnew = xi + h; /* Adjust H to keep XNEW within bounds. */ if (xnew > Xu[i]) { xnew = Xu[i]; if (xnew - xi < hb) { if (xi - Xl[i] > xnew - xi) xnew = fmax( Xl[i], xi - hb ); } h = xnew - xi; } else if (xnew < Xl[i]) { xnew = Xl[i]; if (xi - xnew < hb) { if (Xu[i] - xi > xi - xnew) xnew = fmin( Xu[i], xi + hb ); } h = xnew - xi; } if (h == 0.e0) { G[i] = 0.e0; w[i] = 0.e0; } else { X[i] = xnew; (*dmlcfg)( lxl, x, &fnew, g, &haveg ); *nfvals += 1; /* CALL DMLCFG(LXL, X, FNEW, G, HAVEG) */ G[i] = (fnew - f)/h; w[i] = epsr*(fabs( f ) + fabs( fnew ))/fabs( h ); if (!mode1) { xnew = xi - h; X[i] = xnew; (*dmlcfg)( lxl, x, &fnew, g, &haveg ); *nfvals += 1; /* CALL DMLCFG(LXL, X, FNEW, G, HAVEG) */ xnew = (fnew - f)/(xnew - xi); w[i] += 5.e-4*fabs( xnew - G[i] ); G[i] = .5e0*(G[i] + xnew); } X[i] = xi; } if (dxn == 0.e0) { lnopp = TRUE; if (mode1) w[i] = 0.e0; } else { w[lpp + i] = fmax( 0.e0, fabs( (G[i] - w[lgl + i]) - 10.e0* fabs( w[i] ) )/dxn ) + fac*fabs( f ); tp = 10.e0*fabs( h )*w[lpp + i]; if (mode1) { if (tp > fabs( G[i] )) { if ((kount%10) != 0) { kount = 0; goto L_20; } } w[i] += tp; } else { if (fabs( xnew - G[i] ) > fmin( 10.e0*tp, 1.e-3*fabs( G[i] ) )) w[lpp + i] = -w[lpp + i]; } } if (kdebug != 0) { printf("%3ld%22.14e%10.2e%12.4e%12.4e%11.3e%10.2e\n", i, X[i], X[i] - w[lxl + i], G[i], w[i], h, w[lpp + i]); } w[lxl + i] = xi; w[lgl + i] = G[i]; } return; } /* end of function */ /* ================================================================== */ void /*FUNCTION*/ dmlc21( long ipmore, LOGICAL32 all, long n, long iterc, long nfvals, double f, double tol, double x[], double g[], long nact, long iact[], double par[], double reskt[], double enrmkt) { /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const G = &g[0] - 1; long *const Iact = &iact[0] - 1; double *const Par = &par[0] - 1; double *const Reskt = &reskt[0] - 1; double *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /*>> 1990-07-11 CLL * This subr prints results (first feasible, intermediate, and * final). * ------------------------------------------------------------------ * Arguments * * IPMORE [in, integer] = 0 or 1. 0 means print ITERC, NFVALS, F, X(), * and G(). 1 means also print IACT(), PAR(), and RESKT(). * However see ALL for an exception. * ALL [in, logical] False means G(), PAR(), and RESKT() might not * contain valid data and therefore will not be printed. * N, ITERC, NFVALS, * F, TOL, X(), G(), NACT, IACT(), PAR(), RESKT(), * ENRMKT [in] Euclidian norm of RESKT(). * ------------------------------------------------------------------ */ long int i, k; /* integer I */ /* ------------------------------------------------------------------ */ printf("\n %6ld iterations. %6ld function evals. F = %13.5g\n TOL =%13.5g Norm of RESKT =%13.5g\n", iterc, nfvals, f, tol, enrmkt); printf("\n X ="); for (i = 0; i < n; i+=5){ for (k = i; k < (i < n - 4 ? i+5 : n); k++) printf("%13.5g", x[k]); if (i < n - 4) printf ( "\n " );} if (ipmore > 0) { /* print'('' X ='',5g13.5/(8x,5g13.5))', (X(I),I=1,N) */ if (all) { printf("\n G ="); for (i = 0; i < n; i+=5){ for (k = i; k < (i < n - 4 ? i+5 : n); k++) printf("%13.5g", g[k]); if (i < n - 4) printf ( "\n " );} if (nact == n) { /* print'(a,5g13.5/(8x,5g13.5))',' G =',(G(I),I=1,N) */ printf(" Kuhn-Tucker residual vector is zero.\n"); } else { printf("\n RESKT ="); for (i = 0; i < n; i+=5){ for (k = i; k < (i < n - 4 ? i+5 : n); k++) printf("%13.5g", reskt[k]); if (i < n - 4) printf ( "\n " );} } /* print'('' RESKT ='',5g13.5/(8x,5g13.5))', (RESKT(I),I=1,N) */ } if (nact == 0) { printf(" No active constraints.\n"); } else { printf("\n IACT ="); for (i = 0; i < nact; i+=12){ for (k = i; k < (i < nact - 11 ? i+12 : nact); k++) printf("%5ld", iact[k]); if (i < nact - 11) printf ( "\n " );} if (all) { /* print'(a,12i5/(8x,12i5))',' IACT =',(IACT(I),I=1,NACT) */ printf("\n PAR ="); for (i = 0; i < nact; i+=5){ for (k = i; k < (i < nact - 4 ? i+5 : nact); k++) printf("%13.5g", par[k]); if (i < nact - 4) printf ( "\n " );} printf("\n"); } /* print'('' PAR ='',5g13.5/(8x,5g13.5))', (PAR(I),I=1,NACT) */ } } return; } /* end of function */