trsbox.f 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380
  1. SUBROUTINE TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,
  2. 1 XNEW,D,GNEW,XBDI,S,HS,HRED,DSQ,CRVMIN)
  3. IMPLICIT REAL*8 (A-H,O-Z)
  4. DIMENSION XPT(NPT,*),XOPT(*),GOPT(*),HQ(*),PQ(*),SL(*),SU(*),
  5. 1 XNEW(*),D(*),GNEW(*),XBDI(*),S(*),HS(*),HRED(*)
  6. C
  7. C The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same
  8. C meanings as the corresponding arguments of BOBYQB.
  9. C DELTA is the trust region radius for the present calculation, which
  10. C seeks a small value of the quadratic model within distance DELTA of
  11. C XOPT subject to the bounds on the variables.
  12. C XNEW will be set to a new vector of variables that is approximately
  13. C the one that minimizes the quadratic model within the trust region
  14. C subject to the SL and SU constraints on the variables. It satisfies
  15. C as equations the bounds that become active during the calculation.
  16. C D is the calculated trial step from XOPT, generated iteratively from an
  17. C initial value of zero. Thus XNEW is XOPT+D after the final iteration.
  18. C GNEW holds the gradient of the quadratic model at XOPT+D. It is updated
  19. C when D is updated.
  20. C XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is
  21. C set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the
  22. C I-th variable has become fixed at a bound, the bound being SL(I) or
  23. C SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This
  24. C information is accumulated during the construction of XNEW.
  25. C The arrays S, HS and HRED are also used for working space. They hold the
  26. C current search direction, and the changes in the gradient of Q along S
  27. C and the reduced D, respectively, where the reduced D is the same as D,
  28. C except that the components of the fixed variables are zero.
  29. C DSQ will be set to the square of the length of XNEW-XOPT.
  30. C CRVMIN is set to zero if D reaches the trust region boundary. Otherwise
  31. C it is set to the least curvature of H that occurs in the conjugate
  32. C gradient searches that are not restricted by any constraints. The
  33. C value CRVMIN=-1.0D0 is set, however, if all of these searches are
  34. C constrained.
  35. C
  36. C A version of the truncated conjugate gradient is applied. If a line
  37. C search is restricted by a constraint, then the procedure is restarted,
  38. C the values of the variables that are at their bounds being fixed. If
  39. C the trust region boundary is reached, then further changes may be made
  40. C to D, each one being in the two dimensional space that is spanned
  41. C by the current D and the gradient of Q at XOPT+D, staying on the trust
  42. C region boundary. Termination occurs when the reduction in Q seems to
  43. C be close to the greatest reduction that can be achieved.
  44. C
  45. C Set some constants.
  46. C
  47. HALF=0.5D0
  48. ONE=1.0D0
  49. ONEMIN=-1.0D0
  50. ZERO=0.0D0
  51. C
  52. C The sign of GOPT(I) gives the sign of the change to the I-th variable
  53. C that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether
  54. C or not to fix the I-th variable at one of its bounds initially, with
  55. C NACT being set to the number of fixed variables. D and GNEW are also
  56. C set for the first iteration. DELSQ is the upper bound on the sum of
  57. C squares of the free variables. QRED is the reduction in Q so far.
  58. C
  59. ITERC=0
  60. NACT=0
  61. SQSTP=ZERO
  62. DO 10 I=1,N
  63. XBDI(I)=ZERO
  64. IF (XOPT(I) .LE. SL(I)) THEN
  65. IF (GOPT(I) .GE. ZERO) XBDI(I)=ONEMIN
  66. ELSE IF (XOPT(I) .GE. SU(I)) THEN
  67. IF (GOPT(I) .LE. ZERO) XBDI(I)=ONE
  68. END IF
  69. IF (XBDI(I) .NE. ZERO) NACT=NACT+1
  70. D(I)=ZERO
  71. 10 GNEW(I)=GOPT(I)
  72. DELSQ=DELTA*DELTA
  73. QRED=ZERO
  74. CRVMIN=ONEMIN
  75. C
  76. C Set the next search direction of the conjugate gradient method. It is
  77. C the steepest descent direction initially and when the iterations are
  78. C restarted because a variable has just been fixed by a bound, and of
  79. C course the components of the fixed variables are zero. ITERMAX is an
  80. C upper bound on the indices of the conjugate gradient iterations.
  81. C
  82. 20 BETA=ZERO
  83. 30 STEPSQ=ZERO
  84. DO 40 I=1,N
  85. IF (XBDI(I) .NE. ZERO) THEN
  86. S(I)=ZERO
  87. ELSE IF (BETA .EQ. ZERO) THEN
  88. S(I)=-GNEW(I)
  89. ELSE
  90. S(I)=BETA*S(I)-GNEW(I)
  91. END IF
  92. 40 STEPSQ=STEPSQ+S(I)**2
  93. IF (STEPSQ .EQ. ZERO) GOTO 190
  94. IF (BETA .EQ. ZERO) THEN
  95. GREDSQ=STEPSQ
  96. ITERMAX=ITERC+N-NACT
  97. END IF
  98. IF (GREDSQ*DELSQ .LE. 1.0D-4*QRED*QRED) GO TO 190
  99. C
  100. C Multiply the search direction by the second derivative matrix of Q and
  101. C calculate some scalars for the choice of steplength. Then set BLEN to
  102. C the length of the the step to the trust region boundary and STPLEN to
  103. C the steplength, ignoring the simple bounds.
  104. C
  105. GOTO 210
  106. 50 RESID=DELSQ
  107. DS=ZERO
  108. SHS=ZERO
  109. DO 60 I=1,N
  110. IF (XBDI(I) .EQ. ZERO) THEN
  111. RESID=RESID-D(I)**2
  112. DS=DS+S(I)*D(I)
  113. SHS=SHS+S(I)*HS(I)
  114. END IF
  115. 60 CONTINUE
  116. IF (RESID .LE. ZERO) GOTO 90
  117. TEMP=DSQRT(STEPSQ*RESID+DS*DS)
  118. IF (DS .LT. ZERO) THEN
  119. BLEN=(TEMP-DS)/STEPSQ
  120. ELSE
  121. BLEN=RESID/(TEMP+DS)
  122. END IF
  123. STPLEN=BLEN
  124. IF (SHS .GT. ZERO) THEN
  125. STPLEN=DMIN1(BLEN,GREDSQ/SHS)
  126. END IF
  127. C
  128. C Reduce STPLEN if necessary in order to preserve the simple bounds,
  129. C letting IACT be the index of the new constrained variable.
  130. C
  131. IACT=0
  132. DO 70 I=1,N
  133. IF (S(I) .NE. ZERO) THEN
  134. XSUM=XOPT(I)+D(I)
  135. IF (S(I) .GT. ZERO) THEN
  136. TEMP=(SU(I)-XSUM)/S(I)
  137. ELSE
  138. TEMP=(SL(I)-XSUM)/S(I)
  139. END IF
  140. IF (TEMP .LT. STPLEN) THEN
  141. STPLEN=TEMP
  142. IACT=I
  143. END IF
  144. END IF
  145. 70 CONTINUE
  146. C
  147. C Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q.
  148. C
  149. SDEC=ZERO
  150. IF (STPLEN .GT. ZERO) THEN
  151. ITERC=ITERC+1
  152. TEMP=SHS/STEPSQ
  153. IF (IACT .EQ. 0 .AND. TEMP .GT. ZERO) THEN
  154. CRVMIN=DMIN1(CRVMIN,TEMP)
  155. IF (CRVMIN .EQ. ONEMIN) CRVMIN=TEMP
  156. END IF
  157. GGSAV=GREDSQ
  158. GREDSQ=ZERO
  159. DO 80 I=1,N
  160. GNEW(I)=GNEW(I)+STPLEN*HS(I)
  161. IF (XBDI(I) .EQ. ZERO) GREDSQ=GREDSQ+GNEW(I)**2
  162. 80 D(I)=D(I)+STPLEN*S(I)
  163. SDEC=DMAX1(STPLEN*(GGSAV-HALF*STPLEN*SHS),ZERO)
  164. QRED=QRED+SDEC
  165. END IF
  166. C
  167. C Restart the conjugate gradient method if it has hit a new bound.
  168. C
  169. IF (IACT .GT. 0) THEN
  170. NACT=NACT+1
  171. XBDI(IACT)=ONE
  172. IF (S(IACT) .LT. ZERO) XBDI(IACT)=ONEMIN
  173. DELSQ=DELSQ-D(IACT)**2
  174. IF (DELSQ .LE. ZERO) GOTO 90
  175. GOTO 20
  176. END IF
  177. C
  178. C If STPLEN is less than BLEN, then either apply another conjugate
  179. C gradient iteration or RETURN.
  180. C
  181. IF (STPLEN .LT. BLEN) THEN
  182. IF (ITERC .EQ. ITERMAX) GOTO 190
  183. IF (SDEC .LE. 0.01D0*QRED) GOTO 190
  184. BETA=GREDSQ/GGSAV
  185. GOTO 30
  186. END IF
  187. 90 CRVMIN=ZERO
  188. C
  189. C Prepare for the alternative iteration by calculating some scalars and
  190. C by multiplying the reduced D by the second derivative matrix of Q.
  191. C
  192. 100 IF (NACT .GE. N-1) GOTO 190
  193. DREDSQ=ZERO
  194. DREDG=ZERO
  195. GREDSQ=ZERO
  196. DO 110 I=1,N
  197. IF (XBDI(I) .EQ. ZERO) THEN
  198. DREDSQ=DREDSQ+D(I)**2
  199. DREDG=DREDG+D(I)*GNEW(I)
  200. GREDSQ=GREDSQ+GNEW(I)**2
  201. S(I)=D(I)
  202. ELSE
  203. S(I)=ZERO
  204. END IF
  205. 110 CONTINUE
  206. ITCSAV=ITERC
  207. GOTO 210
  208. C
  209. C Let the search direction S be a linear combination of the reduced D
  210. C and the reduced G that is orthogonal to the reduced D.
  211. C
  212. 120 ITERC=ITERC+1
  213. TEMP=GREDSQ*DREDSQ-DREDG*DREDG
  214. IF (TEMP .LE. 1.0D-4*QRED*QRED) GOTO 190
  215. TEMP=DSQRT(TEMP)
  216. DO 130 I=1,N
  217. IF (XBDI(I) .EQ. ZERO) THEN
  218. S(I)=(DREDG*D(I)-DREDSQ*GNEW(I))/TEMP
  219. ELSE
  220. S(I)=ZERO
  221. END IF
  222. 130 CONTINUE
  223. SREDG=-TEMP
  224. C
  225. C By considering the simple bounds on the variables, calculate an upper
  226. C bound on the tangent of half the angle of the alternative iteration,
  227. C namely ANGBD, except that, if already a free variable has reached a
  228. C bound, there is a branch back to label 100 after fixing that variable.
  229. C
  230. ANGBD=ONE
  231. IACT=0
  232. DO 140 I=1,N
  233. IF (XBDI(I) .EQ. ZERO) THEN
  234. TEMPA=XOPT(I)+D(I)-SL(I)
  235. TEMPB=SU(I)-XOPT(I)-D(I)
  236. IF (TEMPA .LE. ZERO) THEN
  237. NACT=NACT+1
  238. XBDI(I)=ONEMIN
  239. GOTO 100
  240. ELSE IF (TEMPB .LE. ZERO) THEN
  241. NACT=NACT+1
  242. XBDI(I)=ONE
  243. GOTO 100
  244. END IF
  245. RATIO=ONE
  246. SSQ=D(I)**2+S(I)**2
  247. TEMP=SSQ-(XOPT(I)-SL(I))**2
  248. IF (TEMP .GT. ZERO) THEN
  249. TEMP=DSQRT(TEMP)-S(I)
  250. IF (ANGBD*TEMP .GT. TEMPA) THEN
  251. ANGBD=TEMPA/TEMP
  252. IACT=I
  253. XSAV=ONEMIN
  254. END IF
  255. END IF
  256. TEMP=SSQ-(SU(I)-XOPT(I))**2
  257. IF (TEMP .GT. ZERO) THEN
  258. TEMP=DSQRT(TEMP)+S(I)
  259. IF (ANGBD*TEMP .GT. TEMPB) THEN
  260. ANGBD=TEMPB/TEMP
  261. IACT=I
  262. XSAV=ONE
  263. END IF
  264. END IF
  265. END IF
  266. 140 CONTINUE
  267. C
  268. C Calculate HHD and some curvatures for the alternative iteration.
  269. C
  270. GOTO 210
  271. 150 SHS=ZERO
  272. DHS=ZERO
  273. DHD=ZERO
  274. DO 160 I=1,N
  275. IF (XBDI(I) .EQ. ZERO) THEN
  276. SHS=SHS+S(I)*HS(I)
  277. DHS=DHS+D(I)*HS(I)
  278. DHD=DHD+D(I)*HRED(I)
  279. END IF
  280. 160 CONTINUE
  281. C
  282. C Seek the greatest reduction in Q for a range of equally spaced values
  283. C of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of
  284. C the alternative iteration.
  285. C
  286. REDMAX=ZERO
  287. ISAV=0
  288. REDSAV=ZERO
  289. IU=17.0D0*ANGBD+3.1D0
  290. DO 170 I=1,IU
  291. ANGT=ANGBD*DFLOAT(I)/DFLOAT(IU)
  292. STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
  293. TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
  294. REDNEW=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
  295. IF (REDNEW .GT. REDMAX) THEN
  296. REDMAX=REDNEW
  297. ISAV=I
  298. RDPREV=REDSAV
  299. ELSE IF (I .EQ. ISAV+1) THEN
  300. RDNEXT=REDNEW
  301. END IF
  302. 170 REDSAV=REDNEW
  303. C
  304. C Return if the reduction is zero. Otherwise, set the sine and cosine
  305. C of the angle of the alternative iteration, and calculate SDEC.
  306. C
  307. IF (ISAV .EQ. 0) GOTO 190
  308. IF (ISAV .LT. IU) THEN
  309. TEMP=(RDNEXT-RDPREV)/(REDMAX+REDMAX-RDPREV-RDNEXT)
  310. ANGT=ANGBD*(DFLOAT(ISAV)+HALF*TEMP)/DFLOAT(IU)
  311. END IF
  312. CTH=(ONE-ANGT*ANGT)/(ONE+ANGT*ANGT)
  313. STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
  314. TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
  315. SDEC=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
  316. IF (SDEC .LE. ZERO) GOTO 190
  317. C
  318. C Update GNEW, D and HRED. If the angle of the alternative iteration
  319. C is restricted by a bound on a free variable, that variable is fixed
  320. C at the bound.
  321. C
  322. DREDG=ZERO
  323. GREDSQ=ZERO
  324. DO 180 I=1,N
  325. GNEW(I)=GNEW(I)+(CTH-ONE)*HRED(I)+STH*HS(I)
  326. IF (XBDI(I) .EQ. ZERO) THEN
  327. D(I)=CTH*D(I)+STH*S(I)
  328. DREDG=DREDG+D(I)*GNEW(I)
  329. GREDSQ=GREDSQ+GNEW(I)**2
  330. END IF
  331. 180 HRED(I)=CTH*HRED(I)+STH*HS(I)
  332. QRED=QRED+SDEC
  333. IF (IACT .GT. 0 .AND. ISAV .EQ. IU) THEN
  334. NACT=NACT+1
  335. XBDI(IACT)=XSAV
  336. GOTO 100
  337. END IF
  338. C
  339. C If SDEC is sufficiently small, then RETURN after setting XNEW to
  340. C XOPT+D, giving careful attention to the bounds.
  341. C
  342. IF (SDEC .GT. 0.01D0*QRED) GOTO 120
  343. 190 DSQ=ZERO
  344. DO 200 I=1,N
  345. XNEW(I)=DMAX1(DMIN1(XOPT(I)+D(I),SU(I)),SL(I))
  346. IF (XBDI(I) .EQ. ONEMIN) XNEW(I)=SL(I)
  347. IF (XBDI(I) .EQ. ONE) XNEW(I)=SU(I)
  348. D(I)=XNEW(I)-XOPT(I)
  349. 200 DSQ=DSQ+D(I)**2
  350. RETURN
  351. C The following instructions multiply the current S-vector by the second
  352. C derivative matrix of the quadratic model, putting the product in HS.
  353. C They are reached from three different parts of the software above and
  354. C they can be regarded as an external subroutine.
  355. C
  356. 210 IH=0
  357. DO 220 J=1,N
  358. HS(J)=ZERO
  359. DO 220 I=1,J
  360. IH=IH+1
  361. IF (I .LT. J) HS(J)=HS(J)+HQ(IH)*S(I)
  362. 220 HS(I)=HS(I)+HQ(IH)*S(J)
  363. DO 250 K=1,NPT
  364. IF (PQ(K) .NE. ZERO) THEN
  365. TEMP=ZERO
  366. DO 230 J=1,N
  367. 230 TEMP=TEMP+XPT(K,J)*S(J)
  368. TEMP=TEMP*PQ(K)
  369. DO 240 I=1,N
  370. 240 HS(I)=HS(I)+TEMP*XPT(K,I)
  371. END IF
  372. 250 CONTINUE
  373. IF (CRVMIN .NE. ZERO) GOTO 50
  374. IF (ITERC .GT. ITCSAV) GOTO 150
  375. DO 260 I=1,N
  376. 260 HRED(I)=HS(I)
  377. GOTO 120
  378. END