1 SUBROUTINE SMLC01 (SMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,ACC, 2 * IPRINT, NFMAX, IW, LIW, W, LW) 3c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 4c ALL RIGHTS RESERVED. 5c Based on Government Sponsored Research NAS7-03001. 6c File: SMLC.for Minimization with linear constraints. 7C>> 2001-10-05 SMLC Krogh Minor format changes -- 1PE => 1P,E 8C>> 2001-06-18 SMLC Krogh Changed ". LT." to " .LT." 9C>> 1996-07-08 SMLC Krogh Multiple changes for C conversion. 10c>> 1996-03-30 SMLC Krogh Removed MAX from type stmt., added external. 11c>> 1995-11-15 SMLC Krogh Moved formats up for C conversion. 12C>> 1994-11-11 SMLC Krogh Declared all vars. 13c>> 1994-11-02 SMLC Krogh Changes to use M77CON 14c>> 1992-04-27 SMLC CLL 15c>> 1992-04-17 CLL 16c>> 1992-03-27 CLL 17c>> 1992-02-05 CLL Declared all floating point variables. 18c>> 1992-01-16 CLL 19c>> 1991-06-10 CLL & FTK Editing comments. 20c>> 1991-04-30 CLL & FTK 21c>> 1991-04-19 CLL 22c>> 1991-04-15 FTK & CLL Made FIRST an argument of SMLC20 23c>> 1990-09-12 C. L. Lawson and F. T. Krogh, JPL 24c>> 1990-07-25 C. L. Lawson, JPL 25c>> 1990-07-12 C. L. Lawson, JPL 26c>> 1990-04-03 C. L. Lawson, JPL 27c--S replaces "?": ?MLC,?MLC01,?MLC02,?MLCFG,?MLC20,?MLC21,?MLC04,?MLC05 28c--& ?MLC13,?MLC15,?MLC03,?MLC12,?MLC06,?MLC11,?MLC09,?MLC16,?MLC19 29c--& ?MLC07,?MLC18,?MLC08,?MLC14,?MLC10,?MLC17,?ERV1 30c Also calls ERMOR, ERMSG, IERM1, IERV1 31C ------------------------------------------------------------------ 32c This algorithm and the original code is due to 33c M. J. D. Powell, Cambridge Univ., 1989. 34c 1990-04-03 This version adapted for usage at JPL by 35c C. L. Lawson. Instrumented with "c--" lines to automate type 36c conversions. 37c 1990-07-11 CLL Added new argument NFMAX to SMLC01. Added new 38c subroutine SMLC20 to compute gradient by finite differences. 39c Added argument HAVEG in the user-provided subr SMLCFG. 40c Added args XL and XU to subr _MLC09 in order to be able to pass 41c them on to SMLC20. 42c Added reference to [D/R]1MACH in _MLC15. 43c 1990-07-12 CLL Transposed indices in the array A(,). 44c The constraints are now Ax .le. b rather than (A**t)x .le. b. 45c 1990-07-19 CLL Reorganized the arg list of SMLC01. 46c Prev args INFO, NACT, IACT() are embedded in new arg IW(). 47c New item, FVAL, and previous args W() and PAR() are combined 48c into new arg W(). The new item FVAL is the final value of the 49c objective function. 50c 1990-07-24 CLL Setting INFO = 0 initially. It keeps this value 51c till it is changed for some reason. Previously it was set to 52c 4 initially. 53c Changed specification of IPRINT. It now contains 5 print 54c parameters. This gives more flexibility, such as requesting 55c printing only at termination. 56c Moved the main block of intermediate printing into a new subr, 57c SMLC21. This allows cleaner logic for testing when to print. 58c SMLC21 is called both from SMLC02 and from _MLC05. In revising 59c the print control we removed variables NFVALK and ITERP, and 60c added variables NFVPR and TOLPRT. 61c Printing of error messages is now done using JPL MATH77 error 62c printing subroutines. Requires linking of 63c SERV1, ERFIN, ERMOR, ERMSG, IERM1, and IERV1. The calls are in 64c SMLC01, SMLC02, and _MLC15. ERFIN is not called directly. 65c Setting INFO = 9 in SMLC02 when terminate due to solution being 66c determined just by the constraints. Also compute and return 67c the func value in this case. 68c 1990-09-12 CLL and FTK Added arg. W20 for _MLC20. 69c 1991-04-19 CLL Change to have NFVALS count every call to _MLCFG. 70c Move test of IPFREQ to be sure last iteration gets printed if 71c it is a multiple of IPFREQ. Changed usage of NFVPR. 72c 1992 March/April CLL Additions to the code in MINFUN (C05) and 73c LSRCH (C09). The code will now return INFO = 2 in some cases 74c in which it previously returned INFO = 3. It will also finish in 75c fewer iterations in some cases. Added more to error msg on quit- 76c ting with INFO = 3. Added HAVEG into arg list of _MLC05. 77c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 78c-- Begin Mask code changes. 79c Names of subroutines have been changed as follows: 80c Original names New DP names New SP names 81c GETMIN DMLC01 SMLC01 82c MINFLC DMLC02 SMLC02 83c EQCONS DMLC03 SMLC03 84c GETFES DMLC04 SMLC04 85c MINFUN DMLC05 SMLC05 86c CONRES DMLC06 SMLC06 87c GETD DMLC07 SMLC07 88c SDEGEN DMLC08 SMLC08 89c LSRCH DMLC09 SMLC09 90c NEWCON DMLC10 SMLC10 91c SATACT DMLC11 SMLC11 92c ADDCON DMLC12 SMLC12 93c ADJTOL DMLC13 SMLC13 94c DELCON DMLC14 SMLC14 95c INITZU DMLC15 SMLC15 96c KTVEC DMLC16 SMLC16 97c SDIRN DMLC17 SMLC17 98c STEPBD DMLC18 SMLC18 99c ZBFGS DMLC19 SMLC19 100c FGCALC DMLCFG SMLCFG 101c New code by FTK --> DMLC20 SMLC20 102c New code by CLL --> DMLC21 SMLC21 103c-- End Mask code changes. 104 105C 106C This is main entry point to a package of subroutines that calculate 107C the least value of a differentiable function of several variables 108C subject to linear constraints on the values of the variables. 109c ------------------------------------------------------------------ 110c Subroutine Arguments 111c 112C SMLCFG [in] Name of subroutine having an interface of the form 113C SUBROUTINE SMLCFG (N, X, F, G, HAVEG) 114c INTEGER N 115c [DOUBLE PRECISION or REAL] X(N), F, G(N) 116C LOGICAL HAVEG 117c provided by the user to compute values of the objective function, 118c F, and optionally its gradient vector, G(), at the argument value 119c given by X(). The user may choose either to write code to compute 120c G() or not. If SMLCFG computes G() it must also set HAVEG = .true 121c If SMLCFG does not compute G() it must set HAVEG = .false. 122c In this latter case the package will make additional calls to 123c SMLCFG to compute a finite difference approximation to the 124c gradient vector. 125c The vectors X() and G() will be of length N. The subroutine must 126c not change the values of N or X(). 127c SMLCFG may use COMMON block(s) if neccessary to communicate 128c data between the user's main program and SMLCFG. 129C N is the number of variables and must be set by the user. 130C M is the number of linear constraints (excluding simple bounds) and 131C must be set by the user. 132C MEQ is the number of constraints that are equalities and must be set 133C by the user. 134C A(.,.) is a 2-dimensional array whose rows are the gradients of 135C the M constraint functions. Its entries must be set by the user 136C and its dimensions must be at least M by N. 137C LDA is the actual first dimension of the array A that is supplied by 138C the user, so its value may not be less than M. 139C B(.) is a vector of constraint right hand sides that must also be set 140C by the user. Specifically the constraints on the variables X(I) 141C I=1(1)N are 142C A(K,1)*X(1)+...+A(K,N)*X(N) .EQ. B(K) K=1,...,MEQ 143C A(K,1)*X(1)+...+A(K,N)*X(N) .LE. B(K) K=MEQ+1,...,M . 144C Note that the data that define the equality constraints come 145C before the data of the inequalities. 146C XL(.) and XU(.) are vectors whose components must be set to lower and 147C upper bounds on the variables. Choose very large negative and 148C positive entries if a component should be unconstrained, or set 149C XL(I)=XU(I) to freeze the I-th variable. Specifically these 150C simple bounds are 151C XL(I) .LE. X(I) and X(I) .LE. XU(I) I=1,...,N . 152C X(.) is the vector of variables of the optimization calculation. Its 153C initial elements must be set by the user to an estimate of the 154C required solution. The subroutines can usually cope with poor 155C estimates, and there is no need for X(.) to be feasible initially. 156C These variables are adjusted automatically and the values that 157C give the least feasible calculated value of the objective function 158C are available in X(.) on the return from SMLC01. 159C ACC is a tolerance on the first order conditions at the calculated 160C solution of the optimization problem. These first order 161C conditions state that, if X(.) is a solution, then there is a set 162C of active constraints with indices IACT(K) K=1(1)NACT, say, such 163C that X(.) is on the boundaries of these constraints, and the 164C gradient of the objective function can be expressed in the form 165C GRAD(F)=PAR(1)*GRAD(C(IACT(1)))+... 166C ...+PAR(NACT)*GRAD(C(IACT(NACT))) . 167C Here PAR(K) K=1(1)NACT are Lagrange multipliers that are 168C nonpositive 169C for inequality constraints, and GRAD(C(IACT(K))) is the gradient 170C of the IACT(K)-th constraint function, so it is A(IACT(K),.) if 171C IACT(K) .LE. M, and it is minus or plus the J-th coordinate vector 172c if the constraint is the lower or upper bound on X(J) 173c respectively. The normal return from the calculation occurs when 174c X(.) is feasible and the sum of squares of components of the 175c vector RESKT(.) is at most ACC**2, where RESKT(.) is the 176c N-component vector of residuals of the first order condition that 177c is displayed above. Sometimes the package cannot satisfy this 178c condition, because noise in the function values can prevent a 179c change to the variables, no line search being allowed to increase 180c the objective function. 181C IPRINT [in, integer] Contains print control parameters, packed as 182c follows: 183c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST 184c Printing only begins after a first feasible point is found. 185c The items to be printed are selected by IPMORE. The other 186c parameters determine when to print. 187c IPFREQ is zero or positive. If positive, printing will be done on 188c iterations 0, IPFREQ, 2*IPFREQ, etc. 189c IPTOL = 0 or 1. 1 means to print the new TOL value 190c and the standard items each time TOL is changed. 191c IPFRST = 0 or 1. 1 means to print on the first iteration, i.e. on 192c iteration No. 0. 193c IPMORE = 0 or 1. 0 means the items to be printed are ITERC, 194c NFVALS, F, X(1:N), and G(1:N). 1 means to also print 195c IACT(1:NAXT), PAR(1:NACT) and RESKT(1:N). 196c IPLAST = 0 or 1. 1 means to print the final results, and the 197c reason for termination. 198c NFMAX [in, integer] If positive this sets an upper limit on the 199c number of function evaluations. When gradients are estimated by 200c differences, the count includes the function evaluations for this 201c purpose also. If zero there is no upper limit. 202c IW(.) [work, out, integer] Integer work space. Also used to return 203c status information. Its length must be at least LIW. 204c On return: 205c IW(1) contains INFO. Indicates reason for termination. 206c See below for explanation of values. 207c IW(2) contains ITERC. No. of iterations. 208c IW(3) contains NFVALS. No. of function evaluations, not 209c counting extra func evals done to estimate the gradient 210c when the gradient is not computed explicitly. In this 211c latter case the actual No. of func evals will be 212c (K+1)*NFVALS, where K is the number of solution 213c components whose lower and upper bounds are not equal. 214c IW(4) contains NACT. The number of active constraints at the 215c solution point. Will be in the interval [0, N]. 216c IW(5:4+NACT) contains IACT(1:NACT). IACT() is used as work 217c space of length M + 2*N. On return the first NACT 218c locations contain the indices of the constraints that 219c are active at the final point. Indices [1:M] refer to 220c rows of the system Ax .le. b. Indices [M+1:M+N] 221c refer to component lower bounds. Indices [M+N+1:M+2*N] 222c refer to component upper bounds. 223C 224c LIW [in, integer] Dimension of the IW() array. Require 225c LIW .ge. 4 + M + 2*N. 226c 227C W(.) [work, out, float] Working space array of float variables. 228c Also used to return results. 229c Its length must be at least LW. 230c On return: 231c W(1) contains FVAL, the final value of the objective function. 232c W(2:1+N) contains the final gradient vector, GRAD(1:N). 233c W(2+N:1+2*N) contains the Kuhn-Tucker residual vector, RESKT(1:N). 234c W(2+2*N:1+2*N+NACT) contains the Lagrange multipliers, PAR(1:NACT) 235c where NACT has a value in [0,N]. GRAD(), RESKT(), and PAR() 236c are mentioned above in the description of ACC. 237c LW [in, integer] Dimension of the W() array. Require 238c LW .ge. 3 + M + 16*N + N**2. 239c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 240c Return values of INFO 241c 242c INFO indicates the reason for termination. 243c 244C INFO=1 X(.) is feasible and the condition that depends on 245C ACC is satisfied. 246C INFO=2 X(.) is feasible and rounding errors are preventing 247C further progress. 248C INFO=3 X(.) is feasible but the objective function fails to 249C decrease although a decrease is predicted by the current gradient 250C vector. If this return occurs and RESKT(.) has large components 251C then the user's calculation of the gradient of the objective 252C function may be incorrect. One should also question the coding of 253C the gradient when the final rate of convergence is slow. 254C INFO=4 In this case the calculation cannot begin because IA 255C is less than N or because the lower bound on a variable is greater 256C than the upper bound. 257C INFO=5 This value indicates that the equality constraints 258C are inconsistent. These constraints include any components of 259C X(.) that are frozen by setting XL(I)=XU(I). 260C INFO=6 In this case there is an error return because the 261C equality constraints and the bounds on the variables are found to 262C be inconsistent. 263C INFO=7 This value indicates that there is no vector of 264C variables that satisfies all of the constraints. Specifically, 265C when this return or an INFO=6 return occurs, the current active 266C constraints (whose indices are IACT(K) K=1(1)NACT) prevent any 267C change in X(.) that reduces the sum of constraint violations, 268C where only bounds are included in this sum if INFO=6. 269C INFO=8 In this case the limit on the number of calls of 270C subroutine SMLCFG (see below) has been reached, and there would 271C have been further calculation otherwise. 272c INFO = 9 The solution is determined just by the constraints. 273c In this case PAR() and RESKT() are not defined on return. 274c FVAL will be defined on return. G() will only be defined on 275c return if SMLCFG returns HAVEG = .true. 276c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 277C Note 1. The variables N, M, MEQ, LDA, ACC, IPRINT, NFMAX and the 278C elements of the arrays A(.,.), B(.), XL(.) and XU(.) are not 279c altered by the optimization procedure. Their values, and the 280c initial components of X(.) must be set on entry to SMLC01. 281c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 282C Note 2. A paper on the method of calculation and a report on the 283C main features of the computer code are available from the author 284C M.J.D.Powell (D.A.M.T.P., University of Cambridge, Silver Street, 285C Cambridge CB3 9EW, England). 286c ------------------------------------------------------------------ 287 external SMLCFG 288 integer LDA, LIW, LW 289 integer I, IBRES, ID 290 integer IFVAL, IG, IGM, IGS, IIACT, IINFO, IITERC, INACT, INFO 291 integer INFVAL, IPAR, IPRINT, IRESKT, IU, IW(LIW), IW20, IXBIG 292 integer IXS, IZ, IZTG, KTNORM, M, MEQ, N, NFMAX 293 real A(LDA,*),ACC, B(*),TEMP, XL(*),XU(*),X(*),W(LW) 294 parameter(IFVAL = 1, KTNORM = 2, IG = 3) 295 parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5) 296C ------------------------------------------------------------------ 297c Check sizes of LIW and LW 298 INFO = 0 299 if(LIW .lt. 4 + M + 2*N) then 300 INFO = 4 301 call IERM1('SMLC01', INFO, 0,'Dimension LIW is too small.', 302 * 'LIW',LIW, ',') 303 call IERV1('Need',4 + M + 2*N,'.') 304 endif 305c 306 if(LW .lt. 3 + M + N*(16 + N)) then 307 INFO = 4 308 call IERM1('SMLC01', INFO, 0,'Dimension LW is too small.', 309 * 'LW',LW, ',') 310 call IERV1('Need',3 + M + N*(16 + N),'.') 311 endif 312 if(INFO .ne. 0) then 313 IW(1) = INFO 314 return 315 endif 316c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 317C Partition the array W() 318c parameter(IFVAL = 1, KTNORM = 2, IG = 3) 319 IRESKT=IG+N 320 IPAR = IRESKT + N 321 IZ=IPAR + N 322 IU=IZ+N*N 323 IXBIG=IU+N 324 IBRES=IXBIG+N 325 ID=IBRES+M+N+N 326 IZTG=ID+N 327 IGM=IZTG+N 328 IXS=IGM+N 329 IGS=IXS+N 330C Need 4*N+1 locations beginning at IW20. 331 IW20=IGS+N 332C Partition the array IW() 333c parameter(IINFO = 1, IITERC = 2, INFVAL = 3, INACT = 4, IIACT = 5) 334C 335c Zero the RESKT vector so it will be safe to compute its norm on 336c return from SMLC02 even in cases when SMLC02 does not compute 337c RESKT. 338c 339 do 10 I = IRESKT, IRESKT + N -1 340 W(I) = 0.0E0 341 10 continue 342c 343C Call the optimization package. 344C 345 CALL SMLC02 (SMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X,ACC, 346 * IW(IIACT),IW(INACT),W(IPAR),IPRINT, NFMAX, 347 * IW(IINFO), IW(IITERC), IW(INFVAL), 348 * W(IG),W(IZ),W(IU),W(IXBIG),W(IRESKT),W(KTNORM),W(IBRES),W(ID), 349 * W(IZTG),W(IGM),W(IXS),W(IGS), W(IFVAL), W(IW20)) 350c 351 TEMP = 0.0E0 352 do 20 I = IRESKT, IRESKT+N-1 353 TEMP = TEMP + W(I)**2 354 20 continue 355 W(KTNORM) = sqrt(TEMP) 356 RETURN 357 END 358c ================================================================== 359 SUBROUTINE SMLC02 (SMLCFG, N,M,MEQ,A,LDA,B,XL,XU,X, 360 * ACC,IACT,NACT,PAR, IPRINT, NFMAX, INFO, ITERC, NFVALS, 361 * G,Z,U,XBIG,RESKT, ENRMKT, BRES,D,ZTG,GM,XS,GS, FVAL, W20) 362c 363c Main subroutine to control the Powell constrained minimization 364c package. 365c 1990-07-19 CLL Added NFMAX, ITERC, NFVALS, and FVAL to arg list. 366c Also added FVAL to arg list of SMLC05 that is called 367c from this subr. 368c NFMAX provides limit on no. of func evaluations, but zero means 369c no limit. This was previously overloaded onto INFO. 370c ITERC & NFVALS give output of counts of iterations and 371c function evaluations. 372c FVAL provides output of the final objective function value. 373c ------------------------------------------------------------------ 374c Arguments 375c 376c SMLCFG, N,M,MEQ,A(),LDA,B(),XL(),XU(),X(), 377c ACC,IACT(),NACT,PAR(), 378c IPRINT [in, integer] Contains print parameters packed as follows: 379c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST 380c NFMAX [in, integer] Limit on No. of function evaluations. But zero 381c means no limit. 382c INFO [out, integer] Termination status flag. 383c ITERC [out, integer] Count of No. of iterations. Initialized to 0 384c in this subr. Incremented in SMLC05. 385c NFVALS [out, integer] Count of number of function evaluations. 386c G(),Z(),U(),XBIG(), 387c RESKT() 388c ENRMKT Euclidean norm of RESKT(). 389c BRES(),D(), 390c ZTG(),GM(),XS(),GS(), 391c FVAL [out, float] Final (best) value of F. 392c W20 [out, work, float] Working storage for SMLC20. 393c ------------------------------------------------------------------ 394c Descriptions of some of the internal variables. 395c 396c GMODE [integer] Initialized to 0 here. Sent to SMLC05 for eventual 397c use by SMLC20. 398c HAVEG [logical] Arg in call to SMLCFG, not used in this subroutine. 399c NFVPR [integer] Used to avoid printing results from the same 400c function evaluation more than once. Initialized to -1 in this 401c subr. Modified in SMLC05. Used in SMLC02 and DNLC05. 402c ------------------------------------------------------------------ 403 external SMLCFG 404 405 integer GMODE, I, IACT(*), INFO 406 integer IPFREQ, IPFRST, IPLAST, IPMORE, IPRINT 407 integer IPTOL, ITERC, K, LDA, M, MEQ, MEQL, MP, MSAT 408 integer MTOT, N, NACT, NFMAX, NFVALS, NFVPR 409 410 real A(LDA,*),ACC, B(*),BRES(*),D(*),ENRMKT, FVAL 411 real G(*), GM(*), GS(*), PAR(*), RELACC, RESKT(*) 412 real SSQKT, TOL, U(*) 413 real W20(0:*), XBIG(*), XL(*), XU(*), X(*), XS(*) 414 real Z(*), ZTG(*), ZZNORM 415 logical HAVEG 416c ------------------------------------------------------------------ 417C Initialize ZZNORM, ITERC, NFVPR, NFVALS, INFO 418C 419 ZZNORM=-1.0E0 420 W20(0) = -1.0e0 421 W20(1) = -1.0e0 422 GMODE = 0 423 ITERC=0 424 NFVPR = -1 425 NFVALS=0 426 INFO = 1 427c 428c Decompose IPRINT into 5 integers: 429c IPRINT = IPFREQ*100 + IPTOL*8 + IPFRST*4 + IPMORE*2 + IPLAST 430C 431 IPFREQ = IPRINT/100 432 I = IPRINT - IPFREQ * 100 433 IPTOL = I/8 434 I = I - IPTOL * 8 435 IPFRST = I/4 436 I = I - IPFRST * 4 437 IPMORE = I/2 438 IPLAST = I - IPMORE * 2 439c 440C Check the bounds on N, M and MEQ. 441C 442 if(N .lt. 1 .or. M .lt. 0 .or. MEQ .lt. 0 .or. M .lt. MEQ) then 443 INFO = 4 444 call IERM1('SMLC02', INFO, 0,'Bad value for M, N, or MEQ.', 445 * 'M',M, ',') 446 call IERV1('N',N,',') 447 call IERV1('MEQ',MEQ,'.') 448 GOTO 40 449 END IF 450C 451C Initialize RELACC, Z, U and TOL. 452C 453 CALL SMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC) 454 IF (INFO .EQ. 4) go to 40 455 TOL=max(0.01E0,10.0E0*RELACC) 456 if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or. 457 * IPLAST .ne. 0 .or. IPTOL .ne. 0) then 458 print'(/6x,''Beginning subroutine SMLC02, called from SMLC01:'' 459 * /6x,''N ='',i4,'', M ='',i4,'', MEQ ='',i4,'', ACC ='', 460 * g11.3/6x,''RELACC ='',g11.3,'', TOL = '',g11.3)', 461 * N, M, MEQ, ACC, RELACC, TOL 462 endif 463C 464C Add any equality constraints to the active set. 465C 466 IF (MEQ .GT. 0) THEN 467 CALL SMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC,XS, 468 1 GS) 469 IF (INFO .EQ. 5) THEN 470 call ERMSG('SMLC02',INFO,0, 471 * 'Equality constraints are inconsistent.','.') 472 GOTO 40 473 END IF 474 END IF 475 NACT=MEQL 476 MSAT=MEQL 477C 478C Add the bounds to the list of constraints. 479C 480 MTOT=NACT 481 DO 10 I=1,N 482 IF (XL(I) .LT. XU(I)) THEN 483 MTOT=MTOT+2 484 IACT(MTOT-1)=M+I 485 IACT(MTOT)=M+N+I 486 END IF 487 10 CONTINUE 488C 489C Try to satisfy the bound constraints. 490C 491 CALL SMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG, 492 1 RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20) 493 IF (MSAT .LT. MTOT) THEN 494 INFO=6 495 call ERMSG('SMLC02',INFO,0, 496 * 'Equalities and bounds are inconsistent.','.') 497 GOTO 40 498 END IF 499C 500C Add the ordinary inequalities to the list of constraints. 501C 502 IF (M .GT. MEQ) THEN 503 MP=MEQ+1 504 DO 20 K=MP,M 505 MTOT=MTOT+1 506 20 IACT(MTOT)=K 507 END IF 508c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 509 30 continue 510c Begin main loop 511C 512C Correct any constraint violations. 513C 514 CALL SMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z,U,XBIG, 515 1 RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,RESKT,XS,GS,W20) 516 IF (MSAT .LT. MTOT) THEN 517 INFO=7 518 call ERMSG('SMLC02',INFO,0, 519 * 'Constraints are inconsistent.','.') 520 GOTO 40 521 ELSE IF (MEQL .EQ. N) THEN 522c 523c Here the soln is determined just by the constraints. 524c PAR() and RESKT() are not defined on return in this case. 525c FVAL has not yet been evaluated. We evaluate FVAL here in 526c order to return its value. We leave G() not computed if HAVEG 527c is false, however it would be easy to compute G by calling 528c SMLC20 if that seemed preferable. 529c 530C%% (*smlcfg)( n, x, fval, g, &haveg ); 531 call SMLCFG(N, X, FVAL, G, HAVEG) 532 NFVALS = NFVALS+1 533 ENRMKT = 0.0e0 534 INFO = 9 535 536 if(IPFREQ .ne. 0 .or. IPFRST .ne. 0 .or. 537 * IPLAST .ne. 0 .or. IPTOL .ne. 0) then 538 print '(/'' [1] SMLC02..''/ 539 * '' The solution is determined by the equality constraints.'')' 540 call SMLC21(IPMORE, .false., N, ITERC, NFVALS, 541 * FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) 542 endif 543 GOTO 40 544 END IF 545C 546C Minimize the objective function in the case when constraints are 547C treated as degenerate if their residuals are less than TOL. 548C 549c SMLC05 is MINFUN 550 CALL SMLC05 (SMLCFG, N,M,A,LDA,B,XL,XU,X,ACC,IACT,NACT, 551 * PAR, IPFREQ, IPTOL, IPFRST, IPMORE, INFO, G, Z, 552 * U,XBIG,RELACC,ZZNORM,TOL,MEQL,MTOT,ITERC, NFVPR, NFVALS, 553 * NFMAX,RESKT, SSQKT, BRES,D,ZTG,GM,XS,GS, FVAL, GMODE, 554 * HAVEG, W20) 555C 556C Reduce TOL if necessary. 557C 558 IF (TOL .GT. RELACC .AND. NACT .GT. 0) THEN 559 if (NFMAX .le. 0 .or. NFVALS .lt. NFMAX) then 560C SMLC13 is ADJTOL 561 CALL SMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC, 562 * TOL, MEQL) 563 GOTO 30 564 else 565 INFO=8 566 endif 567 END IF 568c 569c Here when INFO = 1, 2, 3, or 8. 570c We treat 1 and 2 as normal conditions and 3 and 8 as errors. 571c 572 if (INFO .EQ. 1) then 573 if (IPLAST .ne. 0) print*, 574 * '[1] SMLC01 quitting because requested accuracy', 575 * ' has been attained.' 576 577 elseif (INFO .EQ. 2) then 578 if (IPLAST .ne. 0) print*, 579 * '[2] SMLC01 quitting due to limitation of', 580 * ' computational precision.' 581 582 elseif (INFO .EQ. 3) then 583 call ERMSG('SMLC02', INFO, 0, 584 * 'F not decreasing, although gradient indicates it should.',',') 585 if(HAVEG) then 586 call ERMOR('Could be due to limitation of precision',',') 587 call ERMOR('or to incorrect code for the gradient.','.') 588 else 589 call ERMOR('Could be due to limitation of precision.','.') 590 endif 591 592 elseif (INFO .EQ. 8) then 593 call IERM1('SMLC02', INFO, 0, 594 * 'Reached specified max number of function evaluations.', 595 * 'Count',NFVALS, ',') 596 call IERV1('NFMAX',NFMAX,'.') 597 endif 598c 599c Test to print results before quitting. 600c 601 ENRMKT = sqrt(SSQKT) 602 if(IPLAST .ne. 0 .and. NFVALS .ne. NFVPR) then 603c print*,'SMLC02.. Debug.. Call C21 on leaving C02.' 604 call SMLC21(IPMORE, .true., N, ITERC, NFVALS, 605 * FVAL, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) 606 endif 607c 608 40 return 609 end 610c ================================================================== 611 SUBROUTINE SMLC03 (N,M,MEQ,A,LDA,B,XU,IACT,MEQL,INFO,Z,U,RELACC, 612 1 AM,CGRAD) 613C 614C Try to add the next equality constraint to the active set. 615c ------------------------------------------------------------------ 616 integer N, M, MEQ, LDA, IACT(*), MEQL, INFO 617 integer I, IZ, J, JM, K, KEQ, NP 618 real A(LDA,*),B(*),XU(*) 619 real RELACC, RHS, SUM, SUMABS, VMULT 620 real Z(*),U(*),AM(*),CGRAD(*) 621c ------------------------------------------------------------------ 622 DO 50 KEQ=1,MEQ 623 IF (MEQL .LT. N) THEN 624 NP=MEQL+1 625 IACT(NP)=KEQ 626 CALL SMLC12 (N,M,A,LDA,IACT,MEQL,Z,U,RELACC,NP,AM,CGRAD) 627 IF (MEQL .EQ. NP) GOTO 50 628 END IF 629C 630C If linear dependence occurs then find the multipliers of the 631C dependence relation and apply them to the right hand sides. 632C 633 SUM=B(KEQ) 634 SUMABS=abs(B(KEQ)) 635 IF (MEQL .GT. 0) THEN 636 DO 10 I=1,N 637 10 AM(I)=A(KEQ,I) 638 K=MEQL 639 20 VMULT=0.0E0 640 IZ=K 641 DO 30 I=1,N 642 VMULT=VMULT+Z(IZ)*AM(I) 643 30 IZ=IZ+N 644 VMULT=VMULT*U(K) 645 J=IACT(K) 646 IF (J .LE. M) THEN 647 DO 40 I=1,N 648 40 AM(I)=AM(I)-VMULT*A(J,I) 649 RHS=B(J) 650 ELSE 651 JM=J-M-N 652 AM(JM)=AM(JM)-VMULT 653 RHS=XU(JM) 654 END IF 655 SUM=SUM-RHS*VMULT 656 SUMABS=SUMABS+abs(RHS*VMULT) 657 K=K-1 658 IF (K .GE. 1) GOTO 20 659 END IF 660C 661C Error return if the constraints are inconsistent. 662C 663 IF (abs(SUM) .GT. RELACC*SUMABS) THEN 664 INFO=5 665 GOTO 60 666 END IF 667 50 CONTINUE 668 60 RETURN 669 END 670c ================================================================== 671 SUBROUTINE SMLC04 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,INFO,G,Z, 672 1 U,XBIG,RELACC,TOL,MEQL,MSAT,MTOT,BRES,D,ZTG,GM,GMNEW,PARNEW, 673 2 CGRAD,W20) 674C 675C Make the correction to X for the active constraints. 676c ------------------------------------------------------------------ 677 integer N, M, LDA, IACT(*), NACT, INFO, MEQL, MSAT, MTOT 678 integer I, ITEST, MSATK, INDXBD 679 real A(LDA,*),B(*),XL(*),XU(*),X(*),PAR(*),G(*) 680 real U(*),XBIG(*),BRES(*),D(*) 681 real GM(*),GMNEW(*),PARNEW(*), CGRAD(*) 682 real W20(0:*), RELACC, STEPCB, SUMRES, SUMRSK, TOL 683 real Z(*), ZTG(*) 684c ------------------------------------------------------------------ 685 INFO=0 686 10 CALL SMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG,RELACC, 687 1 TOL,MEQL) 688 IF (INFO .GT. 0) MSAT=NACT 689 IF (MSAT .EQ. MTOT) GOTO 60 690C 691C Try to correct the infeasibility. 692C 693 20 MSATK=MSAT 694 SUMRSK=0.0E0 695 30 CALL SMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES, 696 1 D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD,GM,GMNEW, 697 2 PARNEW,CGRAD,W20) 698C 699C Include the new constraint in the active set. 700C 701 IF (STEPCB .GT. 0.0E0) THEN 702 DO 40 I=1,N 703 X(I)=X(I)+STEPCB*D(I) 704 40 XBIG(I)=max(XBIG(I),abs(X(I))) 705 CALL SMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,GMNEW, 706 * CGRAD) 707 END IF 708C 709C Test whether to continue the search for feasibility. 710C 711 IF (MSAT .LT. MTOT) THEN 712 IF (STEPCB .EQ. 0.0E0) GOTO 50 713 IF (MSATK .LT. MSAT) GOTO 20 714 IF (SUMRSK .EQ. 0.0E0 .OR. SUMRES .LT. SUMRSK) THEN 715 SUMRSK=SUMRES 716 ITEST=0 717 END IF 718 ITEST=ITEST+1 719 IF (ITEST .LE. 2) GOTO 30 720C 721C Reduce TOL if it may be too large to allow feasibility. 722C 723 50 IF (TOL .GT. RELACC) THEN 724 CALL SMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC, 725 1 TOL,MEQL) 726 GOTO 10 727 END IF 728 END IF 729 60 RETURN 730 END 731c ================================================================== 732 SUBROUTINE SMLC05 (SMLCFG, N,M,A,LDA,B,XL,XU,X,ACC,IACT,NACT, 733 * PAR, IPFREQ, IPTOL, IPFRST, IPMORE, 734 * INFO,G,Z,U,XBIG,RELACC,ZZNORM,TOL,MEQL,MTOT,ITERC,NFVPR,NFVALS, 735 * NFMAX,RESKT, SSQKT, BRES,D,ZTG,GM,XS,GS, FVAL, GMODE, 736 * HAVEG, W20) 737C 738C MINFUN: Solve the minimization problem using the current value of 739c TOL. 740c 1990-07-19 CLL Added FVAL to argument list to return final value 741c of F. 742c 1990-07-23 CLL Removed IPRINT from arg list and added IPFREQ, 743c IPTOL, IPFRST, IPMORE, and NFVPR. NFVPR keeps a record of 744c last function evaluation which printing has been done. 745c 1990-07-23 CLL Added logical variable TOLPRT. 746c ------------------------------------------------------------------ 747c Arguments 748c 749c SMLCFG, N,M,A(),LDA,B(),XL(),XU(),X(),ACC,IACT(),NACT, 750c PAR() 751c IPFREQ [in, integer] Zero or positive. If positive, printing of 752c intermediate results using SMLC21 will be done on iterations 753c 0, IPFREQ, 2*IPFREQ, etc. 754c IPTOL [in, integer] = 0 or 1. 1 means to print the new TOL value 755c and items printed by SMLC21 each time this subr is entered. 756c IPFRST [in, integer] = 0 or 1. 1 means to print using SMLC21 on 757c the first iteration, i.e. iteration No. 0. 758c IPMORE [in, integer] = 0 or 1. Will be passed to SMLC21 to cause 759c printing to include more items when = 1 than when = 0. 760c INFO [inout, integer] Status flag. 761c G(),Z(),U(),XBIG(),RELACC,ZZNORM,TOL,MEQL,MTOT, 762c NFVALS [inout, integer] Count of calls for function evaluation. 763c Only counts calls to SMLCFG when the algorithm directly needs a 764c function value, not calls to SMLCFG that occur in SMLC20 for 765c estimation of the gradient. 766c NFMAX [in, integer] Limit on No. of function evaluations, but zero 767c means no limit. 768c RESKT() 769c SSQKT Sum of Squares of elts of RESKT(). 770c BRES(),D(),ZTG(),GM(),XS(),GS(), 771c FVAL [out, float] On return set to current (best) value of F. 772c GMODE [inout, integer] For use by SMLC20. 773c HAVEG [out, logical] Primarily an internal variable in this subr, 774c but returned so calling subr can use it to select proper error msg 775c when INFO = 3. 776c Set true by SMLCFG if SMLCFG provides the gradient 777c vector and false otherwise. If false, this subr calls SMLC20 to 778c compute a finite difference approximation for the gradient. 779c W20 [out, work, float] Working storage for SMLC20. 780c ------------------------------------------------------------------ 781c Description of some of the internal variables. 782c 783c DOPRNT [logical] When IPTOL = 1 this subr prints TOL on entry and 784c sets DOPRNT = true. Also if IPFRST = 1 and ITERC = 0 this subr 785c sets DOPRNT = true. 786c Having DOPRNT = true causes subsequent printing of intermediate 787c results for the current iteration, after PAR() and RESKT() have 788c been computed. 789c HAVEG [logical] Set true by SMLCFG if SMLCFG provides the gradient 790c vector and false otherwise. If false, this subr calls SMLC20 to 791c compute a finite difference approximation for the gradient. 792c ------------------------------------------------------------------ 793 external R1MACH, SMLCFG 794 integer N, M, LDA, IACT(*), NACT, IPFREQ, IPTOL, IPFRST, IPMORE, 795 * INFO, MEQL, MTOT, ITERC, NFVPR, NFVALS, NFMAX, GMODE 796 integer I, ITERK, K, MSAT, INDXBD 797 real R1MACH 798 real A(LDA,N), ACC, B(*), BRES(*), DDOTG, DIFF 799 real D(*), EPS, F, FPREV, FVAL, G(N), GM(*), GS(*) 800 real PAR(*), RELACC, RELAXF, RESKT(*) 801 real SSQKT, STEP, STEPCB, SUM, TOL 802 real U(*), W20(0:*) 803 real X(*), XBIG(*), XL(N), XS(*), XU(N) 804 real Z(*), ZTG(*), ZZNORM 805 logical DOPRNT, HAVEG 806 save EPS, F 807 data EPS / 0.0e0 / 808c ------------------------------------------------------------------ 809c Initialize the minimization calculation. 810c 811* print'(/''SMLC05.. Debug.. On entry, INFO = '',i6,'', ITERC ='', 812* i6)', INFO, ITERC 813 if(EPS .eq. 0.0e0) then 814 EPS = R1MACH(4) 815 endif 816 MSAT=MTOT 817 ITERK=ITERC 818 IF (NFVALS .EQ. 0 .OR. INFO .EQ. 1) THEN 819c++ CODE for ~.C. is active 820 CALL SMLCFG (N,X,F,G, HAVEG) 821c++ CODE for .C. is inactive 822C%% (*smlcfg)( n, x, &f, g, haveg ); 823C HAVEG = HAVEG 824c++ END 825* print'('' SMLC05 798.. X()='',2g23.15/(19x,2g23.15))', 826* * (X(I),I=1,N) 827* print'('' .. F ='',g23.15)', F 828* print'('' ............ G()='',2g23.15/(19x,2g23.15))', 829* * (G(I),I=1,N) 830 831 NFVALS=NFVALS+1 832 if(.not. HAVEG) CALL SMLC20(SMLCFG,N,X,F,G,XL,XU, 833 * GMODE, NFVALS, W20) 834 END IF 835 FPREV=abs(F+F+1.0E0) 836c Setting DIFF here is not needed for the 837c algorithm, but is convenient when doing 838c debug printing. CLL 3/26/92 839 DIFF = 1.0e0 840 IF (IPTOL .ne. 0) THEN 841 print'(/'' New value of TOL = '',g13.5)', TOL 842 DOPRNT = .true. 843 elseif(IPFRST .ne. 0 .and. ITERC .eq. 0) then 844 DOPRNT = .true. 845 else 846 DOPRNT = .false. 847 END IF 848C 849C Calculate the next search direction. 850C SMLC06 is CONRES 851 10 CALL SMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG,BRES,D, 852 1 ZTG,RELACC,TOL,STEPCB,DDOTG,MEQL,MSAT,MTOT,INDXBD,GM,RESKT,XS, 853 2 GS,W20) 854C 855C Calculate the Kuhn Tucker residual vector. 856C SMLC16 is KTVEC 857 CALL SMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES,RELAXF,MEQL, 858 1 SSQKT,XS,GS) 859c 860c Test for printing results every IPFREQ iterations. 861c 862 if(IPFREQ .ne. 0) then 863 if( mod(ITERC, IPFREQ) .eq. 0) DOPRNT = .true. 864 endif 865c 866 if(DOPRNT) then 867 call SMLC21(IPMORE, .true., N, ITERC, NFVALS, 868 * F, TOL, X, G, NACT, IACT, PAR, RESKT, sqrt(SSQKT)) 869 NFVPR = NFVALS 870 DOPRNT = .false. 871 endif 872C 873C Test for convergence. 874C 875* print'(/'' SMLC05 849..''/'' SMLC05 849..'',3g13.5/ 876* * '' F, FPREV, FPREV-F='',3g13.5/ 877* * '' DIFF='',g13.5/ 878* * '' TOL, RELACC, NACT='',2g13.5,i5)', 879* * ACC**2, SSQKT, DDOTG, F, FPREV, FPREV-F, DIFF, 880* * TOL, RELACC, NACT 881 IF (SSQKT .LE. ACC*ACC) THEN 882 INFO=1 883 GOTO 70 884 END IF 885 IF (DDOTG .GE. 0.0E0) THEN 886 INFO=2 887 GOTO 70 888 END IF 889C 890C Test for termination due to no decrease in F. 891C 892 IF (F .GE. FPREV) THEN 893 IF (TOL .EQ. RELACC .OR. NACT .EQ. 0) THEN 894 IF (DIFF .GT. 0.0E0) GOTO 20 895 END IF 896* print'(/'' SMLC05 869.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F 897 if(abs(DDOTG) .le. EPS * abs(F)) then 898 INFO = 2 899 else 900 INFO=3 901 endif 902 GOTO 70 903 END IF 904 20 continue 905 DIFF=FPREV-F 906 FPREV=F 907C 908C Test that more calls of SMLCFG are allowed. 909C 910 IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) THEN 911 INFO=8 912 GOTO 70 913 END IF 914C 915C Test whether to reduce TOL. 916C 917 IF (TOL .GT. RELACC .AND. ITERC .GT. ITERK .AND. 918 1 0.1E0*RELAXF .GE. max(DIFF,-0.5E0*DDOTG)) then 919* print'('' SMLC05.. Debug.. return to reduce TOL.''/ 920* * '' ........ TOL='',g13.5,'', RELACC='',g13.5,'', INFO='', 921* * i5)', TOL, RELACC, INFO 922 GOTO 70 923 endif 924c 925C Calculate the step along the search direction. 926C 927 ITERC=ITERC+1 928* print'(/'' SMLC05 895.. ITERC = '',i6)', ITERC 929C SMLC09 is LSRCH 930 CALL SMLC09 (SMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB,DDOTG,F, 931 * STEP,NFVALS, NFMAX,BRES, XL, XU, GMODE, W20) 932 IF (STEP .EQ. 0.0E0) THEN 933* print'(/'' SMLC05 918.. DDOTG, EPS*F='',2g13.5)', DDOTG, EPS*F 934 INFO=3 935 if(abs(DDOTG) .le. EPS * abs(F)) then 936 INFO = 2 937 else 938 SUM=0.0E0 939 DO 50 I=1,N 940 50 SUM=SUM+abs(D(I)*GS(I)) 941* print'(/'' SMLC05 890.. DDOTG, SUM='',2g13.5/ 942* * '' RELACC*SUM, DDOTG+RELACC*SUM='',2g13.5)', 943* * DDOTG, SUM, RELACC*SUM, DDOTG+RELACC*SUM 944 IF (DDOTG+RELACC*SUM .GE. 0.0E0) INFO=2 945 endif 946 GOTO 70 947 END IF 948C 949C Revise XBIG. 950C 951 DO 60 I=1,N 952 60 XBIG(I)=max(XBIG(I),abs(X(I))) 953C 954C Revise the second derivative approximation. 955C 956 CALL SMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM) 957C 958C Add a constraint to the active set if it restricts the step. 959C 960 IF (STEP .EQ. STEPCB) THEN 961 K=IACT(INDXBD) 962 IF (K .GT. M) THEN 963 K=K-M 964 IF (K .LE. N) THEN 965 X(K)=XL(K) 966 ELSE 967 X(K-N)=XU(K-N) 968 END IF 969 END IF 970 CALL SMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,XS,GS) 971 END IF 972 GOTO 10 973C 974 70 continue 975 FVAL = F 976 return 977 end 978c ================================================================== 979 SUBROUTINE SMLC06 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,PAR,G,Z,U,XBIG, 980 1 BRES,D,ZTG,RELACC,TOL,STEPCB,SUMRES,MEQL,MSAT,MTOT,INDXBD, 981 2 GM,GMNEW,PARNEW,CGRAD,W20) 982C 983c CONRES 984C Calculate and partition the residuals of the inactive constraints, 985C and set the gradient vector when seeking feasibility. 986C W20 has length 4*N+1. Only locations 1 to N are used here. 987c ------------------------------------------------------------------ 988 integer N,M,LDA,IACT(*),NACT,MEQL,MSAT,MTOT,INDXBD 989 integer I, IDIFF, J, JM, K, KL, MDEG, MSATK 990 real A(LDA,*),B(*), DDOTG, PAR(*) 991 real G(*), U(*),XBIG(*),BRES(*) 992 real D(*), GM(*),GMNEW(*),PARNEW(*), CGRAD(*) 993 real RELACC, RES, RESABS, STEPCB, SUM, SUMRES 994 real TEMP, TOL 995 real W20(0:*), XL(*),XU(*),X(*) 996 real Z(*), ZTG(*) 997c ------------------------------------------------------------------ 998 IDIFF=MTOT-MSAT 999C 1000 IF (IDIFF .GT. 0) THEN 1001 DO 10 I=1,N 1002 10 G(I)=0.0E0 1003 SUMRES=0.0E0 1004 END IF 1005 MSATK=MSAT 1006 MDEG=NACT 1007 MSAT=NACT 1008 KL=MEQL+1 1009 DO 50 K=KL,MTOT 1010 J=IACT(K) 1011C 1012C Calculate the residual of the current constraint. 1013C 1014 IF (J .LE. M) THEN 1015 RES=B(J) 1016 RESABS=abs(B(J)) 1017 DO 20 I=1,N 1018 RES=RES-X(I)*A(J,I) 1019 20 RESABS=RESABS+abs(XBIG(I)*A(J,I)) 1020 ELSE 1021 JM=J-M 1022 IF (JM .LE. N) THEN 1023 RES=X(JM)-XL(JM) 1024 RESABS=abs(XBIG(JM))+abs(XL(JM)) 1025 ELSE 1026 JM=JM-N 1027 RES=XU(JM)-X(JM) 1028 RESABS=abs(XBIG(JM))+abs(XU(JM)) 1029 END IF 1030 END IF 1031 BRES(J)=RES 1032C 1033C Set TEMP to the relative residual. 1034C 1035 TEMP=0.0E0 1036 IF (RESABS .NE. 0.0E0) TEMP=RES/RESABS 1037 IF (K .GT. MSATK .AND. TEMP .LT. 0.0E0) THEN 1038 IF (TEMP+RELACC .GE. 0.0E0) THEN 1039 IF (J .LE. M) THEN 1040 SUM=abs(B(J)) 1041 DO 30 I=1,N 1042 30 SUM=SUM+abs(X(I)*A(J,I)) 1043 ELSE 1044 JM=J-M 1045 IF (JM .LE. N) THEN 1046 SUM=abs(X(JM))+abs(XL(JM)) 1047 ELSE 1048 SUM=abs(X(JM-N))+abs(XU(JM-N)) 1049 END IF 1050 END IF 1051 IF (abs(RES) .LE. SUM*RELACC) TEMP=0.0E0 1052 END IF 1053 END IF 1054C 1055C Place the residual in the appropriate position. 1056C 1057 IF (K .LE. NACT) GOTO 50 1058 IF (K .LE. MSATK .OR. TEMP .GE. 0.0E0) THEN 1059 MSAT=MSAT+1 1060 IF (MSAT .LT. K) THEN 1061 IACT(K)=IACT(MSAT) 1062 END IF 1063 IF (TEMP .GT. TOL) THEN 1064 IACT(MSAT)=J 1065 ELSE 1066 MDEG=MDEG+1 1067 IACT(MSAT)=IACT(MDEG) 1068 IACT(MDEG)=J 1069 END IF 1070C 1071C Update the gradient and SUMRES if the constraint is violated when 1072C seeking feasibility. 1073C 1074 ELSE 1075 IF (J .LE. M) THEN 1076 DO 40 I=1,N 1077 40 G(I)=G(I)+A(J,I) 1078 ELSE 1079 J=J-M 1080 IF (J .LE. N) THEN 1081 G(J)=G(J)-1.0E0 1082 ELSE 1083 G(J-N)=G(J-N)+1.0E0 1084 END IF 1085 END IF 1086 SUMRES=SUMRES+abs(RES) 1087 END IF 1088 50 CONTINUE 1089C 1090C Seek the next search direction unless SMLC06 was called from 1091C SMLC04 [GETFES] and feasibility has been achieved. 1092C 1093 STEPCB=0.0E0 1094 IF (IDIFF .GT. 0 .AND. MSAT .EQ. MTOT) GOTO 60 1095c GETD 1096 CALL SMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC,DDOTG, 1097 * MEQL, MDEG,GM,GMNEW,PARNEW,CGRAD,W20) 1098C 1099C Calculate the (bound on the) step-length due to the constraints. 1100C 1101 IF (DDOTG .LT. 0.0E0) THEN 1102 CALL SMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT, 1103 1 MTOT,INDXBD) 1104 END IF 1105 IF (IDIFF .EQ. 0) SUMRES=DDOTG 1106 60 RETURN 1107 END 1108c ================================================================== 1109 SUBROUTINE SMLC07 (N,M,A,LDA,IACT,NACT,PAR,G,Z,U,D,ZTG,RELACC, 1110 1 DDOTG,MEQL,MDEG,GM,GMNEW,PARNEW,CGRAD,W20) 1111C 1112c GETD 1113C Initialize GM and cycle backwards through the active set. 1114C W20 has length 4*N+1. Only locations 0 to N are used here. 1115c If W20(1) is nonnegative it means [D/S]MLC20 has been called and 1116c has stored into W20(1:N) estimates of the errors in the gradient 1117c values stored in G(1:N). 1118c If W20(1) is negative the user is computing the gradients and 1119c W20(1:N) is not being otherwise used. 1120c This subr stores a value into W20(0). This will only be used if 1121c [D/S]MLC20 is called. 1122c ------------------------------------------------------------------ 1123 integer N, M, LDA, IACT(*), NACT, MEQL, MDEG 1124 integer I, IZ, J, JM, K 1125 real A(LDA,*), CGRAD(*), D(*), DDOTG, DDOTGM 1126 real G(*), GM(*), GMNEW(*) 1127 real PAR(*), PARNEW(*), RELACC, SIZE, TEMP, U(*) 1128 real W20(0:*), Z(*), ZTG(*) 1129c ------------------------------------------------------------------ 1130 10 DO 20 I=1,N 1131 20 GM(I)=G(I) 1132 K=NACT 1133 30 IF (K .GT. 0) THEN 1134C 1135C Set TEMP to the next multiplier, but reduce the active set if 1136C TEMP has an unacceptable sign. 1137C 1138 TEMP=0.0E0 1139 IZ=K 1140 DO 40 I=1,N 1141 TEMP=TEMP+Z(IZ)*GM(I) 1142 40 IZ=IZ+N 1143 TEMP=TEMP*U(K) 1144 IF (K .GT. MEQL .AND. TEMP .GT. 0.0E0) THEN 1145 CALL SMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,K) 1146 GOTO 10 1147 END IF 1148C 1149C Update GM using the multiplier that has just been calculated. 1150C 1151 J=IACT(K) 1152 IF (J .LE. M) THEN 1153 DO 50 I=1,N 1154 50 GM(I)=GM(I)-TEMP*A(J,I) 1155 ELSE 1156 JM=J-M 1157 IF (JM .LE. N) THEN 1158 GM(JM)=GM(JM)+TEMP 1159 ELSE 1160 GM(JM-N)=GM(JM-N)-TEMP 1161 END IF 1162 END IF 1163 PAR(K)=TEMP 1164 K=K-1 1165 GOTO 30 1166 END IF 1167C 1168C Calculate the search direction and DDOTG. 1169C 1170 DDOTG=0.0E0 1171 IF (NACT .LT. N) THEN 1172c SDEGEN 1173 CALL SMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC, 1174 1 DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD) 1175 IF (DDOTGM .LT. 0.0E0) THEN 1176 SIZE = 0.0E0 1177 DO 60 I=1,N 1178 TEMP = D(I) * G(I) 1179 DDOTG=DDOTG + TEMP 1180 if (W20(1) .ge. 0.0E0) then 1181 SIZE = SIZE + abs(D(I)*W20(I)) 1182 else 1183 SIZE = SIZE + abs(TEMP) 1184 endif 1185 60 continue 1186 1187 if(W20(1) .lt. 0.0e0) SIZE = RELACC * SIZE 1188 if (DDOTG .ge. -SIZE) then 1189 DDOTG = 0.0E0 1190 W20(0) = 0.0E0 1191 else 1192 W20(0) = SIZE/abs(DDOTG) 1193 endif 1194 ELSE 1195 W20(0) = 0.0E0 1196 END IF 1197 END IF 1198 RETURN 1199 END 1200c ================================================================== 1201 SUBROUTINE SMLC08 (N,M,A,LDA,IACT,NACT,PAR,Z,U,D,ZTG,GM,RELACC, 1202 1 DDOTGM,MEQL,MDEG,GMNEW,PARNEW,CGRAD) 1203C 1204c SDEGEN 1205C Calculate the search direction and branch if it is not downhill. 1206c ------------------------------------------------------------------ 1207 integer N, M, LDA, IACT(*), NACT, MEQL, MDEG 1208 integer I, IDROP, ITEST, IZ, J, JM, K, KU, MP, NP 1209 real A(LDA,*), CGRAD(*), D(*), DDOTGM, DTEST 1210 real GM(*), GMNEW(*) 1211 real PAR(*), PARNEW(*), U(*) 1212 real RATIO, RELACC, SUM, TEMP, THETA 1213 real Z(*), ZTG(*) 1214c ------------------------------------------------------------------ 1215 MP=MEQL+1 1216 DTEST=0.0E0 1217C 1218 10 CALL SMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM) 1219 IF (DDOTGM .EQ. 0.0E0) GOTO 120 1220C 1221C Branch if there is no need to consider any degenerate constraints. 1222C The test gives termination if two consecutive additions to the 1223C active set fail to increase the predicted new value of F. 1224C 1225 IF (NACT .EQ. MDEG) GOTO 120 1226 NP=NACT+1 1227 SUM=0.0E0 1228 DO 20 J=NP,N 1229 20 SUM=SUM+ZTG(J)**2 1230 IF (DTEST .GT. 0.0E0 .AND. SUM .GE. DTEST) THEN 1231 IF (ITEST .EQ. 1) GOTO 120 1232 ITEST=1 1233 ELSE 1234 DTEST=SUM 1235 ITEST=0 1236 END IF 1237C 1238C Add a constraint to the active set if there are any significant 1239C violations of degenerate constraints. 1240C 1241 K=NACT 1242 CALL SMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,GMNEW,PARNEW, 1243 1 CGRAD) 1244 IF (NACT .EQ. K) GOTO 120 1245 PAR(NACT)=0.0E0 1246C 1247C Calculate the new reduced gradient and Lagrange parameters. 1248C 1249 30 DO 40 I=1,N 1250 40 GMNEW(I)=GM(I) 1251 K=NACT 1252 50 TEMP=0.0E0 1253 IZ=K 1254 DO 60 I=1,N 1255 TEMP=TEMP+Z(IZ)*GMNEW(I) 1256 60 IZ=IZ+N 1257 TEMP=TEMP*U(K) 1258 PARNEW(K)=PAR(K)+TEMP 1259 IF (K .EQ. NACT) PARNEW(K)=min(PARNEW(K),0.0E0) 1260 J=IACT(K) 1261 IF (J .LE. M) THEN 1262 DO 70 I=1,N 1263 70 GMNEW(I)=GMNEW(I)-TEMP*A(J,I) 1264 ELSE 1265 JM=J-M 1266 IF (JM .LE. N) THEN 1267 GMNEW(JM)=GMNEW(JM)+TEMP 1268 ELSE 1269 GMNEW(JM-N)=GMNEW(JM-N)-TEMP 1270 END IF 1271 END IF 1272 K=K-1 1273 IF (K .GT. MEQL) GOTO 50 1274C 1275C Set RATIO for linear interpolation between PAR and PARNEW. 1276C 1277 RATIO=0.0E0 1278 IF (MP .LT. NACT) THEN 1279 KU=NACT-1 1280 DO 80 K=MP,KU 1281 IF (PARNEW(K) .GT. 0.0E0) THEN 1282 RATIO=PARNEW(K)/(PARNEW(K)-PAR(K)) 1283 IDROP=K 1284 END IF 1285 80 CONTINUE 1286 END IF 1287C 1288C Apply the linear interpolation. 1289C 1290 THETA=1.0E0-RATIO 1291 DO 90 K=MP,NACT 1292 90 PAR(K)=min(THETA*PARNEW(K)+RATIO*PAR(K),0.0E0) 1293 DO 100 I=1,N 1294 100 GM(I)=THETA*GMNEW(I)+RATIO*GM(I) 1295C 1296C Drop a constraint if RATIO is positive. 1297C 1298 IF (RATIO .GT. 0.0E0) THEN 1299 CALL SMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) 1300 DO 110 K=IDROP,NACT 1301 110 PAR(K)=PAR(K+1) 1302 GOTO 30 1303 END IF 1304C 1305C Return if there is no freedom for a new search direction. 1306C 1307 IF (NACT .LT. N) GOTO 10 1308 DDOTGM=0.0E0 1309 120 RETURN 1310 END 1311c ================================================================== 1312 SUBROUTINE SMLC09 (SMLCFG, N,X,G,D,XS,GS,RELACC,STEPCB, DDOTG, 1313 * F,STEP, NFVALS,NFMAX,GOPT, XL, XU, GMODE, W20) 1314c LSRCH Line search 1315c 4/18/91 CLL changed handling of ICOUNT and NFVALS. Added GMODE. 1316c ------------------------------------------------------------------ 1317C Subroutine Arguments 1318c 1319c Input: SMLCFG, N, D(), RELACC, STEPCB, DDOTG, NFMAX, XL(), XU() 1320c Inout: X(), G(), F, NFVALS, GMODE, W20() 1321c Out: GS(), STEP 1322c Work space: XS(), GOPT() 1323c 1324c SMLCFG, N,X(),G(),D() 1325c XS() [out] Copy of the initial X(). 1326c GS() [out] Copy of the initial G(). 1327c RELACC,STEPCB, DDOTG, 1328c F,STEP, NFVALS,NFMAX,GOPT(), XL(), XU(), 1329c GMODE [integer,inout] For use by SMLC20. 1330c W20() [float,inout] For use by SMLC20. 1331c ------------------------------------------------------------------ 1332 external SMLCFG 1333 integer N, NFVALS, NFMAX, GMODE 1334 integer I, ICOUNT 1335 real D(*) 1336 real DDOTG, DDOTGB, DGHGH, DGKNOT, DGLOW, DGMID, DGOPT 1337 real F, FBASE, FHGH, FLOW, FOPT 1338 real G(*), GOPT(*), GS(*) 1339 real RATIO, RELACC, RELINT, SBASE, SIZE, STEP, STEPCB 1340 real STPHGH, STPLOW, STPMIN, STPOPT, TEMP 1341 real W20(0:*), X(*), XL(*), XS(*), XU(*) 1342 logical HAVEG, KEEP, NOBACK 1343c ------------------------------------------------------------------ 1344C 1345C Initialization. 1346c 1347* print'(/'' SMLC09 1290.. STEPCB='',g13.5/ 1348* * '' .. D()='',3g13.5/(22x,3g13.5))', 1349* * STEPCB, (D(I),I=1,N) 1350 RELINT=0.9E0 1351 ICOUNT=0 1352 RATIO=-1.0E0 1353 DO 10 I=1,N 1354 XS(I)=X(I) 1355 GS(I)=G(I) 1356 GOPT(I)=G(I) 1357 IF (D(I) .NE. 0.0E0) THEN 1358 TEMP=abs(X(I)/D(I)) 1359 IF (RATIO .LT. 0.0E0 .OR. TEMP .LT. RATIO) RATIO=TEMP 1360 END IF 1361 10 CONTINUE 1362 STEP=min(1.0E0,STEPCB) 1363 STPMIN=max(RELACC*RATIO,1.0E-12*STEP) 1364 STEP=max(STPMIN,STEP) 1365 SBASE=0.0E0 1366 FBASE=F 1367 DDOTGB=DDOTG 1368 STPLOW=0.0E0 1369 FLOW=F 1370 DGLOW=DDOTG 1371 STPHGH=0.0E0 1372 STPOPT=0.0E0 1373 FOPT=F 1374 DGOPT=abs(DDOTG) 1375 NOBACK = .false. 1376* print'('' SMLC09 1318.. RATIO,STPMIN='',2g13.5/ 1377* * '' DDOTG ='',g13.5)', 1378* * RATIO, STPMIN, DDOTG 1379C 1380C Calculate another function and gradient value. 1381C 1382 20 continue 1383 if(NOBACK .and. STEP .lt. STPOPT) then 1384* print*,' SMLC09 1355.. Ending SMLC09 on NOBACK.' 1385 go to 70 1386 endif 1387 NOBACK = .false. 1388 DO 30 I=1,N 1389 30 X(I)=XS(I)+STEP*D(I) 1390C%% (*smlcfg)( n, x, f, g, &haveg ); 1391 call SMLCFG(N, X, F, G, HAVEG) 1392 NFVALS = NFVALS+1 1393 if(.not. HAVEG) CALL SMLC20(SMLCFG,N,X,F,G,XL,XU,GMODE,NFVALS,W20) 1394 ICOUNT=ICOUNT+1 1395* print'(/'' SMLC09 1321.. ICOUNT='',i3/ 1396* * '' .. STEP,STPOPT ='',2g23.15)', ICOUNT, STEP, STPOPT 1397* print'('' .. X()='',2g23.15/(19x,2g23.15))', 1398* * (X(I),I=1,N) 1399* print'('' .. F ='',g23.15)',F 1400* print'('' .. G()='',2g23.15/(19x,2g23.15))', 1401* * (G(I),I=1,N) 1402 DGMID=0.0E0 1403 SIZE = 0.0e0 1404 DO 40 I=1,N 1405 TEMP = D(I)*G(I) 1406 DGMID = DGMID + TEMP 1407 SIZE = SIZE + abs(TEMP) 1408 40 continue 1409 SIZE = SIZE * RELACC 1410* print'(/'' SMLC09 1335.. FOPT - F ='',g13.5/ 1411* * '' .. DGMID ='',g13.5/ 1412* * '' .. SIZE ='',g13.5/ 1413* * '' .. DGOPT - abs(DGMID) ='',g13.5)', 1414* * FOPT - F, DGMID, SIZE, DGOPT - abs(DGMID) 1415 IF (F .LE. FOPT) THEN 1416 if(ICOUNT .eq. 1 .and. DGMID .le. 0.0e0) then 1417* print '('' SMLC09 1381.. Setting NOBACK on new test.'')' 1418 NOBACK = .true. 1419 KEEP = .true. 1420 else 1421 KEEP = F .LT. FOPT .OR. abs(DGMID) .LT. DGOPT 1422 endif 1423 if(KEEP) then 1424 STPOPT=STEP 1425* print'(/'' SMLC09 1384.. New STPOPT='',g13.5)', STPOPT 1426 FOPT=F 1427 DO 50 I=1,N 1428 50 GOPT(I)=G(I) 1429 DGOPT=abs(DGMID) 1430 END IF 1431 END IF 1432 IF (NFMAX .gt. 0 .and. NFVALS .ge. NFMAX) GOTO 70 1433C 1434C Modify the bounds on the steplength or convergence. 1435C 1436 IF (F .GE. FBASE+0.1E0*(STEP-SBASE)*DDOTGB) THEN 1437 IF (STPHGH .GT. 0.0E0 .OR. F .GT. FBASE .OR. DGMID .GT. 1438 1 0.5E0*DDOTG) THEN 1439 STPHGH=STEP 1440 FHGH=F 1441 DGHGH=DGMID 1442 GOTO 60 1443 END IF 1444 SBASE=STEP 1445 FBASE=F 1446 DDOTGB=DGMID 1447 END IF 1448 IF (DGMID .GE. 0.7E0*DDOTGB) GOTO 70 1449 STPLOW=STEP 1450 FLOW=F 1451 DGLOW=DGMID 1452 60 IF (STPHGH .GT. 0.0E0 .AND. STPLOW .GE. RELINT*STPHGH) GOTO 70 1453C 1454C Calculate the next step length or end the iterations. 1455C 1456* print'(/'' SMLC09 1366.. ICOUNT, STEP='',i23,g23.15/ 1457* * '' STPLOW, STPHGH='',2g23.15)', 1458* * ICOUNT, STEP, STPLOW, STPHGH 1459 IF (STPHGH .EQ. 0.0E0) THEN 1460 IF (STEP .EQ. STEPCB) GOTO 70 1461 TEMP=10.0E0 1462 IF (DGMID .GT. 0.9E0*DDOTG) TEMP=DDOTG/(DDOTG-DGMID) 1463 STEP=min(TEMP*STEP,STEPCB) 1464* print'('' SMLC09 1374.. To 20 with STEP='',g23.15)', STEP 1465 GOTO 20 1466 ELSE IF (ICOUNT .EQ. 1 .OR. STPLOW .GT. 0.0E0) THEN 1467 DGKNOT=2.0E0*(FHGH-FLOW)/(STPHGH-STPLOW)-0.5E0*(DGLOW+DGHGH) 1468 IF (DGKNOT .GE. 0.0E0) THEN 1469 RATIO=max(0.1E0,0.5E0*DGLOW/(DGLOW-DGKNOT)) 1470 ELSE 1471 RATIO=(0.5E0*DGHGH-DGKNOT)/(DGHGH-DGKNOT) 1472 END IF 1473 STEP=STPLOW+RATIO*(STPHGH-STPLOW) 1474* print'('' SMLC09 1384.. To 20 with STEP='',g23.15)', STEP 1475 GOTO 20 1476 ELSE 1477 STEP=0.1E0*STEP 1478 if (STEP .GE. STPMIN) then 1479* print'('' SMLC09 1389.. To 20 with STEP='',g23.15)', STEP 1480 GOTO 20 1481 endif 1482 END IF 1483C 1484C Return from subroutine. 1485C 1486 70 continue 1487 IF (STEP .NE. STPOPT) THEN 1488 STEP=STPOPT 1489 F=FOPT 1490 DO 80 I=1,N 1491 X(I)=XS(I)+STEP*D(I) 1492 80 G(I)=GOPT(I) 1493 END IF 1494 RETURN 1495 END 1496c ================================================================== 1497 SUBROUTINE SMLC10 (N,M,A,LDA,IACT,NACT,Z,U,D,RELACC,MDEG,ZZDIAG, 1498 1 GMNEW,CGRAD) 1499c NEWCON 1500c ------------------------------------------------------------------ 1501 integer N, M, LDA, IACT(*), NACT, MDEG 1502 integer I, IADD, IZ, J, JM, JMV, K, KHIGH, NP 1503 real A(LDA,*), CGRAD(*), CVIOL, CVMAX, D(*), GMNEW(*) 1504 real RELACC, SAVABS, SAVSUM, SUM, SUMABS, SUMD, TEMP 1505 real U(*) 1506 real Z(*), ZZDIAG(*) 1507c ------------------------------------------------------------------ 1508C 1509C Initialization. 1510C 1511 NP=NACT+1 1512 KHIGH=MDEG 1513 IZ=0 1514 DO 20 I=1,N 1515 ZZDIAG(I)=0.0E0 1516 DO 10 J=NP,N 1517 10 ZZDIAG(I)=ZZDIAG(I)+Z(IZ+J)**2 1518 20 IZ=IZ+N 1519C 1520C Calculate the scalar products of D with its constraints. 1521C 1522 30 CVMAX=0.0E0 1523 DO 50 K=NP,KHIGH 1524 J=IACT(K) 1525 IF (J .LE. M) THEN 1526 SUM=0.0E0 1527 SUMABS=0.0E0 1528 SUMD=0.0E0 1529 DO 40 I=1,N 1530 TEMP=D(I)*A(J,I) 1531 SUM=SUM+TEMP 1532 SUMABS=SUMABS+abs(TEMP) 1533 40 SUMD=SUMD+ZZDIAG(I)*A(J,I)**2 1534 ELSE 1535 JM=J-M 1536 IF (JM .LE. N) THEN 1537 SUM=-D(JM) 1538 ELSE 1539 JM=JM-N 1540 SUM=D(JM) 1541 END IF 1542 SUMABS=abs(SUM) 1543 SUMD=ZZDIAG(JM) 1544 END IF 1545C 1546C Pick out the most violated constraint, or return if the 1547C violation is negligible. 1548C 1549 IF (SUM .GT. RELACC*SUMABS) THEN 1550 CVIOL=SUM*SUM/SUMD 1551 IF (CVIOL .GT. CVMAX) THEN 1552 CVMAX=CVIOL 1553 IADD=K 1554 SAVSUM=SUM 1555 SAVABS=SUMABS 1556 END IF 1557 END IF 1558 50 CONTINUE 1559 IF (CVMAX .LE. 0.0E0) GOTO 140 1560 IF (NACT .EQ. 0) GOTO 120 1561C 1562C Set GMNEW to the gradient of the most violated constraint. 1563C 1564 J=IACT(IADD) 1565 IF (J .LE. M) THEN 1566 JMV=0 1567 DO 60 I=1,N 1568 60 GMNEW(I)=A(J,I) 1569 ELSE 1570 JMV=J-M 1571 DO 70 I=1,N 1572 70 GMNEW(I)=0.0E0 1573 IF (JMV .LE. N) THEN 1574 GMNEW(JMV)=-1.0E0 1575 ELSE 1576 JMV=JMV-N 1577 GMNEW(JMV)=1.0E0 1578 END IF 1579 END IF 1580C 1581C Modify GMNEW for the next active constraint. 1582C 1583 K=NACT 1584 80 TEMP=0.0E0 1585 IZ=K 1586 DO 90 I=1,N 1587 TEMP=TEMP+Z(IZ)*GMNEW(I) 1588 90 IZ=IZ+N 1589 TEMP=TEMP*U(K) 1590 J=IACT(K) 1591 IF (J .LE. M) THEN 1592 DO 100 I=1,N 1593 100 GMNEW(I)=GMNEW(I)-TEMP*A(J,I) 1594 ELSE 1595 JM=J-M 1596 IF (JM .LE. N) THEN 1597 GMNEW(JM)=GMNEW(JM)+TEMP 1598 ELSE 1599 GMNEW(JM-N)=GMNEW(JM-N)-TEMP 1600 END IF 1601 END IF 1602C 1603C Revise the values of SAVSUM and SAVABS. 1604C 1605 SUM=0.0E0 1606 SUMABS=0.0E0 1607 DO 110 I=1,N 1608 TEMP=D(I)*GMNEW(I) 1609 SUM=SUM+TEMP 1610 110 SUMABS=SUMABS+abs(TEMP) 1611 SAVSUM=min(SAVSUM,SUM) 1612 SAVABS=max(SAVABS,SUMABS) 1613 K=K-1 1614 IF (K .GE. 1) GOTO 80 1615C 1616C Add the new constraint to the active set if the constraint 1617C violation is still significant. 1618C 1619 IF (JMV .GT. 0) D(JMV)=0.0E0 1620 IF (SAVSUM .LE. RELACC*SAVABS) GOTO 130 1621 120 K=NACT 1622 CALL SMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IADD,GMNEW,CGRAD) 1623 IF (NACT .GT. K) GOTO 140 1624C 1625C Seek another constraint violation. 1626C 1627 IADD=NP 1628 130 IF (NP .LT. KHIGH) THEN 1629 K=IACT(KHIGH) 1630 IACT(KHIGH)=IACT(IADD) 1631 IACT(IADD)=K 1632 KHIGH=KHIGH-1 1633 GOTO 30 1634 END IF 1635 140 RETURN 1636 END 1637c ================================================================== 1638 SUBROUTINE SMLC11 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,INFO,Z,U,XBIG, 1639 1 RELACC,TOL,MEQL) 1640c SATACT 1641c ------------------------------------------------------------------ 1642 integer N, M, LDA, IACT(*), NACT, INFO 1643 integer I, IDROP, IZ, J, JX, K, MEQL 1644 real A(LDA,*),B(*) 1645 real RELACC, RES, RESABS, RESBIG, SAVEX, SCALE 1646 real TEMP, TEMPA, TOL, U(*) 1647 real X(*), XBIG(*), XL(*),XU(*), Z(*) 1648c ------------------------------------------------------------------ 1649 IF (NACT .EQ. 0) GOTO 50 1650 DO 30 K=1,NACT 1651C 1652C Calculate the next constraint residual. 1653C 1654 J=IACT(K) 1655 IF (J .LE. M) THEN 1656 RES=B(J) 1657 RESABS=abs(B(J)) 1658 RESBIG=RESABS 1659 DO 10 I=1,N 1660 TEMPA=A(J,I) 1661 TEMP=TEMPA*X(I) 1662 RES=RES-TEMP 1663 RESABS=RESABS+abs(TEMP) 1664 10 RESBIG=RESBIG+abs(TEMPA)*XBIG(I) 1665 ELSE 1666 JX=J-M 1667 IF (JX .LE. N) THEN 1668 RES=X(JX)-XL(JX) 1669 RESABS=abs(X(JX))+abs(XL(JX)) 1670 RESBIG=XBIG(JX)+abs(XL(JX)) 1671 SAVEX=XL(JX) 1672 ELSE 1673 JX=JX-N 1674 RES=XU(JX)-X(JX) 1675 RESABS=abs(X(JX))+abs(XU(JX)) 1676 RESBIG=XBIG(JX)+abs(XU(JX)) 1677 SAVEX=XU(JX) 1678 END IF 1679 END IF 1680C 1681C Shift X if necessary. 1682C 1683 IF (RES .NE. 0.0E0) THEN 1684 TEMP=RES/RESABS 1685 IF (K .LE. MEQL) TEMP=-abs(TEMP) 1686 IF (TOL .EQ. RELACC .OR. TEMP+RELACC .LT. 0.0E0) THEN 1687 INFO=1 1688 SCALE=RES*U(K) 1689 IZ=K 1690 DO 20 I=1,N 1691 X(I)=X(I)+SCALE*Z(IZ) 1692 IZ=IZ+N 1693 20 XBIG(I)=max(XBIG(I),abs(X(I))) 1694 IF (J .GT. M) X(JX)=SAVEX 1695C 1696C Else flag a constraint deletion if necessary. 1697C 1698 ELSE IF (RES/RESBIG .GT. TOL) THEN 1699 IACT(K)=-IACT(K) 1700 END IF 1701 END IF 1702 30 CONTINUE 1703C 1704C Delete any flagged constraints and then return. 1705C 1706 IDROP=NACT 1707 40 IF (IACT(IDROP) .LT. 0) THEN 1708 IACT(IDROP)=-IACT(IDROP) 1709 CALL SMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) 1710 END IF 1711 IDROP=IDROP-1 1712 IF (IDROP .GT. MEQL) GOTO 40 1713 50 RETURN 1714 END 1715c ================================================================== 1716 SUBROUTINE SMLC12 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,INDXBD,ZTC, 1717 1 CGRAD) 1718c ADDCON 1719c ------------------------------------------------------------------ 1720 integer N, M, LDA, IACT(*), NACT, INDXBD 1721 integer I, ICON, INEWBD, IPIV, IZ, IZNBD, J, JP, NP 1722 real A(LDA,*), CGRAD(*) 1723 real RELACC, SUM, SUMABS, TEMP, TEMPA, TEMPB 1724 real U(*), WCOS, WPIV, WSIN 1725 real Z(*), ZTC(*) 1726c ------------------------------------------------------------------ 1727 NP=NACT+1 1728 ICON=IACT(INDXBD) 1729 IACT(INDXBD)=IACT(NP) 1730 IACT(NP)=ICON 1731C 1732C Form ZTC when the new constraint is a bound. 1733C 1734 IF (ICON .GT. M) THEN 1735 INEWBD=ICON-M 1736 IF (INEWBD .LE. N) THEN 1737 TEMP=-1.0E0 1738 ELSE 1739 INEWBD=INEWBD-N 1740 TEMP=1.0E0 1741 END IF 1742 IZNBD=INEWBD*N-N 1743 DO 10 J=1,N 1744 10 ZTC(J)=TEMP*Z(IZNBD+J) 1745C 1746C Else form ZTC for an ordinary constraint. 1747C 1748 ELSE 1749 DO 20 I=1,N 1750 20 CGRAD(I)=A(ICON,I) 1751 DO 30 J=1,N 1752 ZTC(J)=0.0E0 1753 IZ=J 1754 DO 30 I=1,N 1755 ZTC(J)=ZTC(J)+Z(IZ)*CGRAD(I) 1756 30 IZ=IZ+N 1757 END IF 1758C 1759C Find any Givens rotations to apply to the last columns of Z. 1760C 1761 J=N 1762 40 JP=J 1763 J=J-1 1764 IF (J .GT. NACT) THEN 1765 IF (ZTC(JP) .EQ. 0.0E0) GOTO 40 1766 IF (abs(ZTC(JP)) .LE. RELACC*abs(ZTC(J))) THEN 1767 TEMP=abs(ZTC(J)) 1768 ELSE IF (abs(ZTC(J)) .LE. RELACC*abs(ZTC(JP))) THEN 1769 TEMP=abs(ZTC(JP)) 1770 ELSE 1771 TEMP=abs(ZTC(JP))*sqrt(1.0E0+(ZTC(J)/ZTC(JP))**2) 1772 END IF 1773 WCOS=ZTC(J)/TEMP 1774 WSIN=ZTC(JP)/TEMP 1775 ZTC(J)=TEMP 1776C 1777C Apply the rotation when the new constraint is a bound. 1778C 1779 IZ=J 1780 IF (ICON .GT. M) THEN 1781 DO 50 I=1,N 1782 TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) 1783 Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) 1784 Z(IZ+1)=TEMP 1785 50 IZ=IZ+N 1786 Z(IZNBD+JP)=0.0E0 1787C 1788C Else apply the rotation for an ordinary constraint. 1789C 1790 ELSE 1791 WPIV=0.0E0 1792 DO 60 I=1,N 1793 TEMPA=WCOS*Z(IZ+1) 1794 TEMPB=WSIN*Z(IZ) 1795 TEMP=abs(CGRAD(I))*(abs(TEMPA)+abs(TEMPB)) 1796 IF (TEMP .GT. WPIV) THEN 1797 WPIV=TEMP 1798 IPIV=I 1799 END IF 1800 Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) 1801 Z(IZ+1)=TEMPA-TEMPB 1802 60 IZ=IZ+N 1803C 1804C Ensure orthogonality of Z(.,JP) to CGRAD. 1805C 1806 SUM=0.0E0 1807 IZ=JP 1808 DO 70 I=1,N 1809 SUM=SUM+Z(IZ)*CGRAD(I) 1810 70 IZ=IZ+N 1811 IF (SUM .NE. 0.0E0) THEN 1812 IZ=IPIV*N-N+JP 1813 Z(IZ)=Z(IZ)-SUM/CGRAD(IPIV) 1814 END IF 1815 END IF 1816 GO TO 40 1817 END IF 1818C 1819C Test for linear independence in the proposed new active set. 1820C 1821 IF (ZTC(NP) .EQ. 0.0E0) GOTO 90 1822 IF (ICON .LE. M) THEN 1823 SUM=0.0E0 1824 SUMABS=0.0E0 1825 IZ=NP 1826 DO 80 I=1,N 1827 TEMP=Z(IZ)*CGRAD(I) 1828 SUM=SUM+TEMP 1829 SUMABS=SUMABS+abs(TEMP) 1830 80 IZ=IZ+N 1831 IF (abs(SUM) .LE. RELACC*SUMABS) GOTO 90 1832 END IF 1833C 1834C Set the new diagonal element of U() and return. 1835C 1836 U(NP)=1.0E0/ZTC(NP) 1837 NACT=NP 1838 90 RETURN 1839 END 1840c ================================================================== 1841 SUBROUTINE SMLC13 (N,M,A,LDA,B,XL,XU,X,IACT,NACT,XBIG,RELACC,TOL, 1842 1 MEQL) 1843C 1844C ADJTOL: Change the tolerance TOL to a smaller value, if possible. 1845c May also recompute XBIG(). 1846c ------------------------------------------------------------------ 1847 integer N, M, LDA, IACT(*), NACT, MEQL 1848 integer I, J, JM, K, KL 1849 real A(LDA,*), B(*) 1850 real RELACC, RES, RESABS, TOL, VIOL 1851 real X(*), XBIG(*), XL(*), XU(*) 1852c ------------------------------------------------------------------ 1853C Set VIOL to the greatest relative constraint residual of the first 1854C NACT constraints. 1855 VIOL=0.0E0 1856 IF (NACT .GT. MEQL) THEN 1857 KL=MEQL+1 1858 DO 20 K=KL,NACT 1859 J=IACT(K) 1860 IF (J .LE. M) THEN 1861 RES=B(J) 1862 RESABS=abs(B(J)) 1863 DO 10 I=1,N 1864 RES=RES-A(J,I)*X(I) 1865 10 RESABS=RESABS+abs(A(J,I)*XBIG(I)) 1866 ELSE 1867 JM=J-M 1868 IF (JM .LE. N) THEN 1869 RES=X(JM)-XL(JM) 1870 RESABS=XBIG(JM)+abs(XL(JM)) 1871 ELSE 1872 JM=JM-N 1873 RES=XU(JM)-X(JM) 1874 RESABS=XBIG(JM)+abs(XU(JM)) 1875 END IF 1876 END IF 1877 IF (RES .GT. 0.0E0) VIOL=max(VIOL,RES/RESABS) 1878 20 CONTINUE 1879 END IF 1880C 1881C Adjust TOL. 1882C 1883 TOL=0.1E0*min(TOL,VIOL) 1884 IF (TOL .LE. RELACC+RELACC) THEN 1885 TOL=RELACC 1886 DO 30 I=1,N 1887 30 XBIG(I)=abs(X(I)) 1888 END IF 1889 RETURN 1890 END 1891c ================================================================== 1892 SUBROUTINE SMLC14 (N,M,A,LDA,IACT,NACT,Z,U,RELACC,IDROP) 1893C 1894C DELCON: Cycle through the constraint exchanges that are needed. 1895c ------------------------------------------------------------------ 1896 integer N, M, LDA, IACT(*), NACT, IDROP 1897 integer I, IBD, ICON, IPIV, ISAVE, IZ, IZBD, J, JP, NM 1898 real A(LDA,*), DENOM 1899 real RELACC, RJJP, SUM, TEMP, TEMPA, TEMPB 1900 real U(*), UJP, WCOS, WPIV, WSIN 1901 real Z(*) 1902c ------------------------------------------------------------------ 1903 NM=NACT-1 1904 IF (IDROP .EQ. NACT) GOTO 60 1905 ISAVE=IACT(IDROP) 1906C 1907C Cycle through the constraint exchanges that are needed. 1908C 1909 DO 50 J=IDROP,NM 1910 JP=J+1 1911 ICON=IACT(JP) 1912 IACT(J)=ICON 1913C 1914C Calculate the (J,JP) element of R. 1915C 1916 IF (ICON .LE. M) THEN 1917 RJJP=0.0E0 1918 IZ=J 1919 DO 10 I=1,N 1920 RJJP=RJJP+Z(IZ)*A(ICON,I) 1921 10 IZ=IZ+N 1922 ELSE 1923 IBD=ICON-M 1924 IF (IBD .LE. N) THEN 1925 IZBD=IBD*N-N 1926 RJJP=-Z(IZBD+J) 1927 ELSE 1928 IBD=IBD-N 1929 IZBD=IBD*N-N 1930 RJJP=Z(IZBD+J) 1931 END IF 1932 END IF 1933C 1934C Calculate the parameters of the next rotation. 1935C 1936 UJP=U(JP) 1937 TEMP=RJJP*UJP 1938 DENOM=abs(TEMP) 1939 IF (DENOM*RELACC .LT. 1.0E0) DENOM=sqrt(1.0E0+DENOM*DENOM) 1940 WCOS=TEMP/DENOM 1941 WSIN=1.0E0/DENOM 1942C 1943C Rotate Z when a bound constraint is promoted. 1944C 1945 IZ=J 1946 IF (ICON .GT. M) THEN 1947 DO 20 I=1,N 1948 TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) 1949 Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) 1950 Z(IZ+1)=TEMP 1951 20 IZ=IZ+N 1952 Z(IZBD+JP)=0.0E0 1953C 1954C Rotate Z when an ordinary constraint is promoted. 1955C 1956 ELSE 1957 WPIV=0.0E0 1958 DO 30 I=1,N 1959 TEMPA=WCOS*Z(IZ+1) 1960 TEMPB=WSIN*Z(IZ) 1961 TEMP=abs(A(ICON,I))*(abs(TEMPA)+abs(TEMPB)) 1962 IF (TEMP .GT. WPIV) THEN 1963 WPIV=TEMP 1964 IPIV=I 1965 END IF 1966 Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) 1967 Z(IZ+1)=TEMPA-TEMPB 1968 30 IZ=IZ+N 1969C 1970C Ensure orthogonality to promoted constraint. 1971C 1972 SUM=0.0E0 1973 IZ=JP 1974 DO 40 I=1,N 1975 SUM=SUM+Z(IZ)*A(ICON,I) 1976 40 IZ=IZ+N 1977 IF (SUM .NE. 0.0E0) THEN 1978 IZ=IPIV*N-N+JP 1979 Z(IZ)=Z(IZ)-SUM/A(ICON,IPIV) 1980 END IF 1981 END IF 1982C 1983C Set the new diagonal elements of U. 1984C 1985 U(JP)=-DENOM*U(J) 1986 U(J)=UJP/DENOM 1987 50 CONTINUE 1988C 1989C Return. 1990C 1991 IACT(NACT)=ISAVE 1992 60 NACT=NM 1993 RETURN 1994 END 1995c ================================================================== 1996 SUBROUTINE SMLC15 (N,M,XL,XU,X,IACT,MEQL,INFO,Z,U,XBIG,RELACC) 1997c INITZU: Initialize RELACC, Z(), U(). Adjust X(). Set XBIG(). 1998c 1990-07-20 CLL Changed error handling and setting of INFO. 1999c Expect INFO to be set to 1 before entering this subr. This subr 2000c sets INFO = 4 if it finds XL(j) > XU(j) for some j. 2001c ------------------------------------------------------------------ 2002 external R1MACH 2003 integer N, M, IACT(*), MEQL, INFO 2004 integer I, IZ, J, JACT, NN 2005 real R1MACH 2006 real RELACC 2007 real U(*) 2008 real XL(*),XU(*),X(*), XBIG(*), Z(*) 2009c ------------------------------------------------------------------ 2010C 2011C Set RELACC. 2012c July 1990 CLL Commented out following 6 lines and replaced with 2013c reference to R1MACH(4) which is the smallest x such that the 2014c stored value of 1+x is greater than 1. Efforts to determine 2015c such an x with portable code such as the following 6 lines have 2016c generally eventually failed on some new computer. Also using 2017c R1MACH allows adjustments to be made for known deficiencies 2018c in particular computers, for example for the Cray X/MP & Y/MP. 2019c 2020* ZTPAR=100.0E0 2021* RELACC=1.0E0 2022* 10 RELACC=0.5E0*RELACC 2023* TEMPA=ZTPAR+0.5E0*RELACC 2024* TEMPB=ZTPAR+RELACC 2025* IF (ZTPAR .LT. TEMPA .AND. TEMPA .LT. TEMPB) GOTO 10 2026 RELACC = 100.0E0 * R1MACH(4) 2027C 2028C Seek bound inconsistencies and bound equality constraints. 2029C 2030 MEQL=0 2031 DO 20 J=1,N 2032 IF (XL(J) .GT. XU(J)) then 2033 INFO = 4 2034 call IERM1('SMLC15',INFO,0,'Bad bounds: XL(j) > XU(j)', 2035 * 'j',J,',') 2036 call SERV1('XL(j)',XL(J),',') 2037 call SERV1('XU(j)',XU(J),'.') 2038 return 2039 endif 2040 IF (XL(J) .EQ. XU(J)) MEQL=MEQL+1 2041 20 CONTINUE 2042C 2043C Initialize U, Z and XBIG. 2044C 2045 JACT=0 2046 NN=N*N 2047 DO 30 I=1,NN 2048 30 Z(I)=0.0E0 2049 IZ=0 2050 DO 40 I=1,N 2051 IF (XL(I) .EQ. XU(I)) THEN 2052 X(I)=XU(I) 2053 JACT=JACT+1 2054 U(JACT)=1.0E0 2055 IACT(JACT)=I+M+N 2056 J=JACT 2057 ELSE 2058 J=I+MEQL-JACT 2059 END IF 2060 Z(IZ+J)=1.0E0 2061 IZ=IZ+N 2062 40 XBIG(I)=abs(X(I)) 2063 RETURN 2064 END 2065c ================================================================== 2066 SUBROUTINE SMLC16 (N,M,A,LDA,IACT,NACT,PAR,G,RESKT,Z,U,BRES, 2067 * RELAXF, MEQL,SSQKT,PARW,RESKTW) 2068C 2069C Calculate the Lagrange parameters and the residual vector. 2070c 2071c Input: N, M, A(), LDA, IACT, NACT, G(), Z(), U(), BRES(), MEQL 2072c Output: PAR(), RESKT, RELAXF, SSQKT 2073c Work space: PARW(), RESKTW() 2074c ------------------------------------------------------------------ 2075 integer N, M, LDA, IACT(*), NACT, MEQL 2076 integer I, ICASE, IZ, J, JM, K, KK, KL 2077 real A(LDA,*), BRES(*), G(*) 2078 real PAR(*), PARW(*) 2079 real RELAXF, RESKT(*), RESKTW(*) 2080 real SSQKT, SSQKTW, TEMP, U(*), Z(*) 2081c ------------------------------------------------------------------ 2082 DO 10 I=1,N 2083 10 RESKT(I)=G(I) 2084* print'(/'' SMLC16 1993.. G()='',2g23.15/(19x,2g23.15))', 2085* * (G(I),I=1,N) 2086* print'('' .. U()='', 3g23.15/(8x,3g23.15))', (U(I),I=1,3) 2087* print'('' .. Z()='', 2g23.15/(8x,2g23.15))', (Z(I),I=1,16) 2088* print'('' .. A(,)='',2g23.15/(9x,2g23.15))',((A(I,J),J=1,4),I=1,3) 2089 IF (NACT .GT. 0) THEN 2090 ICASE=0 2091 20 DO 50 KK=1,NACT 2092 K=NACT+1-KK 2093 J=IACT(K) 2094 TEMP=0.0E0 2095 IZ=K 2096 DO 30 I=1,N 2097 TEMP=TEMP+Z(IZ)*RESKT(I) 2098* print'('' SMLC16 1991.. Z(),RESKT(),TEMP='',3g13.5)', 2099* * Z(IZ),RESKT(I),TEMP 2100 30 IZ=IZ+N 2101 TEMP=TEMP*U(K) 2102* print'('' SMLC16 1995.. U(K),TEMP='', 2g13.5)', U(K), TEMP 2103 IF (ICASE .EQ. 0) PAR(K)=0.0E0 2104 IF (K .LE. MEQL .OR. PAR(K)+TEMP .LT. 0.0E0) THEN 2105 PAR(K)=PAR(K)+TEMP 2106* print'('' SMLC16 1999.. TEMP,PAR(K)='',2g13.5)', 2107* * TEMP,PAR(K) 2108 ELSE 2109 TEMP=-PAR(K) 2110 PAR(K)=0.0E0 2111 END IF 2112 IF (TEMP .NE. 0.0E0) THEN 2113 IF (J .LE. M) THEN 2114 DO 40 I=1,N 2115 RESKT(I)=RESKT(I)-TEMP*A(J,I) 2116* print'('' SMLC16 2009.. TEMP,A='',2g13.5/ 2117* * '' RESKT(I) ='',g13.5)', 2118* * TEMP, A(J,I), RESKT(I) 2119 40 continue 2120 ELSE 2121 JM=J-M 2122 IF (JM .LE. N) THEN 2123 RESKT(JM)=RESKT(JM)+TEMP 2124* print'('' SMLC16 2018.. TEMP,RESKT(JM)='',2g13.5)', 2125* * TEMP,RESKT(JM) 2126 ELSE 2127 RESKT(JM-N)=RESKT(JM-N)-TEMP 2128* print'('' SMLC16 2022.. TEMP,RESKT(JM-N)='',2g13.5)', 2129* * TEMP,RESKT(JM-N) 2130 END IF 2131 END IF 2132 END IF 2133 50 CONTINUE 2134C 2135C Calculate the sum of squares of the KT residual vector. 2136C 2137 SSQKT=0.0E0 2138 IF (NACT .EQ. N) GOTO 130 2139* print'(/'' SMLC16 2014.. RESKT()='',4g12.4/(23x,4g12.4))', 2140* * (RESKT(I),I=1,N) 2141 DO 60 I=1,N 2142 60 SSQKT=SSQKT+RESKT(I)**2 2143* print'('' SMLC16 2018.. SSQKT='',g13.5)', SSQKT 2144C 2145C Apply iterative refinement to the residual vector. 2146C 2147 IF (ICASE .EQ. 0) THEN 2148 ICASE=1 2149 DO 70 K=1,NACT 2150 70 PARW(K)=PAR(K) 2151 DO 80 I=1,N 2152 80 RESKTW(I)=RESKT(I) 2153 SSQKTW=SSQKT 2154 GOTO 20 2155 END IF 2156C 2157C Undo the iterative refinement if it does not reduce SSQKT. 2158C 2159 IF (SSQKTW .LT. SSQKT) THEN 2160 DO 90 K=1,NACT 2161 90 PAR(K)=PARW(K) 2162 DO 100 I=1,N 2163 100 RESKT(I)=RESKTW(I) 2164 SSQKT=SSQKTW 2165 END IF 2166C 2167C Calculate SSQKT when there are no active constraints. 2168C 2169 ELSE 2170 SSQKT=0.0E0 2171 DO 110 I=1,N 2172 110 SSQKT=SSQKT+G(I)**2 2173 END IF 2174C 2175C Predict the reduction in F if one corrects any positive residuals 2176C of active inequality constraints. 2177C 2178 RELAXF=0.0E0 2179 IF (MEQL .LT. NACT) THEN 2180 KL=MEQL+1 2181 DO 120 K=KL,NACT 2182 J=IACT(K) 2183 IF (BRES(J) .GT. 0.0E0) THEN 2184 RELAXF=RELAXF-PAR(K)*BRES(J) 2185 END IF 2186 120 CONTINUE 2187 END IF 2188 130 RETURN 2189 END 2190c ================================================================== 2191 SUBROUTINE SMLC17 (N,NACT,Z,D,ZTG,GM,RELACC,DDOTGM) 2192c SDIRN: Compute a search direction. 2193c ------------------------------------------------------------------ 2194 integer N, NACT 2195 integer I, IZ, J, NP 2196 real D(*), DDOTGM, GM(*) 2197 real RELACC, SUM, SUMABS, TEMP 2198 real Z(*), ZTG(*) 2199c ------------------------------------------------------------------ 2200 DDOTGM=0.0E0 2201 IF (NACT .GE. N) GOTO 60 2202C 2203C Premultiply GM by the transpose of Z. 2204C 2205 NP=NACT+1 2206 DO 20 J=NP,N 2207 SUM=0.0E0 2208 SUMABS=0.0E0 2209 IZ=J 2210 DO 10 I=1,N 2211 TEMP=Z(IZ)*GM(I) 2212 SUM=SUM+TEMP 2213 SUMABS=SUMABS+abs(TEMP) 2214 10 IZ=IZ+N 2215 IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0E0 2216 20 ZTG(J)=SUM 2217C 2218C Form D by premultiplying ZTG by -Z. 2219C 2220 IZ=0 2221 DO 40 I=1,N 2222 SUM=0.0E0 2223 SUMABS=0.0E0 2224 DO 30 J=NP,N 2225 TEMP=Z(IZ+J)*ZTG(J) 2226 SUM=SUM-TEMP 2227 30 SUMABS=SUMABS+abs(TEMP) 2228 IF (abs(SUM) .LE. RELACC*SUMABS) SUM=0.0E0 2229 D(I)=SUM 2230 40 IZ=IZ+N 2231C 2232C Test that the search direction is downhill. 2233C 2234 SUMABS=0.0E0 2235 DO 50 I=1,N 2236 TEMP=D(I)*GM(I) 2237 DDOTGM=DDOTGM+TEMP 2238 50 SUMABS=SUMABS+abs(TEMP) 2239 IF (DDOTGM+RELACC*SUMABS .GE. 0.0E0) DDOTGM=0.0E0 2240 60 RETURN 2241 END 2242c ================================================================== 2243 SUBROUTINE SMLC18 (N,M,A,LDA,IACT,BRES,D,STEPCB,DDOTG,MDEG,MSAT, 2244 1 MTOT,INDXBD) 2245c STEPBD: 2246C Set steps to constraint boundaries and find the least positive 2247C one. 2248c 2249c>> 1990-04-02 CLL Changes to avoid jumping into scope of Block If. 2250c ------------------------------------------------------------------ 2251 integer N, M, LDA, IACT(*), MDEG, MSAT, MTOT, INDXBD 2252 integer I, IFLAG, J, JM, K, KL 2253 real A(LDA,*),BRES(*),D(*), DDOTG, SP, STEPCB, TEMP 2254c ------------------------------------------------------------------ 2255 IFLAG=0 2256 STEPCB=0.0E0 2257 INDXBD=0 2258 K=MDEG 2259 10 K=K+1 2260 IF (K .gt. MTOT) go to 30 2261c 2262c Use remote code to compute scalar product. 2263 go to 80 2264 20 continue 2265C 2266C Set BRES(J) to indicate the status of the j-th constraint. 2267C 2268 IF (SP*BRES(J) .LE. 0.0E0) THEN 2269 BRES(J)=0.0E0 2270 ELSE 2271 BRES(J)=BRES(J)/SP 2272 IF (STEPCB .EQ. 0.0E0 .OR. BRES(J) .LT. STEPCB) THEN 2273 STEPCB=BRES(J) 2274 INDXBD=K 2275 END IF 2276 END IF 2277 GO TO 10 2278 30 continue 2279C 2280C Try to pass through the boundary of a violated constraint. 2281C 2282 40 continue 2283 IF (INDXBD .le. MSAT) go to 70 2284 IFLAG=1 2285 K=INDXBD 2286c Use remote code to compute scalar product. 2287 go to 80 2288 50 continue 2289 MSAT=MSAT+1 2290 IACT(INDXBD)=IACT(MSAT) 2291 IACT(MSAT)=J 2292 BRES(J)=0.0E0 2293 INDXBD=MSAT 2294 DDOTG=DDOTG-SP 2295 IF (DDOTG .LT. 0.0E0 .AND. MSAT .LT. MTOT) THEN 2296C 2297C Seek the next constraint boundary along the search direction. 2298C 2299 TEMP=0.0E0 2300 KL=MDEG+1 2301 DO 60 K=KL,MTOT 2302 J=IACT(K) 2303 IF (BRES(J) .GT. 0.0E0) THEN 2304 IF (TEMP .EQ. 0.0E0 .OR. BRES(J) .LT. TEMP) THEN 2305 TEMP=BRES(J) 2306 INDXBD=K 2307 END IF 2308 END IF 2309 60 CONTINUE 2310 IF (TEMP .GT. 0.0E0) THEN 2311 STEPCB=TEMP 2312 GOTO 40 2313 END IF 2314 END IF 2315 70 continue 2316 RETURN 2317c ------------------------------------------------------------------ 2318c This is a remote block of code to compute the 2319C scalar product of D with the current constraint normal. 2320C 2321 80 J=IACT(K) 2322 IF (J .LE. M) THEN 2323 SP=0.0E0 2324 DO 90 I=1,N 2325 90 SP=SP+D(I)*A(J,I) 2326 ELSE 2327 JM=J-M 2328 IF (JM .LE. N) THEN 2329 SP=-D(JM) 2330 ELSE 2331 SP=D(JM-N) 2332 END IF 2333 END IF 2334C Return from remote block is selected by IFLAG. 2335 IF (IFLAG .EQ. 0) then 2336 go to 20 2337 ELSE 2338 go to 50 2339 ENDIF 2340 END 2341c ================================================================== 2342 SUBROUTINE SMLC19 (N,X,NACT,G,Z,ZTG,XS,GS,ZZNORM) 2343C 2344C Test if there is sufficient convexity for the update. 2345c ------------------------------------------------------------------ 2346 integer N,NACT 2347 integer I, IZ, K, KP, KM, NP 2348 real DD, DG, G(*), GS(*), SUM, TEMP, WCOS, WSIN 2349 real X(*), XS(*), Z(*),ZTG(*), ZZNORM 2350c ------------------------------------------------------------------ 2351 DD=0.0E0 2352 DG=0.0E0 2353 TEMP=0.0E0 2354 DO 10 I=1,N 2355 XS(I)=X(I)-XS(I) 2356 DD=DD+XS(I)**2 2357 TEMP=TEMP+GS(I)*XS(I) 2358 GS(I)=G(I)-GS(I) 2359 10 DG=DG+GS(I)*XS(I) 2360 IF (DG .LT. 0.1E0*abs(TEMP)) GOTO 90 2361C 2362C Transform the Z matrix. 2363C 2364 K=N 2365 20 KP=K 2366 K=K-1 2367 IF (K .GT. NACT) THEN 2368 IF (ZTG(KP) .EQ. 0.0E0) GOTO 20 2369 TEMP=abs(ZTG(KP))*sqrt(1.0E0+(ZTG(K)/ZTG(KP))**2) 2370 WCOS=ZTG(K)/TEMP 2371 WSIN=ZTG(KP)/TEMP 2372 ZTG(K)=TEMP 2373 IZ=K 2374 DO 30 I=1,N 2375 TEMP=WCOS*Z(IZ+1)-WSIN*Z(IZ) 2376 Z(IZ)=WCOS*Z(IZ)+WSIN*Z(IZ+1) 2377 Z(IZ+1)=TEMP 2378 30 IZ=IZ+N 2379 GOTO 20 2380 END IF 2381C 2382C Update the value of ZZNORM. 2383C 2384 IF (ZZNORM .LT. 0.0E0) THEN 2385 ZZNORM=DD/DG 2386 ELSE 2387 TEMP=sqrt(ZZNORM*DD/DG) 2388 ZZNORM=min(ZZNORM,TEMP) 2389 ZZNORM=max(ZZNORM,0.1E0*TEMP) 2390 END IF 2391C 2392C Complete the updating of Z. 2393C 2394 NP=NACT+1 2395 TEMP=sqrt(DG) 2396 IZ=NP 2397 DO 40 I=1,N 2398 Z(IZ)=XS(I)/TEMP 2399 40 IZ=IZ+N 2400 IF (NP .LT. N) THEN 2401 KM=NP+1 2402 DO 80 K=KM,N 2403 TEMP=0.0E0 2404 IZ=K 2405 DO 50 I=1,N 2406 TEMP=TEMP+GS(I)*Z(IZ) 2407 50 IZ=IZ+N 2408 TEMP=TEMP/DG 2409 SUM=0.0E0 2410 IZ=K 2411 DO 60 I=1,N 2412 Z(IZ)=Z(IZ)-TEMP*XS(I) 2413 SUM=SUM+Z(IZ)**2 2414 60 IZ=IZ+N 2415 IF (SUM .LT. ZZNORM) THEN 2416 TEMP=sqrt(ZZNORM/SUM) 2417 IZ=K 2418 DO 70 I=1,N 2419 Z(IZ)=TEMP*Z(IZ) 2420 70 IZ=IZ+N 2421 END IF 2422 80 CONTINUE 2423 END IF 2424 90 RETURN 2425 END 2426c ================================================================== 2427 SUBROUTINE SMLC20 (SMLCFG, N, X, F, G, XL, XU, GMODE, NFVALS, W) 2428c>> 1991-04-15 FTK & CLL Made FIRST an argument of SMLC20 2429c>> 1990-09-11 Fred T. Krogh, Jet Propulsion Laboratory, Pasadena, CA. 2430c 2431c This subr computes a finite difference approximation to a gradient 2432c vector. Values of F are not computed outside the bounds given by XL() 2433c and XU(). For any I for which XL(I) = XU(I), G(I) will be set to zero 2434c and computation to estimate G(I) will be skipped. 2435c On entry F must have already been evaluated at X(). 2436c The algorithm used is based on the following observations. 2437c 1. Since the gradient vector is being used in an iteration, there is 2438c no point to computing it more accurately than is useful in the 2439c iteration. An extra N function values to get better error 2440c estimates is not likely to pay for itself. 2441c 2. If the delta x used for the difference is in the direction of the 2442c next move, then discretization error in approximating the gradient 2443c is likely to help rather than hurt. 2444c 3. The sign of the change in X(I) for the next move is likely to be 2445c the same as it was on the last move. 2446c 4. Because of 3 and 4, the increment used is larger than would be 2447c suggested from a consideration of round off and discretization 2448c errors. 2449c 5. When the distance between values on successive iterations gets very 2450c small it is useful to have second derivative information in order 2451c to balance the effects of round off and discretization errors. 2452c 2453c *************************** Variables ******************************** 2454c 2455c D The difference between the current X(I) and the value on the 2456c previous iteration. 2457c DXN L2 norm of the move just taken. 2458c SMLCFG [in, subroutine name] Name of user-provided subroutine for 2459c computing function values. When called from this subr, SMLCFG must 2460c compute F, set HAVEG = .false., and not store into G(). 2461c EPSR Save variable = relative machine precision. 2462c F [in, float] On entry F will contain the function value 2463c evaluated at the given X(). 2464c FAC Save variable = sqrt(relative machine precision) 2465c FACX Save variable, 4. * FAC. 2466c G() [out, float] On return G() will contain an approximation of 2467c the N-dimensional gradient of f with respect to x at X() computed 2468c using one-sided finite differences. 2469c GMODE [inout, integer] Has value 0, or 1. 2470c 0 means initial entry. Compute 1-sided differences, assuming no 2471c saved info, and set GMODE to 1. 2472c 1 Compute differences, using saved info. 2473c H Increment in X(I) used for computing gradient. 2474c HB Lower bound desired for value of abs(H). 2475c KOUNT Counts interations. If multiple of 10, two sided differences 2476c are computed. 2477c LGL W(LGL+I) contains the value of G(I) from previous iteration. 2478c LNOPP Value of NOPP for start of next iteration. 2479c LPP W(LPP+I) contains estimate of the second derivative (d/dv)G(I), 2480c where v is the variable consisting of the linear combinations of 2481c X's that gave the last move. 2482c LXL W(LXL+I) contains X() from the previous iteration. 2483c MODE1 .true. if have 1-sided difference, = .false. if have 2-sided. 2484c N [in, integer] Number of components in X(), G(), XL(), & XU(). 2485c NFVALS [inout, integer] Counter of number of calls to SMLCFG. 2486c NOPP Logical variable set .true. if second derivative information is 2487c not available. 2488c TP Temporary storage. 2489c X() [in, work, float] On entry contains the current N-dimensional 2490c parameter vector. This subr will alter components in X() but will 2491c reset X() to its original contents on return. 2492c XI Current base value for X(I) when computing F for a gradient. 2493c XL() [in, float] Lower bounds for components of X(). 2494c XNEW Current value for X(I) in computing F for a gradient. Also 2495c used for temporary storage of 1-sided gradient. 2496c XU() [in, float] Upper bounds for components of X(). 2497c W() [inout, float] Work array of dimension at least 4N+1. 2498c W(0) and W(1) are initialized to -1.0 in [D/S]MLC02. If this 2499c subr is never called, i.e., if the user provides computed 2500c gradients, then W(1) will remain at this initial value. 2501c If this subr is called, it will (generallly) set W(1) positive. 2502c [D/S]MLC07 will (generally) set W(0) positive for subsequent 2503c testing when and if this subr is called. 2504c W(1:N) contains an estimate of the error in G(I) for use in 2505c [D/S]MLC07. 2506c W(N+1 : 4*N) is used as temporary work space in this subr. See the 2507c usage of the indices LGL, LPP, LXL. 2508c ------------------------------------------------------------------ 2509c 2510c ********************** Variable Declarations ************************* 2511c 2512 external R1MACH, SMLCFG 2513 integer GMODE, LGL, LPP, LXL, N 2514 integer I, KDEBUG, KOUNT, NFVALS 2515 real R1MACH 2516 real D, DXN, EPSR 2517 real F, FAC, FACX, FNEW, G(N), H, HB, TP 2518 real X(N), XI, XL(N), XNEW, XU(N), W(0:*) 2519 logical MODE1, NOPP, LNOPP 2520c%% LOGICAL32 haveg; 2521 logical HAVEG 2522 save EPSR, FAC, FACX, KOUNT, LNOPP 2523 save KDEBUG 2524 data KDEBUG / 0 / 2525 data FAC / 0.0E0 /, KOUNT / 0 / 2526 1000 format ('0Gradient #=', I3, ' F=', 1p,e24.17, ' DXN=', 2527 1 1p,e15.8 / ' I', 11x, 'X', 14X, 'DX', 10x, 'G', 10x, 'ERR', 2528 2 9X, 'H', 9X, 'dG/dV') 2529c 2530c ********************** Start of Executable Code ********************** 2531c 2532 if(FAC .eq. 0.0E0) then 2533c Get machine constants 2534 EPSR = R1MACH(4) 2535 FAC = sqrt(EPSR) 2536 FACX = 4.E0*FAC 2537 endif 2538c Set up base locations in W() 2539 LXL = N 2540 LGL = LXL + LXL 2541 LPP = LGL + LXL 2542 DXN = 0.E0 2543 if (GMODE .eq. 0) then 2544 GMODE = 1 2545 KOUNT = 1 2546 NOPP = .true. 2547 else 2548 KOUNT = KOUNT + 1 2549 if (mod(KOUNT,10) .eq. 1) then 2550 if (W(0) .gt. 1.E-5) KOUNT = 0 2551 else 2552 if (W(0) .gt. 0.01E0) KOUNT = 0 2553 end if 2554 do 10 I = 1, LXL 2555 DXN = DXN + (X(I) - W(LXL+I))**2 2556 10 continue 2557 DXN = sqrt(DXN) 2558 NOPP = LNOPP 2559 end if 2560 LNOPP = .false. 2561 if (KDEBUG .ne. 0) then 2562 print 1000, KDEBUG, F, DXN 2563 KDEBUG = KDEBUG + 1 2564 end if 2565 do 40 I = 1, LXL 2566 XI = X(I) 2567 D = XI - W(LXL+I) 2568 if (NOPP) then 2569 HB = FAC * (1.E0 + abs(XI)) 2570 else 2571 HB = FAC * (sqrt(abs(F/W(LPP+I))) + FACX*(abs(X(I))+FACX)) 2572 end if 2573 20 if ((W(LPP+I) .lt. 0.E0) .or. (mod(KOUNT, 10) .eq. 0)) then 2574 if ((XI + HB .le. XU(I)) .and. (XI - HB .ge. XL(I))) then 2575 MODE1 = .false. 2576 H = sign(min(100.E0*HB, min(XU(I)-XI, XI-XL(I))), D) 2577 go to 30 2578 end if 2579 end if 2580 MODE1 = .true. 2581 H = sign(max(HB, abs(1.E-4*D)), D) 2582 30 XNEW = XI + H 2583c Adjust H to keep XNEW within bounds. 2584 if (XNEW .gt. XU(I)) then 2585 XNEW = XU(I) 2586 if (XNEW - XI .lt. HB) then 2587 if (XI - XL(I) .gt. XNEW - XI) XNEW = max(XL(I), XI-HB) 2588 end if 2589 H = XNEW - XI 2590 else if (XNEW .lt. XL(I)) then 2591 XNEW = XL(I) 2592 if (XI - XNEW .lt. HB) then 2593 if (XU(I) - XI .gt. XI - XNEW) XNEW = min(XU(I), XI+HB) 2594 end if 2595 H = XNEW - XI 2596 end if 2597 if(H .eq. 0.E0) then 2598 G(I) = 0.E0 2599 W(I) = 0.E0 2600 else 2601 X(I) = XNEW 2602C%% (*smlcfg)( lxl, x, &fnew, g, &haveg ); 2603 CALL SMLCFG(LXL, X, FNEW, G, HAVEG) 2604 NFVALS = NFVALS+1 2605 G(I) = (FNEW - F) / H 2606 W(I) = EPSR * (abs(F) + abs(FNEW)) / abs(H) 2607 if (.not. MODE1) then 2608 XNEW = XI - H 2609 X(I) = XNEW 2610C%% (*smlcfg)( lxl, x, &fnew, g, &haveg ); 2611 CALL SMLCFG(LXL, X, FNEW, G, HAVEG) 2612 NFVALS = NFVALS+1 2613 XNEW = (FNEW - F) / (XNEW - XI) 2614 W(I) = W(I) + 5.E-4 * abs(XNEW - G(I)) 2615 G(I) = .5E0 * (G(I) + XNEW ) 2616 end if 2617 X(I) = XI 2618 end if 2619 if (DXN .eq. 0.E0) then 2620 LNOPP = .true. 2621 if (MODE1) W(I) = 0.E0 2622 else 2623 W(LPP+I) = max(0.E0, abs((G(I) - W(LGL+I)) - 10.E0 * 2624 1 abs(W(I))) / DXN) + FAC * abs(F) 2625 TP = 10.E0 * abs(H) * W(LPP+I) 2626 if (MODE1) then 2627 if (TP .gt. abs(G(I))) then 2628 if (mod(KOUNT, 10) .ne. 0) then 2629 KOUNT = 0 2630 go to 20 2631 end if 2632 end if 2633 W(I) = W(I) + TP 2634 else 2635 if (abs(XNEW-G(I)) .gt. min(10.E0*TP, 1.E-3 * abs(G(I)))) 2636 1 W(LPP+I) = -W(LPP+I) 2637 end if 2638 end if 2639 if (KDEBUG .ne. 0) then 2640 print '(I3, 1P,E22.14, E10.2, E12.4, E12.4, E11.3, E10.2)', 2641 * I,X(I),X(I)-W(LXL+I),G(I),W(I),H,W(LPP+I) 2642 end if 2643 W(LXL+I) = XI 2644 W(LGL+I) = G(I) 2645 40 continue 2646 return 2647 end 2648c ================================================================== 2649 SUBROUTINE SMLC21(IPMORE, ALL, N, ITERC, NFVALS, 2650 * F, TOL, X, G, NACT, IACT, PAR, RESKT, ENRMKT) 2651c>> 1990-07-11 CLL 2652c This subr prints results (first feasible, intermediate, and 2653c final). 2654c ------------------------------------------------------------------ 2655c Arguments 2656c 2657c IPMORE [in, integer] = 0 or 1. 0 means print ITERC, NFVALS, F, X(), 2658c and G(). 1 means also print IACT(), PAR(), and RESKT(). 2659c However see ALL for an exception. 2660c ALL [in, logical] False means G(), PAR(), and RESKT() might not 2661c contain valid data and therefore will not be printed. 2662c N, ITERC, NFVALS, 2663c F, TOL, X(), G(), NACT, IACT(), PAR(), RESKT(), 2664c ENRMKT [in] Euclidian norm of RESKT(). 2665c ------------------------------------------------------------------ 2666c%% long int i, k; 2667 integer I 2668 integer IPMORE, N, ITERC, NFVALS, NACT, IACT(*) 2669 real ENRMKT, F, X(N), G(N), PAR(*), RESKT(N), TOL 2670 logical ALL 2671c ------------------------------------------------------------------ 2672 print'(/1x,i6,'' iterations. '',i6,'' function evals. F = '', 2673 * g13.5/1x,'' TOL ='',g13.5,'' Norm of RESKT ='',g13.5)', 2674 * ITERC, NFVALS, F, TOL, ENRMKT 2675C%% printf( "\n X =" ); 2676C%% for (i = 0; i < n; i+=5){ 2677C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) 2678C%% printf( "%13.5g", x[k] ); 2679C%% if (i < n - 4) printf ( "\n " );} 2680 print'('' X ='',5g13.5/(8x,5g13.5))', (X(I),I=1,N) 2681 if (IPMORE .gt. 0) then 2682 if(ALL) then 2683C%% printf( "\n G =" ); 2684C%% for (i = 0; i < n; i+=5){ 2685C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) 2686C%% printf( "%13.5g", g[k] ); 2687C%% if (i < n - 4) printf ( "\n " );} 2688 print'(a,5g13.5/(8x,5g13.5))',' G =',(G(I),I=1,N) 2689 if (NACT .eq. N) then 2690 print*,' Kuhn-Tucker residual vector is zero.' 2691 else 2692C%% printf( "\n RESKT =" ); 2693C%% for (i = 0; i < n; i+=5){ 2694C%% for (k = i; k < (i < n - 4 ? i+5 : n); k++) 2695C%% printf( "%13.5g", reskt[k] ); 2696C%% if (i < n - 4) printf ( "\n " );} 2697 print'('' RESKT ='',5g13.5/(8x,5g13.5))', (RESKT(I),I=1,N) 2698 end if 2699 end if 2700 if (NACT .eq. 0) then 2701 print*,' No active constraints.' 2702 else 2703C%% printf( "\n IACT =" ); 2704C%% for (i = 0; i < nact; i+=12){ 2705C%% for (k = i; k < (i < nact - 11 ? i+12 : nact); k++) 2706C%% printf( "%5ld", iact[k] ); 2707C%% if (i < nact - 11) printf ( "\n " );} 2708 print'(a,12i5/(8x,12i5))',' IACT =',(IACT(I),I=1,NACT) 2709 if (ALL) then 2710C%% printf( "\n PAR =" ); 2711C%% for (i = 0; i < nact; i+=5){ 2712C%% for (k = i; k < (i < nact - 4 ? i+5 : nact); k++) 2713C%% printf( "%13.5g", par[k] ); 2714C%% if (i < nact - 4) printf ( "\n " );} 2715C%% printf( "\n" ); 2716 print'('' PAR ='',5g13.5/(8x,5g13.5))', (PAR(I),I=1,NACT) 2717 end if 2718 end if 2719 end if 2720 return 2721 end 2722 2723 2724 2725 2726 2727 2728