The Fortran version of LINCOA is attached. Its purpose is to seek the least value of a function F of several variables subject to general linear inequality constraints on the variables, when derivatives of F are not available. The name LINCOA denotes LINearly Constrained Optimization Algorithm. F is specified by the user through a subroutine called CALFUN. The algorithm is intended to change the variables to values that are close to a local constrained minimum of F. The user, however, should assume responsibility for finding out if the calculations are adequate. It may be helpful to employ several starting points in the space of the variables and to try different values of the parameters NPT and RHOEND. I intend to write a paper that explains briefly the main features of the software. LINCOA is not suitable for very large numbers of variables because no attention is given to any sparsity. A few calculations with 1000 variables, however, have been run successfully overnight, and the performance of LINCOA is satisfactory usually for small numbers of variables. Several calculations of the objective function may be required at points that do not satisfy the linear constraints, especially if an equality constraint is expressed as two inequalities. The attachments in sequence are a suitable Makefile, followed by a main program and a CALFUN routine for the "PtsinTet" problem, in order to provide an example for testing. Then LINCOA and its six auxiliary routines, namely LINCOB, GETACT, PRELIM, QMSTEP, TRSTEP and UPDATE, are given. Finally, the output from the author's computer for the PtsinTet problem is listed. In addition to providing CALFUN, the linear constraints, and an initial vector of variables, the user has to set the values of RHOBEG, RHOEND and NPT. After scaling the individual variables, so that the magnitudes of their expected changes are similar, RHOBEG is the initial steplength for changes to the variables, a reasonable choice being the mesh size of a coarse grid search. Further, RHOEND should be suitable for a search on a very fine grid. Typically, the final vector of variables is within distance 10*RHOEND of a local minimum. The parameter NPT specifies the number of interpolation conditions on each quadratic model, the value NPT=2*N+1 being recommended for a start, where N is the number of variables. The way of calling LINCOA should be clear from the PtsinTet example and from the comments near the beginning of SUBROUTINE LINCOA. There are no restrictions on or charges for the use of the software. I hope that the time and effort I have spent on developing the package will be helpful to much research and to many applications. December 6th, 2013 M.J.D. Powell (mjdp@cam.ac.uk) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Makefile %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% main: calfun.o getact.o lincoa.o lincob.o main.o prelim.o qmstep.o trstep.o update.o ifort -g calfun.o getact.o lincoa.o lincob.o main.o prelim.o qmstep.o trstep.o update.o calfun.o: calfun.f ifort -g -c calfun.f getact.o: getact.f ifort -g -c getact.f lincoa.o: lincoa.f ifort -g -c lincoa.f lincob.o: lincob.f ifort -g -c lincob.f main.o: main.f ifort -g -c main.f prelim.o: prelim.f ifort -g -c prelim.f qmstep.o: qmstep.f ifort -g -c qmstep.f trstep.o: trstep.f ifort -g -c trstep.f update.o: update.f ifort -g -c update.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% main.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Calculate the tetrahedron of least volume that encloses the points C (XP(J),YP(J),ZP(J)), J=1,2,...,NP. Our method requires the origin C to be strictly inside the convex hull of these points. There are C twelve variables that define the four faces of each tetrahedron C that is considered. Each face has the form ALPHA*X+BETA*Y+GAMMA*Z=1, C the variables X(3K-2), X(3K-1) and X(3K) being the values of ALPHA, C BETA and GAMMA for the K-th face, K=1,2,3,4. Let the set T contain C all points in three dimensions that can be reached from the origin C without crossing a face. Because the volume of T may be infinite, C the objective function is the smaller of FMAX and the volume of T, C where FMAX is set to an upper bound on the final volume initially. C There are 4*NP linear constraints on the variables, namely that each C of the given points (XP(J),YP(J),ZP(J)) shall be in T. Let XS = min C XP(J), YS = min YP(J), ZS = min ZP(J) and SS = max XP(J)+YP(J)+ZP(J), C where J runs from 1 to NP. The initial values of the variables are C X(1)=1/XS, X(5)=1/YS, X(9)=1/ZS, X(2)=X(3)=X(4)=X(6)=X(7) =X(8)=0 C and X(10)=X(11)=X(12)=1/SS, which satisfy the linear constraints, C and which provide the bound FMAX=(SS-XS-YS-ZS)**3/6. Other details C of the test calculation are given below, including the choice of C the data points (XP(J),YP(J),ZP(J)), J=1,2,...,NP. The smaller final C value of the objective function in the case NPT=35 shows that the C problem has local minima. C IMPLICIT REAL*8 (A-H,O-Z) COMMON FMAX DIMENSION XP(50),YP(50),ZP(50),A(12,200),B(200),X(12),W(500000) C C Set some constants. C ONE=1.0D0 TWO=2.0D0 ZERO=0.0D0 PI=4.0D0*DATAN(ONE) IA=12 N=12 C C Set the data points. C NP=50 SUMX=ZERO SUMY=ZERO SUMZ=ZERO DO 10 J=1,NP THETA=DFLOAT(J-1)*PI/DFLOAT(NP-1) XP(J)=DCOS(THETA)*DCOS(TWO*THETA) SUMX=SUMX+XP(J) YP(J)=DSIN(THETA)*DCOS(TWO*THETA) SUMY=SUMY+YP(J) ZP(J)=DSIN(TWO*THETA) 10 SUMZ=SUMZ+ZP(J) SUMX=SUMX/DFLOAT(NP) SUMY=SUMY/DFLOAT(NP) SUMZ=SUMZ/DFLOAT(NP) DO 20 J=1,NP XP(J)=XP(J)-SUMX YP(J)=YP(J)-SUMY 20 ZP(J)=ZP(J)-SUMZ C C Set the linear constraints. C M=4*NP DO 30 K=1,M B(K)=ONE DO 30 I=1,N 30 A(I,K)=ZERO DO 40 J=1,NP DO 40 I=1,4 K=4*J+I-4 IW=3*I A(IW-2,K)=XP(J) A(IW-1,K)=YP(J) 40 A(IW,K)=ZP(J) C C Set the initial vector of variables. The JCASE=1,6 loop gives six C different choices of NPT when LINCOA is called. C XS=ZERO YS=ZERO ZS=ZERO SS=ZERO DO 50 J=1,NP XS=DMIN1(XS,XP(J)) YS=DMIN1(YS,YP(J)) ZS=DMIN1(ZS,ZP(J)) 50 SS=DMAX1(SS,XP(J)+YP(J)+ZP(J)) FMAX=(SS-XS-YS-ZS)**3/6.0D0 DO 80 JCASE=1,6 DO 60 I=2,8 60 X(I)=ZERO X(1)=ONE/XS X(5)=ONE/YS X(9)=ONE/ZS X(10)=ONE/SS X(11)=ONE/SS X(12)=ONE/SS C C Call of LINCOA, which provides the printing given at the end of this C note. C NPT=5*JCASE+10 RHOBEG=1.0D0 RHOEND=1.0D-6 IPRINT=1 MAXFUN=10000 PRINT 70, NPT,RHOEND 70 FORMAT (//4X,'Output from LINCOA with NPT =',I4, 1 ' and RHOEND =',1PD12.4) CALL LINCOA (N,NPT,M,A,IA,B,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W) 80 CONTINUE STOP END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% calfun.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE CALFUN (N,X,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON FMAX DIMENSION X(*) ZERO=0.0D0 F=FMAX V12=X(1)*X(5)-X(4)*X(2) V13=X(1)*X(8)-X(7)*X(2) V14=X(1)*X(11)-X(10)*X(2) V23=X(4)*X(8)-X(7)*X(5) V24=X(4)*X(11)-X(10)*X(5) V34=X(7)*X(11)-X(10)*X(8) DEL1=V23*X(12)-V24*X(9)+V34*X(6) IF (DEL1 .LE. ZERO) GOTO 10 DEL2=-V34*X(3)-V13*X(12)+V14*X(9) IF (DEL2 .LE. ZERO) GOTO 10 DEL3=-V14*X(6)+V24*X(3)+V12*X(12) IF (DEL3 .LE. ZERO) GOTO 10 DEL4=-V12*X(9)+V13*X(6)-V23*X(3) IF (DEL4 .LE. ZERO) GOTO 10 TEMP=(DEL1+DEL2+DEL3+DEL4)**3/(DEL1*DEL2*DEL3*DEL4) F=DMIN1(TEMP/6.0D0,FMAX) 10 CONTINUE RETURN END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lincoa.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lincob.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% getact.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% prelim.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% qmstep.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% trstep.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% update.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% PtsinTet output %%%%%%%%%%%%%%%%%%%%%%%%%%%% Output from LINCOA with NPT = 15 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 144 Least value of F = 2.761312530862147D+00 The corresponding X is: -1.233791D+00 -1.231348D+00 -8.870616D-01 1.129888D+00 -1.128285D+00 1.084034D+00 8.431310D-01 7.535391D-01 -8.461760D-01 -6.585848D-01 1.556805D+00 5.853412D-01 Output from LINCOA with NPT = 20 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 170 Least value of F = 2.761312522745138D+00 The corresponding X is: -1.233791D+00 -1.231348D+00 -8.870616D-01 1.129888D+00 -1.128285D+00 1.084034D+00 8.431307D-01 7.535403D-01 -8.461757D-01 -6.585848D-01 1.556805D+00 5.853412D-01 Output from LINCOA with NPT = 25 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 175 Least value of F = 2.761312556579145D+00 The corresponding X is: -8.431310D-01 7.535391D-01 8.461760D-01 1.233791D+00 -1.231348D+00 8.870608D-01 -1.129888D+00 -1.128285D+00 -1.084034D+00 6.585848D-01 1.556805D+00 -5.853412D-01 Output from LINCOA with NPT = 30 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 210 Least value of F = 2.761312522213631D+00 The corresponding X is: -8.431308D-01 7.535397D-01 8.461758D-01 1.233791D+00 -1.231348D+00 8.870616D-01 -1.129888D+00 -1.128285D+00 -1.084034D+00 6.585848D-01 1.556805D+00 -5.853412D-01 Output from LINCOA with NPT = 35 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 131 Least value of F = 2.751575570976407D+00 The corresponding X is: -7.788130D-01 1.062498D+00 7.734507D-01 1.210293D+00 -1.212501D+00 9.512774D-01 -1.157576D+00 -1.159364D+00 -1.048821D+00 7.068259D-01 1.408297D+00 -6.554271D-01 Output from LINCOA with NPT = 40 and RHOEND = 1.0000D-06 At the return from LINCOA Number of function values = 182 Least value of F = 2.761312522228224D+00 The corresponding X is: -8.431308D-01 7.535397D-01 8.461758D-01 1.233791D+00 -1.231348D+00 8.870616D-01 -1.129888D+00 -1.128285D+00 -1.084034D+00 6.585848D-01 1.556805D+00 -5.853412D-01