SUBROUTINE DMLC01 (DMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,ACC, * IPRINT, NFMAX, IW, LIW, W, LW) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c File: DMLC.for Minimization with linear constraints. C>> 2001-10-05 DMLC Krogh Minor format changes -- 1PE => 1P,E C>> 2001-06-18 DMLC Krogh Changed ". LT." to " .LT." C>> 1996-07-08 DMLC Krogh Multiple changes for C conversion. c>> 1996-03-30 DMLC Krogh Removed MAX from type stmt., added external. c>> 1995-11-15 DMLC Krogh Moved formats up for C conversion. C>> 1994-11-11 DMLC Krogh Declared all vars. c>> 1994-11-02 DMLC Krogh Changes to use M77CON c>> 1992-04-27 DMLC CLL c>> 1992-04-17 CLL c>> 1992-03-27 CLL c>> 1992-02-05 CLL Declared all floating point variables. c>> 1992-01-16 CLL c>> 1991-06-10 CLL & FTK Editing comments. c>> 1991-04-30 CLL & FTK c>> 1991-04-19 CLL c>> 1991-04-15 FTK & CLL Made FIRST an argument of DMLC20 c>> 1990-09-12 C. L. Lawson and F. T. Krogh, JPL c>> 1990-07-25 C. L. Lawson, JPL c>> 1990-07-12 C. L. Lawson, JPL c>> 1990-04-03 C. L. Lawson, JPL c--D replaces "?": ?MLC,?MLC01,?MLC02,?MLCFG,?MLC20,?MLC21,?MLC04,?MLC05 c--& ?MLC13,?MLC15,?MLC03,?MLC12,?MLC06,?MLC11,?MLC09,?MLC16,?MLC19 c--& ?MLC07,?MLC18,?MLC08,?MLC14,?MLC10,?MLC17,?ERV1 c Also calls ERMOR, ERMSG, IERM1, IERV1 C ------------------------------------------------------------------ c This algorithm and the original code is due to c M. J. D. Powell, Cambridge Univ., 1989. c 1990-04-03 This version adapted for usage at JPL by c C. L. Lawson. Instrumented with "c--" lines to automate type c conversions. c 1990-07-11 CLL Added new argument NFMAX to DMLC01. Added new c subroutine DMLC20 to compute gradient by finite differences. c Added argument HAVEG in the user-provided subr DMLCFG. c Added args XL and XU to subr _MLC09 in order to be able to pass c them on to DMLC20. c Added reference to [D/R]1MACH in _MLC15. c 1990-07-12 CLL Transposed indices in the array A(,). c The constraints are now Ax .le. b rather than (A**t)x .le. b. c 1990-07-19 CLL Reorganized the arg list of DMLC01. c Prev args INFO, NACT, IACT() are embedded in new arg IW(). c New item, FVAL, and previous args W() and PAR() are combined c into new arg W(). The new item FVAL is the final value of the c objective function. c 1990-07-24 CLL Setting INFO = 0 initially. It keeps this value c till it is changed for some reason. Previously it was set to c 4 initially. c Changed specification of IPRINT. It now contains 5 print c parameters. This gives more flexibility, such as requesting c printing only at termination. c Moved the main block of intermediate printing into a new subr, c DMLC21. This allows cleaner logic for testing when to print. c DMLC21 is called both from DMLC02 and from _MLC05. In revising c the print control we removed variables NFVALK and ITERP, and c added variables NFVPR and TOLPRT. c Printing of error messages is now done using JPL MATH77 error c printing subroutines. Requires linking of c DERV1, ERFIN, ERMOR, ERMSG, IERM1, and IERV1. The calls are in c DMLC01, DMLC02, and _MLC15. ERFIN is not called directly. c Setting INFO = 9 in DMLC02 when terminate due to solution being c determined just by the constraints. Also compute and return c the func value in this case. c 1990-09-12 CLL and FTK Added arg. W20 for _MLC20. c 1991-04-19 CLL Change to have NFVALS count every call to _MLCFG. c Move test of IPFREQ to be sure last iteration gets printed if c it is a multiple of IPFREQ. Changed usage of NFVPR. c 1992 March/April CLL Additions to the code in MINFUN (C05) and c LSRCH (C09). The code will now return INFO = 2 in some cases c in which it previously returned INFO = 3. It will also finish in c fewer iterations in some cases. Added more to error msg on quit- c ting with INFO = 3. Added HAVEG into arg list of _MLC05. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c-- Begin Mask code changes. c Names of subroutines have been changed as follows: c Original names New DP names New SP names c GETMIN DMLC01 SMLC01 c MINFLC DMLC02 SMLC02 c EQCONS DMLC03 SMLC03 c GETFES DMLC04 SMLC04 c MINFUN DMLC05 SMLC05 c CONRES DMLC06 SMLC06 c GETD DMLC07 SMLC07 c SDEGEN DMLC08 SMLC08 c LSRCH DMLC09 SMLC09 c NEWCON DMLC10 SMLC10 c SATACT DMLC11 SMLC11 c ADDCON DMLC12 SMLC12 c ADJTOL DMLC13 SMLC13 c DELCON DMLC14 SMLC14 c INITZU DMLC15 SMLC15 c KTVEC DMLC16 SMLC16 c SDIRN DMLC17 SMLC17 c STEPBD DMLC18 SMLC18 c ZBFGS DMLC19 SMLC19 c FGCALC DMLCFG SMLCFG c New code by FTK --> DMLC20 SMLC20 c New code by CLL --> DMLC21 SMLC21 c-- End Mask code changes. C C This is main entry point to a package of subroutines that calculate C the least value of a differentiable function of several variables C subject to linear constraints on the values of the variables. c ------------------------------------------------------------------ c Subroutine Arguments c C DMLCFG [in] Name of subroutine having an interface of the form C SUBROUTINE DMLCFG (N, X, F, G, HAVEG) c INTEGER N c [DOUBLE PRECISION or REAL] X(N), F, G(N) C LOGICAL HAVEG c provided by the user to compute values of the objective function, c F, and optionally its gradient vector, G(), at the argument value c given by X(). The user may choose either to write code to compute c G() or not. If DMLCFG computes G() it must also set HAVEG = .true c If DMLCFG does not compute G() it must set HAVEG = .false. c In this latter case the package will make additional calls to c DMLCFG to compute a finite difference approximation to the c gradient vector. c The vectors X() and G() will be of length N. The subroutine must c not change the values of N or X(). c DMLCFG may use COMMON block(s) if neccessary to communicate c data between the user's main program and DMLCFG. C N is the number of variables and must be set by the user. C M is the number of linear constraints (excluding simple bounds) and C must be set by the user. C MEQ is the number of constraints that are equalities and must be set C by the user. C A(.,.) is a 2-dimensional array whose rows are the gradients of C the M constraint functions. Its entries must be set by the user C and its dimensions must be at least M by N. C LDA is the actual first dimension of the array A that is supplied by C the user, so its value may not be less than M. C B(.) is a vector of constraint right hand sides that must also be set C by the user. Specifically the constraints on the variables X(I) C I=1(1)N are C A(K,1)*X(1)+...+A(K,N)*X(N) .EQ. B(K) K=1,...,MEQ C A(K,1)*X(1)+...+A(K,N)*X(N) .LE. B(K) K=MEQ+1,...,M . C Note that the data that define the equality constraints come C before the data of the inequalities. C XL(.) and XU(.) are vectors whose components must be set to lower and C upper bounds on the variables. Choose very large negative and C positive entries if a component should be unconstrained, or set C XL(I)=XU(I) to freeze the I-th variable. Specifically these C simple bounds are C XL(I) .LE. X(I) and X(I) .LE. XU(I) I=1,...,N . C X(.) is the vector of variables of the optimization calculation. Its C initial elements must be set by the user to an estimate of the C required solution. The subroutines can usually cope with poor C estimates, and there is no need for X(.) to be feasible initially. C These variables are adjusted automatically and the values that C give the least feasible calculated value of the objective function C are available in X(.) on the return from DMLC01. C ACC is a tolerance on the first order conditions at the calculated C solution of the optimization problem. These first order C conditions state that, if X(.) is a solution, then there is a set C of active constraints with indices IACT(K) K=1(1)NACT, say, such C that X(.) is on the boundaries of these constraints, and the C gradient of the objective function can be expressed in the form C GRAD(F)=PAR(1)*GRAD(C(IACT(1)))+... C ...+PAR(NACT)*GRAD(C(IACT(NACT))) . C Here PAR(K) K=1(1)NACT are Lagrange multipliers that are C nonpositive C for inequality constraints, and GRAD(C(IACT(K))) is the gradient C of the IACT(K)-th constraint function, so it is A(IACT(K),.) if C IACT(K) .LE. M, and it is minus or plus the J-th coordinate vector c if the constraint is the lower or upper bound on X(J) c respectively. The normal return from the calculation occurs when c X(.) is feasible and the sum of squares of components of the c vector RESKT(.) is at most ACC**2, where RESKT(.) is the c N-component vector of residuals of the first order condition that c is displayed above. Sometimes the package cannot satisfy this c condition, because noise in the function values can prevent a c change to the variables, no line search being allowed to increase c the objective function. C IPRINT [in, integer] Contains print control parameters, packed as c follows: c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST c Printing only begins after a first feasible point is found. c The items to be printed are selected by IPMORE. The other c parameters determine when to print. c IPFREQ is zero or positive. If positive, printing will be done on c iterations 0, IPFREQ, 2*IPFREQ, etc. c IPTOL = 0 or 1. 1 means to print the new TOL value c and the standard items each time TOL is changed. c IPFRST = 0 or 1. 1 means to print on the first iteration, i.e. on c iteration No. 0. c IPMORE = 0 or 1. 0 means the items to be printed are ITERC, c NFVALS, F, X(1:N), and G(1:N). 1 means to also print c IACT(1:NAXT), PAR(1:NACT) and RESKT(1:N). c IPLAST = 0 or 1. 1 means to print the final results, and the c reason for termination. c NFMAX [in, integer] If positive this sets an upper limit on the c number of function evaluations. When gradients are estimated by c differences, the count includes the function evaluations for this c purpose also. If zero there is no upper limit. c IW(.) [work, out, integer] Integer work space. Also used to return c status information. Its length must be at least LIW. c On return: c IW(1) contains INFO. Indicates reason for termination. c See below for explanation of values. c IW(2) contains ITERC. No. of iterations. c IW(3) contains NFVALS. No. of function evaluations, not c counting extra func evals done to estimate the gradient c when the gradient is not computed explicitly. In this c latter case the actual No. of func evals will be c (K+1)*NFVALS, where K is the number of solution c components whose lower and upper bounds are not equal. c IW(4) contains NACT. The number of active constraints at the c solution point. Will be in the interval [0, N]. c IW(5:4+NACT) contains IACT(1:NACT). IACT() is used as work c space of length M + 2*N. On return the first NACT c locations contain the indices of the constraints that c are active at the final point. Indices [1:M] refer to c rows of the system Ax .le. b. Indices [M+1:M+N] c refer to component lower bounds. Indices [M+N+1:M+2*N] c refer to component upper bounds. C c LIW [in, integer] Dimension of the IW() array. Require c LIW .ge. 4 + M + 2*N. c C W(.) [work, out, float] Working space array of float variables. c Also used to return results. c Its length must be at least LW. c On return: c W(1) contains FVAL, the final value of the objective function. c W(2:1+N) contains the final gradient vector, GRAD(1:N). c W(2+N:1+2*N) contains the Kuhn-Tucker residual vector, RESKT(1:N). c W(2+2*N:1+2*N+NACT) contains the Lagrange multipliers, PAR(1:NACT) c where NACT has a value in [0,N]. GRAD(), RESKT(), and PAR() c are mentioned above in the description of ACC. c LW [in, integer] Dimension of the W() array. Require c LW .ge. 3 + M + 16*N + N**2. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Return values of INFO c c INFO indicates the reason for termination. c C INFO=1 X(.) is feasible and the condition that depends on C ACC is satisfied. C INFO=2 X(.) is feasible and rounding errors are preventing C further progress. C INFO=3 X(.) is feasible but the objective function fails to C decrease although a decrease is predicted by the current gradient C vector. If this return occurs and RESKT(.) has large components C then the user's calculation of the gradient of the objective C function may be incorrect. One should also question the coding of C the gradient when the final rate of convergence is slow. C INFO=4 In this case the calculation cannot begin because IA C is less than N or because the lower bound on a variable is greater C than the upper bound. C INFO=5 This value indicates that the equality constraints C are inconsistent. These constraints include any components of C X(.) that are frozen by setting XL(I)=XU(I). C INFO=6 In this case there is an error return because the C equality constraints and the bounds on the variables are found to C be inconsistent. C INFO=7 This value indicates that there is no vector of C variables that satisfies all of the constraints. Specifically, C when this return or an INFO=6 return occurs, the current active C constraints (whose indices are IACT(K) K=1(1)NACT) prevent any C change in X(.) that reduces the sum of constraint violations, C where only bounds are included in this sum if INFO=6. C INFO=8 In this case the limit on the number of calls of C subroutine DMLCFG (see below) has been reached, and there would C have been further calculation otherwise. c INFO = 9 The solution is determined just by the constraints. c In this case PAR() and RESKT() are not defined on return. c FVAL will be defined on return. G() will only be defined on c return if DMLCFG returns HAVEG = .true. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Note 1. The variables N, M, MEQ, LDA, ACC, IPRINT, NFMAX and the C elements of the arrays A(.,.), B(.), XL(.) and XU(.) are not c altered by the optimization procedure. Their values, and the c initial components of X(.) must be set on entry to DMLC01. c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Note 2. A paper on the method of calculation and a report on the C main features of the computer code are available from the author C M.J.D.Powell (D.A.M.T.P., University of Cambridge, Silver Street, C Cambridge CB3 9EW, England). c ------------------------------------------------------------------ external DMLCFG integer LDA, LIW, LW integer I, IBRES, ID integer IFVAL, IG, IGM, IGS, IIACT, IINFO, IITERC, INACT, INFO integer INFVAL, IPAR, IPRINT, IRESKT, IU, IW(LIW), IW20, IXBIG integer IXS, IZ, IZTG, KTNORM, M, MEQ, N, NFMAX double precision A(LDA,*),ACC, B(*),TEMP, XL(*),XU(*),X(*),W(LW) parameter(IFVAL = 1, KTNORM = 2, IG = 3) parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5) C ------------------------------------------------------------------ c Check sizes of LIW and LW INFO = 0 if(LIW .lt. 4 + M + 2*N) then INFO = 4 call IERM1('DMLC01', INFO, 0,'Dimension LIW is too small.', * 'LIW',LIW, ',') call IERV1('Need',4 + M + 2*N,'.') endif c if(LW .lt. 3 + M + N*(16 + N)) then INFO = 4 call IERM1('DMLC01', INFO, 0,'Dimension LW is too small.', * 'LW',LW, ',') call IERV1('Need',3 + M + N*(16 + N),'.') endif if(INFO .ne. 0) then IW(1) = INFO return endif c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Partition the array W() c 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 C Need 4*N+1 locations beginning at IW20. IW20=IGS+N C Partition the array IW() c parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5) C c Zero the RESKT vector so it will be safe to compute its norm on c return from DMLC02 even in cases when DMLC02 does not compute c RESKT. c do 10 I = IRESKT, IRESKT + N -1 W(I) = 0.0D0 10 continue c C Call the optimization package. C CALL 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)) c TEMP = 0.0D0 do 20 I = IRESKT, IRESKT+N-1 TEMP = TEMP + W(I)**2 20 continue W(KTNORM) = sqrt(TEMP) RETURN END c ================================================================== SUBROUTINE DMLC02 (DMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X, * ACC,IACT,NACT,PAR, IPRINT, NFMAX, INFO, ITERC, NFVALS, * G,Z,U,XBIG,RESKT, ENRMKT, BRES,D,ZTG,GM,XS,GS, FVAL, W20) c c Main subroutine to control the Powell constrained minimization c package. c 1990-07-19 CLL Added NFMAX, ITERC, NFVALS, and FVAL to arg list. c Also added FVAL to arg list of DMLC05 that is called c from this subr. c NFMAX provides limit on no. of func evaluations, but zero means c no limit. This was previously overloaded onto INFO. c ITERC & NFVALS give output of counts of iterations and c function evaluations. c FVAL provides output of the final objective function value. c ------------------------------------------------------------------ c Arguments c c DMLCFG, N,M,MEQ,A(),LDA,B(),XL(),XU(),X(), c ACC,IACT(),NACT,PAR(), c IPRINT [in, integer] Contains print parameters packed as follows: c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST c NFMAX [in, integer] Limit on No. of function evaluations. But zero c means no limit. c INFO [out, integer] Termination status flag. c ITERC [out, integer] Count of No. of iterations. Initialized to 0 c in this subr. Incremented in DMLC05. c NFVALS [out, integer] Count of number of function evaluations. c G(),Z(),U(),XBIG(), c RESKT() c ENRMKT Euclidean norm of RESKT(). c BRES(),D(), c ZTG(),GM(),XS(),GS(), c FVAL [out, float] Final (best) value of F. c W20 [out, work, float] Working storage for DMLC20. c ------------------------------------------------------------------ c Descriptions of some of the internal variables. c c GMODE [integer] Initialized to 0 here. Sent to DMLC05 for eventual c use by DMLC20. c HAVEG [logical] Arg in call to DMLCFG, not used in this subroutine. c NFVPR [integer] Used to avoid printing results from the same c function evaluation more than once. Initialized to -1 in this c subr. Modified in DMLC05. Used in DMLC02 and DNLC05. c ------------------------------------------------------------------ external DMLCFG integer GMODE, I, IACT(*), INFO integer IPFREQ, IPFRST, IPLAST, IPMORE, IPRINT integer IPTOL, ITERC, K, LDA, M, MEQ, MEQL, MP, MSAT integer MTOT, N, NACT, NFMAX, NFVALS, NFVPR double precision A(LDA,*),ACC, B(*),BRES(*),D(*),ENRMKT, FVAL double precision G(*), GM(*), GS(*), PAR(*), RELACC, RESKT(*) double precision SSQKT, TOL, U(*) double precision W20(0:*), XBIG(*), XL(*), XU(*), X(*), XS(*) double precision Z(*), ZTG(*), ZZNORM logical HAVEG c ------------------------------------------------------------------ C Initialize ZZNORM, ITERC, NFVPR, NFVALS, INFO C ZZNORM=-1.0D0 W20(0) = -1.0d0 W20(1) = -1.0d0 GMODE = 0 ITERC=0 NFVPR = -1 NFVALS=0 INFO = 1 c c Decompose IPRINT into 5 integers: c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST C IPFREQ = IPRINT/100 I = IPRINT - IPFREQ * 100 IPTOL = I/8 I = I - IPTOL * 8 IPFRST = I/4 I = I - IPFRST * 4 IPMORE = I/2 IPLAST = I - IPMORE * 2 c C Check the bounds on N, M and MEQ. C if(N .lt. 1 .or. M .lt. 0 .or. MEQ .lt. 0 .or. M .lt. MEQ) then INFO = 4 call IERM1('DMLC02', INFO, 0,'Bad value for M, N, or MEQ.', * 'M',M, ',') call IERV1('N',N,',') call IERV1('MEQ',MEQ,'.') GOTO 40 END IF C C Initialize RELACC, Z, U and TOL. C CALL DMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC) IF (INFO .EQ. 4) go to 40 TOL=max(0.01D0,10.0D0*RELACC) if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or. * IPLAST .ne. 0 .or. IPTOL .ne. 0) then print'(/6x,''Beginning subroutine DMLC02, called from DMLC01:'' * /6x,''N ='',i4,'', M ='',i4,'', MEQ ='',i4,'', ACC ='', * g11.3/6x,''RELACC ='',g11.3,'', TOL = '',g11.3)', * N, M, MEQ, ACC, RELACC, TOL endif C C Add any equality constraints to the active set. C IF (MEQ .GT. 0) THEN CALL DMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC,XS, 1 GS) IF (INFO .EQ. 5) THEN call ERMSG('DMLC02',INFO,0, * 'Equality constraints are inconsistent.','.') GOTO 40 END IF END IF NACT=MEQL MSAT=MEQL C C Add the bounds to the list of constraints. C MTOT=NACT DO 10 I=1,N IF (XL(I) .LT. XU(I)) THEN MTOT=MTOT+2 IACT(MTOT-1)=M+I IACT(MTOT)=M+N+I END IF 10 CONTINUE C C Try to satisfy the bound constraints. C CALL DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG, 1 RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20) IF (MSAT .LT. MTOT) THEN INFO=6 call ERMSG('DMLC02',INFO,0, * 'Equalities and bounds are inconsistent.','.') GOTO 40 END IF C C Add the ordinary inequalities to the list of constraints. C IF (M .GT. MEQ) THEN MP=MEQ+1 DO 20 K=MP,M MTOT=MTOT+1 20 IACT(MTOT)=K END IF c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 30 continue c Begin main loop C C Correct any constraint violations. C CALL DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG, 1 RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20) IF (MSAT .LT. MTOT) THEN INFO=7 call ERMSG('DMLC02',INFO,0, * 'Constraints are inconsistent.','.') GOTO 40 ELSE IF (MEQL .EQ. N) THEN c c Here the soln is determined just by the constraints. c PAR() and RESKT() are not defined on return in this case. c FVAL has not yet been evaluated. We evaluate FVAL here in c order to return its value. We leave G() not computed if HAVEG c is false, however it would be easy to compute G by calling c DMLC20 if that seemed preferable. c C%% (*dmlcfg)( n, x, fval, g, &haveg ); call DMLCFG(N, X, FVAL, G, HAVEG) NFVALS = NFVALS+1 ENRMKT = 0.0d0 INFO = 9 if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or. * IPLAST .ne. 0 .or. IPTOL .ne. 0) then print '(/'' [1] DMLC02..''/ * '' The solution is determined by the equality constraints.'')' call DMLC21(IPMORE, .false., N, ITERC, NFVALS, * FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) endif GOTO 40 END IF C C Minimize the objective function in the case when constraints are C treated as degenerate if their residuals are less than TOL. C c DMLC05 is MINFUN CALL 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) C C Reduce TOL if necessary. C IF (TOL .GT. RELACC .AND. NACT .GT. 0) THEN if (NFMAX .le. 0 .or. NFVALS .lt. NFMAX) then C DMLC13 is ADJTOL CALL DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC, * TOL, MEQL) GOTO 30 else INFO=8 endif END IF c c Here when INFO = 1, 2, 3, or 8. c We treat 1 and 2 as normal conditions and 3 and 8 as errors. c if (INFO .EQ. 1) then if (IPLAST .ne. 0) print*, * '[1] DMLC01 quitting because requested accuracy', * ' has been attained.' elseif (INFO .EQ. 2) then if (IPLAST .ne. 0) print*, * '[2] DMLC01 quitting due to limitation of', * ' computational precision.' elseif (INFO .EQ. 3) then call ERMSG('DMLC02', INFO, 0, * 'F not decreasing, although gradient indicates it should.',',') if(HAVEG) then call ERMOR('Could be due to limitation of precision',',') call ERMOR('or to incorrect code for the gradient.','.') else call ERMOR('Could be due to limitation of precision.','.') endif elseif (INFO .EQ. 8) then call IERM1('DMLC02', INFO, 0, * 'Reached specified max number of function evaluations.', * 'Count',NFVALS, ',') call IERV1('NFMAX',NFMAX,'.') endif c c Test to print results before quitting. c ENRMKT = sqrt(SSQKT) if(IPLAST .ne. 0 .and. NFVALS .ne. NFVPR) then c print*,'DMLC02.. Debug.. Call C21 on leaving C02.' call DMLC21(IPMORE, .true., N, ITERC, NFVALS, * FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) endif c 40 return end c ================================================================== SUBROUTINE DMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC, 1 AM,CGRAD) C C Try to add the next equality constraint to the active set. c ------------------------------------------------------------------ integer N, M, MEQ, LDA, IACT(*), MEQL, INFO integer I, IZ, J, JM, K, KEQ, NP double precision A(LDA,*),B(*),XU(*) double precision RELACC, RHS, SUM, SUMABS, VMULT double precision Z(*),U(*),AM(*),CGRAD(*) c ------------------------------------------------------------------ DO 50 KEQ=1,MEQ IF (MEQL .LT. N) THEN NP=MEQL+1 IACT(NP)=KEQ CALL DMLC12 (N,M,A,LDA,IACT,MEQL,Z,U,RELACC,NP,AM,CGRAD) IF (MEQL .EQ. NP) GOTO 50 END IF C C If linear dependence occurs then find the multipliers of the C dependence relation and apply them to the right hand sides. C SUM=B(KEQ) SUMABS=abs(B(KEQ)) IF (MEQL .GT. 0) THEN DO 10 I=1,N 10 AM(I)=A(KEQ,I) K=MEQL 20 VMULT=0.0D0 IZ=K DO 30 I=1,N VMULT=VMULT+Z(IZ)*AM(I) 30 IZ=IZ+N VMULT=VMULT*U(K) J=IACT(K) IF (J .LE. M) THEN DO 40 I=1,N 40 AM(I)=AM(I)-VMULT*A(J,I) RHS=B(J) ELSE JM=J-M-N AM(JM)=AM(JM)-VMULT RHS=XU(JM) END IF SUM=SUM-RHS*VMULT SUMABS=SUMABS+abs(RHS*VMULT) K=K-1 IF (K .GE. 1) GOTO 20 END IF C C Error return if the constraints are inconsistent. C IF (abs(SUM) .GT. RELACC*SUMABS) THEN INFO=5 GOTO 60 END IF 50 CONTINUE 60 RETURN END c ================================================================== SUBROUTINE DMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z, 1 U,XBIG,RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,GMNEW,PARNEW, 2 CGRAD,W20) C C Make the correction to X for the active constraints. c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, INFO, MEQL, MSAT, MTOT integer I, ITEST, MSATK, INDXBD double precision A(LDA,*),B(*),XL(*),XU(*),X(*),PAR(*),G(*) double precision U(*),XBIG(*),BRES(*),D(*) double precision GM(*),GMNEW(*),PARNEW(*), CGRAD(*) double precision W20(0:*), RELACC, STEPCB, SUMRES, SUMRSK, TOL double precision Z(*), ZTG(*) c ------------------------------------------------------------------ INFO=0 10 CALL DMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG,RELACC, 1 TOL,MEQL) IF (INFO .GT. 0) MSAT=NACT IF (MSAT .EQ. MTOT) GOTO 60 C C Try to correct the infeasibility. C 20 MSATK=MSAT SUMRSK=0.0D0 30 CALL DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES, 1 D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD,GM,GMNEW, 2 PARNEW,CGRAD,W20) C C Include the new constraint in the active set. C IF (STEPCB .GT. 0.0D0) THEN DO 40 I=1,N X(I)=X(I)+STEPCB*D(I) 40 XBIG(I)=max(XBIG(I),abs(X(I))) CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,GMNEW, * CGRAD) END IF C C Test whether to continue the search for feasibility. C IF (MSAT .LT. MTOT) THEN IF (STEPCB .EQ. 0.0D0) GOTO 50 IF (MSATK .LT. MSAT) GOTO 20 IF (SUMRSK .EQ. 0.0D0 .OR. SUMRES .LT. SUMRSK) THEN SUMRSK=SUMRES ITEST=0 END IF ITEST=ITEST+1 IF (ITEST .LE. 2) GOTO 30 C C Reduce TOL if it may be too large to allow feasibility. C 50 IF (TOL .GT. RELACC) THEN CALL DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC, 1 TOL,MEQL) GOTO 10 END IF END IF 60 RETURN END c ================================================================== SUBROUTINE 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) C C MINFUN: Solve the minimization problem using the current value of c TOL. c 1990-07-19 CLL Added FVAL to argument list to return final value c of F. c 1990-07-23 CLL Removed IPRINT from arg list and added IPFREQ, c IPTOL, IPFRST, IPMORE, and NFVPR. NFVPR keeps a record of c last function evaluation which printing has been done. c 1990-07-23 CLL Added logical variable TOLPRT. c ------------------------------------------------------------------ c Arguments c c DMLCFG, N,M,A(),LDA,B(),XL(),XU(),X(),ACC,IACT(),NACT, c PAR() c IPFREQ [in, integer] Zero or positive. If positive, printing of c intermediate results using DMLC21 will be done on iterations c 0, IPFREQ, 2*IPFREQ, etc. c IPTOL [in, integer] = 0 or 1. 1 means to print the new TOL value c and items printed by DMLC21 each time this subr is entered. c IPFRST [in, integer] = 0 or 1. 1 means to print using DMLC21 on c the first iteration, i.e. iteration No. 0. c IPMORE [in, integer] = 0 or 1. Will be passed to DMLC21 to cause c printing to include more items when = 1 than when = 0. c INFO [inout, integer] Status flag. c G(),Z(),U(),XBIG(),RELACC,ZZNORM,TOL,MEQL,MTOT, c NFVALS [inout, integer] Count of calls for function evaluation. c Only counts calls to DMLCFG when the algorithm directly needs a c function value, not calls to DMLCFG that occur in DMLC20 for c estimation of the gradient. c NFMAX [in, integer] Limit on No. of function evaluations, but zero c means no limit. c RESKT() c SSQKT Sum of Squares of elts of RESKT(). c BRES(),D(),ZTG(),GM(),XS(),GS(), c FVAL [out, float] On return set to current (best) value of F. c GMODE [inout, integer] For use by DMLC20. c HAVEG [out, logical] Primarily an internal variable in this subr, c but returned so calling subr can use it to select proper error msg c when INFO = 3. c Set true by DMLCFG if DMLCFG provides the gradient c vector and false otherwise. If false, this subr calls DMLC20 to c compute a finite difference approximation for the gradient. c W20 [out, work, float] Working storage for DMLC20. c ------------------------------------------------------------------ c Description of some of the internal variables. c c DOPRNT [logical] When IPTOL = 1 this subr prints TOL on entry and c sets DOPRNT = true. Also if IPFRST = 1 and ITERC = 0 this subr c sets DOPRNT = true. c Having DOPRNT = true causes subsequent printing of intermediate c results for the current iteration, after PAR() and RESKT() have c been computed. c HAVEG [logical] Set true by DMLCFG if DMLCFG provides the gradient c vector and false otherwise. If false, this subr calls DMLC20 to c compute a finite difference approximation for the gradient. c ------------------------------------------------------------------ external D1MACH, DMLCFG integer N, M, LDA, IACT(*), NACT, IPFREQ, IPTOL, IPFRST, IPMORE, * INFO, MEQL, MTOT, ITERC, NFVPR, NFVALS, NFMAX, GMODE integer I, ITERK, K, MSAT, INDXBD double precision D1MACH double precision A(LDA,N), ACC, B(*), BRES(*), DDOTG, DIFF double precision D(*), EPS, F, FPREV, FVAL, G(N), GM(*), GS(*) double precision PAR(*), RELACC, RELAXF, RESKT(*) double precision SSQKT, STEP, STEPCB, SUM, TOL double precision U(*), W20(0:*) double precision X(*), XBIG(*), XL(N), XS(*), XU(N) double precision Z(*), ZTG(*), ZZNORM logical DOPRNT, HAVEG save EPS, F data EPS / 0.0d0 / c ------------------------------------------------------------------ c Initialize the minimization calculation. c * print'(/''DMLC05.. Debug.. On entry, INFO = '',i6,'', ITERC ='', * i6)', INFO, ITERC if(EPS .eq. 0.0d0) then EPS = D1MACH(4) endif MSAT=MTOT ITERK=ITERC IF (NFVALS .EQ. 0 .OR. INFO .EQ. 1) THEN c++ CODE for ~.C. is active CALL DMLCFG (N,X,F,G, HAVEG) c++ CODE for .C. is inactive C%% (*dmlcfg)( n, x, &f, g, haveg ); C HAVEG = HAVEG c++ 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=NFVALS+1 if(.not. HAVEG) CALL DMLC20(DMLCFG,N,X,F,G,XL,XU, * GMODE, NFVALS, W20) END IF FPREV=abs(F+F+1.0D0) c Setting DIFF here is not needed for the c algorithm, but is convenient when doing c debug printing. CLL 3/26/92 DIFF = 1.0d0 IF (IPTOL .ne. 0) THEN print'(/'' New value of TOL = '',g13.5)', TOL DOPRNT = .true. elseif(IPFRST .ne. 0 .and. ITERC .eq. 0) then DOPRNT = .true. else DOPRNT = .false. END IF C C Calculate the next search direction. C DMLC06 is CONRES 10 CALL DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES,D, 1 ZTG,RELACC,TOL,STEPCB,DDOTG,MEQL,MSAT,MTOT,INDXBD,GM,RESKT,XS, 2 GS,W20) C C Calculate the Kuhn Tucker residual vector. C DMLC16 is KTVEC CALL DMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES,RELAXF,MEQL, 1 SSQKT,XS,GS) c c Test for printing results every IPFREQ iterations. c if(IPFREQ .ne. 0) then if( mod(ITERC, IPFREQ) .eq. 0) DOPRNT = .true. endif c if(DOPRNT) then call DMLC21(IPMORE, .true., N, ITERC, NFVALS, * F, TOL, X, G, NACT, IACT, PAR, RESKT, sqrt(SSQKT)) NFVPR = NFVALS DOPRNT = .false. endif C C Test for convergence. C * 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 .LE. ACC*ACC) THEN INFO=1 GOTO 70 END IF IF (DDOTG .GE. 0.0D0) THEN INFO=2 GOTO 70 END IF C C Test for termination due to no decrease in F. C IF (F .GE. FPREV) THEN IF (TOL .EQ. RELACC .OR. NACT .EQ. 0) THEN IF (DIFF .GT. 0.0D0) GOTO 20 END IF * print'(/'' DMLC05 869.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F if(abs(DDOTG) .le. EPS * abs(F)) then INFO = 2 else INFO=3 endif GOTO 70 END IF 20 continue DIFF=FPREV-F FPREV=F C C Test that more calls of DMLCFG are allowed. C IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) THEN INFO=8 GOTO 70 END IF C C Test whether to reduce TOL. C IF (TOL .GT. RELACC .AND. ITERC .GT. ITERK .AND. 1 0.1D0*RELAXF .GE. max(DIFF,-0.5D0*DDOTG)) then * print'('' DMLC05.. Debug.. return to reduce TOL.''/ * * '' ........ TOL='',g13.5,'', RELACC='',g13.5,'', INFO='', * * i5)', TOL, RELACC, INFO GOTO 70 endif c C Calculate the step along the search direction. C ITERC=ITERC+1 * print'(/'' DMLC05 895.. ITERC = '',i6)', ITERC C DMLC09 is LSRCH CALL DMLC09 (DMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB,DDOTG,F, * STEP,NFVALS, NFMAX,BRES, XL, XU, GMODE, W20) IF (STEP .EQ. 0.0D0) THEN * print'(/'' DMLC05 918.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F INFO=3 if(abs(DDOTG) .le. EPS * abs(F)) then INFO = 2 else SUM=0.0D0 DO 50 I=1,N 50 SUM=SUM+abs(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 .GE. 0.0D0) INFO=2 endif GOTO 70 END IF C C Revise XBIG. C DO 60 I=1,N 60 XBIG(I)=max(XBIG(I),abs(X(I))) C C Revise the second derivative approximation. C CALL DMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM) C C Add a constraint to the active set if it restricts the step. C IF (STEP .EQ. STEPCB) THEN K=IACT(INDXBD) IF (K .GT. M) THEN K=K-M IF (K .LE. N) THEN X(K)=XL(K) ELSE X(K-N)=XU(K-N) END IF END IF CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,XS,GS) END IF GOTO 10 C 70 continue FVAL = F return end c ================================================================== SUBROUTINE DMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG, 1 BRES,D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD, 2 GM,GMNEW,PARNEW,CGRAD,W20) C c CONRES C Calculate and partition the residuals of the inactive constraints, C and set the gradient vector when seeking feasibility. C W20 has length 4*N+1. Only locations 1 to N are used here. c ------------------------------------------------------------------ integer N,M,LDA,IACT(*),NACT,MEQL,MSAT,MTOT,INDXBD integer I, IDIFF, J, JM, K, KL, MDEG, MSATK double precision A(LDA,*),B(*), DDOTG, PAR(*) double precision G(*), U(*),XBIG(*),BRES(*) double precision D(*), GM(*),GMNEW(*),PARNEW(*), CGRAD(*) double precision RELACC, RES, RESABS, STEPCB, SUM, SUMRES double precision TEMP, TOL double precision W20(0:*), XL(*),XU(*),X(*) double precision Z(*), ZTG(*) c ------------------------------------------------------------------ IDIFF=MTOT-MSAT C IF (IDIFF .GT. 0) THEN DO 10 I=1,N 10 G(I)=0.0D0 SUMRES=0.0D0 END IF MSATK=MSAT MDEG=NACT MSAT=NACT KL=MEQL+1 DO 50 K=KL,MTOT J=IACT(K) C C Calculate the residual of the current constraint. C IF (J .LE. M) THEN RES=B(J) RESABS=abs(B(J)) DO 20 I=1,N RES=RES-X(I)*A(J,I) 20 RESABS=RESABS+abs(XBIG(I)*A(J,I)) ELSE JM=J-M IF (JM .LE. N) THEN RES=X(JM)-XL(JM) RESABS=abs(XBIG(JM))+abs(XL(JM)) ELSE JM=JM-N RES=XU(JM)-X(JM) RESABS=abs(XBIG(JM))+abs(XU(JM)) END IF END IF BRES(J)=RES C C Set TEMP to the relative residual. C TEMP=0.0D0 IF (RESABS .NE. 0.0D0) TEMP=RES/RESABS IF (K .GT. MSATK .AND. TEMP .LT. 0.0D0) THEN IF (TEMP+RELACC .GE. 0.0D0) THEN IF (J .LE. M) THEN SUM=abs(B(J)) DO 30 I=1,N 30 SUM=SUM+abs(X(I)*A(J,I)) ELSE JM=J-M IF (JM .LE. N) THEN SUM=abs(X(JM))+abs(XL(JM)) ELSE SUM=abs(X(JM-N))+abs(XU(JM-N)) END IF END IF IF (abs(RES) .LE. SUM*RELACC) TEMP=0.0D0 END IF END IF C C Place the residual in the appropriate position. C IF (K .LE. NACT) GOTO 50 IF (K .LE. MSATK .OR. TEMP .GE. 0.0D0) THEN MSAT=MSAT+1 IF (MSAT .LT. K) THEN IACT(K)=IACT(MSAT) END IF IF (TEMP .GT. TOL) THEN IACT(MSAT)=J ELSE MDEG=MDEG+1 IACT(MSAT)=IACT(MDEG) IACT(MDEG)=J END IF C C Update the gradient and SUMRES if the constraint is violated when C seeking feasibility. C ELSE IF (J .LE. M) THEN DO 40 I=1,N 40 G(I)=G(I)+A(J,I) ELSE J=J-M IF (J .LE. N) THEN G(J)=G(J)-1.0D0 ELSE G(J-N)=G(J-N)+1.0D0 END IF END IF SUMRES=SUMRES+abs(RES) END IF 50 CONTINUE C C Seek the next search direction unless DMLC06 was called from C DMLC04 [GETFES] and feasibility has been achieved. C STEPCB=0.0D0 IF (IDIFF .GT. 0 .AND. MSAT .EQ. MTOT) GOTO 60 c GETD CALL DMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC,DDOTG, * MEQL, MDEG,GM,GMNEW,PARNEW,CGRAD,W20) C C Calculate the (bound on the) step-length due to the constraints. C IF (DDOTG .LT. 0.0D0) THEN CALL DMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT, 1 MTOT,INDXBD) END IF IF (IDIFF .EQ. 0) SUMRES=DDOTG 60 RETURN END c ================================================================== SUBROUTINE DMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC, 1 DDOTG,MEQL,MDEG,GM,GMNEW,PARNEW,CGRAD,W20) C c GETD C Initialize GM and cycle backwards through the active set. C W20 has length 4*N+1. Only locations 0 to N are used here. c If W20(1) is nonnegative it means [D/S]MLC20 has been called and c has stored into W20(1:N) estimates of the errors in the gradient c values stored in G(1:N). c If W20(1) is negative the user is computing the gradients and c W20(1:N) is not being otherwise used. c This subr stores a value into W20(0). This will only be used if c [D/S]MLC20 is called. c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, MEQL, MDEG integer I, IZ, J, JM, K double precision A(LDA,*), CGRAD(*), D(*), DDOTG, DDOTGM double precision G(*), GM(*), GMNEW(*) double precision PAR(*), PARNEW(*), RELACC, SIZE, TEMP, U(*) double precision W20(0:*), Z(*), ZTG(*) c ------------------------------------------------------------------ 10 DO 20 I=1,N 20 GM(I)=G(I) K=NACT 30 IF (K .GT. 0) THEN C C Set TEMP to the next multiplier, but reduce the active set if C TEMP has an unacceptable sign. C TEMP=0.0D0 IZ=K DO 40 I=1,N TEMP=TEMP+Z(IZ)*GM(I) 40 IZ=IZ+N TEMP=TEMP*U(K) IF (K .GT. MEQL .AND. TEMP .GT. 0.0D0) THEN CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,K) GOTO 10 END IF C C Update GM using the multiplier that has just been calculated. C J=IACT(K) IF (J .LE. M) THEN DO 50 I=1,N 50 GM(I)=GM(I)-TEMP*A(J,I) ELSE JM=J-M IF (JM .LE. N) THEN GM(JM)=GM(JM)+TEMP ELSE GM(JM-N)=GM(JM-N)-TEMP END IF END IF PAR(K)=TEMP K=K-1 GOTO 30 END IF C C Calculate the search direction and DDOTG. C DDOTG=0.0D0 IF (NACT .LT. N) THEN c SDEGEN CALL DMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC, 1 DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD) IF (DDOTGM .LT. 0.0D0) THEN SIZE = 0.0D0 DO 60 I=1,N TEMP = D(I) * G(I) DDOTG=DDOTG + TEMP if (W20(1) .ge. 0.0D0) then SIZE = SIZE + abs(D(I)*W20(I)) else SIZE = SIZE + abs(TEMP) endif 60 continue if(W20(1) .lt. 0.0d0) SIZE = RELACC * SIZE if (DDOTG .ge. -SIZE) then DDOTG = 0.0D0 W20(0) = 0.0D0 else W20(0) = SIZE/abs(DDOTG) endif ELSE W20(0) = 0.0D0 END IF END IF RETURN END c ================================================================== SUBROUTINE DMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC, 1 DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD) C c SDEGEN C Calculate the search direction and branch if it is not downhill. c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, MEQL, MDEG integer I, IDROP, ITEST, IZ, J, JM, K, KU, MP, NP double precision A(LDA,*), CGRAD(*), D(*), DDOTGM, DTEST double precision GM(*), GMNEW(*) double precision PAR(*), PARNEW(*), U(*) double precision RATIO, RELACC, SUM, TEMP, THETA double precision Z(*), ZTG(*) c ------------------------------------------------------------------ MP=MEQL+1 DTEST=0.0D0 C 10 CALL DMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM) IF (DDOTGM .EQ. 0.0D0) GOTO 120 C C Branch if there is no need to consider any degenerate constraints. C The test gives termination if two consecutive additions to the C active set fail to increase the predicted new value of F. C IF (NACT .EQ. MDEG) GOTO 120 NP=NACT+1 SUM=0.0D0 DO 20 J=NP,N 20 SUM=SUM+ZTG(J)**2 IF (DTEST .GT. 0.0D0 .AND. SUM .GE. DTEST) THEN IF (ITEST .EQ. 1) GOTO 120 ITEST=1 ELSE DTEST=SUM ITEST=0 END IF C C Add a constraint to the active set if there are any significant C violations of degenerate constraints. C K=NACT CALL DMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,GMNEW,PARNEW, 1 CGRAD) IF (NACT .EQ. K) GOTO 120 PAR(NACT)=0.0D0 C C Calculate the new reduced gradient and Lagrange parameters. C 30 DO 40 I=1,N 40 GMNEW(I)=GM(I) K=NACT 50 TEMP=0.0D0 IZ=K DO 60 I=1,N TEMP=TEMP+Z(IZ)*GMNEW(I) 60 IZ=IZ+N TEMP=TEMP*U(K) PARNEW(K)=PAR(K)+TEMP IF (K .EQ. NACT) PARNEW(K)=min(PARNEW(K),0.0D0) J=IACT(K) IF (J .LE. M) THEN DO 70 I=1,N 70 GMNEW(I)=GMNEW(I)-TEMP*A(J,I) ELSE JM=J-M IF (JM .LE. N) THEN GMNEW(JM)=GMNEW(JM)+TEMP ELSE GMNEW(JM-N)=GMNEW(JM-N)-TEMP END IF END IF K=K-1 IF (K .GT. MEQL) GOTO 50 C C Set RATIO for linear interpolation between PAR and PARNEW. C RATIO=0.0D0 IF (MP .LT. NACT) THEN KU=NACT-1 DO 80 K=MP,KU IF (PARNEW(K) .GT. 0.0D0) THEN RATIO=PARNEW(K)/(PARNEW(K)-PAR(K)) IDROP=K END IF 80 CONTINUE END IF C C Apply the linear interpolation. C THETA=1.0D0-RATIO DO 90 K=MP,NACT 90 PAR(K)=min(THETA*PARNEW(K)+RATIO*PAR(K),0.0D0) DO 100 I=1,N 100 GM(I)=THETA*GMNEW(I)+RATIO*GM(I) C C Drop a constraint if RATIO is positive. C IF (RATIO .GT. 0.0D0) THEN CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) DO 110 K=IDROP,NACT 110 PAR(K)=PAR(K+1) GOTO 30 END IF C C Return if there is no freedom for a new search direction. C IF (NACT .LT. N) GOTO 10 DDOTGM=0.0D0 120 RETURN END c ================================================================== SUBROUTINE DMLC09 (DMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB, DDOTG, * F,STEP, NFVALS,NFMAX,GOPT, XL, XU, GMODE, W20) c LSRCH Line search c 4/18/91 CLL changed handling of ICOUNT and NFVALS. Added GMODE. c ------------------------------------------------------------------ C Subroutine Arguments c c Input: DMLCFG, N, D(), RELACC, STEPCB, DDOTG, NFMAX, XL(), XU() c Inout: X(), G(), F, NFVALS, GMODE, W20() c Out: GS(), STEP c Work space: XS(), GOPT() c c DMLCFG, N,X(),G(),D() c XS() [out] Copy of the initial X(). c GS() [out] Copy of the initial G(). c RELACC,STEPCB, DDOTG, c F,STEP, NFVALS,NFMAX,GOPT(), XL(), XU(), c GMODE [integer,inout] For use by DMLC20. c W20() [float,inout] For use by DMLC20. c ------------------------------------------------------------------ external DMLCFG integer N, NFVALS, NFMAX, GMODE integer I, ICOUNT double precision D(*) double precision DDOTG, DDOTGB, DGHGH, DGKNOT, DGLOW, DGMID, DGOPT double precision F, FBASE, FHGH, FLOW, FOPT double precision G(*), GOPT(*), GS(*) double precision RATIO, RELACC, RELINT, SBASE, SIZE, STEP, STEPCB double precision STPHGH, STPLOW, STPMIN, STPOPT, TEMP double precision W20(0:*), X(*), XL(*), XS(*), XU(*) logical HAVEG, KEEP, NOBACK c ------------------------------------------------------------------ C C Initialization. c * print'(/'' DMLC09 1290.. STEPCB='',g13.5/ * * '' .. D()='',3g13.5/(22x,3g13.5))', * * STEPCB, (D(I),I=1,N) RELINT=0.9D0 ICOUNT=0 RATIO=-1.0D0 DO 10 I=1,N XS(I)=X(I) GS(I)=G(I) GOPT(I)=G(I) IF (D(I) .NE. 0.0D0) THEN TEMP=abs(X(I)/D(I)) IF (RATIO .LT. 0.0D0 .OR. TEMP .LT. RATIO) RATIO=TEMP END IF 10 CONTINUE STEP=min(1.0D0,STEPCB) STPMIN=max(RELACC*RATIO,1.0D-12*STEP) STEP=max(STPMIN,STEP) SBASE=0.0D0 FBASE=F DDOTGB=DDOTG STPLOW=0.0D0 FLOW=F DGLOW=DDOTG STPHGH=0.0D0 STPOPT=0.0D0 FOPT=F DGOPT=abs(DDOTG) NOBACK = .false. * print'('' DMLC09 1318.. RATIO,STPMIN='',2g13.5/ * * '' DDOTG ='',g13.5)', * * RATIO, STPMIN, DDOTG C C Calculate another function and gradient value. C 20 continue if(NOBACK .and. STEP .lt. STPOPT) then * print*,' DMLC09 1355.. Ending DMLC09 on NOBACK.' go to 70 endif NOBACK = .false. DO 30 I=1,N 30 X(I)=XS(I)+STEP*D(I) C%% (*dmlcfg)( n, x, f, g, &haveg ); call DMLCFG(N, X, F, G, HAVEG) NFVALS = NFVALS+1 if(.not. HAVEG) CALL DMLC20(DMLCFG,N,X,F,G,XL,XU,GMODE,NFVALS,W20) ICOUNT=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.0D0 SIZE = 0.0d0 DO 40 I=1,N TEMP = D(I)*G(I) DGMID = DGMID + TEMP SIZE = SIZE + abs(TEMP) 40 continue SIZE = 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 .LE. FOPT) THEN if(ICOUNT .eq. 1 .and. DGMID .le. 0.0d0) then * print '('' DMLC09 1381.. Setting NOBACK on new test.'')' NOBACK = .true. KEEP = .true. else KEEP = F .LT. FOPT .OR. abs(DGMID) .LT. DGOPT endif if(KEEP) then STPOPT=STEP * print'(/'' DMLC09 1384.. New STPOPT='',g13.5)', STPOPT FOPT=F DO 50 I=1,N 50 GOPT(I)=G(I) DGOPT=abs(DGMID) END IF END IF IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) GOTO 70 C C Modify the bounds on the steplength or convergence. C IF (F .GE. FBASE+0.1D0*(STEP-SBASE)*DDOTGB) THEN IF (STPHGH .GT. 0.0D0 .OR. F .GT. FBASE .OR. DGMID .GT. 1 0.5D0*DDOTG) THEN STPHGH=STEP FHGH=F DGHGH=DGMID GOTO 60 END IF SBASE=STEP FBASE=F DDOTGB=DGMID END IF IF (DGMID .GE. 0.7D0*DDOTGB) GOTO 70 STPLOW=STEP FLOW=F DGLOW=DGMID 60 IF (STPHGH .GT. 0.0D0 .AND. STPLOW .GE. RELINT*STPHGH) GOTO 70 C C Calculate the next step length or end the iterations. C * print'(/'' DMLC09 1366.. ICOUNT, STEP='',i23,g23.15/ * * '' STPLOW, STPHGH='',2g23.15)', * * ICOUNT, STEP, STPLOW, STPHGH IF (STPHGH .EQ. 0.0D0) THEN IF (STEP .EQ. STEPCB) GOTO 70 TEMP=10.0D0 IF (DGMID .GT. 0.9D0*DDOTG) TEMP=DDOTG/(DDOTG-DGMID) STEP=min(TEMP*STEP,STEPCB) * print'('' DMLC09 1374.. To 20 with STEP='',g23.15)', STEP GOTO 20 ELSE IF (ICOUNT .EQ. 1 .OR. STPLOW .GT. 0.0D0) THEN DGKNOT=2.0D0*(FHGH-FLOW)/(STPHGH-STPLOW)-0.5D0*(DGLOW+DGHGH) IF (DGKNOT .GE. 0.0D0) THEN RATIO=max(0.1D0,0.5D0*DGLOW/(DGLOW-DGKNOT)) ELSE RATIO=(0.5D0*DGHGH-DGKNOT)/(DGHGH-DGKNOT) END IF STEP=STPLOW+RATIO*(STPHGH-STPLOW) * print'('' DMLC09 1384.. To 20 with STEP='',g23.15)', STEP GOTO 20 ELSE STEP=0.1D0*STEP if (STEP .GE. STPMIN) then * print'('' DMLC09 1389.. To 20 with STEP='',g23.15)', STEP GOTO 20 endif END IF C C Return from subroutine. C 70 continue IF (STEP .NE. STPOPT) THEN STEP=STPOPT F=FOPT DO 80 I=1,N X(I)=XS(I)+STEP*D(I) 80 G(I)=GOPT(I) END IF RETURN END c ================================================================== SUBROUTINE DMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,ZZDIAG, 1 GMNEW,CGRAD) c NEWCON c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, MDEG integer I, IADD, IZ, J, JM, JMV, K, KHIGH, NP double precision A(LDA,*), CGRAD(*), CVIOL, CVMAX, D(*), GMNEW(*) double precision RELACC, SAVABS, SAVSUM, SUM, SUMABS, SUMD, TEMP double precision U(*) double precision Z(*), ZZDIAG(*) c ------------------------------------------------------------------ C C Initialization. C NP=NACT+1 KHIGH=MDEG IZ=0 DO 20 I=1,N ZZDIAG(I)=0.0D0 DO 10 J=NP,N 10 ZZDIAG(I)=ZZDIAG(I)+Z(IZ+J)**2 20 IZ=IZ+N C C Calculate the scalar products of D with its constraints. C 30 CVMAX=0.0D0 DO 50 K=NP,KHIGH J=IACT(K) IF (J .LE. M) THEN SUM=0.0D0 SUMABS=0.0D0 SUMD=0.0D0 DO 40 I=1,N TEMP=D(I)*A(J,I) SUM=SUM+TEMP SUMABS=SUMABS+abs(TEMP) 40 SUMD=SUMD+ZZDIAG(I)*A(J,I)**2 ELSE JM=J-M IF (JM .LE. N) THEN SUM=-D(JM) ELSE JM=JM-N SUM=D(JM) END IF SUMABS=abs(SUM) SUMD=ZZDIAG(JM) END IF C C Pick out the most violated constraint, or return if the C violation is negligible. C IF (SUM .GT. RELACC*SUMABS) THEN CVIOL=SUM*SUM/SUMD IF (CVIOL .GT. CVMAX) THEN CVMAX=CVIOL IADD=K SAVSUM=SUM SAVABS=SUMABS END IF END IF 50 CONTINUE IF (CVMAX .LE. 0.0D0) GOTO 140 IF (NACT .EQ. 0) GOTO 120 C C Set GMNEW to the gradient of the most violated constraint. C J=IACT(IADD) IF (J .LE. M) THEN JMV=0 DO 60 I=1,N 60 GMNEW(I)=A(J,I) ELSE JMV=J-M DO 70 I=1,N 70 GMNEW(I)=0.0D0 IF (JMV .LE. N) THEN GMNEW(JMV)=-1.0D0 ELSE JMV=JMV-N GMNEW(JMV)=1.0D0 END IF END IF C C Modify GMNEW for the next active constraint. C K=NACT 80 TEMP=0.0D0 IZ=K DO 90 I=1,N TEMP=TEMP+Z(IZ)*GMNEW(I) 90 IZ=IZ+N TEMP=TEMP*U(K) J=IACT(K) IF (J .LE. M) THEN DO 100 I=1,N 100 GMNEW(I)=GMNEW(I)-TEMP*A(J,I) ELSE JM=J-M IF (JM .LE. N) THEN GMNEW(JM)=GMNEW(JM)+TEMP ELSE GMNEW(JM-N)=GMNEW(JM-N)-TEMP END IF END IF C C Revise the values of SAVSUM and SAVABS. C SUM=0.0D0 SUMABS=0.0D0 DO 110 I=1,N TEMP=D(I)*GMNEW(I) SUM=SUM+TEMP 110 SUMABS=SUMABS+abs(TEMP) SAVSUM=min(SAVSUM,SUM) SAVABS=max(SAVABS,SUMABS) K=K-1 IF (K .GE. 1) GOTO 80 C C Add the new constraint to the active set if the constraint C violation is still significant. C IF (JMV .GT. 0) D(JMV)=0.0D0 IF (SAVSUM .LE. RELACC*SAVABS) GOTO 130 120 K=NACT CALL DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IADD,GMNEW,CGRAD) IF (NACT .GT. K) GOTO 140 C C Seek another constraint violation. C IADD=NP 130 IF (NP .LT. KHIGH) THEN K=IACT(KHIGH) IACT(KHIGH)=IACT(IADD) IACT(IADD)=K KHIGH=KHIGH-1 GOTO 30 END IF 140 RETURN END c ================================================================== SUBROUTINE DMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG, 1 RELACC,TOL,MEQL) c SATACT c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, INFO integer I, IDROP, IZ, J, JX, K, MEQL double precision A(LDA,*),B(*) double precision RELACC, RES, RESABS, RESBIG, SAVEX, SCALE double precision TEMP, TEMPA, TOL, U(*) double precision X(*), XBIG(*), XL(*),XU(*), Z(*) c ------------------------------------------------------------------ IF (NACT .EQ. 0) GOTO 50 DO 30 K=1,NACT C C Calculate the next constraint residual. C J=IACT(K) IF (J .LE. M) THEN RES=B(J) RESABS=abs(B(J)) RESBIG=RESABS DO 10 I=1,N TEMPA=A(J,I) TEMP=TEMPA*X(I) RES=RES-TEMP RESABS=RESABS+abs(TEMP) 10 RESBIG=RESBIG+abs(TEMPA)*XBIG(I) ELSE JX=J-M IF (JX .LE. N) THEN RES=X(JX)-XL(JX) RESABS=abs(X(JX))+abs(XL(JX)) RESBIG=XBIG(JX)+abs(XL(JX)) SAVEX=XL(JX) ELSE JX=JX-N RES=XU(JX)-X(JX) RESABS=abs(X(JX))+abs(XU(JX)) RESBIG=XBIG(JX)+abs(XU(JX)) SAVEX=XU(JX) END IF END IF C C Shift X if necessary. C IF (RES .NE. 0.0D0) THEN TEMP=RES/RESABS IF (K .LE. MEQL) TEMP=-abs(TEMP) IF (TOL .EQ. RELACC .OR. TEMP+RELACC .LT. 0.0D0) THEN INFO=1 SCALE=RES*U(K) IZ=K DO 20 I=1,N X(I)=X(I)+SCALE*Z(IZ) IZ=IZ+N 20 XBIG(I)=max(XBIG(I),abs(X(I))) IF (J .GT. M) X(JX)=SAVEX C C Else flag a constraint deletion if necessary. C ELSE IF (RES/RESBIG .GT. TOL) THEN IACT(K)=-IACT(K) END IF END IF 30 CONTINUE C C Delete any flagged constraints and then return. C IDROP=NACT 40 IF (IACT(IDROP) .LT. 0) THEN IACT(IDROP)=-IACT(IDROP) CALL DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) END IF IDROP=IDROP-1 IF (IDROP .GT. MEQL) GOTO 40 50 RETURN END c ================================================================== SUBROUTINE DMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,ZTC, 1 CGRAD) c ADDCON c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, INDXBD integer I, ICON, INEWBD, IPIV, IZ, IZNBD, J, JP, NP double precision A(LDA,*), CGRAD(*) double precision RELACC, SUM, SUMABS, TEMP, TEMPA, TEMPB double precision U(*), WCOS, WPIV, WSIN double precision Z(*), ZTC(*) c ------------------------------------------------------------------ NP=NACT+1 ICON=IACT(INDXBD) IACT(INDXBD)=IACT(NP) IACT(NP)=ICON C C Form ZTC when the new constraint is a bound. C IF (ICON .GT. M) THEN INEWBD=ICON-M IF (INEWBD .LE. N) THEN TEMP=-1.0D0 ELSE INEWBD=INEWBD-N TEMP=1.0D0 END IF IZNBD=INEWBD*N-N DO 10 J=1,N 10 ZTC(J)=TEMP*Z(IZNBD+J) C C Else form ZTC for an ordinary constraint. C ELSE DO 20 I=1,N 20 CGRAD(I)=A(ICON,I) DO 30 J=1,N ZTC(J)=0.0D0 IZ=J DO 30 I=1,N ZTC(J)=ZTC(J)+Z(IZ)*CGRAD(I) 30 IZ=IZ+N END IF C C Find any Givens rotations to apply to the last columns of Z. C J=N 40 JP=J J=J-1 IF (J .GT. NACT) THEN IF (ZTC(JP) .EQ. 0.0D0) GOTO 40 IF (abs(ZTC(JP)) .LE. RELACC*abs(ZTC(J))) THEN TEMP=abs(ZTC(J)) ELSE IF (abs(ZTC(J)) .LE. RELACC*abs(ZTC(JP))) THEN TEMP=abs(ZTC(JP)) ELSE TEMP=abs(ZTC(JP))*sqrt(1.0D0+(ZTC(J)/ZTC(JP))**2) END IF WCOS=ZTC(J)/TEMP WSIN=ZTC(JP)/TEMP ZTC(J)=TEMP C C Apply the rotation when the new constraint is a bound. C IZ=J IF (ICON .GT. M) THEN DO 50 I=1,N TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) Z(IZ+1)=TEMP 50 IZ=IZ+N Z(IZNBD+JP)=0.0D0 C C Else apply the rotation for an ordinary constraint. C ELSE WPIV=0.0D0 DO 60 I=1,N TEMPA=WCOS*Z(IZ+1) TEMPB=WSIN*Z(IZ) TEMP=abs(CGRAD(I))*(abs(TEMPA)+abs(TEMPB)) IF (TEMP .GT. WPIV) THEN WPIV=TEMP IPIV=I END IF Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) Z(IZ+1)=TEMPA-TEMPB 60 IZ=IZ+N C C Ensure orthogonality of Z(.,JP) to CGRAD. C SUM=0.0D0 IZ=JP DO 70 I=1,N SUM=SUM+Z(IZ)*CGRAD(I) 70 IZ=IZ+N IF (SUM .NE. 0.0D0) THEN IZ=IPIV*N-N+JP Z(IZ)=Z(IZ)-SUM/CGRAD(IPIV) END IF END IF GO TO 40 END IF C C Test for linear independence in the proposed new active set. C IF (ZTC(NP) .EQ. 0.0D0) GOTO 90 IF (ICON .LE. M) THEN SUM=0.0D0 SUMABS=0.0D0 IZ=NP DO 80 I=1,N TEMP=Z(IZ)*CGRAD(I) SUM=SUM+TEMP SUMABS=SUMABS+abs(TEMP) 80 IZ=IZ+N IF (abs(SUM) .LE. RELACC*SUMABS) GOTO 90 END IF C C Set the new diagonal element of U() and return. C U(NP)=1.0D0/ZTC(NP) NACT=NP 90 RETURN END c ================================================================== SUBROUTINE DMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC,TOL, 1 MEQL) C C ADJTOL: Change the tolerance TOL to a smaller value, if possible. c May also recompute XBIG(). c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, MEQL integer I, J, JM, K, KL double precision A(LDA,*), B(*) double precision RELACC, RES, RESABS, TOL, VIOL double precision X(*), XBIG(*), XL(*), XU(*) c ------------------------------------------------------------------ C Set VIOL to the greatest relative constraint residual of the first C NACT constraints. VIOL=0.0D0 IF (NACT .GT. MEQL) THEN KL=MEQL+1 DO 20 K=KL,NACT J=IACT(K) IF (J .LE. M) THEN RES=B(J) RESABS=abs(B(J)) DO 10 I=1,N RES=RES-A(J,I)*X(I) 10 RESABS=RESABS+abs(A(J,I)*XBIG(I)) ELSE JM=J-M IF (JM .LE. N) THEN RES=X(JM)-XL(JM) RESABS=XBIG(JM)+abs(XL(JM)) ELSE JM=JM-N RES=XU(JM)-X(JM) RESABS=XBIG(JM)+abs(XU(JM)) END IF END IF IF (RES .GT. 0.0D0) VIOL=max(VIOL,RES/RESABS) 20 CONTINUE END IF C C Adjust TOL. C TOL=0.1D0*min(TOL,VIOL) IF (TOL .LE. RELACC+RELACC) THEN TOL=RELACC DO 30 I=1,N 30 XBIG(I)=abs(X(I)) END IF RETURN END c ================================================================== SUBROUTINE DMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) C C DELCON: Cycle through the constraint exchanges that are needed. c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, IDROP integer I, IBD, ICON, IPIV, ISAVE, IZ, IZBD, J, JP, NM double precision A(LDA,*), DENOM double precision RELACC, RJJP, SUM, TEMP, TEMPA, TEMPB double precision U(*), UJP, WCOS, WPIV, WSIN double precision Z(*) c ------------------------------------------------------------------ NM=NACT-1 IF (IDROP .EQ. NACT) GOTO 60 ISAVE=IACT(IDROP) C C Cycle through the constraint exchanges that are needed. C DO 50 J=IDROP,NM JP=J+1 ICON=IACT(JP) IACT(J)=ICON C C Calculate the (J,JP) element of R. C IF (ICON .LE. M) THEN RJJP=0.0D0 IZ=J DO 10 I=1,N RJJP=RJJP+Z(IZ)*A(ICON,I) 10 IZ=IZ+N ELSE IBD=ICON-M IF (IBD .LE. N) THEN IZBD=IBD*N-N RJJP=-Z(IZBD+J) ELSE IBD=IBD-N IZBD=IBD*N-N RJJP=Z(IZBD+J) END IF END IF C C Calculate the parameters of the next rotation. C UJP=U(JP) TEMP=RJJP*UJP DENOM=abs(TEMP) IF (DENOM*RELACC .LT. 1.0D0) DENOM=sqrt(1.0D0+DENOM*DENOM) WCOS=TEMP/DENOM WSIN=1.0D0/DENOM C C Rotate Z when a bound constraint is promoted. C IZ=J IF (ICON .GT. M) THEN DO 20 I=1,N TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) Z(IZ+1)=TEMP 20 IZ=IZ+N Z(IZBD+JP)=0.0D0 C C Rotate Z when an ordinary constraint is promoted. C ELSE WPIV=0.0D0 DO 30 I=1,N TEMPA=WCOS*Z(IZ+1) TEMPB=WSIN*Z(IZ) TEMP=abs(A(ICON,I))*(abs(TEMPA)+abs(TEMPB)) IF (TEMP .GT. WPIV) THEN WPIV=TEMP IPIV=I END IF Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) Z(IZ+1)=TEMPA-TEMPB 30 IZ=IZ+N C C Ensure orthogonality to promoted constraint. C SUM=0.0D0 IZ=JP DO 40 I=1,N SUM=SUM+Z(IZ)*A(ICON,I) 40 IZ=IZ+N IF (SUM .NE. 0.0D0) THEN IZ=IPIV*N-N+JP Z(IZ)=Z(IZ)-SUM/A(ICON,IPIV) END IF END IF C C Set the new diagonal elements of U. C U(JP)=-DENOM*U(J) U(J)=UJP/DENOM 50 CONTINUE C C Return. C IACT(NACT)=ISAVE 60 NACT=NM RETURN END c ================================================================== SUBROUTINE DMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC) c INITZU: Initialize RELACC, Z(), U(). Adjust X(). Set XBIG(). c 1990-07-20 CLL Changed error handling and setting of INFO. c Expect INFO to be set to 1 before entering this subr. This subr c sets INFO = 4 if it finds XL(j) > XU(j) for some j. c ------------------------------------------------------------------ external D1MACH integer N, M, IACT(*), MEQL, INFO integer I, IZ, J, JACT, NN double precision D1MACH double precision RELACC double precision U(*) double precision XL(*),XU(*),X(*), XBIG(*), Z(*) c ------------------------------------------------------------------ C C Set RELACC. c July 1990 CLL Commented out following 6 lines and replaced with c reference to D1MACH(4) which is the smallest x such that the c stored value of 1+x is greater than 1. Efforts to determine c such an x with portable code such as the following 6 lines have c generally eventually failed on some new computer. Also using c D1MACH allows adjustments to be made for known deficiencies c in particular computers, for example for the Cray X/MP & Y/MP. c * 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.0D0 * D1MACH(4) C C Seek bound inconsistencies and bound equality constraints. C MEQL=0 DO 20 J=1,N IF (XL(J) .GT. XU(J)) then INFO = 4 call IERM1('DMLC15',INFO,0,'Bad bounds: XL(j) > XU(j)', * 'j',J,',') call DERV1('XL(j)',XL(J),',') call DERV1('XU(j)',XU(J),'.') return endif IF (XL(J) .EQ. XU(J)) MEQL=MEQL+1 20 CONTINUE C C Initialize U, Z and XBIG. C JACT=0 NN=N*N DO 30 I=1,NN 30 Z(I)=0.0D0 IZ=0 DO 40 I=1,N IF (XL(I) .EQ. XU(I)) THEN X(I)=XU(I) JACT=JACT+1 U(JACT)=1.0D0 IACT(JACT)=I+M+N J=JACT ELSE J=I+MEQL-JACT END IF Z(IZ+J)=1.0D0 IZ=IZ+N 40 XBIG(I)=abs(X(I)) RETURN END c ================================================================== SUBROUTINE DMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES, * RELAXF, MEQL,SSQKT,PARW,RESKTW) C C Calculate the Lagrange parameters and the residual vector. c c Input: N, M, A(), LDA, IACT, NACT, G(), Z(), U(), BRES(), MEQL c Output: PAR(), RESKT, RELAXF, SSQKT c Work space: PARW(), RESKTW() c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), NACT, MEQL integer I, ICASE, IZ, J, JM, K, KK, KL double precision A(LDA,*), BRES(*), G(*) double precision PAR(*), PARW(*) double precision RELAXF, RESKT(*), RESKTW(*) double precision SSQKT, SSQKTW, TEMP, U(*), Z(*) c ------------------------------------------------------------------ DO 10 I=1,N 10 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 .GT. 0) THEN ICASE=0 20 DO 50 KK=1,NACT K=NACT+1-KK J=IACT(K) TEMP=0.0D0 IZ=K DO 30 I=1,N TEMP=TEMP+Z(IZ)*RESKT(I) * print'('' DMLC16 1991.. Z(),RESKT(),TEMP='',3g13.5)', * * Z(IZ),RESKT(I),TEMP 30 IZ=IZ+N TEMP=TEMP*U(K) * print'('' DMLC16 1995.. U(K),TEMP='', 2g13.5)', U(K), TEMP IF (ICASE .EQ. 0) PAR(K)=0.0D0 IF (K .LE. MEQL .OR. PAR(K)+TEMP .LT. 0.0D0) THEN PAR(K)=PAR(K)+TEMP * print'('' DMLC16 1999.. TEMP,PAR(K)='',2g13.5)', * * TEMP,PAR(K) ELSE TEMP=-PAR(K) PAR(K)=0.0D0 END IF IF (TEMP .NE. 0.0D0) THEN IF (J .LE. M) THEN DO 40 I=1,N RESKT(I)=RESKT(I)-TEMP*A(J,I) * print'('' DMLC16 2009.. TEMP,A='',2g13.5/ * * '' RESKT(I) ='',g13.5)', * * TEMP, A(J,I), RESKT(I) 40 continue ELSE JM=J-M IF (JM .LE. N) THEN RESKT(JM)=RESKT(JM)+TEMP * print'('' DMLC16 2018.. TEMP,RESKT(JM)='',2g13.5)', * * TEMP,RESKT(JM) ELSE RESKT(JM-N)=RESKT(JM-N)-TEMP * print'('' DMLC16 2022.. TEMP,RESKT(JM-N)='',2g13.5)', * * TEMP,RESKT(JM-N) END IF END IF END IF 50 CONTINUE C C Calculate the sum of squares of the KT residual vector. C SSQKT=0.0D0 IF (NACT .EQ. N) GOTO 130 * print'(/'' DMLC16 2014.. RESKT()='',4g12.4/(23x,4g12.4))', * * (RESKT(I),I=1,N) DO 60 I=1,N 60 SSQKT=SSQKT+RESKT(I)**2 * print'('' DMLC16 2018.. SSQKT='',g13.5)', SSQKT C C Apply iterative refinement to the residual vector. C IF (ICASE .EQ. 0) THEN ICASE=1 DO 70 K=1,NACT 70 PARW(K)=PAR(K) DO 80 I=1,N 80 RESKTW(I)=RESKT(I) SSQKTW=SSQKT GOTO 20 END IF C C Undo the iterative refinement if it does not reduce SSQKT. C IF (SSQKTW .LT. SSQKT) THEN DO 90 K=1,NACT 90 PAR(K)=PARW(K) DO 100 I=1,N 100 RESKT(I)=RESKTW(I) SSQKT=SSQKTW END IF C C Calculate SSQKT when there are no active constraints. C ELSE SSQKT=0.0D0 DO 110 I=1,N 110 SSQKT=SSQKT+G(I)**2 END IF C C Predict the reduction in F if one corrects any positive residuals C of active inequality constraints. C RELAXF=0.0D0 IF (MEQL .LT. NACT) THEN KL=MEQL+1 DO 120 K=KL,NACT J=IACT(K) IF (BRES(J) .GT. 0.0D0) THEN RELAXF=RELAXF-PAR(K)*BRES(J) END IF 120 CONTINUE END IF 130 RETURN END c ================================================================== SUBROUTINE DMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM) c SDIRN: Compute a search direction. c ------------------------------------------------------------------ integer N, NACT integer I, IZ, J, NP double precision D(*), DDOTGM, GM(*) double precision RELACC, SUM, SUMABS, TEMP double precision Z(*), ZTG(*) c ------------------------------------------------------------------ DDOTGM=0.0D0 IF (NACT .GE. N) GOTO 60 C C Premultiply GM by the transpose of Z. C NP=NACT+1 DO 20 J=NP,N SUM=0.0D0 SUMABS=0.0D0 IZ=J DO 10 I=1,N TEMP=Z(IZ)*GM(I) SUM=SUM+TEMP SUMABS=SUMABS+abs(TEMP) 10 IZ=IZ+N IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0D0 20 ZTG(J)=SUM C C Form D by premultiplying ZTG by -Z. C IZ=0 DO 40 I=1,N SUM=0.0D0 SUMABS=0.0D0 DO 30 J=NP,N TEMP=Z(IZ+J)*ZTG(J) SUM=SUM-TEMP 30 SUMABS=SUMABS+abs(TEMP) IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0D0 D(I)=SUM 40 IZ=IZ+N C C Test that the search direction is downhill. C SUMABS=0.0D0 DO 50 I=1,N TEMP=D(I)*GM(I) DDOTGM=DDOTGM+TEMP 50 SUMABS=SUMABS+abs(TEMP) IF (DDOTGM+RELACC*SUMABS .GE. 0.0D0) DDOTGM=0.0D0 60 RETURN END c ================================================================== SUBROUTINE DMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT, 1 MTOT,INDXBD) c STEPBD: C Set steps to constraint boundaries and find the least positive C one. c c>> 1990-04-02 CLL Changes to avoid jumping into scope of Block If. c ------------------------------------------------------------------ integer N, M, LDA, IACT(*), MDEG, MSAT, MTOT, INDXBD integer I, IFLAG, J, JM, K, KL double precision A(LDA,*),BRES(*),D(*), DDOTG, SP, STEPCB, TEMP c ------------------------------------------------------------------ IFLAG=0 STEPCB=0.0D0 INDXBD=0 K=MDEG 10 K=K+1 IF (K .gt. MTOT) go to 30 c c Use remote code to compute scalar product. go to 80 20 continue C C Set BRES(J) to indicate the status of the j-th constraint. C IF (SP*BRES(J) .LE. 0.0D0) THEN BRES(J)=0.0D0 ELSE BRES(J)=BRES(J)/SP IF (STEPCB .EQ. 0.0D0 .OR. BRES(J) .LT. STEPCB) THEN STEPCB=BRES(J) INDXBD=K END IF END IF GO TO 10 30 continue C C Try to pass through the boundary of a violated constraint. C 40 continue IF (INDXBD .le. MSAT) go to 70 IFLAG=1 K=INDXBD c Use remote code to compute scalar product. go to 80 50 continue MSAT=MSAT+1 IACT(INDXBD)=IACT(MSAT) IACT(MSAT)=J BRES(J)=0.0D0 INDXBD=MSAT DDOTG=DDOTG-SP IF (DDOTG .LT. 0.0D0 .AND. MSAT .LT. MTOT) THEN C C Seek the next constraint boundary along the search direction. C TEMP=0.0D0 KL=MDEG+1 DO 60 K=KL,MTOT J=IACT(K) IF (BRES(J) .GT. 0.0D0) THEN IF (TEMP .EQ. 0.0D0 .OR. BRES(J) .LT. TEMP) THEN TEMP=BRES(J) INDXBD=K END IF END IF 60 CONTINUE IF (TEMP .GT. 0.0D0) THEN STEPCB=TEMP GOTO 40 END IF END IF 70 continue RETURN c ------------------------------------------------------------------ c This is a remote block of code to compute the C scalar product of D with the current constraint normal. C 80 J=IACT(K) IF (J .LE. M) THEN SP=0.0D0 DO 90 I=1,N 90 SP=SP+D(I)*A(J,I) ELSE JM=J-M IF (JM .LE. N) THEN SP=-D(JM) ELSE SP=D(JM-N) END IF END IF C Return from remote block is selected by IFLAG. IF (IFLAG .EQ. 0) then go to 20 ELSE go to 50 ENDIF END c ================================================================== SUBROUTINE DMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM) C C Test if there is sufficient convexity for the update. c ------------------------------------------------------------------ integer N,NACT integer I, IZ, K, KP, KM, NP double precision DD, DG, G(*), GS(*), SUM, TEMP, WCOS, WSIN double precision X(*), XS(*), Z(*),ZTG(*), ZZNORM c ------------------------------------------------------------------ DD=0.0D0 DG=0.0D0 TEMP=0.0D0 DO 10 I=1,N XS(I)=X(I)-XS(I) DD=DD+XS(I)**2 TEMP=TEMP+GS(I)*XS(I) GS(I)=G(I)-GS(I) 10 DG=DG+GS(I)*XS(I) IF (DG .LT. 0.1D0*abs(TEMP)) GOTO 90 C C Transform the Z matrix. C K=N 20 KP=K K=K-1 IF (K .GT. NACT) THEN IF (ZTG(KP) .EQ. 0.0D0) GOTO 20 TEMP=abs(ZTG(KP))*sqrt(1.0D0+(ZTG(K)/ZTG(KP))**2) WCOS=ZTG(K)/TEMP WSIN=ZTG(KP)/TEMP ZTG(K)=TEMP IZ=K DO 30 I=1,N TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) Z(IZ+1)=TEMP 30 IZ=IZ+N GOTO 20 END IF C C Update the value of ZZNORM. C IF (ZZNORM .LT. 0.0D0) THEN ZZNORM=DD/DG ELSE TEMP=sqrt(ZZNORM*DD/DG) ZZNORM=min(ZZNORM,TEMP) ZZNORM=max(ZZNORM,0.1D0*TEMP) END IF C C Complete the updating of Z. C NP=NACT+1 TEMP=sqrt(DG) IZ=NP DO 40 I=1,N Z(IZ)=XS(I)/TEMP 40 IZ=IZ+N IF (NP .LT. N) THEN KM=NP+1 DO 80 K=KM,N TEMP=0.0D0 IZ=K DO 50 I=1,N TEMP=TEMP+GS(I)*Z(IZ) 50 IZ=IZ+N TEMP=TEMP/DG SUM=0.0D0 IZ=K DO 60 I=1,N Z(IZ)=Z(IZ)-TEMP*XS(I) SUM=SUM+Z(IZ)**2 60 IZ=IZ+N IF (SUM .LT. ZZNORM) THEN TEMP=sqrt(ZZNORM/SUM) IZ=K DO 70 I=1,N Z(IZ)=TEMP*Z(IZ) 70 IZ=IZ+N END IF 80 CONTINUE END IF 90 RETURN END c ================================================================== SUBROUTINE DMLC20 (DMLCFG, N, X, F, G, XL, XU, GMODE, NFVALS, W) c>> 1991-04-15 FTK & CLL Made FIRST an argument of DMLC20 c>> 1990-09-11 Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA. c c This subr computes a finite difference approximation to a gradient c vector. Values of F are not computed outside the bounds given by XL() c and XU(). For any I for which XL(I) = XU(I), G(I) will be set to zero c and computation to estimate G(I) will be skipped. c On entry F must have already been evaluated at X(). c The algorithm used is based on the following observations. c 1. Since the gradient vector is being used in an iteration, there is c no point to computing it more accurately than is useful in the c iteration. An extra N function values to get better error c estimates is not likely to pay for itself. c 2. If the delta x used for the difference is in the direction of the c next move, then discretization error in approximating the gradient c is likely to help rather than hurt. c 3. The sign of the change in X(I) for the next move is likely to be c the same as it was on the last move. c 4. Because of 3 and 4, the increment used is larger than would be c suggested from a consideration of round off and discretization c errors. c 5. When the distance between values on successive iterations gets very c small it is useful to have second derivative information in order c to balance the effects of round off and discretization errors. c c *************************** Variables ******************************** c c D The difference between the current X(I) and the value on the c previous iteration. c DXN L2 norm of the move just taken. c DMLCFG [in, subroutine name] Name of user-provided subroutine for c computing function values. When called from this subr, DMLCFG must c compute F, set HAVEG = .false., and not store into G(). c EPSR Save variable = relative machine precision. c F [in, float] On entry F will contain the function value c evaluated at the given X(). c FAC Save variable = sqrt(relative machine precision) c FACX Save variable, 4. * FAC. c G() [out, float] On return G() will contain an approximation of c the N-dimensional gradient of f with respect to x at X() computed c using one-sided finite differences. c GMODE [inout, integer] Has value 0, or 1. c 0 means initial entry. Compute 1-sided differences, assuming no c saved info, and set GMODE to 1. c 1 Compute differences, using saved info. c H Increment in X(I) used for computing gradient. c HB Lower bound desired for value of abs(H). c KOUNT Counts interations. If multiple of 10, two sided differences c are computed. c LGL W(LGL+I) contains the value of G(I) from previous iteration. c LNOPP Value of NOPP for start of next iteration. c LPP W(LPP+I) contains estimate of the second derivative (d/dv)G(I), c where v is the variable consisting of the linear combinations of c X's that gave the last move. c LXL W(LXL+I) contains X() from the previous iteration. c MODE1 .true. if have 1-sided difference, = .false. if have 2-sided. c N [in, integer] Number of components in X(), G(), XL(), & XU(). c NFVALS [inout, integer] Counter of number of calls to DMLCFG. c NOPP Logical variable set .true. if second derivative information is c not available. c TP Temporary storage. c X() [in, work, float] On entry contains the current N-dimensional c parameter vector. This subr will alter components in X() but will c reset X() to its original contents on return. c XI Current base value for X(I) when computing F for a gradient. c XL() [in, float] Lower bounds for components of X(). c XNEW Current value for X(I) in computing F for a gradient. Also c used for temporary storage of 1-sided gradient. c XU() [in, float] Upper bounds for components of X(). c W() [inout, float] Work array of dimension at least 4N+1. c W(0) and W(1) are initialized to -1.0 in [D/S]MLC02. If this c subr is never called, i.e., if the user provides computed c gradients, then W(1) will remain at this initial value. c If this subr is called, it will (generallly) set W(1) positive. c [D/S]MLC07 will (generally) set W(0) positive for subsequent c testing when and if this subr is called. c W(1:N) contains an estimate of the error in G(I) for use in c [D/S]MLC07. c W(N+1 : 4*N) is used as temporary work space in this subr. See the c usage of the indices LGL, LPP, LXL. c ------------------------------------------------------------------ c c ********************** Variable Declarations ************************* c external D1MACH, DMLCFG integer GMODE, LGL, LPP, LXL, N integer I, KDEBUG, KOUNT, NFVALS double precision D1MACH double precision D, DXN, EPSR double precision F, FAC, FACX, FNEW, G(N), H, HB, TP double precision X(N), XI, XL(N), XNEW, XU(N), W(0:*) logical MODE1, NOPP, LNOPP c%% LOGICAL32 haveg; logical HAVEG save EPSR, FAC, FACX, KOUNT, LNOPP save KDEBUG data KDEBUG / 0 / data FAC / 0.0D0 /, KOUNT / 0 / 1000 format ('0Gradient #=', I3, ' F=', 1p,e24.17, ' DXN=', 1 1p,e15.8 / ' I', 11x, 'X', 14X, 'DX', 10x, 'G', 10x, 'ERR', 2 9X, 'H', 9X, 'dG/dV') c c ********************** Start of Executable Code ********************** c if(FAC .eq. 0.0D0) then c Get machine constants EPSR = D1MACH(4) FAC = sqrt(EPSR) FACX = 4.D0*FAC endif c Set up base locations in W() LXL = N LGL = LXL + LXL LPP = LGL + LXL DXN = 0.D0 if (GMODE .eq. 0) then GMODE = 1 KOUNT = 1 NOPP = .true. else KOUNT = KOUNT + 1 if (mod(KOUNT,10) .eq. 1) then if (W(0) .gt. 1.D-5) KOUNT = 0 else if (W(0) .gt. 0.01D0) KOUNT = 0 end if do 10 I = 1, LXL DXN = DXN + (X(I) - W(LXL+I))**2 10 continue DXN = sqrt(DXN) NOPP = LNOPP end if LNOPP = .false. if (KDEBUG .ne. 0) then print 1000, KDEBUG, F, DXN KDEBUG = KDEBUG + 1 end if do 40 I = 1, LXL XI = X(I) D = XI - W(LXL+I) if (NOPP) then HB = FAC * (1.D0 + abs(XI)) else HB = FAC * (sqrt(abs(F/W(LPP+I))) + FACX*(abs(X(I))+FACX)) end if 20 if ((W(LPP+I) .lt. 0.D0) .or. (mod(KOUNT, 10) .eq. 0)) then if ((XI + HB .le. XU(I)) .and. (XI - HB .ge. XL(I))) then MODE1 = .false. H = sign(min(100.D0*HB, min(XU(I)-XI, XI-XL(I))), D) go to 30 end if end if MODE1 = .true. H = sign(max(HB, abs(1.D-4*D)), D) 30 XNEW = XI + H c Adjust H to keep XNEW within bounds. if (XNEW .gt. XU(I)) then XNEW = XU(I) if (XNEW - XI .lt. HB) then if (XI - XL(I) .gt. XNEW - XI) XNEW = max(XL(I), XI-HB) end if H = XNEW - XI else if (XNEW .lt. XL(I)) then XNEW = XL(I) if (XI - XNEW .lt. HB) then if (XU(I) - XI .gt. XI - XNEW) XNEW = min(XU(I), XI+HB) end if H = XNEW - XI end if if(H .eq. 0.D0) then G(I) = 0.D0 W(I) = 0.D0 else X(I) = XNEW C%% (*dmlcfg)( lxl, x, &fnew, g, &haveg ); CALL DMLCFG(LXL, X, FNEW, G, HAVEG) NFVALS = NFVALS+1 G(I) = (FNEW - F) / H W(I) = EPSR * (abs(F) + abs(FNEW)) / abs(H) if (.not. MODE1) then XNEW = XI - H X(I) = XNEW C%% (*dmlcfg)( lxl, x, &fnew, g, &haveg ); CALL DMLCFG(LXL, X, FNEW, G, HAVEG) NFVALS = NFVALS+1 XNEW = (FNEW - F) / (XNEW - XI) W(I) = W(I) + 5.D-4 * abs(XNEW - G(I)) G(I) = .5D0 * (G(I) + XNEW ) end if X(I) = XI end if if (DXN .eq. 0.D0) then LNOPP = .true. if (MODE1) W(I) = 0.D0 else W(LPP+I) = max(0.D0, abs((G(I) - W(LGL+I)) - 10.D0 * 1 abs(W(I))) / DXN) + FAC * abs(F) TP = 10.D0 * abs(H) * W(LPP+I) if (MODE1) then if (TP .gt. abs(G(I))) then if (mod(KOUNT, 10) .ne. 0) then KOUNT = 0 go to 20 end if end if W(I) = W(I) + TP else if (abs(XNEW-G(I)) .gt. min(10.D0*TP, 1.D-3 * abs(G(I)))) 1 W(LPP+I) = -W(LPP+I) end if end if if (KDEBUG .ne. 0) then print '(I3, 1P,E22.14, E10.2, E12.4, E12.4, E11.3, E10.2)', * I,X(I),X(I)-W(LXL+I),G(I),W(I),H,W(LPP+I) end if W(LXL+I) = XI W(LGL+I) = G(I) 40 continue return end c ================================================================== SUBROUTINE DMLC21(IPMORE, ALL, N, ITERC, NFVALS, * F, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) c>> 1990-07-11 CLL c This subr prints results (first feasible, intermediate, and c final). c ------------------------------------------------------------------ c Arguments c c IPMORE [in, integer] = 0 or 1. 0 means print ITERC, NFVALS, F, X(), c and G(). 1 means also print IACT(), PAR(), and RESKT(). c However see ALL for an exception. c ALL [in, logical] False means G(), PAR(), and RESKT() might not c contain valid data and therefore will not be printed. c N, ITERC, NFVALS, c F, TOL, X(), G(), NACT, IACT(), PAR(), RESKT(), c ENRMKT [in] Euclidian norm of RESKT(). c ------------------------------------------------------------------ c%% long int i, k; integer I integer IPMORE, N, ITERC, NFVALS, NACT, IACT(*) double precision ENRMKT, F, X(N), G(N), PAR(*), RESKT(N), TOL logical ALL c ------------------------------------------------------------------ print'(/1x,i6,'' iterations. '',i6,'' function evals. F = '', * g13.5/1x,'' TOL ='',g13.5,'' Norm of RESKT ='',g13.5)', * ITERC, NFVALS, F, TOL, ENRMKT C%% printf( "\n X =" ); C%% for (i = 0; i < n; i+=5){ C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) C%% printf( "%13.5g", x[k] ); C%% if (i < n - 4) printf ( "\n " );} print'('' X ='',5g13.5/(8x,5g13.5))', (X(I),I=1,N) if (IPMORE .gt. 0) then if(ALL) then C%% printf( "\n G =" ); C%% for (i = 0; i < n; i+=5){ C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) C%% printf( "%13.5g", g[k] ); C%% if (i < n - 4) printf ( "\n " );} print'(a,5g13.5/(8x,5g13.5))',' G =',(G(I),I=1,N) if (NACT .eq. N) then print*,' Kuhn-Tucker residual vector is zero.' else C%% printf( "\n RESKT =" ); C%% for (i = 0; i < n; i+=5){ C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) C%% printf( "%13.5g", reskt[k] ); C%% if (i < n - 4) printf ( "\n " );} print'('' RESKT ='',5g13.5/(8x,5g13.5))', (RESKT(I),I=1,N) end if end if if (NACT .eq. 0) then print*,' No active constraints.' else C%% printf( "\n IACT =" ); C%% for (i = 0; i < nact; i+=12){ C%% for (k = i; k < (i < nact - 11 ? i+12 : nact); k++) C%% printf( "%5ld", iact[k] ); C%% if (i < nact - 11) printf ( "\n " );} print'(a,12i5/(8x,12i5))',' IACT =',(IACT(I),I=1,NACT) if (ALL) then C%% printf( "\n PAR =" ); C%% for (i = 0; i < nact; i+=5){ C%% for (k = i; k < (i < nact - 4 ? i+5 : nact); k++) C%% printf( "%13.5g", par[k] ); C%% if (i < nact - 4) printf ( "\n " );} C%% printf( "\n" ); print'('' PAR ='',5g13.5/(8x,5g13.5))', (PAR(I),I=1,NACT) end if end if end if return end .