SUBROUTINE LINCOA (N,NPT,M,A,IA,B,X,RHOBEG,RHOEND,IPRINT, 1 MAXFUN,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(IA,*),B(*),X(*),W(*) C C This subroutine seeks the least value of a function of many variables, C subject to general linear inequality constraints, by a trust region C method that forms quadratic models by interpolation. Usually there C is much freedom in each new model after satisfying the interpolation C conditions, which is taken up by minimizing the Frobenius norm of C the change to the second derivative matrix of the model. One new C function value is calculated on each iteration, usually at a point C where the current model predicts a reduction in the least value so C far of the objective function subject to the linear constraints. C Alternatively, a new vector of variables may be chosen to replace C an interpolation point that may be too far away for reliability, and C then the new point does not have to satisfy the linear constraints. C The arguments of the subroutine are as follows. C C N must be set to the number of variables and must be at least two. C NPT must be set to the number of interpolation conditions, which is C required to be in the interval [N+2,(N+1)(N+2)/2]. Typical choices C of the author are NPT=N+6 and NPT=2*N+1. Larger values tend to be C highly inefficent when the number of variables is substantial, due C to the amount of work and extra difficulty of adjusting more points. C M must be set to the number of linear inequality constraints. C A is a matrix whose columns are the constraint gradients, which are C required to be nonzero. C IA is the first dimension of the array A, which must be at least N. C B is the vector of right hand sides of the constraints, the J-th C constraint being that the scalar product of A(.,J) with X(.) is at C most B(J). The initial vector X(.) is made feasible by increasing C the value of B(J) if necessary. C X is the vector of variables. Initial values of X(1),X(2),...,X(N) C must be supplied. If they do not satisfy the constraints, then B C is increased as mentioned above. X contains on return the variables C that have given the least calculated F subject to the constraints. C RHOBEG and RHOEND must be set to the initial and final values of a C trust region radius, so both must be positive with RHOEND<=RHOBEG. C Typically, RHOBEG should be about one tenth of the greatest expected C change to a variable, and RHOEND should indicate the accuracy that C is required in the final values of the variables. C The value of IPRINT should be set to 0, 1, 2 or 3, which controls the C amount of printing. Specifically, there is no output if IPRINT=0 and C there is output only at the return if IPRINT=1. Otherwise, the best C feasible vector of variables so far and the corresponding value of C the objective function are printed whenever RHO is reduced, where C RHO is the current lower bound on the trust region radius. Further, C each new value of F with its variables are output if IPRINT=3. C MAXFUN must be set to an upper bound on the number of calls of CALFUN, C its value being at least NPT+1. C W is an array used for working space. Its length must be at least C M*(2+N) + NPT*(4+N+NPT) + N*(9+3*N) + MAX [ M+3*N, 2*M+N, 2*NPT ]. C On return, W(1) is set to the final value of F, and W(2) is set to C the total number of function evaluations plus 0.5. C C SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set C F to the value of the objective function for the variables X(1), C X(2),...,X(N). The value of the argument F is positive when CALFUN C is called if and only if the current X satisfies the constraints C to working accuracy. C C Check that N, NPT and MAXFUN are acceptable. C ZERO=0.0D0 SMALLX=1.0D-6*RHOEND NP=N+1 NPTM=NPT-NP IF (N .LE. 1) THEN PRINT 10 10 FORMAT (/4X,'Return from LINCOA because N is less than 2.') GOTO 80 END IF IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN PRINT 20 20 FORMAT (/4X,'Return from LINCOA because NPT is not in', 1 ' the required interval.') GOTO 80 END IF IF (MAXFUN .LE. NPT) THEN PRINT 30 30 FORMAT (/4X,'Return from LINCOA because MAXFUN is less', 1 ' than NPT+1.') GOTO 80 END IF C C Normalize the constraints, and copy the resultant constraint matrix C and right hand sides into working space, after increasing the right C hand sides if necessary so that the starting point is feasible. C IAMAT=MAX0(M+3*N,2*M+N,2*NPT)+1 IB=IAMAT+M*N IFLAG=0 IF (M .GT. 0) THEN IW=IAMAT-1 DO 60 J=1,M SUM=ZERO TEMP=ZERO DO 40 I=1,N SUM=SUM+A(I,J)*X(I) 40 TEMP=TEMP+A(I,J)**2 IF (TEMP .EQ. ZERO) THEN PRINT 50 50 FORMAT (/4X,'Return from LINCOA because the gradient of', 1 ' a constraint is zero.') GOTO 80 END IF TEMP=DSQRT(TEMP) IF (SUM-B(J) .GT. SMALLX*TEMP) IFLAG=1 W(IB+J-1)=DMAX1(B(J),SUM)/TEMP DO 60 I=1,N IW=IW+1 60 W(IW)=A(I,J)/TEMP END IF IF (IFLAG .EQ. 1) THEN IF (IPRINT .GT. 0) PRINT 70 70 FORMAT (/4X,'LINCOA has made the initial X feasible by', 1 ' increasing part(s) of B.') END IF C C Partition the working space array, so that different parts of it can be C treated separately by the subroutine that performs the main calculation. C NDIM=NPT+N IXB=IB+M IXP=IXB+N IFV=IXP+N*NPT IXS=IFV+NPT IXO=IXS+N IGO=IXO+N IHQ=IGO+N IPQ=IHQ+(N*NP)/2 IBMAT=IPQ+NPT IZMAT=IBMAT+NDIM*N ISTP=IZMAT+NPT*NPTM ISP=ISTP+N IXN=ISP+NPT+NPT IAC=IXN+N IRC=IAC+N IQF=IRC+M IRF=IQF+N*N IPQW=IRF+(N*NP)/2 C C The above settings provide a partition of W for subroutine LINCOB. C CALL LINCOB (N,NPT,M,W(IAMAT),W(IB),X,RHOBEG,RHOEND,IPRINT, 1 MAXFUN,W(IXB),W(IXP),W(IFV),W(IXS),W(IXO),W(IGO),W(IHQ), 2 W(IPQ),W(IBMAT),W(IZMAT),NDIM,W(ISTP),W(ISP),W(IXN),W(IAC), 3 W(IRC),W(IQF),W(IRF),W(IPQW),W) 80 RETURN END SUBROUTINE GETACT (N,M,AMAT,B,NACT,IACT,QFAC,RFAC,SNORM, 1 RESNEW,RESACT,G,DW,VLAM,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMAT(N,*),B(*),IACT(*),QFAC(N,*),RFAC(*), 1 RESNEW(*),RESACT(*),G(*),DW(*),VLAM(*),W(*) C C N, M, AMAT, B, NACT, IACT, QFAC and RFAC are the same as the terms C with these names in SUBROUTINE LINCOB. The current values must be C set on entry. NACT, IACT, QFAC and RFAC are kept up to date when C GETACT changes the current active set. C SNORM, RESNEW, RESACT, G and DW are the same as the terms with these C names in SUBROUTINE TRSTEP. The elements of RESNEW and RESACT are C also kept up to date. C VLAM and W are used for working space, the vector VLAM being reserved C for the Lagrange multipliers of the calculation. Their lengths must C be at least N. C The main purpose of GETACT is to pick the current active set. It is C defined by the property that the projection of -G into the space C orthogonal to the active constraint normals is as large as possible, C subject to this projected steepest descent direction moving no closer C to the boundary of every constraint whose current residual is at most C 0.2*SNORM. On return, the settings in NACT, IACT, QFAC and RFAC are C all appropriate to this choice of active set. C Occasionally this projected direction is zero, and then the final value C of W(1) is set to zero. Otherwise, the direction itself is returned C in DW, and W(1) is set to the square of the length of the direction. C C Set some constants and a temporary VLAM. C ONE=1.0D0 TINY=1.0D-60 ZERO=0.0D0 TDEL=0.2D0*SNORM DDSAV=ZERO DO 10 I=1,N DDSAV=DDSAV+G(I)**2 10 VLAM(I)=ZERO DDSAV=DDSAV+DDSAV C C Set the initial QFAC to the identity matrix in the case NACT=0. C IF (NACT .EQ. 0) THEN DO 30 I=1,N DO 20 J=1,N 20 QFAC(I,J)=ZERO 30 QFAC(I,I)=ONE GOTO 100 END IF C C Remove any constraints from the initial active set whose residuals C exceed TDEL. C IFLAG=1 IC=NACT 40 IF (RESACT(IC) .GT. TDEL) GOTO 800 50 IC=IC-1 IF (IC .GT. 0) GOTO 40 C C Remove any constraints from the initial active set whose Lagrange C multipliers are nonnegative, and set the surviving multipliers. C IFLAG=2 60 IF (NACT .EQ. 0) GOTO 100 IC=NACT 70 TEMP=ZERO DO 80 I=1,N 80 TEMP=TEMP+QFAC(I,IC)*G(I) IDIAG=(IC*IC+IC)/2 IF (IC .LT. NACT) THEN JW=IDIAG+IC DO 90 J=IC+1,NACT TEMP=TEMP-RFAC(JW)*VLAM(J) 90 JW=JW+J END IF IF (TEMP .GE. ZERO) GOTO 800 VLAM(IC)=TEMP/RFAC(IDIAG) IC=IC-1 IF (IC .GT. 0) GOTO 70 C C Set the new search direction D. Terminate if the 2-norm of D is zero C or does not decrease, or if NACT=N holds. The situation NACT=N C occurs for sufficiently large SNORM if the origin is in the convex C hull of the constraint gradients. C 100 IF (NACT .EQ. N) GOTO 290 DO 110 J=NACT+1,N W(J)=ZERO DO 110 I=1,N 110 W(J)=W(J)+QFAC(I,J)*G(I) DD=ZERO DO 130 I=1,N DW(I)=ZERO DO 120 J=NACT+1,N 120 DW(I)=DW(I)-W(J)*QFAC(I,J) 130 DD=DD+DW(I)**2 IF (DD .GE. DDSAV) GOTO 290 IF (DD .EQ. ZERO) GOTO 300 DDSAV=DD DNORM=DSQRT(DD) C C Pick the next integer L or terminate, a positive value of L being C the index of the most violated constraint. The purpose of CTOL C below is to estimate whether a positive value of VIOLMX may be C due to computer rounding errors. C L=0 IF (M .GT. 0) THEN TEST=DNORM/SNORM VIOLMX=ZERO DO 150 J=1,M IF (RESNEW(J) .GT. ZERO .AND. RESNEW(J) .LE. TDEL) THEN SUM=ZERO DO 140 I=1,N 140 SUM=SUM+AMAT(I,J)*DW(I) IF (SUM .GT. TEST*RESNEW(J)) THEN IF (SUM .GT. VIOLMX) THEN L=J VIOLMX=SUM END IF END IF END IF 150 CONTINUE CTOL=ZERO TEMP=0.01D0*DNORM IF (VIOLMX .GT. ZERO .AND. VIOLMX .LT. TEMP) THEN IF (NACT .GT. 0) THEN DO 170 K=1,NACT J=IACT(K) SUM=ZERO DO 160 I=1,N 160 SUM=SUM+DW(I)*AMAT(I,J) 170 CTOL=DMAX1(CTOL,DABS(SUM)) END IF END IF END IF W(1)=ONE IF (L .EQ. 0) GOTO 300 IF (VIOLMX .LE. 10.0D0*CTOL) GOTO 300 C C Apply Givens rotations to the last (N-NACT) columns of QFAC so that C the first (NACT+1) columns of QFAC are the ones required for the C addition of the L-th constraint, and add the appropriate column C to RFAC. C NACTP=NACT+1 IDIAG=(NACTP*NACTP-NACTP)/2 RDIAG=ZERO DO 200 J=N,1,-1 SPROD=ZERO DO 180 I=1,N 180 SPROD=SPROD+QFAC(I,J)*AMAT(I,L) IF (J .LE. NACT) THEN RFAC(IDIAG+J)=SPROD ELSE IF (DABS(RDIAG) .LE. 1.0D-20*DABS(SPROD)) THEN RDIAG=SPROD ELSE TEMP=DSQRT(SPROD*SPROD+RDIAG*RDIAG) COSV=SPROD/TEMP SINV=RDIAG/TEMP RDIAG=TEMP DO 190 I=1,N TEMP=COSV*QFAC(I,J)+SINV*QFAC(I,J+1) QFAC(I,J+1)=-SINV*QFAC(I,J)+COSV*QFAC(I,J+1) 190 QFAC(I,J)=TEMP END IF END IF 200 CONTINUE IF (RDIAG .LT. ZERO) THEN DO 210 I=1,N 210 QFAC(I,NACTP)=-QFAC(I,NACTP) END IF RFAC(IDIAG+NACTP)=DABS(RDIAG) NACT=NACTP IACT(NACT)=L RESACT(NACT)=RESNEW(L) VLAM(NACT)=ZERO RESNEW(L)=ZERO C C Set the components of the vector VMU in W. C 220 W(NACT)=ONE/RFAC((NACT*NACT+NACT)/2)**2 IF (NACT .GT. 1) THEN DO 240 I=NACT-1,1,-1 IDIAG=(I*I+I)/2 JW=IDIAG+I SUM=ZERO DO 230 J=I+1,NACT SUM=SUM-RFAC(JW)*W(J) 230 JW=JW+J 240 W(I)=SUM/RFAC(IDIAG) END IF C C Calculate the multiple of VMU to subtract from VLAM, and update VLAM. C VMULT=VIOLMX IC=0 J=1 250 IF (J .LT. NACT) THEN IF (VLAM(J) .GE. VMULT*W(J)) THEN IC=J VMULT=VLAM(J)/W(J) END IF J=J+1 GOTO 250 END IF DO 260 J=1,NACT 260 VLAM(J)=VLAM(J)-VMULT*W(J) IF (IC .GT. 0) VLAM(IC)=ZERO VIOLMX=DMAX1(VIOLMX-VMULT,ZERO) IF (IC .EQ. 0) VIOLMX=ZERO C C Reduce the active set if necessary, so that all components of the C new VLAM are negative, with resetting of the residuals of the C constraints that become inactive. C IFLAG=3 IC=NACT 270 IF (VLAM(IC) .LT. ZERO) GOTO 280 RESNEW(IACT(IC))=DMAX1(RESACT(IC),TINY) GOTO 800 280 IC=IC-1 IF (IC .GT. 0) GOTO 270 C C Calculate the next VMU if VIOLMX is positive. Return if NACT=N holds, C as then the active constraints imply D=0. Otherwise, go to label C 100, to calculate the new D and to test for termination. C IF (VIOLMX .GT. ZERO) GOTO 220 IF (NACT .LT. N) GOTO 100 290 DD=ZERO 300 W(1)=DD RETURN C C These instructions rearrange the active constraints so that the new C value of IACT(NACT) is the old value of IACT(IC). A sequence of C Givens rotations is applied to the current QFAC and RFAC. Then NACT C is reduced by one. C 800 RESNEW(IACT(IC))=DMAX1(RESACT(IC),TINY) JC=IC 810 IF (JC .LT. NACT) THEN JCP=JC+1 IDIAG=JC*JCP/2 JW=IDIAG+JCP TEMP=DSQRT(RFAC(JW-1)**2+RFAC(JW)**2) CVAL=RFAC(JW)/TEMP SVAL=RFAC(JW-1)/TEMP RFAC(JW-1)=SVAL*RFAC(IDIAG) RFAC(JW)=CVAL*RFAC(IDIAG) RFAC(IDIAG)=TEMP IF (JCP .LT. NACT) THEN DO 820 J=JCP+1,NACT TEMP=SVAL*RFAC(JW+JC)+CVAL*RFAC(JW+JCP) RFAC(JW+JCP)=CVAL*RFAC(JW+JC)-SVAL*RFAC(JW+JCP) RFAC(JW+JC)=TEMP 820 JW=JW+J END IF JDIAG=IDIAG-JC DO 830 I=1,N IF (I .LT. JC) THEN TEMP=RFAC(IDIAG+I) RFAC(IDIAG+I)=RFAC(JDIAG+I) RFAC(JDIAG+I)=TEMP END IF TEMP=SVAL*QFAC(I,JC)+CVAL*QFAC(I,JCP) QFAC(I,JCP)=CVAL*QFAC(I,JC)-SVAL*QFAC(I,JCP) 830 QFAC(I,JC)=TEMP IACT(JC)=IACT(JCP) RESACT(JC)=RESACT(JCP) VLAM(JC)=VLAM(JCP) JC=JCP GOTO 810 END IF NACT=NACT-1 GOTO (50,60,280),IFLAG END SUBROUTINE LINCOB (N,NPT,M,AMAT,B,X,RHOBEG,RHOEND,IPRINT, 1 MAXFUN,XBASE,XPT,FVAL,XSAV,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM, 2 STEP,SP,XNEW,IACT,RESCON,QFAC,RFAC,PQW,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMAT(N,*),B(*),X(*),XBASE(*),XPT(NPT,*),FVAL(*), 1 XSAV(*),XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*), 2 ZMAT(NPT,*),STEP(*),SP(*),XNEW(*),IACT(*),RESCON(*), 3 QFAC(N,*),RFAC(*),PQW(*),W(*) C C The arguments N, NPT, M, X, RHOBEG, RHOEND, IPRINT and MAXFUN are C identical to the corresponding arguments in SUBROUTINE LINCOA. C AMAT is a matrix whose columns are the constraint gradients, scaled C so that they have unit length. C B contains on entry the right hand sides of the constraints, scaled C as above, but later B is modified for variables relative to XBASE. C XBASE holds a shift of origin that should reduce the contributions C from rounding errors to values of the model and Lagrange functions. C XPT contains the interpolation point coordinates relative to XBASE. C FVAL holds the values of F at the interpolation points. C XSAV holds the best feasible vector of variables so far, without any C shift of origin. C XOPT is set to XSAV-XBASE, which is the displacement from XBASE of C the feasible vector of variables that provides the least calculated C F so far, this vector being the current trust region centre. C GOPT holds the gradient of the quadratic model at XSAV = XBASE+XOPT. C HQ holds the explicit second derivatives of the quadratic model. C PQ contains the parameters of the implicit second derivatives of the C quadratic model. C BMAT holds the last N columns of the big inverse matrix H. C ZMAT holds the factorization of the leading NPT by NPT submatrix C of H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, C where the elements of DZ are plus or minus one, as specified by IDZ. C NDIM is the first dimension of BMAT and has the value NPT+N. C STEP is employed for trial steps from XOPT. It is also used for working C space when XBASE is shifted and in PRELIM. C SP is reserved for the scalar products XOPT^T XPT(K,.), K=1,2,...,NPT, C followed by STEP^T XPT(K,.), K=1,2,...,NPT. C XNEW is the displacement from XBASE of the vector of variables for C the current calculation of F, except that SUBROUTINE TRSTEP uses it C for working space. C IACT is an integer array for the indices of the active constraints. C RESCON holds useful information about the constraint residuals. Every C nonnegative RESCON(J) is the residual of the J-th constraint at the C current trust region centre. Otherwise, if RESCON(J) is negative, the C J-th constraint holds as a strict inequality at the trust region C centre, its residual being at least |RESCON(J)|; further, the value C of |RESCON(J)| is at least the current trust region radius DELTA. C QFAC is the orthogonal part of the QR factorization of the matrix of C active constraint gradients, these gradients being ordered in C accordance with IACT. When NACT is less than N, columns are added C to QFAC to complete an N by N orthogonal matrix, which is important C for keeping calculated steps sufficiently close to the boundaries C of the active constraints. C RFAC is the upper triangular part of this QR factorization, beginning C with the first diagonal element, followed by the two elements in the C upper triangular part of the second column and so on. C PQW is used for working space, mainly for storing second derivative C coefficients of quadratic functions. Its length is NPT+N. C The array W is also used for working space. The required number of C elements, namely MAX[M+3*N,2*M+N,2*NPT], is set in LINCOA. C C Set some constants. C HALF=0.5D0 ONE=1.0D0 TENTH=0.1D0 ZERO=0.0D0 NP=N+1 NH=(N*NP)/2 NPTM=NPT-NP C C Set the elements of XBASE, XPT, FVAL, XSAV, XOPT, GOPT, HQ, PQ, BMAT, C ZMAT and SP for the first iteration. An important feature is that, C if the interpolation point XPT(K,.) is not feasible, where K is any C integer from [1,NPT], then a change is made to XPT(K,.) if necessary C so that the constraint violation is at least 0.2*RHOBEG. Also KOPT C is set so that XPT(KOPT,.) is the initial trust region centre. C CALL PRELIM (N,NPT,M,AMAT,B,X,RHOBEG,IPRINT,XBASE,XPT,FVAL, 1 XSAV,XOPT,GOPT,KOPT,HQ,PQ,BMAT,ZMAT,IDZ,NDIM,SP,RESCON, 2 STEP,PQW,W) C C Begin the iterative procedure. C NF=NPT FOPT=FVAL(KOPT) RHO=RHOBEG DELTA=RHO IFEAS=0 NACT=0 ITEST=3 10 KNEW=0 NVALA=0 NVALB=0 C C Shift XBASE if XOPT may be too far from XBASE. First make the changes C to BMAT that do not depend on ZMAT. C 20 FSAVE=FOPT XOPTSQ=ZERO DO 30 I=1,N 30 XOPTSQ=XOPTSQ+XOPT(I)**2 IF (XOPTSQ .GE. 1.0D4*DELTA*DELTA) THEN QOPTSQ=0.25D0*XOPTSQ DO 50 K=1,NPT SUM=ZERO DO 40 I=1,N 40 SUM=SUM+XPT(K,I)*XOPT(I) SUM=SUM-HALF*XOPTSQ W(NPT+K)=SUM SP(K)=ZERO DO 50 I=1,N XPT(K,I)=XPT(K,I)-HALF*XOPT(I) STEP(I)=BMAT(K,I) W(I)=SUM*XPT(K,I)+QOPTSQ*XOPT(I) IP=NPT+I DO 50 J=1,I 50 BMAT(IP,J)=BMAT(IP,J)+STEP(I)*W(J)+W(I)*STEP(J) C C Then the revisions of BMAT that depend on ZMAT are calculated. C DO 90 K=1,NPTM SUMZ=ZERO DO 60 I=1,NPT SUMZ=SUMZ+ZMAT(I,K) 60 W(I)=W(NPT+I)*ZMAT(I,K) DO 80 J=1,N SUM=QOPTSQ*SUMZ*XOPT(J) DO 70 I=1,NPT 70 SUM=SUM+W(I)*XPT(I,J) STEP(J)=SUM IF (K .LT. IDZ) SUM=-SUM DO 80 I=1,NPT 80 BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K) DO 90 I=1,N IP=I+NPT TEMP=STEP(I) IF (K .LT. IDZ) TEMP=-TEMP DO 90 J=1,I 90 BMAT(IP,J)=BMAT(IP,J)+TEMP*STEP(J) C C Update the right hand sides of the constraints. C IF (M .GT. 0) THEN DO 110 J=1,M TEMP=ZERO DO 100 I=1,N 100 TEMP=TEMP+AMAT(I,J)*XOPT(I) 110 B(J)=B(J)-TEMP END IF C C The following instructions complete the shift of XBASE, including the C changes to the parameters of the quadratic model. C IH=0 DO 130 J=1,N W(J)=ZERO DO 120 K=1,NPT W(J)=W(J)+PQ(K)*XPT(K,J) 120 XPT(K,J)=XPT(K,J)-HALF*XOPT(J) DO 130 I=1,J IH=IH+1 HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J) 130 BMAT(NPT+I,J)=BMAT(NPT+J,I) DO 140 J=1,N XBASE(J)=XBASE(J)+XOPT(J) XOPT(J)=ZERO 140 XPT(KOPT,J)=ZERO END IF C C In the case KNEW=0, generate the next trust region step by calling C TRSTEP, where SNORM is the current trust region radius initially. C The final value of SNORM is the length of the calculated step, C except that SNORM is zero on return if the projected gradient is C unsuitable for starting the conjugate gradient iterations. C DELSAV=DELTA KSAVE=KNEW IF (KNEW .EQ. 0) THEN SNORM=DELTA DO 150 I=1,N 150 XNEW(I)=GOPT(I) CALL TRSTEP (N,NPT,M,AMAT,B,XPT,HQ,PQ,NACT,IACT,RESCON, 1 QFAC,RFAC,SNORM,STEP,XNEW,W,W(M+1),PQW,PQW(NP),W(M+NP)) C C A trust region step is applied whenever its length, namely SNORM, is at C least HALF*DELTA. It is also applied if its length is at least 0.1999 C times DELTA and if a line search of TRSTEP has caused a change to the C active set. Otherwise there is a branch below to label 530 or 560. C TEMP=HALF*DELTA IF (XNEW(1) .GE. HALF) TEMP=0.1999D0*DELTA IF (SNORM .LE. TEMP) THEN DELTA=HALF*DELTA IF (DELTA .LE. 1.4D0*RHO) DELTA=RHO NVALA=NVALA+1 NVALB=NVALB+1 TEMP=SNORM/RHO IF (DELSAV .GT. RHO) TEMP=ONE IF (TEMP .GE. HALF) NVALA=ZERO IF (TEMP .GE. TENTH) NVALB=ZERO IF (DELSAV .GT. RHO) GOTO 530 IF (NVALA .LT. 5 .AND. NVALB .LT. 3) GOTO 530 IF (SNORM .GT. ZERO) KSAVE=-1 GOTO 560 END IF NVALA=ZERO NVALB=ZERO C C Alternatively, KNEW is positive. Then the model step is calculated C within a trust region of radius DEL, after setting the gradient at C XBASE and the second derivative parameters of the KNEW-th Lagrange C function in W(1) to W(N) and in PQW(1) to PQW(NPT), respectively. C ELSE DEL=DMAX1(TENTH*DELTA,RHO) DO 160 I=1,N 160 W(I)=BMAT(KNEW,I) DO 170 K=1,NPT 170 PQW(K)=ZERO DO 180 J=1,NPTM TEMP=ZMAT(KNEW,J) IF (J .LT. IDZ) TEMP=-TEMP DO 180 K=1,NPT 180 PQW(K)=PQW(K)+TEMP*ZMAT(K,J) CALL QMSTEP (N,NPT,M,AMAT,B,XPT,XOPT,NACT,IACT,RESCON, 1 QFAC,KOPT,KNEW,DEL,STEP,W,PQW,W(NP),W(NP+M),IFEAS) END IF C C Set VQUAD to the change to the quadratic model when the move STEP is C made from XOPT. If STEP is a trust region step, then VQUAD should be C negative. If it is nonnegative due to rounding errors in this case, C there is a branch to label 530 to try to improve the model. C VQUAD=ZERO IH=0 DO 190 J=1,N VQUAD=VQUAD+STEP(J)*GOPT(J) DO 190 I=1,J IH=IH+1 TEMP=STEP(I)*STEP(J) IF (I .EQ. J) TEMP=HALF*TEMP 190 VQUAD=VQUAD+TEMP*HQ(IH) DO 210 K=1,NPT TEMP=ZERO DO 200 J=1,N TEMP=TEMP+XPT(K,J)*STEP(J) 200 SP(NPT+K)=TEMP 210 VQUAD=VQUAD+HALF*PQ(K)*TEMP*TEMP IF (KSAVE .EQ. 0 .AND. VQUAD .GE. ZERO) GOTO 530 C C Calculate the next value of the objective function. The difference C between the actual new value of F and the value predicted by the C model is recorded in DIFF. C 220 NF=NF+1 IF (NF .GT. MAXFUN) THEN NF=NF-1 IF (IPRINT .GT. 0) PRINT 230 230 FORMAT (/4X,'Return from LINCOA because CALFUN has been', 1 ' called MAXFUN times.') GOTO 600 END IF XDIFF=ZERO DO 240 I=1,N XNEW(I)=XOPT(I)+STEP(I) X(I)=XBASE(I)+XNEW(I) 240 XDIFF=XDIFF+(X(I)-XSAV(I))**2 XDIFF=DSQRT(XDIFF) IF (KSAVE .EQ. -1) XDIFF=RHO IF (XDIFF .LE. TENTH*RHO .OR. XDIFF .GE. DELTA+DELTA) THEN IFEAS=0 IF (IPRINT .GT. 0) PRINT 250 250 FORMAT (/4X,'Return from LINCOA because rounding errors', 1 ' prevent reasonable changes to X.') GOTO 600 END IF IF (KSAVE .LE. 0) IFEAS=1 F=DFLOAT(IFEAS) CALL CALFUN (N,X,F) IF (IPRINT .EQ. 3) THEN PRINT 260, NF,F,(X(I),I=1,N) 260 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, 1 ' The corresponding X is:'/(2X,5D15.6)) END IF IF (KSAVE .EQ. -1) GOTO 600 DIFF=F-FOPT-VQUAD C C If X is feasible, then set DFFALT to the difference between the new C value of F and the value predicted by the alternative model. C IF (IFEAS .EQ. 1 .AND. ITEST .LT. 3) THEN DO 270 K=1,NPT PQW(K)=ZERO 270 W(K)=FVAL(K)-FVAL(KOPT) DO 290 J=1,NPTM SUM=ZERO DO 280 I=1,NPT 280 SUM=SUM+W(I)*ZMAT(I,J) IF (J .LT. IDZ) SUM=-SUM DO 290 K=1,NPT 290 PQW(K)=PQW(K)+SUM*ZMAT(K,J) VQALT=ZERO DO 310 K=1,NPT SUM=ZERO DO 300 J=1,N 300 SUM=SUM+BMAT(K,J)*STEP(J) VQALT=VQALT+SUM*W(K) 310 VQALT=VQALT+PQW(K)*SP(NPT+K)*(HALF*SP(NPT+K)+SP(K)) DFFALT=F-FOPT-VQALT END IF IF (ITEST .EQ. 3) THEN DFFALT=DIFF ITEST=0 END IF C C Pick the next value of DELTA after a trust region step. C IF (KSAVE .EQ. 0) THEN RATIO=(F-FOPT)/VQUAD IF (RATIO .LE. TENTH) THEN DELTA=HALF*DELTA ELSE IF (RATIO .LE. 0.7D0) THEN DELTA=DMAX1(HALF*DELTA,SNORM) ELSE TEMP=DSQRT(2.0D0)*DELTA DELTA=DMAX1(HALF*DELTA,SNORM+SNORM) DELTA=DMIN1(DELTA,TEMP) END IF IF (DELTA .LE. 1.4D0*RHO) DELTA=RHO END IF C C Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point C can be moved. If STEP is a trust region step, then KNEW is zero at C present, but a positive value is picked by subroutine UPDATE. C CALL UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP,KOPT, 1 KNEW,PQW,W) IF (KNEW .EQ. 0) THEN IF (IPRINT .GT. 0) PRINT 320 320 FORMAT (/4X,'Return from LINCOA because the denominator' 1 ' of the updating formula is zero.') GOTO 600 END IF C C If ITEST is increased to 3, then the next quadratic model is the C one whose second derivative matrix is least subject to the new C interpolation conditions. Otherwise the new model is constructed C by the symmetric Broyden method in the usual way. C IF (IFEAS .EQ. 1) THEN ITEST=ITEST+1 IF (DABS(DFFALT) .GE. TENTH*DABS(DIFF)) ITEST=0 END IF C C Update the second derivatives of the model by the symmetric Broyden C method, using PQW for the second derivative parameters of the new C KNEW-th Lagrange function. The contribution from the old parameter C PQ(KNEW) is included in the second derivative matrix HQ. W is used C later for the gradient of the new KNEW-th Lagrange function. C IF (ITEST .LT. 3) THEN DO 330 K=1,NPT 330 PQW(K)=ZERO DO 350 J=1,NPTM TEMP=ZMAT(KNEW,J) IF (TEMP .NE. ZERO) THEN IF (J .LT. IDZ) TEMP=-TEMP DO 340 K=1,NPT 340 PQW(K)=PQW(K)+TEMP*ZMAT(K,J) END IF 350 CONTINUE IH=0 DO 360 I=1,N W(I)=BMAT(KNEW,I) TEMP=PQ(KNEW)*XPT(KNEW,I) DO 360 J=1,I IH=IH+1 360 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J) PQ(KNEW)=ZERO DO 370 K=1,NPT 370 PQ(K)=PQ(K)+DIFF*PQW(K) END IF C C Include the new interpolation point with the corresponding updates of C SP. Also make the changes of the symmetric Broyden method to GOPT at C the old XOPT if ITEST is less than 3. C FVAL(KNEW)=F SP(KNEW)=SP(KOPT)+SP(NPT+KOPT) SSQ=ZERO DO 380 I=1,N XPT(KNEW,I)=XNEW(I) 380 SSQ=SSQ+STEP(I)**2 SP(NPT+KNEW)=SP(NPT+KOPT)+SSQ IF (ITEST .LT. 3) THEN DO 390 K=1,NPT TEMP=PQW(K)*SP(K) DO 390 I=1,N 390 W(I)=W(I)+TEMP*XPT(K,I) DO 400 I=1,N 400 GOPT(I)=GOPT(I)+DIFF*W(I) END IF C C Update FOPT, XSAV, XOPT, KOPT, RESCON and SP if the new F is the C least calculated value so far with a feasible vector of variables. C IF (F .LT. FOPT .AND. IFEAS .EQ. 1) THEN FOPT=F DO 410 J=1,N XSAV(J)=X(J) 410 XOPT(J)=XNEW(J) KOPT=KNEW SNORM=DSQRT(SSQ) DO 430 J=1,M IF (RESCON(J) .GE. DELTA+SNORM) THEN RESCON(J)=SNORM-RESCON(J) ELSE RESCON(J)=RESCON(J)+SNORM IF (RESCON(J)+DELTA .GT. ZERO) THEN TEMP=B(J) DO 420 I=1,N 420 TEMP=TEMP-XOPT(I)*AMAT(I,J) TEMP=DMAX1(TEMP,ZERO) IF (TEMP .GE. DELTA) TEMP=-TEMP RESCON(J)=TEMP END IF END IF 430 CONTINUE DO 440 K=1,NPT 440 SP(K)=SP(K)+SP(NPT+K) C C Also revise GOPT when symmetric Broyden updating is applied. C IF (ITEST .LT. 3) THEN IH=0 DO 450 J=1,N DO 450 I=1,J IH=IH+1 IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*STEP(I) 450 GOPT(I)=GOPT(I)+HQ(IH)*STEP(J) DO 460 K=1,NPT TEMP=PQ(K)*SP(NPT+K) DO 460 I=1,N 460 GOPT(I)=GOPT(I)+TEMP*XPT(K,I) END IF END IF C C Replace the current model by the least Frobenius norm interpolant if C this interpolant gives substantial reductions in the predictions C of values of F at feasible points. C IF (ITEST .EQ. 3) THEN DO 470 K=1,NPT PQ(K)=ZERO 470 W(K)=FVAL(K)-FVAL(KOPT) DO 490 J=1,NPTM SUM=ZERO DO 480 I=1,NPT 480 SUM=SUM+W(I)*ZMAT(I,J) IF (J .LT. IDZ) SUM=-SUM DO 490 K=1,NPT 490 PQ(K)=PQ(K)+SUM*ZMAT(K,J) DO 500 J=1,N GOPT(J)=ZERO DO 500 I=1,NPT 500 GOPT(J)=GOPT(J)+W(I)*BMAT(I,J) DO 510 K=1,NPT TEMP=PQ(K)*SP(K) DO 510 I=1,N 510 GOPT(I)=GOPT(I)+TEMP*XPT(K,I) DO 520 IH=1,NH 520 HQ(IH)=ZERO END IF C C If a trust region step has provided a sufficient decrease in F, then C branch for another trust region calculation. Every iteration that C takes a model step is followed by an attempt to take a trust region C step. C KNEW=0 IF (KSAVE .GT. 0) GOTO 20 IF (RATIO .GE. TENTH) GOTO 20 C C Alternatively, find out if the interpolation points are close enough C to the best point so far. C 530 DISTSQ=DMAX1(DELTA*DELTA,4.0D0*RHO*RHO) DO 550 K=1,NPT SUM=ZERO DO 540 J=1,N 540 SUM=SUM+(XPT(K,J)-XOPT(J))**2 IF (SUM .GT. DISTSQ) THEN KNEW=K DISTSQ=SUM END IF 550 CONTINUE C C If KNEW is positive, then branch back for the next iteration, which C will generate a "model step". Otherwise, if the current iteration C has reduced F, or if DELTA was above its lower bound when the last C trust region step was calculated, then try a "trust region" step C instead. C IF (KNEW .GT. 0) GOTO 20 KNEW=0 IF (FOPT .LT. FSAVE) GOTO 20 IF (DELSAV .GT. RHO) GOTO 20 C C The calculations with the current value of RHO are complete. C Pick the next value of RHO. C 560 IF (RHO .GT. RHOEND) THEN DELTA=HALF*RHO IF (RHO .GT. 250.0D0*RHOEND) THEN RHO=TENTH*RHO ELSE IF (RHO .LE. 16.0D0*RHOEND) THEN RHO=RHOEND ELSE RHO=DSQRT(RHO*RHOEND) END IF DELTA=DMAX1(DELTA,RHO) IF (IPRINT .GE. 2) THEN IF (IPRINT .GE. 3) PRINT 570 570 FORMAT (5X) PRINT 580, RHO,NF 580 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of', 1 ' function values =',I6) PRINT 590, FOPT,(XBASE(I)+XOPT(I),I=1,N) 590 FORMAT (4X,'Least value of F =',1PD23.15,9X, 1 'The corresponding X is:'/(2X,5D15.6)) END IF GOTO 10 END IF C C Return from the calculation, after branching to label 220 for another C Newton-Raphson step if it has not been tried before. C IF (KSAVE .EQ. -1) GOTO 220 600 IF (FOPT .LE. F .OR. IFEAS .EQ. 0) THEN DO 610 I=1,N 610 X(I)=XSAV(I) F=FOPT END IF IF (IPRINT .GE. 1) THEN PRINT 620, NF 620 FORMAT (/4X,'At the return from LINCOA',5X, 1 'Number of function values =',I6) PRINT 590, F,(X(I),I=1,N) END IF W(1)=F W(2)=DFLOAT(NF)+HALF RETURN END SUBROUTINE PRELIM (N,NPT,M,AMAT,B,X,RHOBEG,IPRINT,XBASE, 1 XPT,FVAL,XSAV,XOPT,GOPT,KOPT,HQ,PQ,BMAT,ZMAT,IDZ,NDIM, 2 SP,RESCON,STEP,PQW,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMAT(N,*),B(*),X(*),XBASE(*),XPT(NPT,*),FVAL(*), 1 XSAV(*),XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*), 2 SP(*),RESCON(*),STEP(*),PQW(*),W(*) C C The arguments N, NPT, M, AMAT, B, X, RHOBEG, IPRINT, XBASE, XPT, FVAL, C XSAV, XOPT, GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SP and RESCON are the C same as the corresponding arguments in SUBROUTINE LINCOB. C KOPT is set to the integer such that XPT(KOPT,.) is the initial trust C region centre. C IDZ is going to be set to one, so that every element of Diag(DZ) is C one in the product ZMAT times Diag(DZ) times ZMAT^T, which is the C factorization of the leading NPT by NPT submatrix of H. C STEP, PQW and W are used for working space, the arrays STEP and PQW C being taken from LINCOB. The length of W must be at least N+NPT. C C SUBROUTINE PRELIM provides the elements of XBASE, XPT, BMAT and ZMAT C for the first iteration, an important feature being that, if any of C of the columns of XPT is an infeasible point, then the largest of C the constraint violations there is at least 0.2*RHOBEG. It also sets C the initial elements of FVAL, XOPT, GOPT, HQ, PQ, SP and RESCON. C C Set some constants. C HALF=0.5D0 ONE=1.0D0 ZERO=0.0D0 NPTM=NPT-N-1 RHOSQ=RHOBEG*RHOBEG RECIP=ONE/RHOSQ RECIQ=DSQRT(HALF)/RHOSQ TEST=0.2D0*RHOBEG IDZ=1 KBASE=1 C C Set the initial elements of XPT, BMAT, SP and ZMAT to zero. C DO 20 J=1,N XBASE(J)=X(J) DO 10 K=1,NPT 10 XPT(K,J)=ZERO DO 20 I=1,NDIM 20 BMAT(I,J)=ZERO DO 30 K=1,NPT SP(K)=ZERO DO 30 J=1,NPT-N-1 30 ZMAT(K,J)=ZERO C C Set the nonzero coordinates of XPT(K,.), K=1,2,...,min[2*N+1,NPT], C but they may be altered later to make a constraint violation C sufficiently large. The initial nonzero elements of BMAT and of C the first min[N,NPT-N-1] columns of ZMAT are set also. C DO 40 J=1,N XPT(J+1,J)=RHOBEG IF (J .LT. NPT-N) THEN JP=N+J+1 XPT(JP,J)=-RHOBEG BMAT(J+1,J)=HALF/RHOBEG BMAT(JP,J)=-HALF/RHOBEG ZMAT(1,J)=-RECIQ-RECIQ ZMAT(J+1,J)=RECIQ ZMAT(JP,J)=RECIQ ELSE BMAT(1,J)=-ONE/RHOBEG BMAT(J+1,J)=ONE/RHOBEG BMAT(NPT+J,J)=-HALF*RHOSQ END IF 40 CONTINUE C C Set the remaining initial nonzero elements of XPT and ZMAT when the C number of interpolation points exceeds 2*N+1. C IF (NPT .GT. 2*N+1) THEN DO 50 K=N+1,NPT-N-1 ITEMP=(K-1)/N IPT=K-ITEMP*N JPT=IPT+ITEMP IF (JPT .GT. N) JPT=JPT-N XPT(N+K+1,IPT)=RHOBEG XPT(N+K+1,JPT)=RHOBEG ZMAT(1,K)=RECIP ZMAT(IPT+1,K)=-RECIP ZMAT(JPT+1,K)=-RECIP 50 ZMAT(N+K+1,K)=RECIP END IF C C Update the constraint right hand sides to allow for the shift XBASE. C IF (M .GT. 0) THEN DO 70 J=1,M TEMP=ZERO DO 60 I=1,N 60 TEMP=TEMP+AMAT(I,J)*XBASE(I) 70 B(J)=B(J)-TEMP END IF C C Go through the initial points, shifting every infeasible point if C necessary so that its constraint violation is at least 0.2*RHOBEG. C DO 150 NF=1,NPT FEAS=ONE BIGV=ZERO J=0 80 J=J+1 IF (J .LE. M .AND. NF .GE. 2) THEN RESID=-B(J) DO 90 I=1,N 90 RESID=RESID+XPT(NF,I)*AMAT(I,J) IF (RESID .LE. BIGV) GOTO 80 BIGV=RESID JSAV=J IF (RESID .LE. TEST) THEN FEAS=-ONE GOTO 80 END IF FEAS=ZERO END IF IF (FEAS .LT. ZERO) THEN DO 100 I=1,N 100 STEP(I)=XPT(NF,I)+(TEST-BIGV)*AMAT(I,JSAV) DO 110 K=1,NPT SP(NPT+K)=ZERO DO 110 J=1,N 110 SP(NPT+K)=SP(NPT+K)+XPT(K,J)*STEP(J) CALL UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP, 1 KBASE,NF,PQW,W) DO 120 I=1,N 120 XPT(NF,I)=STEP(I) END IF C C Calculate the objective function at the current interpolation point, C and set KOPT to the index of the first trust region centre. C DO 130 J=1,N 130 X(J)=XBASE(J)+XPT(NF,J) F=FEAS CALL CALFUN (N,X,F) IF (IPRINT .EQ. 3) THEN PRINT 140, NF,F,(X(I),I=1,N) 140 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, 1 ' The corresponding X is:'/(2X,5D15.6)) END IF IF (NF .EQ. 1) THEN KOPT=1 ELSE IF (F .LT. FVAL(KOPT) .AND. FEAS .GT. ZERO) THEN KOPT=NF END IF 150 FVAL(NF)=F C C Set PQ for the first quadratic model. C DO 160 J=1,NPTM W(J)=ZERO DO 160 K=1,NPT 160 W(J)=W(J)+ZMAT(K,J)*FVAL(K) DO 170 K=1,NPT PQ(K)=ZERO DO 170 J=1,NPTM 170 PQ(K)=PQ(K)+ZMAT(K,J)*W(J) C C Set XOPT, SP, GOPT and HQ for the first quadratic model. C DO 180 J=1,N XOPT(J)=XPT(KOPT,J) XSAV(J)=XBASE(J)+XOPT(J) 180 GOPT(J)=ZERO DO 200 K=1,NPT SP(K)=ZERO DO 190 J=1,N 190 SP(K)=SP(K)+XPT(K,J)*XOPT(J) TEMP=PQ(K)*SP(K) DO 200 J=1,N 200 GOPT(J)=GOPT(J)+FVAL(K)*BMAT(K,J)+TEMP*XPT(K,J) DO 210 I=1,(N*N+N)/2 210 HQ(I)=ZERO C C Set the initial elements of RESCON. C DO 230 J=1,M TEMP=B(J) DO 220 I=1,N 220 TEMP=TEMP-XOPT(I)*AMAT(I,J) TEMP=DMAX1(TEMP,ZERO) IF (TEMP .GE. RHOBEG) TEMP=-TEMP 230 RESCON(J)=TEMP RETURN END SUBROUTINE QMSTEP (N,NPT,M,AMAT,B,XPT,XOPT,NACT,IACT, 1 RESCON,QFAC,KOPT,KNEW,DEL,STEP,GL,PQW,RSTAT,W,IFEAS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMAT(N,*),B(*),XPT(NPT,*),XOPT(*),IACT(*), 1 RESCON(*),QFAC(N,*),STEP(*),GL(*),PQW(*),RSTAT(*),W(*) C C N, NPT, M, AMAT, B, XPT, XOPT, NACT, IACT, RESCON, QFAC, KOPT are the C same as the terms with these names in SUBROUTINE LINCOB. C KNEW is the index of the interpolation point that is going to be moved. C DEL is the current restriction on the length of STEP, which is never C greater than the current trust region radius DELTA. C STEP will be set to the required step from XOPT to the new point. C GL must be set on entry to the gradient of LFUNC at XBASE, where LFUNC C is the KNEW-th Lagrange function. It is used also for some other C gradients of LFUNC. C PQW provides the second derivative parameters of LFUNC. C RSTAT and W are used for working space. Their lengths must be at least C M and N, respectively. RSTAT(J) is set to -1.0, 0.0, or 1.0 if the C J-th constraint is irrelevant, active, or both inactive and relevant, C respectively. C IFEAS will be set to 0 or 1 if XOPT+STEP is infeasible or feasible. C C STEP is chosen to provide a relatively large value of the modulus of C LFUNC(XOPT+STEP), subject to ||STEP|| .LE. DEL. A projected STEP is C calculated too, within the trust region, that does not alter the C residuals of the active constraints. The projected step is preferred C if its value of | LFUNC(XOPT+STEP) | is at least one fifth of the C original one, but the greatest violation of a linear constraint must C be at least 0.2*DEL, in order to keep the interpolation points apart. C The remedy when the maximum constraint violation is too small is to C restore the original step, which is perturbed if necessary so that C its maximum constraint violation becomes 0.2*DEL. C C Set some constants. C HALF=0.5D0 ONE=1.0D0 TENTH=0.1D0 ZERO=0.0D0 TEST=0.2D0*DEL C C Replace GL by the gradient of LFUNC at the trust region centre, and C set the elements of RSTAT. C DO 20 K=1,NPT TEMP=ZERO DO 10 J=1,N 10 TEMP=TEMP+XPT(K,J)*XOPT(J) TEMP=PQW(K)*TEMP DO 20 I=1,N 20 GL(I)=GL(I)+TEMP*XPT(K,I) IF (M .GT. 0) THEN DO 30 J=1,M RSTAT(J)=ONE 30 IF (DABS(RESCON(J)) .GE. DEL) RSTAT(J)=-ONE DO 40 K=1,NACT 40 RSTAT(IACT(K))=ZERO END IF C C Find the greatest modulus of LFUNC on a line through XOPT and C another interpolation point within the trust region. C IFLAG=0 VBIG=ZERO DO 60 K=1,NPT IF (K .EQ. KOPT) GOTO 60 SS=ZERO SP=ZERO DO 50 I=1,N TEMP=XPT(K,I)-XOPT(I) SS=SS+TEMP*TEMP 50 SP=SP+GL(I)*TEMP STP=-DEL/DSQRT(SS) IF (K .EQ. KNEW) THEN IF (SP*(SP-ONE) .LT. ZERO) STP=-STP VLAG=DABS(STP*SP)+STP*STP*DABS(SP-ONE) ELSE VLAG=DABS(STP*(ONE-STP)*SP) END IF IF (VLAG .GT. VBIG) THEN KSAV=K STPSAV=STP VBIG=VLAG END IF 60 CONTINUE C C Set STEP to the move that gives the greatest modulus calculated above. C This move may be replaced by a steepest ascent step from XOPT. C GG=ZERO DO 70 I=1,N GG=GG+GL(I)**2 70 STEP(I)=STPSAV*(XPT(KSAV,I)-XOPT(I)) VGRAD=DEL*DSQRT(GG) IF (VGRAD .LE. TENTH*VBIG) GOTO 220 C C Make the replacement if it provides a larger value of VBIG. C GHG=ZERO DO 90 K=1,NPT TEMP=ZERO DO 80 J=1,N 80 TEMP=TEMP+XPT(K,J)*GL(J) 90 GHG=GHG+PQW(K)*TEMP*TEMP VNEW=VGRAD+DABS(HALF*DEL*DEL*GHG/GG) IF (VNEW .GT. VBIG) THEN VBIG=VNEW STP=DEL/DSQRT(GG) IF (GHG .LT. ZERO) STP=-STP DO 100 I=1,N 100 STEP(I)=STP*GL(I) END IF IF (NACT .EQ. 0 .OR. NACT .EQ. N) GOTO 220 C C Overwrite GL by its projection. Then set VNEW to the greatest C value of |LFUNC| on the projected gradient from XOPT subject to C the trust region bound. If VNEW is sufficiently large, then STEP C may be changed to a move along the projected gradient. C DO 110 K=NACT+1,N W(K)=ZERO DO 110 I=1,N 110 W(K)=W(K)+GL(I)*QFAC(I,K) GG=ZERO DO 130 I=1,N GL(I)=ZERO DO 120 K=NACT+1,N 120 GL(I)=GL(I)+QFAC(I,K)*W(K) 130 GG=GG+GL(I)**2 VGRAD=DEL*DSQRT(GG) IF (VGRAD .LE. TENTH*VBIG) GOTO 220 GHG=ZERO DO 150 K=1,NPT TEMP=ZERO DO 140 J=1,N 140 TEMP=TEMP+XPT(K,J)*GL(J) 150 GHG=GHG+PQW(K)*TEMP*TEMP VNEW=VGRAD+DABS(HALF*DEL*DEL*GHG/GG) C C Set W to the possible move along the projected gradient. C STP=DEL/DSQRT(GG) IF (GHG .LT. ZERO) STP=-STP WW=ZERO DO 160 I=1,N W(I)=STP*GL(I) 160 WW=WW+W(I)**2 C C Set STEP to W if W gives a sufficiently large value of the modulus C of the Lagrange function, and if W either preserves feasibility C or gives a constraint violation of at least 0.2*DEL. The purpose C of CTOL below is to provide a check on feasibility that includes C a tolerance for contributions from computer rounding errors. C IF (VNEW/VBIG .GE. 0.2D0) THEN IFEAS=1 BIGV=ZERO J=0 170 J=J+1 IF (J .LE. M) THEN IF (RSTAT(J) .EQ. ONE) THEN TEMP=-RESCON(J) DO 180 I=1,N 180 TEMP=TEMP+W(I)*AMAT(I,J) BIGV=DMAX1(BIGV,TEMP) END IF IF (BIGV .LT. TEST) GOTO 170 IFEAS=0 END IF CTOL=ZERO TEMP=0.01D0*DSQRT(WW) IF (BIGV .GT. ZERO .AND. BIGV .LT. TEMP) THEN DO 200 K=1,NACT J=IACT(K) SUM=ZERO DO 190 I=1,N 190 SUM=SUM+W(I)*AMAT(I,J) 200 CTOL=DMAX1(CTOL,DABS(SUM)) END IF IF (BIGV .LE. 10.0D0*CTOL .OR. BIGV .GE. TEST) THEN DO 210 I=1,N 210 STEP(I)=W(I) GOTO 260 END IF END IF C C Calculate the greatest constraint violation at XOPT+STEP with STEP at C its original value. Modify STEP if this violation is unacceptable. C 220 IFEAS=1 BIGV=ZERO RESMAX=ZERO J=0 230 J=J+1 IF (J .LE. M) THEN IF (RSTAT(J) .LT. ZERO) GOTO 230 TEMP=-RESCON(J) DO 240 I=1,N 240 TEMP=TEMP+STEP(I)*AMAT(I,J) RESMAX=DMAX1(RESMAX,TEMP) IF (TEMP .LT. TEST) THEN IF (TEMP .LE. BIGV) GOTO 230 BIGV=TEMP JSAV=J IFEAS=-1 GOTO 230 END IF IFEAS=0 END IF IF (IFEAS .EQ. -1) THEN DO 250 I=1,N 250 STEP(I)=STEP(I)+(TEST-BIGV)*AMAT(I,JSAV) IFEAS=0 END IF C C Return the calculated STEP and the value of IFEAS. C 260 RETURN END SUBROUTINE TRSTEP (N,NPT,M,AMAT,B,XPT,HQ,PQ,NACT,IACT,RESCON, 1 QFAC,RFAC,SNORM,STEP,G,RESNEW,RESACT,D,DW,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AMAT(N,*),B(*),XPT(NPT,*),HQ(*),PQ(*),IACT(*), 1 RESCON(*),QFAC(N,*),RFAC(*),STEP(*),G(*),RESNEW(*),RESACT(*), 2 D(*),DW(*),W(*) C C N, NPT, M, AMAT, B, XPT, HQ, PQ, NACT, IACT, RESCON, QFAC and RFAC C are the same as the terms with these names in LINCOB. If RESCON(J) C is negative, then |RESCON(J)| must be no less than the trust region C radius, so that the J-th constraint can be ignored. C SNORM is set to the trust region radius DELTA initially. On the C return, however, it is the length of the calculated STEP, which is C set to zero if the constraints do not allow a long enough step. C STEP is the total calculated step so far from the trust region centre, C its final value being given by the sequence of CG iterations, which C terminate if the trust region boundary is reached. C G must be set on entry to the gradient of the quadratic model at the C trust region centre. It is used as working space, however, and is C always the gradient of the model at the current STEP, except that C on return the value of G(1) is set to ONE instead of to ZERO if C and only if GETACT is called more than once. C RESNEW, RESACT, D, DW and W are used for working space. A negative C value of RESNEW(J) indicates that the J-th constraint does not C restrict the CG steps of the current trust region calculation, a C zero value of RESNEW(J) indicates that the J-th constraint is active, C and otherwise RESNEW(J) is set to the greater of TINY and the actual C residual of the J-th constraint for the current STEP. RESACT holds C the residuals of the active constraints, which may be positive. C D is the search direction of each line search. DW is either another C search direction or the change in gradient along D. The length of W C must be at least MAX[M,2*N]. C C Set some numbers for the conjugate gradient iterations. C HALF=0.5D0 ONE=1.0D0 TINY=1.0D-60 ZERO=0.0D0 CTEST=0.01D0 SNSQ=SNORM*SNORM C C Set the initial elements of RESNEW, RESACT and STEP. C IF (M .GT. 0) THEN DO 10 J=1,M RESNEW(J)=RESCON(J) IF (RESCON(J) .GE. SNORM) THEN RESNEW(J)=-ONE ELSE IF (RESCON(J) .GE. ZERO) THEN RESNEW(J)=DMAX1(RESNEW(J),TINY) END IF 10 CONTINUE IF (NACT .GT. 0) THEN DO 20 K=1,NACT RESACT(K)=RESCON(IACT(K)) 20 RESNEW(IACT(K))=ZERO END IF END IF DO 30 I=1,N 30 STEP(I)=ZERO SS=ZERO REDUCT=ZERO NCALL=0 C C GETACT picks the active set for the current STEP. It also sets DW to C the vector closest to -G that is orthogonal to the normals of the C active constraints. DW is scaled to have length 0.2*SNORM, as then C a move of DW from STEP is allowed by the linear constraints. C 40 NCALL=NCALL+1 CALL GETACT (N,M,AMAT,B,NACT,IACT,QFAC,RFAC,SNORM,RESNEW, 1 RESACT,G,DW,W,W(N+1)) IF (W(N+1) .EQ. ZERO) GOTO 320 SCALE=0.2D0*SNORM/DSQRT(W(N+1)) DO 50 I=1,N 50 DW(I)=SCALE*DW(I) C C If the modulus of the residual of an active constraint is substantial, C then set D to the shortest move from STEP to the boundaries of the C active constraints. C RESMAX=ZERO IF (NACT .GT. 0) THEN DO 60 K=1,NACT 60 RESMAX=DMAX1(RESMAX,RESACT(K)) END IF GAMMA=ZERO IF (RESMAX .GT. 1.0D-4*SNORM) THEN IR=0 DO 80 K=1,NACT TEMP=RESACT(K) IF (K .GE. 2) THEN DO 70 I=1,K-1 IR=IR+1 70 TEMP=TEMP-RFAC(IR)*W(I) END IF IR=IR+1 80 W(K)=TEMP/RFAC(IR) DO 90 I=1,N D(I)=ZERO DO 90 K=1,NACT 90 D(I)=D(I)+W(K)*QFAC(I,K) C C The vector D that has just been calculated is also the shortest move C from STEP+DW to the boundaries of the active constraints. Set GAMMA C to the greatest steplength of this move that satisfies the trust C region bound. C RHS=SNSQ DS=ZERO DD=ZERO DO 100 I=1,N SUM=STEP(I)+DW(I) RHS=RHS-SUM*SUM DS=DS+D(I)*SUM 100 DD=DD+D(I)**2 IF (RHS .GT. ZERO) THEN TEMP=DSQRT(DS*DS+DD*RHS) IF (DS .LE. ZERO) THEN GAMMA=(TEMP-DS)/DD ELSE GAMMA=RHS/(TEMP+DS) END IF END IF C C Reduce the steplength GAMMA if necessary so that the move along D C also satisfies the linear constraints. C J=0 110 IF (GAMMA .GT. ZERO) THEN J=J+1 IF (RESNEW(J) .GT. ZERO) THEN AD=ZERO ADW=ZERO DO 120 I=1,N AD=AD+AMAT(I,J)*D(I) 120 ADW=ADW+AMAT(I,J)*DW(I) IF (AD .GT. ZERO) THEN TEMP=DMAX1((RESNEW(J)-ADW)/AD,ZERO) GAMMA=DMIN1(GAMMA,TEMP) END IF END IF IF (J .LT. M) GOTO 110 END IF GAMMA=DMIN1(GAMMA,ONE) END IF C C Set the next direction for seeking a reduction in the model function C subject to the trust region bound and the linear constraints. C IF (GAMMA .LE. ZERO) THEN DO 130 I=1,N 130 D(I)=DW(I) ICOUNT=NACT ELSE DO 140 I=1,N 140 D(I)=DW(I)+GAMMA*D(I) ICOUNT=NACT-1 END IF ALPBD=ONE C C Set ALPHA to the steplength from STEP along D to the trust region C boundary. Return if the first derivative term of this step is C sufficiently small or if no further progress is possible. C 150 ICOUNT=ICOUNT+1 RHS=SNSQ-SS IF (RHS .LE. ZERO) GOTO 320 DG=ZERO DS=ZERO DD=ZERO DO 160 I=1,N DG=DG+D(I)*G(I) DS=DS+D(I)*STEP(I) 160 DD=DD+D(I)**2 IF (DG .GE. ZERO) GOTO 320 TEMP=DSQRT(RHS*DD+DS*DS) IF (DS .LE. ZERO) THEN ALPHA=(TEMP-DS)/DD ELSE ALPHA=RHS/(TEMP+DS) END IF IF (-ALPHA*DG .LE. CTEST*REDUCT) GOTO 320 C C Set DW to the change in gradient along D. C IH=0 DO 170 J=1,N DW(J)=ZERO DO 170 I=1,J IH=IH+1 IF (I .LT. J) DW(J)=DW(J)+HQ(IH)*D(I) 170 DW(I)=DW(I)+HQ(IH)*D(J) DO 190 K=1,NPT TEMP=ZERO DO 180 J=1,N 180 TEMP=TEMP+XPT(K,J)*D(J) TEMP=PQ(K)*TEMP DO 190 I=1,N 190 DW(I)=DW(I)+TEMP*XPT(K,I) C C Set DGD to the curvature of the model along D. Then reduce ALPHA if C necessary to the value that minimizes the model. C DGD=ZERO DO 200 I=1,N 200 DGD=DGD+D(I)*DW(I) ALPHT=ALPHA IF (DG+ALPHA*DGD .GT. ZERO) THEN ALPHA=-DG/DGD END IF C C Make a further reduction in ALPHA if necessary to preserve feasibility, C and put some scalar products of D with constraint gradients in W. C ALPHM=ALPHA JSAV=0 IF (M .GT. 0) THEN DO 220 J=1,M AD=ZERO IF (RESNEW(J) .GT. ZERO) THEN DO 210 I=1,N 210 AD=AD+AMAT(I,J)*D(I) IF (ALPHA*AD .GT. RESNEW(J)) THEN ALPHA=RESNEW(J)/AD JSAV=J END IF END IF 220 W(J)=AD END IF ALPHA=DMAX1(ALPHA,ALPBD) ALPHA=DMIN1(ALPHA,ALPHM) IF (ICOUNT .EQ. NACT) ALPHA=DMIN1(ALPHA,ONE) C C Update STEP, G, RESNEW, RESACT and REDUCT. C SS=ZERO DO 230 I=1,N STEP(I)=STEP(I)+ALPHA*D(I) SS=SS+STEP(I)**2 230 G(I)=G(I)+ALPHA*DW(I) IF (M .GT. 0) THEN DO 240 J=1,M IF (RESNEW(J) .GT. ZERO) THEN RESNEW(J)=DMAX1(RESNEW(J)-ALPHA*W(J),TINY) END IF 240 CONTINUE END IF IF (ICOUNT .EQ. NACT .AND. NACT .GT. 0) THEN DO 250 K=1,NACT 250 RESACT(K)=(ONE-GAMMA)*RESACT(K) END IF REDUCT=REDUCT-ALPHA*(DG+HALF*ALPHA*DGD) C C Test for termination. Branch to label 40 if there is a new active C constraint and if the distance from STEP to the trust region C boundary is at least 0.2*SNORM. C IF (ALPHA .EQ. ALPHT) GOTO 320 TEMP=-ALPHM*(DG+HALF*ALPHM*DGD) IF (TEMP .LE. CTEST*REDUCT) GOTO 320 IF (JSAV .GT. 0) THEN IF (SS .LE. 0.64D0*SNSQ) GOTO 40 GOTO 320 END IF IF (ICOUNT .EQ. N) GOTO 320 C C Calculate the next search direction, which is conjugate to the C previous one except in the case ICOUNT=NACT. C IF (NACT .GT. 0) THEN DO 260 J=NACT+1,N W(J)=ZERO DO 260 I=1,N 260 W(J)=W(J)+G(I)*QFAC(I,J) DO 280 I=1,N TEMP=ZERO DO 270 J=NACT+1,N 270 TEMP=TEMP+QFAC(I,J)*W(J) 280 W(N+I)=TEMP ELSE DO 290 I=1,N 290 W(N+I)=G(I) END IF IF (ICOUNT .EQ. NACT) THEN BETA=ZERO ELSE WGD=ZERO DO 300 I=1,N 300 WGD=WGD+W(N+I)*DW(I) BETA=WGD/DGD END IF DO 310 I=1,N 310 D(I)=-W(N+I)+BETA*D(I) ALPBD=ZERO GOTO 150 C C Return from the subroutine. C 320 SNORM=ZERO IF (REDUCT .GT. ZERO) SNORM=DSQRT(SS) G(1)=ZERO IF (NCALL .GT. 1) G(1)=ONE RETURN END SUBROUTINE UPDATE (N,NPT,XPT,BMAT,ZMAT,IDZ,NDIM,SP,STEP, 1 KOPT,KNEW,VLAG,W) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),SP(*),STEP(*), 2 VLAG(*),W(*) C C The arguments N, NPT, XPT, BMAT, ZMAT, IDZ, NDIM ,SP and STEP are C identical to the corresponding arguments in SUBROUTINE LINCOB. C KOPT is such that XPT(KOPT,.) is the current trust region centre. C KNEW on exit is usually positive, and then it is the index of an C interpolation point to be moved to the position XPT(KOPT,.)+STEP(.). C It is set on entry either to its final value or to 0. In the latter C case, the final value of KNEW is chosen to maximize the denominator C of the matrix updating formula times a weighting factor. C VLAG and W are used for working space, the first NPT+N elements of C both of these vectors being required. C C The arrays BMAT and ZMAT with IDZ are updated, the new matrices being C the ones that are suitable after the shift of the KNEW-th point to C the new position XPT(KOPT,.)+STEP(.). A return with KNEW set to zero C occurs if the calculation fails due to a zero denominator in the C updating formula, which should never happen. C C Set some constants. C HALF=0.5D0 ONE=1.0D0 ZERO=0.0D0 NPTM=NPT-N-1 C C Calculate VLAG and BETA for the current choice of STEP. The first NPT C elements of VLAG are set to the values of the Lagrange functions at C XPT(KOPT,.)+STEP(.). The first NPT components of W_check are held C in W, where W_check is defined in a paper on the updating method. C DO 20 K=1,NPT W(K)=SP(NPT+K)*(HALF*SP(NPT+K)+SP(K)) SUM=ZERO DO 10 J=1,N 10 SUM=SUM+BMAT(K,J)*STEP(J) 20 VLAG(K)=SUM BETA=ZERO DO 40 K=1,NPTM SUM=ZERO DO 30 I=1,NPT 30 SUM=SUM+ZMAT(I,K)*W(I) IF (K .LT. IDZ) THEN BETA=BETA+SUM*SUM SUM=-SUM ELSE BETA=BETA-SUM*SUM END IF DO 40 I=1,NPT 40 VLAG(I)=VLAG(I)+SUM*ZMAT(I,K) BSUM=ZERO DX=ZERO SSQ=ZERO DO 70 J=1,N SUM=ZERO DO 50 I=1,NPT 50 SUM=SUM+W(I)*BMAT(I,J) BSUM=BSUM+SUM*STEP(J) JP=NPT+J DO 60 K=1,N 60 SUM=SUM+BMAT(JP,K)*STEP(K) VLAG(JP)=SUM BSUM=BSUM+SUM*STEP(J) DX=DX+STEP(J)*XPT(KOPT,J) 70 SSQ=SSQ+STEP(J)**2 BETA=DX*DX+SSQ*(SP(KOPT)+DX+DX+HALF*SSQ)+BETA-BSUM VLAG(KOPT)=VLAG(KOPT)+ONE C C If KNEW is zero initially, then pick the index of the interpolation C point to be deleted, by maximizing the absolute value of the C denominator of the updating formula times a weighting factor. C C IF (KNEW .EQ. 0) THEN DENMAX=ZERO DO 100 K=1,NPT HDIAG=ZERO DO 80 J=1,NPTM TEMP=ONE IF (J .LT. IDZ) TEMP=-ONE 80 HDIAG=HDIAG+TEMP*ZMAT(K,J)**2 DENABS=DABS(BETA*HDIAG+VLAG(K)**2) DISTSQ=ZERO DO 90 J=1,N 90 DISTSQ=DISTSQ+(XPT(K,J)-XPT(KOPT,J))**2 TEMP=DENABS*DISTSQ*DISTSQ IF (TEMP .GT. DENMAX) THEN DENMAX=TEMP KNEW=K END IF 100 CONTINUE END IF C C Apply the rotations that put zeros in the KNEW-th row of ZMAT. C JL=1 IF (NPTM .GE. 2) THEN DO 120 J=2,NPTM IF (J .EQ. IDZ) THEN JL=IDZ ELSE IF (ZMAT(KNEW,J) .NE. ZERO) THEN TEMP=DSQRT(ZMAT(KNEW,JL)**2+ZMAT(KNEW,J)**2) TEMPA=ZMAT(KNEW,JL)/TEMP TEMPB=ZMAT(KNEW,J)/TEMP DO 110 I=1,NPT TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J) ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL) 110 ZMAT(I,JL)=TEMP ZMAT(KNEW,J)=ZERO END IF 120 CONTINUE END IF C C Put the first NPT components of the KNEW-th column of the Z Z^T matrix C into W, and calculate the parameters of the updating formula. C TEMPA=ZMAT(KNEW,1) IF (IDZ .GE. 2) TEMPA=-TEMPA IF (JL .GT. 1) TEMPB=ZMAT(KNEW,JL) DO 130 I=1,NPT W(I)=TEMPA*ZMAT(I,1) IF (JL .GT. 1) W(I)=W(I)+TEMPB*ZMAT(I,JL) 130 CONTINUE ALPHA=W(KNEW) TAU=VLAG(KNEW) TAUSQ=TAU*TAU DENOM=ALPHA*BETA+TAUSQ VLAG(KNEW)=VLAG(KNEW)-ONE IF (DENOM .EQ. ZERO) THEN KNEW=0 GOTO 180 END IF SQRTDN=DSQRT(DABS(DENOM)) C C Complete the updating of ZMAT when there is only one nonzero element C in the KNEW-th row of the new matrix ZMAT. IFLAG is set to one when C the value of IDZ is going to be reduced. C IFLAG=0 IF (JL .EQ. 1) THEN TEMPA=TAU/SQRTDN TEMPB=ZMAT(KNEW,1)/SQRTDN DO 140 I=1,NPT 140 ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I) IF (DENOM .LT. ZERO) THEN IF (IDZ .EQ. 1) THEN IDZ=2 ELSE IFLAG=1 END IF END IF ELSE C C Complete the updating of ZMAT in the alternative case. C JA=1 IF (BETA .GE. ZERO) JA=JL JB=JL+1-JA TEMP=ZMAT(KNEW,JB)/DENOM TEMPA=TEMP*BETA TEMPB=TEMP*TAU TEMP=ZMAT(KNEW,JA) SCALA=ONE/DSQRT(DABS(BETA)*TEMP*TEMP+TAUSQ) SCALB=SCALA*SQRTDN DO 150 I=1,NPT ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I)) 150 ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I)) IF (DENOM .LE. ZERO) THEN IF (BETA .LT. ZERO) THEN IDZ=IDZ+1 ELSE IFLAG=1 END IF END IF END IF C C Reduce IDZ when the diagonal part of the ZMAT times Diag(DZ) times C ZMAT^T factorization gains another positive element. Then exchange C the first and IDZ-th columns of ZMAT. C IF (IFLAG .EQ. 1) THEN IDZ=IDZ-1 DO 160 I=1,NPT TEMP=ZMAT(I,1) ZMAT(I,1)=ZMAT(I,IDZ) 160 ZMAT(I,IDZ)=TEMP END IF C C Finally, update the matrix BMAT. C DO 170 J=1,N JP=NPT+J W(JP)=BMAT(KNEW,J) TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM DO 170 I=1,JP BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I) IF (I .GT. NPT) BMAT(JP,I-NPT)=BMAT(I,J) 170 CONTINUE 180 RETURN END