123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279 |
- 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
|