C
C Test problem for BOBYQA, the objective function being the sum of
C the reciprocals of all pairwise distances between the points P_I,
C I=1,2,...,M in two dimensions, where M=N/2 and where the components
C of P_I are X(2*I-1) and X(2*I). Thus each vector X of N variables
C defines the M points P_I. The initial X gives equally spaced points
C on a circle. Four different choices of the pairs (N,NPT) are tried,
C namely (10,16), (10,21), (20,26) and (20,41). Convergence to a local
C minimum that is not global occurs in both the N=10 cases. The details
C of the results are highly sensitive to computer rounding errors. The
C choice IPRINT=2 provides the current X and optimal F so far whenever
C RHO is reduced. The bound constraints of the problem require every
C component of X to be in the interval [-1,1].
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(100),XL(100),XU(100),W(500000)
TWOPI=8.0D0*DATAN(1.0D0)
BDL=-1.0D0
BDU=1.0D0
IPRINT=2
MAXFUN=500000
RHOBEG=1.0D-1
RHOEND=1.0D-6
M=5
10 N=2*M
DO 20 I=1,N
XL(I)=BDL
20 XU(I)=BDU
DO 50 JCASE=1,2
NPT=N+6
IF (JCASE .EQ. 2) NPT=2*N+1
PRINT 30, M,N,NPT
30 FORMAT (//5X,'2D output with M =',I4,', N =',I4,
1 ' and NPT =',I4)
DO 40 J=1,M
TEMP=DFLOAT(J)*TWOPI/DFLOAT(M)
X(2*J-1)=DCOS(TEMP)
40 X(2*J)=DSIN(TEMP)
CALL BOBYQA (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN,W)
50 CONTINUE
M=M+M
IF (M .LE. 10) GOTO 10
STOP
END
SUBROUTINE CALFUN (N,X,F)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(*)
F=0.0D0
DO 10 I=4,N,2
DO 10 J=2,I-2,2
TEMP=(X(I-1)-X(J-1))**2+(X(I)-X(J))**2
TEMP=DMAX1(TEMP,1.0D-6)
10 F=F+1.0D0/DSQRT(TEMP)
RETURN
END
SUBROUTINE BOBYQA (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,
1 MAXFUN,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(*),XL(*),XU(*),W(*)
C
C This subroutine seeks the least value of a function of many variables,
C by applying a trust region method that forms quadratic models by
C interpolation. There is usually some freedom in the interpolation
C conditions, which is taken up by minimizing the Frobenius norm of
C the change to the second derivative of the model, beginning with the
C zero matrix. The values of the variables are constrained by upper and
C lower bounds. 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 is the number of interpolation conditions. Its value must be in
C the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not
C recommended.
C Initial values of the variables must be set in X(1),X(2),...,X(N). They
C will be changed to the values that give the least calculated F.
C For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper
C bounds, respectively, on X(I). The construction of quadratic models
C requires XL(I) to be strictly less than XU(I) for each I. Further,
C the contribution to a model from changes to the I-th variable is
C damaged severely by rounding errors if XU(I)-XL(I) is too small.
C RHOBEG and RHOEND must be set to the initial and final values of a trust
C region radius, so both must be positive with RHOEND no greater than
C RHOBEG. Typically, RHOBEG should be about one tenth of the greatest
C expected change to a variable, while RHOEND should indicate the
C accuracy that is required in the final values of the variables. An
C error return occurs if any of the differences XU(I)-XL(I), I=1,...,N,
C is less than 2*RHOBEG.
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, each new
C value of RHO is printed, with the best vector of variables so far and
C the corresponding value of the objective function. Further, each new
C 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 The array W will be used for working space. Its length must be at least
C (NPT+5)*(NPT+N)+3*N*(N+5)/2.
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 current values of the
C variables X(1),X(2),...,X(N), which are generated automatically in a
C way that satisfies the bounds given in XL and XU.
C
C Return if the value of NPT is unacceptable.
C
NP=N+1
IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN
PRINT 10
10 FORMAT (/4X,'Return from BOBYQA because NPT is not in',
1 ' the required interval')
GO TO 40
END IF
C
C Partition the working space array, so that different parts of it can
C be treated separately during the calculation of BOBYQB. The partition
C requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the
C space that is taken by the last array in the argument list of BOBYQB.
C
NDIM=NPT+N
IXB=1
IXP=IXB+N
IFV=IXP+N*NPT
IXO=IFV+NPT
IGO=IXO+N
IHQ=IGO+N
IPQ=IHQ+(N*NP)/2
IBMAT=IPQ+NPT
IZMAT=IBMAT+NDIM*N
ISL=IZMAT+NPT*(NPT-NP)
ISU=ISL+N
IXN=ISU+N
IXA=IXN+N
ID=IXA+N
IVL=ID+N
IW=IVL+NDIM
C
C Return if there is insufficient space between the bounds. Modify the
C initial X if necessary in order to avoid conflicts between the bounds
C and the construction of the first quadratic model. The lower and upper
C bounds on moves from the updated X are set now, in the ISL and ISU
C partitions of W, in order to provide useful and exact information about
C components of X that become within distance RHOBEG from their bounds.
C
ZERO=0.0D0
DO 30 J=1,N
TEMP=XU(J)-XL(J)
IF (TEMP .LT. RHOBEG+RHOBEG) THEN
PRINT 20
20 FORMAT (/4X,'Return from BOBYQA because one of the',
1 ' differences XU(I)-XL(I)'/6X,' is less than 2*RHOBEG.')
GO TO 40
END IF
JSL=ISL+J-1
JSU=JSL+N
W(JSL)=XL(J)-X(J)
W(JSU)=XU(J)-X(J)
IF (W(JSL) .GE. -RHOBEG) THEN
IF (W(JSL) .GE. ZERO) THEN
X(J)=XL(J)
W(JSL)=ZERO
W(JSU)=TEMP
ELSE
X(J)=XL(J)+RHOBEG
W(JSL)=-RHOBEG
W(JSU)=DMAX1(XU(J)-X(J),RHOBEG)
END IF
ELSE IF (W(JSU) .LE. RHOBEG) THEN
IF (W(JSU) .LE. ZERO) THEN
X(J)=XU(J)
W(JSL)=-TEMP
W(JSU)=ZERO
ELSE
X(J)=XU(J)-RHOBEG
W(JSL)=DMIN1(XL(J)-X(J),-RHOBEG)
W(JSU)=RHOBEG
END IF
END IF
30 CONTINUE
C
C Make the call of BOBYQB.
C
CALL BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB),
1 W(IXP),W(IFV),W(IXO),W(IGO),W(IHQ),W(IPQ),W(IBMAT),W(IZMAT),
2 NDIM,W(ISL),W(ISU),W(IXN),W(IXA),W(ID),W(IVL),W(IW))
40 RETURN
END
SUBROUTINE BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,
1 MAXFUN,XBASE,XPT,FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,
2 SL,SU,XNEW,XALT,D,VLAG,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(*),XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),
1 XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),
2 SL(*),SU(*),XNEW(*),XALT(*),D(*),VLAG(*),W(*)
C
C The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
C are identical to the corresponding arguments in SUBROUTINE BOBYQA.
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 is a two-dimensional array that holds the coordinates of the
C interpolation points relative to XBASE.
C FVAL holds the values of F at the interpolation points.
C XOPT is set to the displacement from XBASE of the trust region centre.
C GOPT holds the gradient of the quadratic model at 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 H.
C ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
C this factorization being ZMAT times ZMAT^T, which provides both the
C correct rank and positive semi-definiteness.
C NDIM is the first dimension of BMAT and has the value NPT+N.
C SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
C All the components of every XOPT are going to satisfy the bounds
C SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
C XOPT is on a constraint boundary.
C XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the
C vector of variables for the next call of CALFUN. XNEW also satisfies
C the SL and SU constraints in the way that has just been mentioned.
C XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
C in order to increase the denominator in the updating of UPDATE.
C D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
C VLAG contains the values of the Lagrange functions at a new point X.
C They are part of a product that requires VLAG to be of length NDIM.
C W is a one-dimensional array that is used for working space. Its length
C must be at least 3*NDIM = 3*(NPT+N).
C
C Set some constants.
C
HALF=0.5D0
ONE=1.0D0
TEN=10.0D0
TENTH=0.1D0
TWO=2.0D0
ZERO=0.0D0
NP=N+1
NPTM=NPT-NP
NH=(N*NP)/2
C
C The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
C BMAT and ZMAT for the first iteration, with the corresponding values of
C of NF and KOPT, which are the number of calls of CALFUN so far and the
C index of the interpolation point at the trust region centre. Then the
C initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
C less than NPT. GOPT will be updated if KOPT is different from KBASE.
C
CALL PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,XPT,
1 FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT)
XOPTSQ=ZERO
DO 10 I=1,N
XOPT(I)=XPT(KOPT,I)
10 XOPTSQ=XOPTSQ+XOPT(I)**2
FSAVE=FVAL(1)
IF (NF .LT. NPT) THEN
IF (IPRINT .GT. 0) PRINT 390
GOTO 720
END IF
KBASE=1
C
C Complete the settings that are required for the iterative procedure.
C
RHO=RHOBEG
DELTA=RHO
NRESC=NF
NTRITS=0
DIFFA=ZERO
DIFFB=ZERO
ITEST=0
NFSAV=NF
C
C Update GOPT if necessary before the first iteration and after each
C call of RESCUE that makes a call of CALFUN.
C
20 IF (KOPT .NE. KBASE) THEN
IH=0
DO 30 J=1,N
DO 30 I=1,J
IH=IH+1
IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*XOPT(I)
30 GOPT(I)=GOPT(I)+HQ(IH)*XOPT(J)
IF (NF .GT. NPT) THEN
DO 50 K=1,NPT
TEMP=ZERO
DO 40 J=1,N
40 TEMP=TEMP+XPT(K,J)*XOPT(J)
TEMP=PQ(K)*TEMP
DO 50 I=1,N
50 GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
END IF
END IF
C
C Generate the next point in the trust region that provides a small value
C of the quadratic model subject to the constraints on the variables.
C The integer NTRITS is set to the number "trust region" iterations that
C have occurred since the last "alternative" iteration. If the length
C of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
C label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
C
60 CALL TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,XNEW,D,
1 W,W(NP),W(NP+N),W(NP+2*N),W(NP+3*N),DSQ,CRVMIN)
DNORM=DMIN1(DELTA,DSQRT(DSQ))
IF (DNORM .LT. HALF*RHO) THEN
NTRITS=-1
DISTSQ=(TEN*RHO)**2
IF (NF .LE. NFSAV+2) GOTO 650
C
C The following choice between labels 650 and 680 depends on whether or
C not our work with the current RHO seems to be complete. Either RHO is
C decreased or termination occurs if the errors in the quadratic model at
C the last three interpolation points compare favourably with predictions
C of likely improvements to the model within distance HALF*RHO of XOPT.
C
ERRBIG=DMAX1(DIFFA,DIFFB,DIFFC)
FRHOSQ=0.125D0*RHO*RHO
IF (CRVMIN .GT. ZERO .AND. ERRBIG .GT. FRHOSQ*CRVMIN)
1 GOTO 650
BDTOL=ERRBIG/RHO
DO 80 J=1,N
BDTEST=BDTOL
IF (XNEW(J) .EQ. SL(J)) BDTEST=W(J)
IF (XNEW(J) .EQ. SU(J)) BDTEST=-W(J)
IF (BDTEST .LT. BDTOL) THEN
CURV=HQ((J+J*J)/2)
DO 70 K=1,NPT
70 CURV=CURV+PQ(K)*XPT(K,J)**2
BDTEST=BDTEST+HALF*CURV*RHO
IF (BDTEST .LT. BDTOL) GOTO 650
END IF
80 CONTINUE
GOTO 680
END IF
NTRITS=NTRITS+1
C
C Severe cancellation is likely to occur if XOPT is too far from XBASE.
C If the following test holds, then XBASE is shifted so that XOPT becomes
C zero. The appropriate changes are made to BMAT and to the second
C derivatives of the current model, beginning with the changes to BMAT
C that do not depend on ZMAT. VLAG is used temporarily for working space.
C
90 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN
FRACSQ=0.25D0*XOPTSQ
SUMPQ=ZERO
DO 110 K=1,NPT
SUMPQ=SUMPQ+PQ(K)
SUM=-HALF*XOPTSQ
DO 100 I=1,N
100 SUM=SUM+XPT(K,I)*XOPT(I)
W(NPT+K)=SUM
TEMP=FRACSQ-HALF*SUM
DO 110 I=1,N
W(I)=BMAT(K,I)
VLAG(I)=SUM*XPT(K,I)+TEMP*XOPT(I)
IP=NPT+I
DO 110 J=1,I
110 BMAT(IP,J)=BMAT(IP,J)+W(I)*VLAG(J)+VLAG(I)*W(J)
C
C Then the revisions of BMAT that depend on ZMAT are calculated.
C
DO 150 JJ=1,NPTM
SUMZ=ZERO
SUMW=ZERO
DO 120 K=1,NPT
SUMZ=SUMZ+ZMAT(K,JJ)
VLAG(K)=W(NPT+K)*ZMAT(K,JJ)
120 SUMW=SUMW+VLAG(K)
DO 140 J=1,N
SUM=(FRACSQ*SUMZ-HALF*SUMW)*XOPT(J)
DO 130 K=1,NPT
130 SUM=SUM+VLAG(K)*XPT(K,J)
W(J)=SUM
DO 140 K=1,NPT
140 BMAT(K,J)=BMAT(K,J)+SUM*ZMAT(K,JJ)
DO 150 I=1,N
IP=I+NPT
TEMP=W(I)
DO 150 J=1,I
150 BMAT(IP,J)=BMAT(IP,J)+TEMP*W(J)
C
C The following instructions complete the shift, including the changes
C to the second derivative parameters of the quadratic model.
C
IH=0
DO 170 J=1,N
W(J)=-HALF*SUMPQ*XOPT(J)
DO 160 K=1,NPT
W(J)=W(J)+PQ(K)*XPT(K,J)
160 XPT(K,J)=XPT(K,J)-XOPT(J)
DO 170 I=1,J
IH=IH+1
HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
170 BMAT(NPT+I,J)=BMAT(NPT+J,I)
DO 180 I=1,N
XBASE(I)=XBASE(I)+XOPT(I)
XNEW(I)=XNEW(I)-XOPT(I)
SL(I)=SL(I)-XOPT(I)
SU(I)=SU(I)-XOPT(I)
180 XOPT(I)=ZERO
XOPTSQ=ZERO
END IF
IF (NTRITS .EQ. 0) GOTO 210
GOTO 230
C
C XBASE is also moved to XOPT by a call of RESCUE. This calculation is
C more expensive than the previous shift, because new matrices BMAT and
C ZMAT are generated from scratch, which may include the replacement of
C interpolation points whose positions seem to be causing near linear
C dependence in the interpolation conditions. Therefore RESCUE is called
C only if rounding errors have reduced by at least a factor of two the
C denominator of the formula for updating the H matrix. It provides a
C useful safeguard, but is not invoked in most applications of BOBYQA.
C
190 NFSAV=NF
KBASE=KOPT
CALL RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,FVAL,
1 XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,KOPT,
2 VLAG,W,W(N+NP),W(NDIM+NP))
C
C XOPT is updated now in case the branch below to label 720 is taken.
C Any updating of GOPT occurs after the branch below to label 20, which
C leads to a trust region iteration as does the branch to label 60.
C
XOPTSQ=ZERO
IF (KOPT .NE. KBASE) THEN
DO 200 I=1,N
XOPT(I)=XPT(KOPT,I)
200 XOPTSQ=XOPTSQ+XOPT(I)**2
END IF
IF (NF .LT. 0) THEN
NF=MAXFUN
IF (IPRINT .GT. 0) PRINT 390
GOTO 720
END IF
NRESC=NF
IF (NFSAV .LT. NF) THEN
NFSAV=NF
GOTO 20
END IF
IF (NTRITS .GT. 0) GOTO 60
C
C Pick two alternative vectors of variables, relative to XBASE, that
C are suitable as new positions of the KNEW-th interpolation point.
C Firstly, XNEW is set to the point on a line through XOPT and another
C interpolation point that minimizes the predicted value of the next
C denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
C and SU bounds. Secondly, XALT is set to the best feasible point on
C a constrained version of the Cauchy step of the KNEW-th Lagrange
C function, the corresponding value of the square of this function
C being returned in CAUCHY. The choice between these alternatives is
C going to be made when the denominator is calculated.
C
210 CALL ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
1 KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,W,W(NP),W(NDIM+1))
DO 220 I=1,N
220 D(I)=XNEW(I)-XOPT(I)
C
C Calculate VLAG and BETA for the current choice of D. The scalar
C product of D with XPT(K,.) is going to be held in W(NPT+K) for
C use when VQUAD is calculated.
C
230 DO 250 K=1,NPT
SUMA=ZERO
SUMB=ZERO
SUM=ZERO
DO 240 J=1,N
SUMA=SUMA+XPT(K,J)*D(J)
SUMB=SUMB+XPT(K,J)*XOPT(J)
240 SUM=SUM+BMAT(K,J)*D(J)
W(K)=SUMA*(HALF*SUMA+SUMB)
VLAG(K)=SUM
250 W(NPT+K)=SUMA
BETA=ZERO
DO 270 JJ=1,NPTM
SUM=ZERO
DO 260 K=1,NPT
260 SUM=SUM+ZMAT(K,JJ)*W(K)
BETA=BETA-SUM*SUM
DO 270 K=1,NPT
270 VLAG(K)=VLAG(K)+SUM*ZMAT(K,JJ)
DSQ=ZERO
BSUM=ZERO
DX=ZERO
DO 300 J=1,N
DSQ=DSQ+D(J)**2
SUM=ZERO
DO 280 K=1,NPT
280 SUM=SUM+W(K)*BMAT(K,J)
BSUM=BSUM+SUM*D(J)
JP=NPT+J
DO 290 I=1,N
290 SUM=SUM+BMAT(JP,I)*D(I)
VLAG(JP)=SUM
BSUM=BSUM+SUM*D(J)
300 DX=DX+D(J)*XOPT(J)
BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
VLAG(KOPT)=VLAG(KOPT)+ONE
C
C If NTRITS is zero, the denominator may be increased by replacing
C the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
C rounding errors have damaged the chosen denominator.
C
IF (NTRITS .EQ. 0) THEN
DENOM=VLAG(KNEW)**2+ALPHA*BETA
IF (DENOM .LT. CAUCHY .AND. CAUCHY .GT. ZERO) THEN
DO 310 I=1,N
XNEW(I)=XALT(I)
310 D(I)=XNEW(I)-XOPT(I)
CAUCHY=ZERO
GO TO 230
END IF
IF (DENOM .LE. HALF*VLAG(KNEW)**2) THEN
IF (NF .GT. NRESC) GOTO 190
IF (IPRINT .GT. 0) PRINT 320
320 FORMAT (/5X,'Return from BOBYQA because of much',
1 ' cancellation in a denominator.')
GOTO 720
END IF
C
C Alternatively, if NTRITS is positive, then set KNEW to the index of
C the next interpolation point to be deleted to make room for a trust
C region step. Again RESCUE may be called if rounding errors have damaged
C the chosen denominator, which is the reason for attempting to select
C KNEW before calculating the next value of the objective function.
C
ELSE
DELSQ=DELTA*DELTA
SCADEN=ZERO
BIGLSQ=ZERO
KNEW=0
DO 350 K=1,NPT
IF (K .EQ. KOPT) GOTO 350
HDIAG=ZERO
DO 330 JJ=1,NPTM
330 HDIAG=HDIAG+ZMAT(K,JJ)**2
DEN=BETA*HDIAG+VLAG(K)**2
DISTSQ=ZERO
DO 340 J=1,N
340 DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
IF (TEMP*DEN .GT. SCADEN) THEN
SCADEN=TEMP*DEN
KNEW=K
DENOM=DEN
END IF
BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
350 CONTINUE
IF (SCADEN .LE. HALF*BIGLSQ) THEN
IF (NF .GT. NRESC) GOTO 190
IF (IPRINT .GT. 0) PRINT 320
GOTO 720
END IF
END IF
C
C Put the variables for the next calculation of the objective function
C in XNEW, with any adjustments for the bounds.
C
C
C Calculate the value of the objective function at XBASE+XNEW, unless
C the limit on the number of calculations of F has been reached.
C
360 DO 380 I=1,N
X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XNEW(I)),XU(I))
IF (XNEW(I) .EQ. SL(I)) X(I)=XL(I)
IF (XNEW(I) .EQ. SU(I)) X(I)=XU(I)
380 CONTINUE
IF (NF .GE. MAXFUN) THEN
IF (IPRINT .GT. 0) PRINT 390
390 FORMAT (/4X,'Return from BOBYQA because CALFUN has been',
1 ' called MAXFUN times.')
GOTO 720
END IF
NF=NF+1
CALL CALFUN (N,X,F)
IF (IPRINT .EQ. 3) THEN
PRINT 400, NF,F,(X(I),I=1,N)
400 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
1 ' The corresponding X is:'/(2X,5D15.6))
END IF
IF (NTRITS .EQ. -1) THEN
FSAVE=F
GOTO 720
END IF
C
C Use the quadratic model to predict the change in F due to the step D,
C and set DIFF to the error of this prediction.
C
FOPT=FVAL(KOPT)
VQUAD=ZERO
IH=0
DO 410 J=1,N
VQUAD=VQUAD+D(J)*GOPT(J)
DO 410 I=1,J
IH=IH+1
TEMP=D(I)*D(J)
IF (I .EQ. J) TEMP=HALF*TEMP
410 VQUAD=VQUAD+HQ(IH)*TEMP
DO 420 K=1,NPT
420 VQUAD=VQUAD+HALF*PQ(K)*W(NPT+K)**2
DIFF=F-FOPT-VQUAD
DIFFC=DIFFB
DIFFB=DIFFA
DIFFA=DABS(DIFF)
IF (DNORM .GT. RHO) NFSAV=NF
C
C Pick the next value of DELTA after a trust region step.
C
IF (NTRITS .GT. 0) THEN
IF (VQUAD .GE. ZERO) THEN
IF (IPRINT .GT. 0) PRINT 430
430 FORMAT (/4X,'Return from BOBYQA because a trust',
1 ' region step has failed to reduce Q.')
GOTO 720
END IF
RATIO=(F-FOPT)/VQUAD
IF (RATIO .LE. TENTH) THEN
DELTA=DMIN1(HALF*DELTA,DNORM)
ELSE IF (RATIO. LE. 0.7D0) THEN
DELTA=DMAX1(HALF*DELTA,DNORM)
ELSE
DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
END IF
IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
C
C Recalculate KNEW and DENOM if the new F is less than FOPT.
C
IF (F .LT. FOPT) THEN
KSAV=KNEW
DENSAV=DENOM
DELSQ=DELTA*DELTA
SCADEN=ZERO
BIGLSQ=ZERO
KNEW=0
DO 460 K=1,NPT
HDIAG=ZERO
DO 440 JJ=1,NPTM
440 HDIAG=HDIAG+ZMAT(K,JJ)**2
DEN=BETA*HDIAG+VLAG(K)**2
DISTSQ=ZERO
DO 450 J=1,N
450 DISTSQ=DISTSQ+(XPT(K,J)-XNEW(J))**2
TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
IF (TEMP*DEN .GT. SCADEN) THEN
SCADEN=TEMP*DEN
KNEW=K
DENOM=DEN
END IF
460 BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
IF (SCADEN .LE. HALF*BIGLSQ) THEN
KNEW=KSAV
DENOM=DENSAV
END IF
END IF
END IF
C
C Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
C moved. Also update the second derivative terms of the model.
C
CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
IH=0
PQOLD=PQ(KNEW)
PQ(KNEW)=ZERO
DO 470 I=1,N
TEMP=PQOLD*XPT(KNEW,I)
DO 470 J=1,I
IH=IH+1
470 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
DO 480 JJ=1,NPTM
TEMP=DIFF*ZMAT(KNEW,JJ)
DO 480 K=1,NPT
480 PQ(K)=PQ(K)+TEMP*ZMAT(K,JJ)
C
C Include the new interpolation point, and make the changes to GOPT at
C the old XOPT that are caused by the updating of the quadratic model.
C
FVAL(KNEW)=F
DO 490 I=1,N
XPT(KNEW,I)=XNEW(I)
490 W(I)=BMAT(KNEW,I)
DO 520 K=1,NPT
SUMA=ZERO
DO 500 JJ=1,NPTM
500 SUMA=SUMA+ZMAT(KNEW,JJ)*ZMAT(K,JJ)
SUMB=ZERO
DO 510 J=1,N
510 SUMB=SUMB+XPT(K,J)*XOPT(J)
TEMP=SUMA*SUMB
DO 520 I=1,N
520 W(I)=W(I)+TEMP*XPT(K,I)
DO 530 I=1,N
530 GOPT(I)=GOPT(I)+DIFF*W(I)
C
C Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
C
IF (F .LT. FOPT) THEN
KOPT=KNEW
XOPTSQ=ZERO
IH=0
DO 540 J=1,N
XOPT(J)=XNEW(J)
XOPTSQ=XOPTSQ+XOPT(J)**2
DO 540 I=1,J
IH=IH+1
IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*D(I)
540 GOPT(I)=GOPT(I)+HQ(IH)*D(J)
DO 560 K=1,NPT
TEMP=ZERO
DO 550 J=1,N
550 TEMP=TEMP+XPT(K,J)*D(J)
TEMP=PQ(K)*TEMP
DO 560 I=1,N
560 GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
END IF
C
C Calculate the parameters of the least Frobenius norm interpolant to
C the current data, the gradient of this interpolant at XOPT being put
C into VLAG(NPT+I), I=1,2,...,N.
C
IF (NTRITS .GT. 0) THEN
DO 570 K=1,NPT
VLAG(K)=FVAL(K)-FVAL(KOPT)
570 W(K)=ZERO
DO 590 J=1,NPTM
SUM=ZERO
DO 580 K=1,NPT
580 SUM=SUM+ZMAT(K,J)*VLAG(K)
DO 590 K=1,NPT
590 W(K)=W(K)+SUM*ZMAT(K,J)
DO 610 K=1,NPT
SUM=ZERO
DO 600 J=1,N
600 SUM=SUM+XPT(K,J)*XOPT(J)
W(K+NPT)=W(K)
610 W(K)=SUM*W(K)
GQSQ=ZERO
GISQ=ZERO
DO 630 I=1,N
SUM=ZERO
DO 620 K=1,NPT
620 SUM=SUM+BMAT(K,I)*VLAG(K)+XPT(K,I)*W(K)
IF (XOPT(I) .EQ. SL(I)) THEN
GQSQ=GQSQ+DMIN1(ZERO,GOPT(I))**2
GISQ=GISQ+DMIN1(ZERO,SUM)**2
ELSE IF (XOPT(I) .EQ. SU(I)) THEN
GQSQ=GQSQ+DMAX1(ZERO,GOPT(I))**2
GISQ=GISQ+DMAX1(ZERO,SUM)**2
ELSE
GQSQ=GQSQ+GOPT(I)**2
GISQ=GISQ+SUM*SUM
END IF
630 VLAG(NPT+I)=SUM
C
C Test whether to replace the new quadratic model by the least Frobenius
C norm interpolant, making the replacement if the test is satisfied.
C
ITEST=ITEST+1
IF (GQSQ .LT. TEN*GISQ) ITEST=0
IF (ITEST .GE. 3) THEN
DO 640 I=1,MAX0(NPT,NH)
IF (I .LE. N) GOPT(I)=VLAG(NPT+I)
IF (I .LE. NPT) PQ(I)=W(NPT+I)
IF (I .LE. NH) HQ(I)=ZERO
ITEST=0
640 CONTINUE
END IF
END IF
C
C If a trust region step has provided a sufficient decrease in F, then
C branch for another trust region calculation. The case NTRITS=0 occurs
C when the new interpolation point was reached by an alternative step.
C
IF (NTRITS .EQ. 0) GOTO 60
IF (F .LE. FOPT+TENTH*VQUAD) GOTO 60
C
C Alternatively, find out if the interpolation points are close enough
C to the best point so far.
C
DISTSQ=DMAX1((TWO*DELTA)**2,(TEN*RHO)**2)
650 KNEW=0
DO 670 K=1,NPT
SUM=ZERO
DO 660 J=1,N
660 SUM=SUM+(XPT(K,J)-XOPT(J))**2
IF (SUM .GT. DISTSQ) THEN
KNEW=K
DISTSQ=SUM
END IF
670 CONTINUE
C
C If KNEW is positive, then ALTMOV finds alternative new positions for
C the KNEW-th interpolation point within distance ADELT of XOPT. It is
C reached via label 90. Otherwise, there is a branch to label 60 for
C another trust region iteration, unless the calculations with the
C current RHO are complete.
C
IF (KNEW .GT. 0) THEN
DIST=DSQRT(DISTSQ)
IF (NTRITS .EQ. -1) THEN
DELTA=DMIN1(TENTH*DELTA,HALF*DIST)
IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
END IF
NTRITS=0
ADELT=DMAX1(DMIN1(TENTH*DIST,DELTA),RHO)
DSQ=ADELT*ADELT
GOTO 90
END IF
IF (NTRITS .EQ. -1) GOTO 680
IF (RATIO .GT. ZERO) GOTO 60
IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 60
C
C The calculations with the current value of RHO are complete. Pick the
C next values of RHO and DELTA.
C
680 IF (RHO .GT. RHOEND) THEN
DELTA=HALF*RHO
RATIO=RHO/RHOEND
IF (RATIO .LE. 16.0D0) THEN
RHO=RHOEND
ELSE IF (RATIO .LE. 250.0D0) THEN
RHO=DSQRT(RATIO)*RHOEND
ELSE
RHO=TENTH*RHO
END IF
DELTA=DMAX1(DELTA,RHO)
IF (IPRINT .GE. 2) THEN
IF (IPRINT .GE. 3) PRINT 690
690 FORMAT (5X)
PRINT 700, RHO,NF
700 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of',
1 ' function values =',I6)
PRINT 710, FVAL(KOPT),(XBASE(I)+XOPT(I),I=1,N)
710 FORMAT (4X,'Least value of F =',1PD23.15,9X,
1 'The corresponding X is:'/(2X,5D15.6))
END IF
NTRITS=0
NFSAV=NF
GOTO 60
END IF
C
C Return from the calculation, after another Newton-Raphson step, if
C it is too short to have been tried before.
C
IF (NTRITS .EQ. -1) GOTO 360
720 IF (FVAL(KOPT) .LE. FSAVE) THEN
DO 730 I=1,N
X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XOPT(I)),XU(I))
IF (XOPT(I) .EQ. SL(I)) X(I)=XL(I)
IF (XOPT(I) .EQ. SU(I)) X(I)=XU(I)
730 CONTINUE
F=FVAL(KOPT)
END IF
IF (IPRINT .GE. 1) THEN
PRINT 740, NF
740 FORMAT (/4X,'At the return from BOBYQA',5X,
1 'Number of function values =',I6)
PRINT 710, F,(X(I),I=1,N)
END IF
RETURN
END
SUBROUTINE ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
1 KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,GLAG,HCOL,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XPT(NPT,*),XOPT(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),
1 SU(*),XNEW(*),XALT(*),GLAG(*),HCOL(*),W(*)
C
C The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have
C the same meanings as the corresponding arguments of BOBYQB.
C KOPT is the index of the optimal interpolation point.
C KNEW is the index of the interpolation point that is going to be moved.
C ADELT is the current trust region bound.
C XNEW will be set to a suitable new position for the interpolation point
C XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region
C bounds and it should provide a large denominator in the next call of
C UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the
C straight lines through XOPT and another interpolation point.
C XALT also provides a large value of the modulus of the KNEW-th Lagrange
C function subject to the constraints that have been mentioned, its main
C difference from XNEW being that XALT-XOPT is a constrained version of
C the Cauchy step within the trust region. An exception is that XALT is
C not calculated if all components of GLAG (see below) are zero.
C ALPHA will be set to the KNEW-th diagonal element of the H matrix.
C CAUCHY will be set to the square of the KNEW-th Lagrange function at
C the step XALT-XOPT from XOPT for the vector XALT that is returned,
C except that CAUCHY is set to zero if XALT is not calculated.
C GLAG is a working space vector of length N for the gradient of the
C KNEW-th Lagrange function at XOPT.
C HCOL is a working space vector of length NPT for the second derivative
C coefficients of the KNEW-th Lagrange function.
C W is a working space vector of length 2N that is going to hold the
C constrained Cauchy step from XOPT of the Lagrange function, followed
C by the downhill version of XALT when the uphill step is calculated.
C
C Set the first NPT components of W to the leading elements of the
C KNEW-th column of the H matrix.
C
HALF=0.5D0
ONE=1.0D0
ZERO=0.0D0
CONST=ONE+DSQRT(2.0D0)
DO 10 K=1,NPT
10 HCOL(K)=ZERO
DO 20 J=1,NPT-N-1
TEMP=ZMAT(KNEW,J)
DO 20 K=1,NPT
20 HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J)
ALPHA=HCOL(KNEW)
HA=HALF*ALPHA
C
C Calculate the gradient of the KNEW-th Lagrange function at XOPT.
C
DO 30 I=1,N
30 GLAG(I)=BMAT(KNEW,I)
DO 50 K=1,NPT
TEMP=ZERO
DO 40 J=1,N
40 TEMP=TEMP+XPT(K,J)*XOPT(J)
TEMP=HCOL(K)*TEMP
DO 50 I=1,N
50 GLAG(I)=GLAG(I)+TEMP*XPT(K,I)
C
C Search for a large denominator along the straight lines through XOPT
C and another interpolation point. SLBD and SUBD will be lower and upper
C bounds on the step along each of these lines in turn. PREDSQ will be
C set to the square of the predicted denominator for each line. PRESAV
C will be set to the largest admissible value of PREDSQ that occurs.
C
PRESAV=ZERO
DO 80 K=1,NPT
IF (K .EQ. KOPT) GOTO 80
DDERIV=ZERO
DISTSQ=ZERO
DO 60 I=1,N
TEMP=XPT(K,I)-XOPT(I)
DDERIV=DDERIV+GLAG(I)*TEMP
60 DISTSQ=DISTSQ+TEMP*TEMP
SUBD=ADELT/DSQRT(DISTSQ)
SLBD=-SUBD
ILBD=0
IUBD=0
SUMIN=DMIN1(ONE,SUBD)
C
C Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
C
DO 70 I=1,N
TEMP=XPT(K,I)-XOPT(I)
IF (TEMP .GT. ZERO) THEN
IF (SLBD*TEMP .LT. SL(I)-XOPT(I)) THEN
SLBD=(SL(I)-XOPT(I))/TEMP
ILBD=-I
END IF
IF (SUBD*TEMP .GT. SU(I)-XOPT(I)) THEN
SUBD=DMAX1(SUMIN,(SU(I)-XOPT(I))/TEMP)
IUBD=I
END IF
ELSE IF (TEMP .LT. ZERO) THEN
IF (SLBD*TEMP .GT. SU(I)-XOPT(I)) THEN
SLBD=(SU(I)-XOPT(I))/TEMP
ILBD=I
END IF
IF (SUBD*TEMP .LT. SL(I)-XOPT(I)) THEN
SUBD=DMAX1(SUMIN,(SL(I)-XOPT(I))/TEMP)
IUBD=-I
END IF
END IF
70 CONTINUE
C
C Seek a large modulus of the KNEW-th Lagrange function when the index
C of the other interpolation point on the line through XOPT is KNEW.
C
IF (K .EQ. KNEW) THEN
DIFF=DDERIV-ONE
STEP=SLBD
VLAG=SLBD*(DDERIV-SLBD*DIFF)
ISBD=ILBD
TEMP=SUBD*(DDERIV-SUBD*DIFF)
IF (DABS(TEMP) .GT. DABS(VLAG)) THEN
STEP=SUBD
VLAG=TEMP
ISBD=IUBD
END IF
TEMPD=HALF*DDERIV
TEMPA=TEMPD-DIFF*SLBD
TEMPB=TEMPD-DIFF*SUBD
IF (TEMPA*TEMPB .LT. ZERO) THEN
TEMP=TEMPD*TEMPD/DIFF
IF (DABS(TEMP) .GT. DABS(VLAG)) THEN
STEP=TEMPD/DIFF
VLAG=TEMP
ISBD=0
END IF
END IF
C
C Search along each of the other lines through XOPT and another point.
C
ELSE
STEP=SLBD
VLAG=SLBD*(ONE-SLBD)
ISBD=ILBD
TEMP=SUBD*(ONE-SUBD)
IF (DABS(TEMP) .GT. DABS(VLAG)) THEN
STEP=SUBD
VLAG=TEMP
ISBD=IUBD
END IF
IF (SUBD .GT. HALF) THEN
IF (DABS(VLAG) .LT. 0.25D0) THEN
STEP=HALF
VLAG=0.25D0
ISBD=0
END IF
END IF
VLAG=VLAG*DDERIV
END IF
C
C Calculate PREDSQ for the current line search and maintain PRESAV.
C
TEMP=STEP*(ONE-STEP)*DISTSQ
PREDSQ=VLAG*VLAG*(VLAG*VLAG+HA*TEMP*TEMP)
IF (PREDSQ .GT. PRESAV) THEN
PRESAV=PREDSQ
KSAV=K
STPSAV=STEP
IBDSAV=ISBD
END IF
80 CONTINUE
C
C Construct XNEW in a way that satisfies the bound constraints exactly.
C
DO 90 I=1,N
TEMP=XOPT(I)+STPSAV*(XPT(KSAV,I)-XOPT(I))
90 XNEW(I)=DMAX1(SL(I),DMIN1(SU(I),TEMP))
IF (IBDSAV .LT. 0) XNEW(-IBDSAV)=SL(-IBDSAV)
IF (IBDSAV .GT. 0) XNEW(IBDSAV)=SU(IBDSAV)
C
C Prepare for the iterative method that assembles the constrained Cauchy
C step in W. The sum of squares of the fixed components of W is formed in
C WFIXSQ, and the free components of W are set to BIGSTP.
C
BIGSTP=ADELT+ADELT
IFLAG=0
100 WFIXSQ=ZERO
GGFREE=ZERO
DO 110 I=1,N
W(I)=ZERO
TEMPA=DMIN1(XOPT(I)-SL(I),GLAG(I))
TEMPB=DMAX1(XOPT(I)-SU(I),GLAG(I))
IF (TEMPA .GT. ZERO .OR. TEMPB .LT. ZERO) THEN
W(I)=BIGSTP
GGFREE=GGFREE+GLAG(I)**2
END IF
110 CONTINUE
IF (GGFREE .EQ. ZERO) THEN
CAUCHY=ZERO
GOTO 200
END IF
C
C Investigate whether more components of W can be fixed.
C
120 TEMP=ADELT*ADELT-WFIXSQ
IF (TEMP .GT. ZERO) THEN
WSQSAV=WFIXSQ
STEP=DSQRT(TEMP/GGFREE)
GGFREE=ZERO
DO 130 I=1,N
IF (W(I) .EQ. BIGSTP) THEN
TEMP=XOPT(I)-STEP*GLAG(I)
IF (TEMP .LE. SL(I)) THEN
W(I)=SL(I)-XOPT(I)
WFIXSQ=WFIXSQ+W(I)**2
ELSE IF (TEMP .GE. SU(I)) THEN
W(I)=SU(I)-XOPT(I)
WFIXSQ=WFIXSQ+W(I)**2
ELSE
GGFREE=GGFREE+GLAG(I)**2
END IF
END IF
130 CONTINUE
IF (WFIXSQ .GT. WSQSAV .AND. GGFREE .GT. ZERO) GOTO 120
END IF
C
C Set the remaining free components of W and all components of XALT,
C except that W may be scaled later.
C
GW=ZERO
DO 140 I=1,N
IF (W(I) .EQ. BIGSTP) THEN
W(I)=-STEP*GLAG(I)
XALT(I)=DMAX1(SL(I),DMIN1(SU(I),XOPT(I)+W(I)))
ELSE IF (W(I) .EQ. ZERO) THEN
XALT(I)=XOPT(I)
ELSE IF (GLAG(I) .GT. ZERO) THEN
XALT(I)=SL(I)
ELSE
XALT(I)=SU(I)
END IF
140 GW=GW+GLAG(I)*W(I)
C
C Set CURV to the curvature of the KNEW-th Lagrange function along W.
C Scale W by a factor less than one if that can reduce the modulus of
C the Lagrange function at XOPT+W. Set CAUCHY to the final value of
C the square of this function.
C
CURV=ZERO
DO 160 K=1,NPT
TEMP=ZERO
DO 150 J=1,N
150 TEMP=TEMP+XPT(K,J)*W(J)
160 CURV=CURV+HCOL(K)*TEMP*TEMP
IF (IFLAG .EQ. 1) CURV=-CURV
IF (CURV .GT. -GW .AND. CURV .LT. -CONST*GW) THEN
SCALE=-GW/CURV
DO 170 I=1,N
TEMP=XOPT(I)+SCALE*W(I)
170 XALT(I)=DMAX1(SL(I),DMIN1(SU(I),TEMP))
CAUCHY=(HALF*GW*SCALE)**2
ELSE
CAUCHY=(GW+HALF*CURV)**2
END IF
C
C If IFLAG is zero, then XALT is calculated as before after reversing
C the sign of GLAG. Thus two XALT vectors become available. The one that
C is chosen is the one that gives the larger value of CAUCHY.
C
IF (IFLAG .EQ. 0) THEN
DO 180 I=1,N
GLAG(I)=-GLAG(I)
180 W(N+I)=XALT(I)
CSAVE=CAUCHY
IFLAG=1
GOTO 100
END IF
IF (CSAVE .GT. CAUCHY) THEN
DO 190 I=1,N
190 XALT(I)=W(N+I)
CAUCHY=CSAVE
END IF
200 RETURN
END
SUBROUTINE PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,
1 XPT,FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(*),XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),GOPT(*),
1 HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*)
C
C The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the
C same as the corresponding arguments in SUBROUTINE BOBYQA.
C The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU
C are the same as the corresponding arguments in BOBYQB, the elements
C of SL and SU being set in BOBYQA.
C GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but
C it is set by PRELIM to the gradient of the quadratic model at XBASE.
C If XOPT is nonzero, BOBYQB will change it to its usual value later.
C NF is maintaned as the number of calls of CALFUN so far.
C KOPT will be such that the least calculated value of F so far is at
C the point XPT(KOPT,.)+XBASE in the space of the variables.
C
C SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
C BMAT and ZMAT for the first iteration, and it maintains the values of
C NF and KOPT. The vector X is also changed by PRELIM.
C
C Set some constants.
C
HALF=0.5D0
ONE=1.0D0
TWO=2.0D0
ZERO=0.0D0
RHOSQ=RHOBEG*RHOBEG
RECIP=ONE/RHOSQ
NP=N+1
C
C Set XBASE to the initial vector of variables, and set the initial
C elements of XPT, BMAT, HQ, PQ 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 IH=1,(N*NP)/2
30 HQ(IH)=ZERO
DO 40 K=1,NPT
PQ(K)=ZERO
DO 40 J=1,NPT-NP
40 ZMAT(K,J)=ZERO
C
C Begin the initialization procedure. NF becomes one more than the number
C of function values so far. The coordinates of the displacement of the
C next initial interpolation point from XBASE are set in XPT(NF+1,.).
C
NF=0
50 NFM=NF
NFX=NF-N
NF=NF+1
IF (NFM .LE. 2*N) THEN
IF (NFM .GE. 1 .AND. NFM .LE. N) THEN
STEPA=RHOBEG
IF (SU(NFM) .EQ. ZERO) STEPA=-STEPA
XPT(NF,NFM)=STEPA
ELSE IF (NFM .GT. N) THEN
STEPA=XPT(NF-N,NFX)
STEPB=-RHOBEG
IF (SL(NFX) .EQ. ZERO) STEPB=DMIN1(TWO*RHOBEG,SU(NFX))
IF (SU(NFX) .EQ. ZERO) STEPB=DMAX1(-TWO*RHOBEG,SL(NFX))
XPT(NF,NFX)=STEPB
END IF
ELSE
ITEMP=(NFM-NP)/N
JPT=NFM-ITEMP*N-N
IPT=JPT+ITEMP
IF (IPT .GT. N) THEN
ITEMP=JPT
JPT=IPT-N
IPT=ITEMP
END IF
XPT(NF,IPT)=XPT(IPT+1,IPT)
XPT(NF,JPT)=XPT(JPT+1,JPT)
END IF
C
C Calculate the next value of F. The least function value so far and
C its index are required.
C
DO 60 J=1,N
X(J)=DMIN1(DMAX1(XL(J),XBASE(J)+XPT(NF,J)),XU(J))
IF (XPT(NF,J) .EQ. SL(J)) X(J)=XL(J)
IF (XPT(NF,J) .EQ. SU(J)) X(J)=XU(J)
60 CONTINUE
CALL CALFUN (N,X,F)
IF (IPRINT .EQ. 3) THEN
PRINT 70, NF,F,(X(I),I=1,N)
70 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
1 ' The corresponding X is:'/(2X,5D15.6))
END IF
FVAL(NF)=F
IF (NF .EQ. 1) THEN
FBEG=F
KOPT=1
ELSE IF (F .LT. FVAL(KOPT)) THEN
KOPT=NF
END IF
C
C Set the nonzero initial elements of BMAT and the quadratic model in the
C cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions
C of the NF-th and (NF-N)-th interpolation points may be switched, in
C order that the function value at the first of them contributes to the
C off-diagonal second derivative terms of the initial quadratic model.
C
IF (NF .LE. 2*N+1) THEN
IF (NF .GE. 2 .AND. NF .LE. N+1) THEN
GOPT(NFM)=(F-FBEG)/STEPA
IF (NPT .LT. NF+N) THEN
BMAT(1,NFM)=-ONE/STEPA
BMAT(NF,NFM)=ONE/STEPA
BMAT(NPT+NFM,NFM)=-HALF*RHOSQ
END IF
ELSE IF (NF .GE. N+2) THEN
IH=(NFX*(NFX+1))/2
TEMP=(F-FBEG)/STEPB
DIFF=STEPB-STEPA
HQ(IH)=TWO*(TEMP-GOPT(NFX))/DIFF
GOPT(NFX)=(GOPT(NFX)*STEPB-TEMP*STEPA)/DIFF
IF (STEPA*STEPB .LT. ZERO) THEN
IF (F .LT. FVAL(NF-N)) THEN
FVAL(NF)=FVAL(NF-N)
FVAL(NF-N)=F
IF (KOPT .EQ. NF) KOPT=NF-N
XPT(NF-N,NFX)=STEPB
XPT(NF,NFX)=STEPA
END IF
END IF
BMAT(1,NFX)=-(STEPA+STEPB)/(STEPA*STEPB)
BMAT(NF,NFX)=-HALF/XPT(NF-N,NFX)
BMAT(NF-N,NFX)=-BMAT(1,NFX)-BMAT(NF,NFX)
ZMAT(1,NFX)=DSQRT(TWO)/(STEPA*STEPB)
ZMAT(NF,NFX)=DSQRT(HALF)/RHOSQ
ZMAT(NF-N,NFX)=-ZMAT(1,NFX)-ZMAT(NF,NFX)
END IF
C
C Set the off-diagonal second derivatives of the Lagrange functions and
C the initial quadratic model.
C
ELSE
IH=(IPT*(IPT-1))/2+JPT
ZMAT(1,NFX)=RECIP
ZMAT(NF,NFX)=RECIP
ZMAT(IPT+1,NFX)=-RECIP
ZMAT(JPT+1,NFX)=-RECIP
TEMP=XPT(NF,IPT)*XPT(NF,JPT)
HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/TEMP
END IF
IF (NF .LT. NPT .AND. NF .LT. MAXFUN) GOTO 50
RETURN
END
SUBROUTINE RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,
1 FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,
2 KOPT,VLAG,PTSAUX,PTSID,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),XOPT(*),
1 GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*),
2 VLAG(*),PTSAUX(2,*),PTSID(*),W(*)
C
C The arguments N, NPT, XL, XU, IPRINT, MAXFUN, XBASE, XPT, FVAL, XOPT,
C GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as
C the corresponding arguments of BOBYQB on the entry to RESCUE.
C NF is maintained as the number of calls of CALFUN so far, except that
C NF is set to -1 if the value of MAXFUN prevents further progress.
C KOPT is maintained so that FVAL(KOPT) is the least calculated function
C value. Its correct value must be given on entry. It is updated if a
C new least function value is found, but the corresponding changes to
C XOPT and GOPT have to be made later by the calling program.
C DELTA is the current trust region radius.
C VLAG is a working space vector that will be used for the values of the
C provisional Lagrange functions at each of the interpolation points.
C They are part of a product that requires VLAG to be of length NDIM.
C PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and
C PTSAUX(2,J) specify the two positions of provisional interpolation
C points when a nonzero step is taken along e_J (the J-th coordinate
C direction) through XBASE+XOPT, as specified below. Usually these
C steps have length DELTA, but other lengths are chosen if necessary
C in order to satisfy the given bounds on the variables.
C PTSID is also a working space array. It has NPT components that denote
C provisional new positions of the original interpolation points, in
C case changes are needed to restore the linear independence of the
C interpolation conditions. The K-th point is a candidate for change
C if and only if PTSID(K) is nonzero. In this case let p and q be the
C integer parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p
C and q are both positive, the step from XBASE+XOPT to the new K-th
C interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise
C the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or
C p=0, respectively.
C The first NDIM+NPT elements of the array W are used for working space.
C The final elements of BMAT and ZMAT are set in a well-conditioned way
C to the values that are appropriate for the new interpolation points.
C The elements of GOPT, HQ and PQ are also revised to the values that are
C appropriate to the final quadratic model.
C
C Set some constants.
C
HALF=0.5D0
ONE=1.0D0
ZERO=0.0D0
NP=N+1
SFRAC=HALF/DFLOAT(NP)
NPTM=NPT-NP
C
C Shift the interpolation points so that XOPT becomes the origin, and set
C the elements of ZMAT to zero. The value of SUMPQ is required in the
C updating of HQ below. The squares of the distances from XOPT to the
C other interpolation points are set at the end of W. Increments of WINC
C may be added later to these squares to balance the consideration of
C the choice of point that is going to become current.
C
SUMPQ=ZERO
WINC=ZERO
DO 20 K=1,NPT
DISTSQ=ZERO
DO 10 J=1,N
XPT(K,J)=XPT(K,J)-XOPT(J)
10 DISTSQ=DISTSQ+XPT(K,J)**2
SUMPQ=SUMPQ+PQ(K)
W(NDIM+K)=DISTSQ
WINC=DMAX1(WINC,DISTSQ)
DO 20 J=1,NPTM
20 ZMAT(K,J)=ZERO
C
C Update HQ so that HQ and PQ define the second derivatives of the model
C after XBASE has been shifted to the trust region centre.
C
IH=0
DO 40 J=1,N
W(J)=HALF*SUMPQ*XOPT(J)
DO 30 K=1,NPT
30 W(J)=W(J)+PQ(K)*XPT(K,J)
DO 40 I=1,J
IH=IH+1
40 HQ(IH)=HQ(IH)+W(I)*XOPT(J)+W(J)*XOPT(I)
C
C Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and
C also set the elements of PTSAUX.
C
DO 50 J=1,N
XBASE(J)=XBASE(J)+XOPT(J)
SL(J)=SL(J)-XOPT(J)
SU(J)=SU(J)-XOPT(J)
XOPT(J)=ZERO
PTSAUX(1,J)=DMIN1(DELTA,SU(J))
PTSAUX(2,J)=DMAX1(-DELTA,SL(J))
IF (PTSAUX(1,J)+PTSAUX(2,J) .LT. ZERO) THEN
TEMP=PTSAUX(1,J)
PTSAUX(1,J)=PTSAUX(2,J)
PTSAUX(2,J)=TEMP
END IF
IF (DABS(PTSAUX(2,J)) .LT. HALF*DABS(PTSAUX(1,J))) THEN
PTSAUX(2,J)=HALF*PTSAUX(1,J)
END IF
DO 50 I=1,NDIM
50 BMAT(I,J)=ZERO
FBASE=FVAL(KOPT)
C
C Set the identifiers of the artificial interpolation points that are
C along a coordinate direction from XOPT, and set the corresponding
C nonzero elements of BMAT and ZMAT.
C
PTSID(1)=SFRAC
DO 60 J=1,N
JP=J+1
JPN=JP+N
PTSID(JP)=DFLOAT(J)+SFRAC
IF (JPN .LE. NPT) THEN
PTSID(JPN)=DFLOAT(J)/DFLOAT(NP)+SFRAC
TEMP=ONE/(PTSAUX(1,J)-PTSAUX(2,J))
BMAT(JP,J)=-TEMP+ONE/PTSAUX(1,J)
BMAT(JPN,J)=TEMP+ONE/PTSAUX(2,J)
BMAT(1,J)=-BMAT(JP,J)-BMAT(JPN,J)
ZMAT(1,J)=DSQRT(2.0D0)/DABS(PTSAUX(1,J)*PTSAUX(2,J))
ZMAT(JP,J)=ZMAT(1,J)*PTSAUX(2,J)*TEMP
ZMAT(JPN,J)=-ZMAT(1,J)*PTSAUX(1,J)*TEMP
ELSE
BMAT(1,J)=-ONE/PTSAUX(1,J)
BMAT(JP,J)=ONE/PTSAUX(1,J)
BMAT(J+NPT,J)=-HALF*PTSAUX(1,J)**2
END IF
60 CONTINUE
C
C Set any remaining identifiers with their nonzero elements of ZMAT.
C
IF (NPT .GE. N+NP) THEN
DO 70 K=2*NP,NPT
IW=(DFLOAT(K-NP)-HALF)/DFLOAT(N)
IP=K-NP-IW*N
IQ=IP+IW
IF (IQ .GT. N) IQ=IQ-N
PTSID(K)=DFLOAT(IP)+DFLOAT(IQ)/DFLOAT(NP)+SFRAC
TEMP=ONE/(PTSAUX(1,IP)*PTSAUX(1,IQ))
ZMAT(1,K-NP)=TEMP
ZMAT(IP+1,K-NP)=-TEMP
ZMAT(IQ+1,K-NP)=-TEMP
70 ZMAT(K,K-NP)=TEMP
END IF
NREM=NPT
KOLD=1
KNEW=KOPT
C
C Reorder the provisional points in the way that exchanges PTSID(KOLD)
C with PTSID(KNEW).
C
80 DO 90 J=1,N
TEMP=BMAT(KOLD,J)
BMAT(KOLD,J)=BMAT(KNEW,J)
90 BMAT(KNEW,J)=TEMP
DO 100 J=1,NPTM
TEMP=ZMAT(KOLD,J)
ZMAT(KOLD,J)=ZMAT(KNEW,J)
100 ZMAT(KNEW,J)=TEMP
PTSID(KOLD)=PTSID(KNEW)
PTSID(KNEW)=ZERO
W(NDIM+KNEW)=ZERO
NREM=NREM-1
IF (KNEW .NE. KOPT) THEN
TEMP=VLAG(KOLD)
VLAG(KOLD)=VLAG(KNEW)
VLAG(KNEW)=TEMP
C
C Update the BMAT and ZMAT matrices so that the status of the KNEW-th
C interpolation point can be changed from provisional to original. The
C branch to label 350 occurs if all the original points are reinstated.
C The nonnegative values of W(NDIM+K) are required in the search below.
C
CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
IF (NREM .EQ. 0) GOTO 350
DO 110 K=1,NPT
110 W(NDIM+K)=DABS(W(NDIM+K))
END IF
C
C Pick the index KNEW of an original interpolation point that has not
C yet replaced one of the provisional interpolation points, giving
C attention to the closeness to XOPT and to previous tries with KNEW.
C
120 DSQMIN=ZERO
DO 130 K=1,NPT
IF (W(NDIM+K) .GT. ZERO) THEN
IF (DSQMIN .EQ. ZERO .OR. W(NDIM+K) .LT. DSQMIN) THEN
KNEW=K
DSQMIN=W(NDIM+K)
END IF
END IF
130 CONTINUE
IF (DSQMIN .EQ. ZERO) GOTO 260
C
C Form the W-vector of the chosen original interpolation point.
C
DO 140 J=1,N
140 W(NPT+J)=XPT(KNEW,J)
DO 160 K=1,NPT
SUM=ZERO
IF (K .EQ. KOPT) THEN
CONTINUE
ELSE IF (PTSID(K) .EQ. ZERO) THEN
DO 150 J=1,N
150 SUM=SUM+W(NPT+J)*XPT(K,J)
ELSE
IP=PTSID(K)
IF (IP .GT. 0) SUM=W(NPT+IP)*PTSAUX(1,IP)
IQ=DFLOAT(NP)*PTSID(K)-DFLOAT(IP*NP)
IF (IQ .GT. 0) THEN
IW=1
IF (IP .EQ. 0) IW=2
SUM=SUM+W(NPT+IQ)*PTSAUX(IW,IQ)
END IF
END IF
160 W(K)=HALF*SUM*SUM
C
C Calculate VLAG and BETA for the required updating of the H matrix if
C XPT(KNEW,.) is reinstated in the set of interpolation points.
C
DO 180 K=1,NPT
SUM=ZERO
DO 170 J=1,N
170 SUM=SUM+BMAT(K,J)*W(NPT+J)
180 VLAG(K)=SUM
BETA=ZERO
DO 200 J=1,NPTM
SUM=ZERO
DO 190 K=1,NPT
190 SUM=SUM+ZMAT(K,J)*W(K)
BETA=BETA-SUM*SUM
DO 200 K=1,NPT
200 VLAG(K)=VLAG(K)+SUM*ZMAT(K,J)
BSUM=ZERO
DISTSQ=ZERO
DO 230 J=1,N
SUM=ZERO
DO 210 K=1,NPT
210 SUM=SUM+BMAT(K,J)*W(K)
JP=J+NPT
BSUM=BSUM+SUM*W(JP)
DO 220 IP=NPT+1,NDIM
220 SUM=SUM+BMAT(IP,J)*W(IP)
BSUM=BSUM+SUM*W(JP)
VLAG(JP)=SUM
230 DISTSQ=DISTSQ+XPT(KNEW,J)**2
BETA=HALF*DISTSQ*DISTSQ+BETA-BSUM
VLAG(KOPT)=VLAG(KOPT)+ONE
C
C KOLD is set to the index of the provisional interpolation point that is
C going to be deleted to make way for the KNEW-th original interpolation
C point. The choice of KOLD is governed by the avoidance of a small value
C of the denominator in the updating calculation of UPDATE.
C
DENOM=ZERO
VLMXSQ=ZERO
DO 250 K=1,NPT
IF (PTSID(K) .NE. ZERO) THEN
HDIAG=ZERO
DO 240 J=1,NPTM
240 HDIAG=HDIAG+ZMAT(K,J)**2
DEN=BETA*HDIAG+VLAG(K)**2
IF (DEN .GT. DENOM) THEN
KOLD=K
DENOM=DEN
END IF
END IF
250 VLMXSQ=DMAX1(VLMXSQ,VLAG(K)**2)
IF (DENOM .LE. 1.0D-2*VLMXSQ) THEN
W(NDIM+KNEW)=-W(NDIM+KNEW)-WINC
GOTO 120
END IF
GOTO 80
C
C When label 260 is reached, all the final positions of the interpolation
C points have been chosen although any changes have not been included yet
C in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart
C from the shift of XBASE, the updating of the quadratic model remains to
C be done. The following cycle through the new interpolation points begins
C by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero,
C except that a RETURN occurs if MAXFUN prohibits another value of F.
C
260 DO 340 KPT=1,NPT
IF (PTSID(KPT) .EQ. ZERO) GOTO 340
IF (NF .GE. MAXFUN) THEN
NF=-1
GOTO 350
END IF
IH=0
DO 270 J=1,N
W(J)=XPT(KPT,J)
XPT(KPT,J)=ZERO
TEMP=PQ(KPT)*W(J)
DO 270 I=1,J
IH=IH+1
270 HQ(IH)=HQ(IH)+TEMP*W(I)
PQ(KPT)=ZERO
IP=PTSID(KPT)
IQ=DFLOAT(NP)*PTSID(KPT)-DFLOAT(IP*NP)
IF (IP .GT. 0) THEN
XP=PTSAUX(1,IP)
XPT(KPT,IP)=XP
END IF
IF (IQ .GT. 0) THEN
XQ=PTSAUX(1,IQ)
IF (IP .EQ. 0) XQ=PTSAUX(2,IQ)
XPT(KPT,IQ)=XQ
END IF
C
C Set VQUAD to the value of the current model at the new point.
C
VQUAD=FBASE
IF (IP .GT. 0) THEN
IHP=(IP+IP*IP)/2
VQUAD=VQUAD+XP*(GOPT(IP)+HALF*XP*HQ(IHP))
END IF
IF (IQ .GT. 0) THEN
IHQ=(IQ+IQ*IQ)/2
VQUAD=VQUAD+XQ*(GOPT(IQ)+HALF*XQ*HQ(IHQ))
IF (IP .GT. 0) THEN
IW=MAX0(IHP,IHQ)-IABS(IP-IQ)
VQUAD=VQUAD+XP*XQ*HQ(IW)
END IF
END IF
DO 280 K=1,NPT
TEMP=ZERO
IF (IP .GT. 0) TEMP=TEMP+XP*XPT(K,IP)
IF (IQ .GT. 0) TEMP=TEMP+XQ*XPT(K,IQ)
280 VQUAD=VQUAD+HALF*PQ(K)*TEMP*TEMP
C
C Calculate F at the new interpolation point, and set DIFF to the factor
C that is going to multiply the KPT-th Lagrange function when the model
C is updated to provide interpolation to the new function value.
C
DO 290 I=1,N
W(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XPT(KPT,I)),XU(I))
IF (XPT(KPT,I) .EQ. SL(I)) W(I)=XL(I)
IF (XPT(KPT,I) .EQ. SU(I)) W(I)=XU(I)
290 CONTINUE
NF=NF+1
CALL CALFUN (N,W,F)
IF (IPRINT .EQ. 3) THEN
PRINT 300, NF,F,(W(I),I=1,N)
300 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
1 ' The corresponding X is:'/(2X,5D15.6))
END IF
FVAL(KPT)=F
IF (F .LT. FVAL(KOPT)) KOPT=KPT
DIFF=F-VQUAD
C
C Update the quadratic model. The RETURN from the subroutine occurs when
C all the new interpolation points are included in the model.
C
DO 310 I=1,N
310 GOPT(I)=GOPT(I)+DIFF*BMAT(KPT,I)
DO 330 K=1,NPT
SUM=ZERO
DO 320 J=1,NPTM
320 SUM=SUM+ZMAT(K,J)*ZMAT(KPT,J)
TEMP=DIFF*SUM
IF (PTSID(K) .EQ. ZERO) THEN
PQ(K)=PQ(K)+TEMP
ELSE
IP=PTSID(K)
IQ=DFLOAT(NP)*PTSID(K)-DFLOAT(IP*NP)
IHQ=(IQ*IQ+IQ)/2
IF (IP .EQ. 0) THEN
HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(2,IQ)**2
ELSE
IHP=(IP*IP+IP)/2
HQ(IHP)=HQ(IHP)+TEMP*PTSAUX(1,IP)**2
IF (IQ .GT. 0) THEN
HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(1,IQ)**2
IW=MAX0(IHP,IHQ)-IABS(IQ-IP)
HQ(IW)=HQ(IW)+TEMP*PTSAUX(1,IP)*PTSAUX(1,IQ)
END IF
END IF
END IF
330 CONTINUE
PTSID(KPT)=ZERO
340 CONTINUE
350 RETURN
END
SUBROUTINE TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,
1 XNEW,D,GNEW,XBDI,S,HS,HRED,DSQ,CRVMIN)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XPT(NPT,*),XOPT(*),GOPT(*),HQ(*),PQ(*),SL(*),SU(*),
1 XNEW(*),D(*),GNEW(*),XBDI(*),S(*),HS(*),HRED(*)
C
C The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same
C meanings as the corresponding arguments of BOBYQB.
C DELTA is the trust region radius for the present calculation, which
C seeks a small value of the quadratic model within distance DELTA of
C XOPT subject to the bounds on the variables.
C XNEW will be set to a new vector of variables that is approximately
C the one that minimizes the quadratic model within the trust region
C subject to the SL and SU constraints on the variables. It satisfies
C as equations the bounds that become active during the calculation.
C D is the calculated trial step from XOPT, generated iteratively from an
C initial value of zero. Thus XNEW is XOPT+D after the final iteration.
C GNEW holds the gradient of the quadratic model at XOPT+D. It is updated
C when D is updated.
C XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is
C set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the
C I-th variable has become fixed at a bound, the bound being SL(I) or
C SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This
C information is accumulated during the construction of XNEW.
C The arrays S, HS and HRED are also used for working space. They hold the
C current search direction, and the changes in the gradient of Q along S
C and the reduced D, respectively, where the reduced D is the same as D,
C except that the components of the fixed variables are zero.
C DSQ will be set to the square of the length of XNEW-XOPT.
C CRVMIN is set to zero if D reaches the trust region boundary. Otherwise
C it is set to the least curvature of H that occurs in the conjugate
C gradient searches that are not restricted by any constraints. The
C value CRVMIN=-1.0D0 is set, however, if all of these searches are
C constrained.
C
C A version of the truncated conjugate gradient is applied. If a line
C search is restricted by a constraint, then the procedure is restarted,
C the values of the variables that are at their bounds being fixed. If
C the trust region boundary is reached, then further changes may be made
C to D, each one being in the two dimensional space that is spanned
C by the current D and the gradient of Q at XOPT+D, staying on the trust
C region boundary. Termination occurs when the reduction in Q seems to
C be close to the greatest reduction that can be achieved.
C
C Set some constants.
C
HALF=0.5D0
ONE=1.0D0
ONEMIN=-1.0D0
ZERO=0.0D0
C
C The sign of GOPT(I) gives the sign of the change to the I-th variable
C that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether
C or not to fix the I-th variable at one of its bounds initially, with
C NACT being set to the number of fixed variables. D and GNEW are also
C set for the first iteration. DELSQ is the upper bound on the sum of
C squares of the free variables. QRED is the reduction in Q so far.
C
ITERC=0
NACT=0
SQSTP=ZERO
DO 10 I=1,N
XBDI(I)=ZERO
IF (XOPT(I) .LE. SL(I)) THEN
IF (GOPT(I) .GE. ZERO) XBDI(I)=ONEMIN
ELSE IF (XOPT(I) .GE. SU(I)) THEN
IF (GOPT(I) .LE. ZERO) XBDI(I)=ONE
END IF
IF (XBDI(I) .NE. ZERO) NACT=NACT+1
D(I)=ZERO
10 GNEW(I)=GOPT(I)
DELSQ=DELTA*DELTA
QRED=ZERO
CRVMIN=ONEMIN
C
C Set the next search direction of the conjugate gradient method. It is
C the steepest descent direction initially and when the iterations are
C restarted because a variable has just been fixed by a bound, and of
C course the components of the fixed variables are zero. ITERMAX is an
C upper bound on the indices of the conjugate gradient iterations.
C
20 BETA=ZERO
30 STEPSQ=ZERO
DO 40 I=1,N
IF (XBDI(I) .NE. ZERO) THEN
S(I)=ZERO
ELSE IF (BETA .EQ. ZERO) THEN
S(I)=-GNEW(I)
ELSE
S(I)=BETA*S(I)-GNEW(I)
END IF
40 STEPSQ=STEPSQ+S(I)**2
IF (STEPSQ .EQ. ZERO) GOTO 190
IF (BETA .EQ. ZERO) THEN
GREDSQ=STEPSQ
ITERMAX=ITERC+N-NACT
END IF
IF (GREDSQ*DELSQ .LE. 1.0D-4*QRED*QRED) GO TO 190
C
C Multiply the search direction by the second derivative matrix of Q and
C calculate some scalars for the choice of steplength. Then set BLEN to
C the length of the the step to the trust region boundary and STPLEN to
C the steplength, ignoring the simple bounds.
C
GOTO 210
50 RESID=DELSQ
DS=ZERO
SHS=ZERO
DO 60 I=1,N
IF (XBDI(I) .EQ. ZERO) THEN
RESID=RESID-D(I)**2
DS=DS+S(I)*D(I)
SHS=SHS+S(I)*HS(I)
END IF
60 CONTINUE
IF (RESID .LE. ZERO) GOTO 90
TEMP=DSQRT(STEPSQ*RESID+DS*DS)
IF (DS .LT. ZERO) THEN
BLEN=(TEMP-DS)/STEPSQ
ELSE
BLEN=RESID/(TEMP+DS)
END IF
STPLEN=BLEN
IF (SHS .GT. ZERO) THEN
STPLEN=DMIN1(BLEN,GREDSQ/SHS)
END IF
C
C Reduce STPLEN if necessary in order to preserve the simple bounds,
C letting IACT be the index of the new constrained variable.
C
IACT=0
DO 70 I=1,N
IF (S(I) .NE. ZERO) THEN
XSUM=XOPT(I)+D(I)
IF (S(I) .GT. ZERO) THEN
TEMP=(SU(I)-XSUM)/S(I)
ELSE
TEMP=(SL(I)-XSUM)/S(I)
END IF
IF (TEMP .LT. STPLEN) THEN
STPLEN=TEMP
IACT=I
END IF
END IF
70 CONTINUE
C
C Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q.
C
SDEC=ZERO
IF (STPLEN .GT. ZERO) THEN
ITERC=ITERC+1
TEMP=SHS/STEPSQ
IF (IACT .EQ. 0 .AND. TEMP .GT. ZERO) THEN
CRVMIN=DMIN1(CRVMIN,TEMP)
IF (CRVMIN .EQ. ONEMIN) CRVMIN=TEMP
END IF
GGSAV=GREDSQ
GREDSQ=ZERO
DO 80 I=1,N
GNEW(I)=GNEW(I)+STPLEN*HS(I)
IF (XBDI(I) .EQ. ZERO) GREDSQ=GREDSQ+GNEW(I)**2
80 D(I)=D(I)+STPLEN*S(I)
SDEC=DMAX1(STPLEN*(GGSAV-HALF*STPLEN*SHS),ZERO)
QRED=QRED+SDEC
END IF
C
C Restart the conjugate gradient method if it has hit a new bound.
C
IF (IACT .GT. 0) THEN
NACT=NACT+1
XBDI(IACT)=ONE
IF (S(IACT) .LT. ZERO) XBDI(IACT)=ONEMIN
DELSQ=DELSQ-D(IACT)**2
IF (DELSQ .LE. ZERO) GOTO 90
GOTO 20
END IF
C
C If STPLEN is less than BLEN, then either apply another conjugate
C gradient iteration or RETURN.
C
IF (STPLEN .LT. BLEN) THEN
IF (ITERC .EQ. ITERMAX) GOTO 190
IF (SDEC .LE. 0.01D0*QRED) GOTO 190
BETA=GREDSQ/GGSAV
GOTO 30
END IF
90 CRVMIN=ZERO
C
C Prepare for the alternative iteration by calculating some scalars
C and by multiplying the reduced D by the second derivative matrix of
C Q, where S holds the reduced D in the call of GGMULT.
C
100 IF (NACT .GE. N-1) GOTO 190
DREDSQ=ZERO
DREDG=ZERO
GREDSQ=ZERO
DO 110 I=1,N
IF (XBDI(I) .EQ. ZERO) THEN
DREDSQ=DREDSQ+D(I)**2
DREDG=DREDG+D(I)*GNEW(I)
GREDSQ=GREDSQ+GNEW(I)**2
S(I)=D(I)
ELSE
S(I)=ZERO
END IF
110 CONTINUE
ITCSAV=ITERC
GOTO 210
C
C Let the search direction S be a linear combination of the reduced D
C and the reduced G that is orthogonal to the reduced D.
C
120 ITERC=ITERC+1
TEMP=GREDSQ*DREDSQ-DREDG*DREDG
IF (TEMP .LE. 1.0D-4*QRED*QRED) GOTO 190
TEMP=DSQRT(TEMP)
DO 130 I=1,N
IF (XBDI(I) .EQ. ZERO) THEN
S(I)=(DREDG*D(I)-DREDSQ*GNEW(I))/TEMP
ELSE
S(I)=ZERO
END IF
130 CONTINUE
SREDG=-TEMP
C
C By considering the simple bounds on the variables, calculate an upper
C bound on the tangent of half the angle of the alternative iteration,
C namely ANGBD, except that, if already a free variable has reached a
C bound, there is a branch back to label 100 after fixing that variable.
C
ANGBD=ONE
IACT=0
DO 140 I=1,N
IF (XBDI(I) .EQ. ZERO) THEN
TEMPA=XOPT(I)+D(I)-SL(I)
TEMPB=SU(I)-XOPT(I)-D(I)
IF (TEMPA .LE. ZERO) THEN
NACT=NACT+1
XBDI(I)=ONEMIN
GOTO 100
ELSE IF (TEMPB .LE. ZERO) THEN
NACT=NACT+1
XBDI(I)=ONE
GOTO 100
END IF
RATIO=ONE
SSQ=D(I)**2+S(I)**2
TEMP=SSQ-(XOPT(I)-SL(I))**2
IF (TEMP .GT. ZERO) THEN
TEMP=DSQRT(TEMP)-S(I)
IF (ANGBD*TEMP .GT. TEMPA) THEN
ANGBD=TEMPA/TEMP
IACT=I
XSAV=ONEMIN
END IF
END IF
TEMP=SSQ-(SU(I)-XOPT(I))**2
IF (TEMP .GT. ZERO) THEN
TEMP=DSQRT(TEMP)+S(I)
IF (ANGBD*TEMP .GT. TEMPB) THEN
ANGBD=TEMPB/TEMP
IACT=I
XSAV=ONE
END IF
END IF
END IF
140 CONTINUE
C
C Calculate HHD and some curvatures for the alternative iteration.
C
GOTO 210
150 SHS=ZERO
DHS=ZERO
DHD=ZERO
DO 160 I=1,N
IF (XBDI(I) .EQ. ZERO) THEN
SHS=SHS+S(I)*HS(I)
DHS=DHS+D(I)*HS(I)
DHD=DHD+D(I)*HRED(I)
END IF
160 CONTINUE
C
C Seek the greatest reduction in Q for a range of equally spaced values
C of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of
C the alternative iteration.
C
REDMAX=ZERO
ISAV=0
REDSAV=ZERO
IU=17.0D0*ANGBD+3.1D0
DO 170 I=1,IU
ANGT=ANGBD*DFLOAT(I)/DFLOAT(IU)
STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
REDNEW=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
IF (REDNEW .GT. REDMAX) THEN
REDMAX=REDNEW
ISAV=I
RDPREV=REDSAV
ELSE IF (I .EQ. ISAV+1) THEN
RDNEXT=REDNEW
END IF
170 REDSAV=REDNEW
C
C Return if the reduction is zero. Otherwise, set the sine and cosine
C of the angle of the alternative iteration, and calculate SDEC.
C
IF (ISAV .EQ. 0) GOTO 190
IF (ISAV .LT. IU) THEN
TEMP=(RDNEXT-RDPREV)/(REDMAX+REDMAX-RDPREV-RDNEXT)
ANGT=ANGBD*(DFLOAT(ISAV)+HALF*TEMP)/DFLOAT(IU)
END IF
CTH=(ONE-ANGT*ANGT)/(ONE+ANGT*ANGT)
STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
SDEC=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
IF (SDEC .LE. ZERO) GOTO 190
C
C Update GNEW, D and HRED. If the angle of the alternative iteration
C is restricted by a bound on a free variable, that variable is fixed
C at the bound.
C
DREDG=ZERO
GREDSQ=ZERO
DO 180 I=1,N
GNEW(I)=GNEW(I)+(CTH-ONE)*HRED(I)+STH*HS(I)
IF (XBDI(I) .EQ. ZERO) THEN
D(I)=CTH*D(I)+STH*S(I)
DREDG=DREDG+D(I)*GNEW(I)
GREDSQ=GREDSQ+GNEW(I)**2
END IF
180 HRED(I)=CTH*HRED(I)+STH*HS(I)
QRED=QRED+SDEC
IF (IACT .GT. 0 .AND. ISAV .EQ. IU) THEN
NACT=NACT+1
XBDI(IACT)=XSAV
GOTO 100
END IF
C
C If SDEC is sufficiently small, then RETURN after setting XNEW to
C XOPT+D, giving careful attention to the bounds.
C
IF (SDEC .GT. 0.01D0*QRED) GOTO 120
190 DSQ=ZERO
DO 200 I=1,N
XNEW(I)=DMAX1(DMIN1(XOPT(I)+D(I),SU(I)),SL(I))
IF (XBDI(I) .EQ. ONEMIN) XNEW(I)=SL(I)
IF (XBDI(I) .EQ. ONE) XNEW(I)=SU(I)
D(I)=XNEW(I)-XOPT(I)
200 DSQ=DSQ+D(I)**2
RETURN
C The following instructions multiply the current S-vector by the second
C derivative matrix of the quadratic model, putting the product in HS.
C They are reached from three different parts of the software above and
C they can be regarded as an external subroutine.
C
210 IH=0
DO 220 J=1,N
HS(J)=ZERO
DO 220 I=1,J
IH=IH+1
IF (I .LT. J) HS(J)=HS(J)+HQ(IH)*S(I)
220 HS(I)=HS(I)+HQ(IH)*S(J)
DO 250 K=1,NPT
IF (PQ(K) .NE. ZERO) THEN
TEMP=ZERO
DO 230 J=1,N
230 TEMP=TEMP+XPT(K,J)*S(J)
TEMP=TEMP*PQ(K)
DO 240 I=1,N
240 HS(I)=HS(I)+TEMP*XPT(K,I)
END IF
250 CONTINUE
IF (CRVMIN .NE. ZERO) GOTO 50
IF (ITERC .GT. ITCSAV) GOTO 150
DO 260 I=1,N
260 HRED(I)=HS(I)
GOTO 120
END
SUBROUTINE UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,
1 KNEW,W)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*)
C
C The arrays BMAT and ZMAT are updated, as required by the new position
C of the interpolation point that has the index KNEW. The vector VLAG has
C N+NPT components, set on entry to the first NPT and last N components
C of the product Hw in equation (4.11) of the Powell (2006) paper on
C NEWUOA. Further, BETA is set on entry to the value of the parameter
C with that name, and DENOM is set to the denominator of the updating
C formula. Elements of ZMAT may be treated as zero if their moduli are
C at most ZTEST. The first NDIM elements of W are used for working space.
C
C Set some constants.
C
ONE=1.0D0
ZERO=0.0D0
NPTM=NPT-N-1
ZTEST=ZERO
DO 10 K=1,NPT
DO 10 J=1,NPTM
10 ZTEST=DMAX1(ZTEST,DABS(ZMAT(K,J)))
ZTEST=1.0D-20*ZTEST
C
C Apply the rotations that put zeros in the KNEW-th row of ZMAT.
C
JL=1
DO 30 J=2,NPTM
IF (DABS(ZMAT(KNEW,J)) .GT. ZTEST) THEN
TEMP=DSQRT(ZMAT(KNEW,1)**2+ZMAT(KNEW,J)**2)
TEMPA=ZMAT(KNEW,1)/TEMP
TEMPB=ZMAT(KNEW,J)/TEMP
DO 20 I=1,NPT
TEMP=TEMPA*ZMAT(I,1)+TEMPB*ZMAT(I,J)
ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,1)
20 ZMAT(I,1)=TEMP
END IF
ZMAT(KNEW,J)=ZERO
30 CONTINUE
C
C Put the first NPT components of the KNEW-th column of HLAG into W,
C and calculate the parameters of the updating formula.
C
DO 40 I=1,NPT
W(I)=ZMAT(KNEW,1)*ZMAT(I,1)
40 CONTINUE
ALPHA=W(KNEW)
TAU=VLAG(KNEW)
VLAG(KNEW)=VLAG(KNEW)-ONE
C
C Complete the updating of ZMAT.
C
TEMP=DSQRT(DENOM)
TEMPB=ZMAT(KNEW,1)/TEMP
TEMPA=TAU/TEMP
DO 50 I=1,NPT
50 ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
C
C Finally, update the matrix BMAT.
C
DO 60 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 60 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)
60 CONTINUE
RETURN
END
Submit an article or tip