rescue.f 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385
  1. SUBROUTINE RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,
  2. 1 FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,
  3. 2 KOPT,VLAG,PTSAUX,PTSID,W)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. DIMENSION XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),XOPT(*),
  6. 1 GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*),
  7. 2 VLAG(*),PTSAUX(2,*),PTSID(*),W(*)
  8. C
  9. C The arguments N, NPT, XL, XU, IPRINT, MAXFUN, XBASE, XPT, FVAL, XOPT,
  10. C GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as
  11. C the corresponding arguments of BOBYQB on the entry to RESCUE.
  12. C NF is maintained as the number of calls of CALFUN so far, except that
  13. C NF is set to -1 if the value of MAXFUN prevents further progress.
  14. C KOPT is maintained so that FVAL(KOPT) is the least calculated function
  15. C value. Its correct value must be given on entry. It is updated if a
  16. C new least function value is found, but the corresponding changes to
  17. C XOPT and GOPT have to be made later by the calling program.
  18. C DELTA is the current trust region radius.
  19. C VLAG is a working space vector that will be used for the values of the
  20. C provisional Lagrange functions at each of the interpolation points.
  21. C They are part of a product that requires VLAG to be of length NDIM.
  22. C PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and
  23. C PTSAUX(2,J) specify the two positions of provisional interpolation
  24. C points when a nonzero step is taken along e_J (the J-th coordinate
  25. C direction) through XBASE+XOPT, as specified below. Usually these
  26. C steps have length DELTA, but other lengths are chosen if necessary
  27. C in order to satisfy the given bounds on the variables.
  28. C PTSID is also a working space array. It has NPT components that denote
  29. C provisional new positions of the original interpolation points, in
  30. C case changes are needed to restore the linear independence of the
  31. C interpolation conditions. The K-th point is a candidate for change
  32. C if and only if PTSID(K) is nonzero. In this case let p and q be the
  33. C integer parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p
  34. C and q are both positive, the step from XBASE+XOPT to the new K-th
  35. C interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise
  36. C the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or
  37. C p=0, respectively.
  38. C The first NDIM+NPT elements of the array W are used for working space.
  39. C The final elements of BMAT and ZMAT are set in a well-conditioned way
  40. C to the values that are appropriate for the new interpolation points.
  41. C The elements of GOPT, HQ and PQ are also revised to the values that are
  42. C appropriate to the final quadratic model.
  43. C
  44. C Set some constants.
  45. C
  46. HALF=0.5D0
  47. ONE=1.0D0
  48. ZERO=0.0D0
  49. NP=N+1
  50. SFRAC=HALF/DFLOAT(NP)
  51. NPTM=NPT-NP
  52. C
  53. C Shift the interpolation points so that XOPT becomes the origin, and set
  54. C the elements of ZMAT to zero. The value of SUMPQ is required in the
  55. C updating of HQ below. The squares of the distances from XOPT to the
  56. C other interpolation points are set at the end of W. Increments of WINC
  57. C may be added later to these squares to balance the consideration of
  58. C the choice of point that is going to become current.
  59. C
  60. SUMPQ=ZERO
  61. WINC=ZERO
  62. DO 20 K=1,NPT
  63. DISTSQ=ZERO
  64. DO 10 J=1,N
  65. XPT(K,J)=XPT(K,J)-XOPT(J)
  66. 10 DISTSQ=DISTSQ+XPT(K,J)**2
  67. SUMPQ=SUMPQ+PQ(K)
  68. W(NDIM+K)=DISTSQ
  69. WINC=DMAX1(WINC,DISTSQ)
  70. DO 20 J=1,NPTM
  71. 20 ZMAT(K,J)=ZERO
  72. C
  73. C Update HQ so that HQ and PQ define the second derivatives of the model
  74. C after XBASE has been shifted to the trust region centre.
  75. C
  76. IH=0
  77. DO 40 J=1,N
  78. W(J)=HALF*SUMPQ*XOPT(J)
  79. DO 30 K=1,NPT
  80. 30 W(J)=W(J)+PQ(K)*XPT(K,J)
  81. DO 40 I=1,J
  82. IH=IH+1
  83. 40 HQ(IH)=HQ(IH)+W(I)*XOPT(J)+W(J)*XOPT(I)
  84. C
  85. C Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and
  86. C also set the elements of PTSAUX.
  87. C
  88. DO 50 J=1,N
  89. XBASE(J)=XBASE(J)+XOPT(J)
  90. SL(J)=SL(J)-XOPT(J)
  91. SU(J)=SU(J)-XOPT(J)
  92. XOPT(J)=ZERO
  93. PTSAUX(1,J)=DMIN1(DELTA,SU(J))
  94. PTSAUX(2,J)=DMAX1(-DELTA,SL(J))
  95. IF (PTSAUX(1,J)+PTSAUX(2,J) .LT. ZERO) THEN
  96. TEMP=PTSAUX(1,J)
  97. PTSAUX(1,J)=PTSAUX(2,J)
  98. PTSAUX(2,J)=TEMP
  99. END IF
  100. IF (DABS(PTSAUX(2,J)) .LT. HALF*DABS(PTSAUX(1,J))) THEN
  101. PTSAUX(2,J)=HALF*PTSAUX(1,J)
  102. END IF
  103. DO 50 I=1,NDIM
  104. 50 BMAT(I,J)=ZERO
  105. FBASE=FVAL(KOPT)
  106. C
  107. C Set the identifiers of the artificial interpolation points that are
  108. C along a coordinate direction from XOPT, and set the corresponding
  109. C nonzero elements of BMAT and ZMAT.
  110. C
  111. PTSID(1)=SFRAC
  112. DO 60 J=1,N
  113. JP=J+1
  114. JPN=JP+N
  115. PTSID(JP)=DFLOAT(J)+SFRAC
  116. IF (JPN .LE. NPT) THEN
  117. PTSID(JPN)=DFLOAT(J)/DFLOAT(NP)+SFRAC
  118. TEMP=ONE/(PTSAUX(1,J)-PTSAUX(2,J))
  119. BMAT(JP,J)=-TEMP+ONE/PTSAUX(1,J)
  120. BMAT(JPN,J)=TEMP+ONE/PTSAUX(2,J)
  121. BMAT(1,J)=-BMAT(JP,J)-BMAT(JPN,J)
  122. ZMAT(1,J)=DSQRT(2.0D0)/DABS(PTSAUX(1,J)*PTSAUX(2,J))
  123. ZMAT(JP,J)=ZMAT(1,J)*PTSAUX(2,J)*TEMP
  124. ZMAT(JPN,J)=-ZMAT(1,J)*PTSAUX(1,J)*TEMP
  125. ELSE
  126. BMAT(1,J)=-ONE/PTSAUX(1,J)
  127. BMAT(JP,J)=ONE/PTSAUX(1,J)
  128. BMAT(J+NPT,J)=-HALF*PTSAUX(1,J)**2
  129. END IF
  130. 60 CONTINUE
  131. C
  132. C Set any remaining identifiers with their nonzero elements of ZMAT.
  133. C
  134. IF (NPT .GE. N+NP) THEN
  135. DO 70 K=2*NP,NPT
  136. IW=(DFLOAT(K-NP)-HALF)/DFLOAT(N)
  137. IP=K-NP-IW*N
  138. IQ=IP+IW
  139. IF (IQ .GT. N) IQ=IQ-N
  140. PTSID(K)=DFLOAT(IP)+DFLOAT(IQ)/DFLOAT(NP)+SFRAC
  141. TEMP=ONE/(PTSAUX(1,IP)*PTSAUX(1,IQ))
  142. ZMAT(1,K-NP)=TEMP
  143. ZMAT(IP+1,K-NP)=-TEMP
  144. ZMAT(IQ+1,K-NP)=-TEMP
  145. 70 ZMAT(K,K-NP)=TEMP
  146. END IF
  147. NREM=NPT
  148. KOLD=1
  149. KNEW=KOPT
  150. C
  151. C Reorder the provisional points in the way that exchanges PTSID(KOLD)
  152. C with PTSID(KNEW).
  153. C
  154. 80 DO 90 J=1,N
  155. TEMP=BMAT(KOLD,J)
  156. BMAT(KOLD,J)=BMAT(KNEW,J)
  157. 90 BMAT(KNEW,J)=TEMP
  158. DO 100 J=1,NPTM
  159. TEMP=ZMAT(KOLD,J)
  160. ZMAT(KOLD,J)=ZMAT(KNEW,J)
  161. 100 ZMAT(KNEW,J)=TEMP
  162. PTSID(KOLD)=PTSID(KNEW)
  163. PTSID(KNEW)=ZERO
  164. W(NDIM+KNEW)=ZERO
  165. NREM=NREM-1
  166. IF (KNEW .NE. KOPT) THEN
  167. TEMP=VLAG(KOLD)
  168. VLAG(KOLD)=VLAG(KNEW)
  169. VLAG(KNEW)=TEMP
  170. C
  171. C Update the BMAT and ZMAT matrices so that the status of the KNEW-th
  172. C interpolation point can be changed from provisional to original. The
  173. C branch to label 350 occurs if all the original points are reinstated.
  174. C The nonnegative values of W(NDIM+K) are required in the search below.
  175. C
  176. CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
  177. IF (NREM .EQ. 0) GOTO 350
  178. DO 110 K=1,NPT
  179. 110 W(NDIM+K)=DABS(W(NDIM+K))
  180. END IF
  181. C
  182. C Pick the index KNEW of an original interpolation point that has not
  183. C yet replaced one of the provisional interpolation points, giving
  184. C attention to the closeness to XOPT and to previous tries with KNEW.
  185. C
  186. 120 DSQMIN=ZERO
  187. DO 130 K=1,NPT
  188. IF (W(NDIM+K) .GT. ZERO) THEN
  189. IF (DSQMIN .EQ. ZERO .OR. W(NDIM+K) .LT. DSQMIN) THEN
  190. KNEW=K
  191. DSQMIN=W(NDIM+K)
  192. END IF
  193. END IF
  194. 130 CONTINUE
  195. IF (DSQMIN .EQ. ZERO) GOTO 260
  196. C
  197. C Form the W-vector of the chosen original interpolation point.
  198. C
  199. DO 140 J=1,N
  200. 140 W(NPT+J)=XPT(KNEW,J)
  201. DO 160 K=1,NPT
  202. SUM=ZERO
  203. IF (K .EQ. KOPT) THEN
  204. CONTINUE
  205. ELSE IF (PTSID(K) .EQ. ZERO) THEN
  206. DO 150 J=1,N
  207. 150 SUM=SUM+W(NPT+J)*XPT(K,J)
  208. ELSE
  209. IP=PTSID(K)
  210. IF (IP .GT. 0) SUM=W(NPT+IP)*PTSAUX(1,IP)
  211. IQ=DFLOAT(NP)*PTSID(K)-DFLOAT(IP*NP)
  212. IF (IQ .GT. 0) THEN
  213. IW=1
  214. IF (IP .EQ. 0) IW=2
  215. SUM=SUM+W(NPT+IQ)*PTSAUX(IW,IQ)
  216. END IF
  217. END IF
  218. 160 W(K)=HALF*SUM*SUM
  219. C
  220. C Calculate VLAG and BETA for the required updating of the H matrix if
  221. C XPT(KNEW,.) is reinstated in the set of interpolation points.
  222. C
  223. DO 180 K=1,NPT
  224. SUM=ZERO
  225. DO 170 J=1,N
  226. 170 SUM=SUM+BMAT(K,J)*W(NPT+J)
  227. 180 VLAG(K)=SUM
  228. BETA=ZERO
  229. DO 200 J=1,NPTM
  230. SUM=ZERO
  231. DO 190 K=1,NPT
  232. 190 SUM=SUM+ZMAT(K,J)*W(K)
  233. BETA=BETA-SUM*SUM
  234. DO 200 K=1,NPT
  235. 200 VLAG(K)=VLAG(K)+SUM*ZMAT(K,J)
  236. BSUM=ZERO
  237. DISTSQ=ZERO
  238. DO 230 J=1,N
  239. SUM=ZERO
  240. DO 210 K=1,NPT
  241. 210 SUM=SUM+BMAT(K,J)*W(K)
  242. JP=J+NPT
  243. BSUM=BSUM+SUM*W(JP)
  244. DO 220 IP=NPT+1,NDIM
  245. 220 SUM=SUM+BMAT(IP,J)*W(IP)
  246. BSUM=BSUM+SUM*W(JP)
  247. VLAG(JP)=SUM
  248. 230 DISTSQ=DISTSQ+XPT(KNEW,J)**2
  249. BETA=HALF*DISTSQ*DISTSQ+BETA-BSUM
  250. VLAG(KOPT)=VLAG(KOPT)+ONE
  251. C
  252. C KOLD is set to the index of the provisional interpolation point that is
  253. C going to be deleted to make way for the KNEW-th original interpolation
  254. C point. The choice of KOLD is governed by the avoidance of a small value
  255. C of the denominator in the updating calculation of UPDATE.
  256. C
  257. DENOM=ZERO
  258. VLMXSQ=ZERO
  259. DO 250 K=1,NPT
  260. IF (PTSID(K) .NE. ZERO) THEN
  261. HDIAG=ZERO
  262. DO 240 J=1,NPTM
  263. 240 HDIAG=HDIAG+ZMAT(K,J)**2
  264. DEN=BETA*HDIAG+VLAG(K)**2
  265. IF (DEN .GT. DENOM) THEN
  266. KOLD=K
  267. DENOM=DEN
  268. END IF
  269. END IF
  270. 250 VLMXSQ=DMAX1(VLMXSQ,VLAG(K)**2)
  271. IF (DENOM .LE. 1.0D-2*VLMXSQ) THEN
  272. W(NDIM+KNEW)=-W(NDIM+KNEW)-WINC
  273. GOTO 120
  274. END IF
  275. GOTO 80
  276. C
  277. C When label 260 is reached, all the final positions of the interpolation
  278. C points have been chosen although any changes have not been included yet
  279. C in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart
  280. C from the shift of XBASE, the updating of the quadratic model remains to
  281. C be done. The following cycle through the new interpolation points begins
  282. C by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero,
  283. C except that a RETURN occurs if MAXFUN prohibits another value of F.
  284. C
  285. 260 DO 340 KPT=1,NPT
  286. IF (PTSID(KPT) .EQ. ZERO) GOTO 340
  287. IF (NF .GE. MAXFUN) THEN
  288. NF=-1
  289. GOTO 350
  290. END IF
  291. IH=0
  292. DO 270 J=1,N
  293. W(J)=XPT(KPT,J)
  294. XPT(KPT,J)=ZERO
  295. TEMP=PQ(KPT)*W(J)
  296. DO 270 I=1,J
  297. IH=IH+1
  298. 270 HQ(IH)=HQ(IH)+TEMP*W(I)
  299. PQ(KPT)=ZERO
  300. IP=PTSID(KPT)
  301. IQ=DFLOAT(NP)*PTSID(KPT)-DFLOAT(IP*NP)
  302. IF (IP .GT. 0) THEN
  303. XP=PTSAUX(1,IP)
  304. XPT(KPT,IP)=XP
  305. END IF
  306. IF (IQ .GT. 0) THEN
  307. XQ=PTSAUX(1,IQ)
  308. IF (IP .EQ. 0) XQ=PTSAUX(2,IQ)
  309. XPT(KPT,IQ)=XQ
  310. END IF
  311. C
  312. C Set VQUAD to the value of the current model at the new point.
  313. C
  314. VQUAD=FBASE
  315. IF (IP .GT. 0) THEN
  316. IHP=(IP+IP*IP)/2
  317. VQUAD=VQUAD+XP*(GOPT(IP)+HALF*XP*HQ(IHP))
  318. END IF
  319. IF (IQ .GT. 0) THEN
  320. IHQ=(IQ+IQ*IQ)/2
  321. VQUAD=VQUAD+XQ*(GOPT(IQ)+HALF*XQ*HQ(IHQ))
  322. IF (IP .GT. 0) THEN
  323. IW=MAX0(IHP,IHQ)-IABS(IP-IQ)
  324. VQUAD=VQUAD+XP*XQ*HQ(IW)
  325. END IF
  326. END IF
  327. DO 280 K=1,NPT
  328. TEMP=ZERO
  329. IF (IP .GT. 0) TEMP=TEMP+XP*XPT(K,IP)
  330. IF (IQ .GT. 0) TEMP=TEMP+XQ*XPT(K,IQ)
  331. 280 VQUAD=VQUAD+HALF*PQ(K)*TEMP*TEMP
  332. C
  333. C Calculate F at the new interpolation point, and set DIFF to the factor
  334. C that is going to multiply the KPT-th Lagrange function when the model
  335. C is updated to provide interpolation to the new function value.
  336. C
  337. DO 290 I=1,N
  338. W(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XPT(KPT,I)),XU(I))
  339. IF (XPT(KPT,I) .EQ. SL(I)) W(I)=XL(I)
  340. IF (XPT(KPT,I) .EQ. SU(I)) W(I)=XU(I)
  341. 290 CONTINUE
  342. NF=NF+1
  343. CALL CALFUN (N,W,F)
  344. IF (IPRINT .EQ. 3) THEN
  345. PRINT 300, NF,F,(W(I),I=1,N)
  346. 300 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
  347. 1 ' The corresponding X is:'/(2X,5D15.6))
  348. END IF
  349. FVAL(KPT)=F
  350. IF (F .LT. FVAL(KOPT)) KOPT=KPT
  351. DIFF=F-VQUAD
  352. C
  353. C Update the quadratic model. The RETURN from the subroutine occurs when
  354. C all the new interpolation points are included in the model.
  355. C
  356. DO 310 I=1,N
  357. 310 GOPT(I)=GOPT(I)+DIFF*BMAT(KPT,I)
  358. DO 330 K=1,NPT
  359. SUM=ZERO
  360. DO 320 J=1,NPTM
  361. 320 SUM=SUM+ZMAT(K,J)*ZMAT(KPT,J)
  362. TEMP=DIFF*SUM
  363. IF (PTSID(K) .EQ. ZERO) THEN
  364. PQ(K)=PQ(K)+TEMP
  365. ELSE
  366. IP=PTSID(K)
  367. IQ=DFLOAT(NP)*PTSID(K)-DFLOAT(IP*NP)
  368. IHQ=(IQ*IQ+IQ)/2
  369. IF (IP .EQ. 0) THEN
  370. HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(2,IQ)**2
  371. ELSE
  372. IHP=(IP*IP+IP)/2
  373. HQ(IHP)=HQ(IHP)+TEMP*PTSAUX(1,IP)**2
  374. IF (IQ .GT. 0) THEN
  375. HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(1,IQ)**2
  376. IW=MAX0(IHP,IHQ)-IABS(IQ-IP)
  377. HQ(IW)=HQ(IW)+TEMP*PTSAUX(1,IP)*PTSAUX(1,IQ)
  378. END IF
  379. END IF
  380. END IF
  381. 330 CONTINUE
  382. PTSID(KPT)=ZERO
  383. 340 CONTINUE
  384. 350 RETURN
  385. END