123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380 |
- 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 and
- C by multiplying the reduced D by the second derivative matrix of Q.
- 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
|