bobyqb.f 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670
  1. SUBROUTINE BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,
  2. 1 MAXFUN,XBASE,XPT,FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,
  3. 2 SL,SU,XNEW,XALT,D,VLAG,W)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. DIMENSION X(*),XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),
  6. 1 XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),
  7. 2 SL(*),SU(*),XNEW(*),XALT(*),D(*),VLAG(*),W(*)
  8. C
  9. C The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
  10. C are identical to the corresponding arguments in SUBROUTINE BOBYQA.
  11. C XBASE holds a shift of origin that should reduce the contributions
  12. C from rounding errors to values of the model and Lagrange functions.
  13. C XPT is a two-dimensional array that holds the coordinates of the
  14. C interpolation points relative to XBASE.
  15. C FVAL holds the values of F at the interpolation points.
  16. C XOPT is set to the displacement from XBASE of the trust region centre.
  17. C GOPT holds the gradient of the quadratic model at XBASE+XOPT.
  18. C HQ holds the explicit second derivatives of the quadratic model.
  19. C PQ contains the parameters of the implicit second derivatives of the
  20. C quadratic model.
  21. C BMAT holds the last N columns of H.
  22. C ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
  23. C this factorization being ZMAT times ZMAT^T, which provides both the
  24. C correct rank and positive semi-definiteness.
  25. C NDIM is the first dimension of BMAT and has the value NPT+N.
  26. C SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
  27. C All the components of every XOPT are going to satisfy the bounds
  28. C SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
  29. C XOPT is on a constraint boundary.
  30. C XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the
  31. C vector of variables for the next call of CALFUN. XNEW also satisfies
  32. C the SL and SU constraints in the way that has just been mentioned.
  33. C XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
  34. C in order to increase the denominator in the updating of UPDATE.
  35. C D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
  36. C VLAG contains the values of the Lagrange functions at a new point X.
  37. C They are part of a product that requires VLAG to be of length NDIM.
  38. C W is a one-dimensional array that is used for working space. Its length
  39. C must be at least 3*NDIM = 3*(NPT+N).
  40. C
  41. C Set some constants.
  42. C
  43. HALF=0.5D0
  44. ONE=1.0D0
  45. TEN=10.0D0
  46. TENTH=0.1D0
  47. TWO=2.0D0
  48. ZERO=0.0D0
  49. NP=N+1
  50. NPTM=NPT-NP
  51. NH=(N*NP)/2
  52. C
  53. C The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
  54. C BMAT and ZMAT for the first iteration, with the corresponding values of
  55. C of NF and KOPT, which are the number of calls of CALFUN so far and the
  56. C index of the interpolation point at the trust region centre. Then the
  57. C initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
  58. C less than NPT. GOPT will be updated if KOPT is different from KBASE.
  59. C
  60. CALL PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,XPT,
  61. 1 FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT)
  62. XOPTSQ=ZERO
  63. DO 10 I=1,N
  64. XOPT(I)=XPT(KOPT,I)
  65. 10 XOPTSQ=XOPTSQ+XOPT(I)**2
  66. FSAVE=FVAL(1)
  67. IF (NF .LT. NPT) THEN
  68. IF (IPRINT .GT. 0) PRINT 390
  69. GOTO 720
  70. END IF
  71. KBASE=1
  72. C
  73. C Complete the settings that are required for the iterative procedure.
  74. C
  75. RHO=RHOBEG
  76. DELTA=RHO
  77. NRESC=NF
  78. NTRITS=0
  79. DIFFA=ZERO
  80. DIFFB=ZERO
  81. ITEST=0
  82. NFSAV=NF
  83. C
  84. C Update GOPT if necessary before the first iteration and after each
  85. C call of RESCUE that makes a call of CALFUN.
  86. C
  87. 20 IF (KOPT .NE. KBASE) THEN
  88. IH=0
  89. DO 30 J=1,N
  90. DO 30 I=1,J
  91. IH=IH+1
  92. IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*XOPT(I)
  93. 30 GOPT(I)=GOPT(I)+HQ(IH)*XOPT(J)
  94. IF (NF .GT. NPT) THEN
  95. DO 50 K=1,NPT
  96. TEMP=ZERO
  97. DO 40 J=1,N
  98. 40 TEMP=TEMP+XPT(K,J)*XOPT(J)
  99. TEMP=PQ(K)*TEMP
  100. DO 50 I=1,N
  101. 50 GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
  102. END IF
  103. END IF
  104. C
  105. C Generate the next point in the trust region that provides a small value
  106. C of the quadratic model subject to the constraints on the variables.
  107. C The integer NTRITS is set to the number "trust region" iterations that
  108. C have occurred since the last "alternative" iteration. If the length
  109. C of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
  110. C label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
  111. C
  112. 60 CALL TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,XNEW,D,
  113. 1 W,W(NP),W(NP+N),W(NP+2*N),W(NP+3*N),DSQ,CRVMIN)
  114. DNORM=DMIN1(DELTA,DSQRT(DSQ))
  115. IF (DNORM .LT. HALF*RHO) THEN
  116. NTRITS=-1
  117. DISTSQ=(TEN*RHO)**2
  118. IF (NF .LE. NFSAV+2) GOTO 650
  119. C
  120. C The following choice between labels 650 and 680 depends on whether or
  121. C not our work with the current RHO seems to be complete. Either RHO is
  122. C decreased or termination occurs if the errors in the quadratic model at
  123. C the last three interpolation points compare favourably with predictions
  124. C of likely improvements to the model within distance HALF*RHO of XOPT.
  125. C
  126. ERRBIG=DMAX1(DIFFA,DIFFB,DIFFC)
  127. FRHOSQ=0.125D0*RHO*RHO
  128. IF (CRVMIN .GT. ZERO .AND. ERRBIG .GT. FRHOSQ*CRVMIN)
  129. 1 GOTO 650
  130. BDTOL=ERRBIG/RHO
  131. DO 80 J=1,N
  132. BDTEST=BDTOL
  133. IF (XNEW(J) .EQ. SL(J)) BDTEST=W(J)
  134. IF (XNEW(J) .EQ. SU(J)) BDTEST=-W(J)
  135. IF (BDTEST .LT. BDTOL) THEN
  136. CURV=HQ((J+J*J)/2)
  137. DO 70 K=1,NPT
  138. 70 CURV=CURV+PQ(K)*XPT(K,J)**2
  139. BDTEST=BDTEST+HALF*CURV*RHO
  140. IF (BDTEST .LT. BDTOL) GOTO 650
  141. END IF
  142. 80 CONTINUE
  143. GOTO 680
  144. END IF
  145. NTRITS=NTRITS+1
  146. C
  147. C Severe cancellation is likely to occur if XOPT is too far from XBASE.
  148. C If the following test holds, then XBASE is shifted so that XOPT becomes
  149. C zero. The appropriate changes are made to BMAT and to the second
  150. C derivatives of the current model, beginning with the changes to BMAT
  151. C that do not depend on ZMAT. VLAG is used temporarily for working space.
  152. C
  153. 90 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN
  154. FRACSQ=0.25D0*XOPTSQ
  155. SUMPQ=ZERO
  156. DO 110 K=1,NPT
  157. SUMPQ=SUMPQ+PQ(K)
  158. SUM=-HALF*XOPTSQ
  159. DO 100 I=1,N
  160. 100 SUM=SUM+XPT(K,I)*XOPT(I)
  161. W(NPT+K)=SUM
  162. TEMP=FRACSQ-HALF*SUM
  163. DO 110 I=1,N
  164. W(I)=BMAT(K,I)
  165. VLAG(I)=SUM*XPT(K,I)+TEMP*XOPT(I)
  166. IP=NPT+I
  167. DO 110 J=1,I
  168. 110 BMAT(IP,J)=BMAT(IP,J)+W(I)*VLAG(J)+VLAG(I)*W(J)
  169. C
  170. C Then the revisions of BMAT that depend on ZMAT are calculated.
  171. C
  172. DO 150 JJ=1,NPTM
  173. SUMZ=ZERO
  174. SUMW=ZERO
  175. DO 120 K=1,NPT
  176. SUMZ=SUMZ+ZMAT(K,JJ)
  177. VLAG(K)=W(NPT+K)*ZMAT(K,JJ)
  178. 120 SUMW=SUMW+VLAG(K)
  179. DO 140 J=1,N
  180. SUM=(FRACSQ*SUMZ-HALF*SUMW)*XOPT(J)
  181. DO 130 K=1,NPT
  182. 130 SUM=SUM+VLAG(K)*XPT(K,J)
  183. W(J)=SUM
  184. DO 140 K=1,NPT
  185. 140 BMAT(K,J)=BMAT(K,J)+SUM*ZMAT(K,JJ)
  186. DO 150 I=1,N
  187. IP=I+NPT
  188. TEMP=W(I)
  189. DO 150 J=1,I
  190. 150 BMAT(IP,J)=BMAT(IP,J)+TEMP*W(J)
  191. C
  192. C The following instructions complete the shift, including the changes
  193. C to the second derivative parameters of the quadratic model.
  194. C
  195. IH=0
  196. DO 170 J=1,N
  197. W(J)=-HALF*SUMPQ*XOPT(J)
  198. DO 160 K=1,NPT
  199. W(J)=W(J)+PQ(K)*XPT(K,J)
  200. 160 XPT(K,J)=XPT(K,J)-XOPT(J)
  201. DO 170 I=1,J
  202. IH=IH+1
  203. HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
  204. 170 BMAT(NPT+I,J)=BMAT(NPT+J,I)
  205. DO 180 I=1,N
  206. XBASE(I)=XBASE(I)+XOPT(I)
  207. XNEW(I)=XNEW(I)-XOPT(I)
  208. SL(I)=SL(I)-XOPT(I)
  209. SU(I)=SU(I)-XOPT(I)
  210. 180 XOPT(I)=ZERO
  211. XOPTSQ=ZERO
  212. END IF
  213. IF (NTRITS .EQ. 0) GOTO 210
  214. GOTO 230
  215. C
  216. C XBASE is also moved to XOPT by a call of RESCUE. This calculation is
  217. C more expensive than the previous shift, because new matrices BMAT and
  218. C ZMAT are generated from scratch, which may include the replacement of
  219. C interpolation points whose positions seem to be causing near linear
  220. C dependence in the interpolation conditions. Therefore RESCUE is called
  221. C only if rounding errors have reduced by at least a factor of two the
  222. C denominator of the formula for updating the H matrix. It provides a
  223. C useful safeguard, but is not invoked in most applications of BOBYQA.
  224. C
  225. 190 NFSAV=NF
  226. KBASE=KOPT
  227. CALL RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,FVAL,
  228. 1 XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,KOPT,
  229. 2 VLAG,W,W(N+NP),W(NDIM+NP))
  230. C
  231. C XOPT is updated now in case the branch below to label 720 is taken.
  232. C Any updating of GOPT occurs after the branch below to label 20, which
  233. C leads to a trust region iteration as does the branch to label 60.
  234. C
  235. XOPTSQ=ZERO
  236. IF (KOPT .NE. KBASE) THEN
  237. DO 200 I=1,N
  238. XOPT(I)=XPT(KOPT,I)
  239. 200 XOPTSQ=XOPTSQ+XOPT(I)**2
  240. END IF
  241. IF (NF .LT. 0) THEN
  242. NF=MAXFUN
  243. IF (IPRINT .GT. 0) PRINT 390
  244. GOTO 720
  245. END IF
  246. NRESC=NF
  247. IF (NFSAV .LT. NF) THEN
  248. NFSAV=NF
  249. GOTO 20
  250. END IF
  251. IF (NTRITS .GT. 0) GOTO 60
  252. C
  253. C Pick two alternative vectors of variables, relative to XBASE, that
  254. C are suitable as new positions of the KNEW-th interpolation point.
  255. C Firstly, XNEW is set to the point on a line through XOPT and another
  256. C interpolation point that minimizes the predicted value of the next
  257. C denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
  258. C and SU bounds. Secondly, XALT is set to the best feasible point on
  259. C a constrained version of the Cauchy step of the KNEW-th Lagrange
  260. C function, the corresponding value of the square of this function
  261. C being returned in CAUCHY. The choice between these alternatives is
  262. C going to be made when the denominator is calculated.
  263. C
  264. 210 CALL ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
  265. 1 KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,W,W(NP),W(NDIM+1))
  266. DO 220 I=1,N
  267. 220 D(I)=XNEW(I)-XOPT(I)
  268. C
  269. C Calculate VLAG and BETA for the current choice of D. The scalar
  270. C product of D with XPT(K,.) is going to be held in W(NPT+K) for
  271. C use when VQUAD is calculated.
  272. C
  273. 230 DO 250 K=1,NPT
  274. SUMA=ZERO
  275. SUMB=ZERO
  276. SUM=ZERO
  277. DO 240 J=1,N
  278. SUMA=SUMA+XPT(K,J)*D(J)
  279. SUMB=SUMB+XPT(K,J)*XOPT(J)
  280. 240 SUM=SUM+BMAT(K,J)*D(J)
  281. W(K)=SUMA*(HALF*SUMA+SUMB)
  282. VLAG(K)=SUM
  283. 250 W(NPT+K)=SUMA
  284. BETA=ZERO
  285. DO 270 JJ=1,NPTM
  286. SUM=ZERO
  287. DO 260 K=1,NPT
  288. 260 SUM=SUM+ZMAT(K,JJ)*W(K)
  289. BETA=BETA-SUM*SUM
  290. DO 270 K=1,NPT
  291. 270 VLAG(K)=VLAG(K)+SUM*ZMAT(K,JJ)
  292. DSQ=ZERO
  293. BSUM=ZERO
  294. DX=ZERO
  295. DO 300 J=1,N
  296. DSQ=DSQ+D(J)**2
  297. SUM=ZERO
  298. DO 280 K=1,NPT
  299. 280 SUM=SUM+W(K)*BMAT(K,J)
  300. BSUM=BSUM+SUM*D(J)
  301. JP=NPT+J
  302. DO 290 I=1,N
  303. 290 SUM=SUM+BMAT(JP,I)*D(I)
  304. VLAG(JP)=SUM
  305. BSUM=BSUM+SUM*D(J)
  306. 300 DX=DX+D(J)*XOPT(J)
  307. BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
  308. VLAG(KOPT)=VLAG(KOPT)+ONE
  309. C
  310. C If NTRITS is zero, the denominator may be increased by replacing
  311. C the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
  312. C rounding errors have damaged the chosen denominator.
  313. C
  314. IF (NTRITS .EQ. 0) THEN
  315. DENOM=VLAG(KNEW)**2+ALPHA*BETA
  316. IF (DENOM .LT. CAUCHY .AND. CAUCHY .GT. ZERO) THEN
  317. DO 310 I=1,N
  318. XNEW(I)=XALT(I)
  319. 310 D(I)=XNEW(I)-XOPT(I)
  320. CAUCHY=ZERO
  321. GO TO 230
  322. END IF
  323. IF (DENOM .LE. HALF*VLAG(KNEW)**2) THEN
  324. IF (NF .GT. NRESC) GOTO 190
  325. IF (IPRINT .GT. 0) PRINT 320
  326. 320 FORMAT (/5X,'Return from BOBYQA because of much',
  327. 1 ' cancellation in a denominator.')
  328. GOTO 720
  329. END IF
  330. C
  331. C Alternatively, if NTRITS is positive, then set KNEW to the index of
  332. C the next interpolation point to be deleted to make room for a trust
  333. C region step. Again RESCUE may be called if rounding errors have damaged
  334. C the chosen denominator, which is the reason for attempting to select
  335. C KNEW before calculating the next value of the objective function.
  336. C
  337. ELSE
  338. DELSQ=DELTA*DELTA
  339. SCADEN=ZERO
  340. BIGLSQ=ZERO
  341. KNEW=0
  342. DO 350 K=1,NPT
  343. IF (K .EQ. KOPT) GOTO 350
  344. HDIAG=ZERO
  345. DO 330 JJ=1,NPTM
  346. 330 HDIAG=HDIAG+ZMAT(K,JJ)**2
  347. DEN=BETA*HDIAG+VLAG(K)**2
  348. DISTSQ=ZERO
  349. DO 340 J=1,N
  350. 340 DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
  351. TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
  352. IF (TEMP*DEN .GT. SCADEN) THEN
  353. SCADEN=TEMP*DEN
  354. KNEW=K
  355. DENOM=DEN
  356. END IF
  357. BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
  358. 350 CONTINUE
  359. IF (SCADEN .LE. HALF*BIGLSQ) THEN
  360. IF (NF .GT. NRESC) GOTO 190
  361. IF (IPRINT .GT. 0) PRINT 320
  362. GOTO 720
  363. END IF
  364. END IF
  365. C
  366. C Put the variables for the next calculation of the objective function
  367. C in XNEW, with any adjustments for the bounds.
  368. C
  369. C
  370. C Calculate the value of the objective function at XBASE+XNEW, unless
  371. C the limit on the number of calculations of F has been reached.
  372. C
  373. 360 DO 380 I=1,N
  374. X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XNEW(I)),XU(I))
  375. IF (XNEW(I) .EQ. SL(I)) X(I)=XL(I)
  376. IF (XNEW(I) .EQ. SU(I)) X(I)=XU(I)
  377. 380 CONTINUE
  378. IF (NF .GE. MAXFUN) THEN
  379. IF (IPRINT .GT. 0) PRINT 390
  380. 390 FORMAT (/4X,'Return from BOBYQA because CALFUN has been',
  381. 1 ' called MAXFUN times.')
  382. GOTO 720
  383. END IF
  384. NF=NF+1
  385. CALL CALFUN (N,X,F)
  386. IF (IPRINT .EQ. 3) THEN
  387. PRINT 400, NF,F,(X(I),I=1,N)
  388. 400 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
  389. 1 ' The corresponding X is:'/(2X,5D15.6))
  390. END IF
  391. IF (NTRITS .EQ. -1) THEN
  392. FSAVE=F
  393. GOTO 720
  394. END IF
  395. C
  396. C Use the quadratic model to predict the change in F due to the step D,
  397. C and set DIFF to the error of this prediction.
  398. C
  399. FOPT=FVAL(KOPT)
  400. VQUAD=ZERO
  401. IH=0
  402. DO 410 J=1,N
  403. VQUAD=VQUAD+D(J)*GOPT(J)
  404. DO 410 I=1,J
  405. IH=IH+1
  406. TEMP=D(I)*D(J)
  407. IF (I .EQ. J) TEMP=HALF*TEMP
  408. 410 VQUAD=VQUAD+HQ(IH)*TEMP
  409. DO 420 K=1,NPT
  410. 420 VQUAD=VQUAD+HALF*PQ(K)*W(NPT+K)**2
  411. DIFF=F-FOPT-VQUAD
  412. DIFFC=DIFFB
  413. DIFFB=DIFFA
  414. DIFFA=DABS(DIFF)
  415. IF (DNORM .GT. RHO) NFSAV=NF
  416. C
  417. C Pick the next value of DELTA after a trust region step.
  418. C
  419. IF (NTRITS .GT. 0) THEN
  420. IF (VQUAD .GE. ZERO) THEN
  421. IF (IPRINT .GT. 0) PRINT 430
  422. 430 FORMAT (/4X,'Return from BOBYQA because a trust',
  423. 1 ' region step has failed to reduce Q.')
  424. GOTO 720
  425. END IF
  426. RATIO=(F-FOPT)/VQUAD
  427. IF (RATIO .LE. TENTH) THEN
  428. DELTA=DMIN1(HALF*DELTA,DNORM)
  429. ELSE IF (RATIO .LE. 0.7D0) THEN
  430. DELTA=DMAX1(HALF*DELTA,DNORM)
  431. ELSE
  432. DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
  433. END IF
  434. IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
  435. C
  436. C Recalculate KNEW and DENOM if the new F is less than FOPT.
  437. C
  438. IF (F .LT. FOPT) THEN
  439. KSAV=KNEW
  440. DENSAV=DENOM
  441. DELSQ=DELTA*DELTA
  442. SCADEN=ZERO
  443. BIGLSQ=ZERO
  444. KNEW=0
  445. DO 460 K=1,NPT
  446. HDIAG=ZERO
  447. DO 440 JJ=1,NPTM
  448. 440 HDIAG=HDIAG+ZMAT(K,JJ)**2
  449. DEN=BETA*HDIAG+VLAG(K)**2
  450. DISTSQ=ZERO
  451. DO 450 J=1,N
  452. 450 DISTSQ=DISTSQ+(XPT(K,J)-XNEW(J))**2
  453. TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
  454. IF (TEMP*DEN .GT. SCADEN) THEN
  455. SCADEN=TEMP*DEN
  456. KNEW=K
  457. DENOM=DEN
  458. END IF
  459. 460 BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
  460. IF (SCADEN .LE. HALF*BIGLSQ) THEN
  461. KNEW=KSAV
  462. DENOM=DENSAV
  463. END IF
  464. END IF
  465. END IF
  466. C
  467. C Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
  468. C moved. Also update the second derivative terms of the model.
  469. C
  470. CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
  471. IH=0
  472. PQOLD=PQ(KNEW)
  473. PQ(KNEW)=ZERO
  474. DO 470 I=1,N
  475. TEMP=PQOLD*XPT(KNEW,I)
  476. DO 470 J=1,I
  477. IH=IH+1
  478. 470 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
  479. DO 480 JJ=1,NPTM
  480. TEMP=DIFF*ZMAT(KNEW,JJ)
  481. DO 480 K=1,NPT
  482. 480 PQ(K)=PQ(K)+TEMP*ZMAT(K,JJ)
  483. C
  484. C Include the new interpolation point, and make the changes to GOPT at
  485. C the old XOPT that are caused by the updating of the quadratic model.
  486. C
  487. FVAL(KNEW)=F
  488. DO 490 I=1,N
  489. XPT(KNEW,I)=XNEW(I)
  490. 490 W(I)=BMAT(KNEW,I)
  491. DO 520 K=1,NPT
  492. SUMA=ZERO
  493. DO 500 JJ=1,NPTM
  494. 500 SUMA=SUMA+ZMAT(KNEW,JJ)*ZMAT(K,JJ)
  495. SUMB=ZERO
  496. DO 510 J=1,N
  497. 510 SUMB=SUMB+XPT(K,J)*XOPT(J)
  498. TEMP=SUMA*SUMB
  499. DO 520 I=1,N
  500. 520 W(I)=W(I)+TEMP*XPT(K,I)
  501. DO 530 I=1,N
  502. 530 GOPT(I)=GOPT(I)+DIFF*W(I)
  503. C
  504. C Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
  505. C
  506. IF (F .LT. FOPT) THEN
  507. KOPT=KNEW
  508. XOPTSQ=ZERO
  509. IH=0
  510. DO 540 J=1,N
  511. XOPT(J)=XNEW(J)
  512. XOPTSQ=XOPTSQ+XOPT(J)**2
  513. DO 540 I=1,J
  514. IH=IH+1
  515. IF (I .LT. J) GOPT(J)=GOPT(J)+HQ(IH)*D(I)
  516. 540 GOPT(I)=GOPT(I)+HQ(IH)*D(J)
  517. DO 560 K=1,NPT
  518. TEMP=ZERO
  519. DO 550 J=1,N
  520. 550 TEMP=TEMP+XPT(K,J)*D(J)
  521. TEMP=PQ(K)*TEMP
  522. DO 560 I=1,N
  523. 560 GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
  524. END IF
  525. C
  526. C Calculate the parameters of the least Frobenius norm interpolant to
  527. C the current data, the gradient of this interpolant at XOPT being put
  528. C into VLAG(NPT+I), I=1,2,...,N.
  529. C
  530. IF (NTRITS .GT. 0) THEN
  531. DO 570 K=1,NPT
  532. VLAG(K)=FVAL(K)-FVAL(KOPT)
  533. 570 W(K)=ZERO
  534. DO 590 J=1,NPTM
  535. SUM=ZERO
  536. DO 580 K=1,NPT
  537. 580 SUM=SUM+ZMAT(K,J)*VLAG(K)
  538. DO 590 K=1,NPT
  539. 590 W(K)=W(K)+SUM*ZMAT(K,J)
  540. DO 610 K=1,NPT
  541. SUM=ZERO
  542. DO 600 J=1,N
  543. 600 SUM=SUM+XPT(K,J)*XOPT(J)
  544. W(K+NPT)=W(K)
  545. 610 W(K)=SUM*W(K)
  546. GQSQ=ZERO
  547. GISQ=ZERO
  548. DO 630 I=1,N
  549. SUM=ZERO
  550. DO 620 K=1,NPT
  551. 620 SUM=SUM+BMAT(K,I)*VLAG(K)+XPT(K,I)*W(K)
  552. IF (XOPT(I) .EQ. SL(I)) THEN
  553. GQSQ=GQSQ+DMIN1(ZERO,GOPT(I))**2
  554. GISQ=GISQ+DMIN1(ZERO,SUM)**2
  555. ELSE IF (XOPT(I) .EQ. SU(I)) THEN
  556. GQSQ=GQSQ+DMAX1(ZERO,GOPT(I))**2
  557. GISQ=GISQ+DMAX1(ZERO,SUM)**2
  558. ELSE
  559. GQSQ=GQSQ+GOPT(I)**2
  560. GISQ=GISQ+SUM*SUM
  561. END IF
  562. 630 VLAG(NPT+I)=SUM
  563. C
  564. C Test whether to replace the new quadratic model by the least Frobenius
  565. C norm interpolant, making the replacement if the test is satisfied.
  566. C
  567. ITEST=ITEST+1
  568. IF (GQSQ .LT. TEN*GISQ) ITEST=0
  569. IF (ITEST .GE. 3) THEN
  570. DO 640 I=1,MAX0(NPT,NH)
  571. IF (I .LE. N) GOPT(I)=VLAG(NPT+I)
  572. IF (I .LE. NPT) PQ(I)=W(NPT+I)
  573. IF (I .LE. NH) HQ(I)=ZERO
  574. ITEST=0
  575. 640 CONTINUE
  576. END IF
  577. END IF
  578. C
  579. C If a trust region step has provided a sufficient decrease in F, then
  580. C branch for another trust region calculation. The case NTRITS=0 occurs
  581. C when the new interpolation point was reached by an alternative step.
  582. C
  583. IF (NTRITS .EQ. 0) GOTO 60
  584. IF (F .LE. FOPT+TENTH*VQUAD) GOTO 60
  585. C
  586. C Alternatively, find out if the interpolation points are close enough
  587. C to the best point so far.
  588. C
  589. DISTSQ=DMAX1((TWO*DELTA)**2,(TEN*RHO)**2)
  590. 650 KNEW=0
  591. DO 670 K=1,NPT
  592. SUM=ZERO
  593. DO 660 J=1,N
  594. 660 SUM=SUM+(XPT(K,J)-XOPT(J))**2
  595. IF (SUM .GT. DISTSQ) THEN
  596. KNEW=K
  597. DISTSQ=SUM
  598. END IF
  599. 670 CONTINUE
  600. C
  601. C If KNEW is positive, then ALTMOV finds alternative new positions for
  602. C the KNEW-th interpolation point within distance ADELT of XOPT. It is
  603. C reached via label 90. Otherwise, there is a branch to label 60 for
  604. C another trust region iteration, unless the calculations with the
  605. C current RHO are complete.
  606. C
  607. IF (KNEW .GT. 0) THEN
  608. DIST=DSQRT(DISTSQ)
  609. IF (NTRITS .EQ. -1) THEN
  610. DELTA=DMIN1(TENTH*DELTA,HALF*DIST)
  611. IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
  612. END IF
  613. NTRITS=0
  614. ADELT=DMAX1(DMIN1(TENTH*DIST,DELTA),RHO)
  615. DSQ=ADELT*ADELT
  616. GOTO 90
  617. END IF
  618. IF (NTRITS .EQ. -1) GOTO 680
  619. IF (RATIO .GT. ZERO) GOTO 60
  620. IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 60
  621. C
  622. C The calculations with the current value of RHO are complete. Pick the
  623. C next values of RHO and DELTA.
  624. C
  625. 680 IF (RHO .GT. RHOEND) THEN
  626. DELTA=HALF*RHO
  627. RATIO=RHO/RHOEND
  628. IF (RATIO .LE. 16.0D0) THEN
  629. RHO=RHOEND
  630. ELSE IF (RATIO .LE. 250.0D0) THEN
  631. RHO=DSQRT(RATIO)*RHOEND
  632. ELSE
  633. RHO=TENTH*RHO
  634. END IF
  635. DELTA=DMAX1(DELTA,RHO)
  636. IF (IPRINT .GE. 2) THEN
  637. IF (IPRINT .GE. 3) PRINT 690
  638. 690 FORMAT (5X)
  639. PRINT 700, RHO,NF
  640. 700 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of',
  641. 1 ' function values =',I6)
  642. PRINT 710, FVAL(KOPT),(XBASE(I)+XOPT(I),I=1,N)
  643. 710 FORMAT (4X,'Least value of F =',1PD23.15,9X,
  644. 1 'The corresponding X is:'/(2X,5D15.6))
  645. END IF
  646. NTRITS=0
  647. NFSAV=NF
  648. GOTO 60
  649. END IF
  650. C
  651. C Return from the calculation, after another Newton-Raphson step, if
  652. C it is too short to have been tried before.
  653. C
  654. IF (NTRITS .EQ. -1) GOTO 360
  655. 720 IF (FVAL(KOPT) .LE. FSAVE) THEN
  656. DO 730 I=1,N
  657. X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XOPT(I)),XU(I))
  658. IF (XOPT(I) .EQ. SL(I)) X(I)=XL(I)
  659. IF (XOPT(I) .EQ. SU(I)) X(I)=XU(I)
  660. 730 CONTINUE
  661. F=FVAL(KOPT)
  662. END IF
  663. IF (IPRINT .GE. 1) THEN
  664. PRINT 740, NF
  665. 740 FORMAT (/4X,'At the return from BOBYQA',5X,
  666. 1 'Number of function values =',I6)
  667. PRINT 710, F,(X(I),I=1,N)
  668. END IF
  669. RETURN
  670. END