Line data Source code
1 : module mod_bobyqa
2 : use const_def, only: dp
3 :
4 : contains
5 :
6 3 : subroutine do_BOBYQA (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN,W,CALFUN,max_valid_value)
7 : implicit REAL(dp) (A-H,O-Z)
8 : integer :: N, NPT, IPRINT, MAXFUN
9 : dimension X(:),XL(:),XU(:),W(*)
10 : interface
11 : #include "num_bobyqa_proc.dek"
12 : end interface
13 : real(dp), intent(in) :: max_valid_value
14 :
15 : ! This subroutine seeks the least value of a function of many variables,
16 : ! by applying a trust region method that forms quadratic models by
17 : ! interpolation. There is usually some freedom in the interpolation
18 : ! conditions, which is taken up by minimizing the Frobenius norm of
19 : ! the change to the second derivative of the model, beginning with the
20 : ! zero matrix. The values of the variables are constrained by upper and
21 : ! lower bounds. The arguments of the subroutine are as follows.
22 : !
23 : ! N must be set to the number of variables and must be at least two.
24 : ! NPT is the number of interpolation conditions. Its value must be in
25 : ! the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not
26 : ! recommended.
27 : ! Initial values of the variables must be set in X(1),X(2),...,X(N). They
28 : ! will be changed to the values that give the least calculated F.
29 : ! For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper
30 : ! bounds, respectively, on X(I). The construction of quadratic models
31 : ! requires XL(I) to be strictly less than XU(I) for each I. Further,
32 : ! the contribution to a model from changes to the I-th variable is
33 : ! damaged severely by rounding errors if XU(I)-XL(I) is too small.
34 : ! RHOBEG and RHOEND must be set to the initial and final values of a trust
35 : ! region radius, so both must be positive with RHOEND no greater than
36 : ! RHOBEG. Typically, RHOBEG should be about one tenth of the greatest
37 : ! expected change to a variable, while RHOEND should indicate the
38 : ! accuracy that is required in the final values of the variables. An
39 : ! error return occurs if any of the differences XU(I)-XL(I), I=1,...,N,
40 : ! is less than 2*RHOBEG.
41 : ! The value of IPRINT should be set to 0, 1, 2 or 3, which controls the
42 : ! amount of printing. Specifically, there is no output if IPRINT=0 and
43 : ! there is output only at the return if IPRINT=1. Otherwise, each new
44 : ! value of RHO is printed, with the best vector of variables so far and
45 : ! the corresponding value of the objective function. Further, each new
46 : ! value of F with its variables are output if IPRINT=3.
47 : ! MAXFUN must be set to an upper bound on the number of calls of CALFUN.
48 : ! The array W will be used for working space. Its length must be at least
49 : ! (NPT+5)*(NPT+N)+3*N*(N+5)/2.
50 : !
51 : ! subroutine CALFUN (N,X,F) has to be provided by the user. It must set
52 : ! F to the value of the objective function for the current values of the
53 : ! variables X(1),X(2),...,X(N), which are generated automatically in a
54 : ! way that satisfies the bounds given in XL and XU.
55 : !
56 : ! Return if the value of NPT is unacceptable.
57 : !
58 3 : NP=N+1
59 3 : if (NPT < N+2 .or. NPT > ((N+2)*NP)/2) then
60 0 : PRINT 10
61 : 10 FORMAT (/4X,'Return from BOBYQA because NPT is not in the required interval')
62 0 : GOTO 40
63 : end if
64 : !
65 : ! Partition the working space array, so that different parts of it can
66 : ! be treated separately during the calculation of BOBYQB. The partition
67 : ! requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the
68 : ! space that is taken by the last array in the argument list of BOBYQB.
69 : !
70 3 : NDIM=NPT+N
71 3 : IXB=1
72 3 : IXP=IXB+N
73 3 : IFV=IXP+N*NPT
74 3 : IXO=IFV+NPT
75 3 : IGO=IXO+N
76 3 : IHQ=IGO+N
77 3 : IPQ=IHQ+(N*NP)/2
78 3 : IBMAT=IPQ+NPT
79 3 : IZMAT=IBMAT+NDIM*N
80 3 : ISL=IZMAT+NPT*(NPT-NP)
81 3 : ISU=ISL+N
82 3 : IXN=ISU+N
83 3 : IXA=IXN+N
84 3 : ID=IXA+N
85 3 : IVL=ID+N
86 3 : IW=IVL+NDIM
87 : !
88 : ! Return if there is insufficient space between the bounds. Modify the
89 : ! initial X if necessary in order to avoid conflicts between the bounds
90 : ! and the construction of the first quadratic model. The lower and upper
91 : ! bounds on moves from the updated X are set now, in the ISL and ISU
92 : ! partitions of W, in order to provide useful and exact information about
93 : ! components of X that become within distance RHOBEG from their bounds.
94 : !
95 6 : ZERO=0.0D0
96 15 : do J=1,N
97 3 : TEMP=XU(J)-XL(J)
98 12 : if (TEMP < RHOBEG+RHOBEG) then
99 0 : PRINT 20
100 : 20 FORMAT (/4X,'Return from BOBYQA because one of the differences XU(I)-XL(I)'/6X,' is less than 2*RHOBEG.')
101 0 : GOTO 40
102 : end if
103 12 : JSL=ISL+J-1
104 12 : JSU=JSL+N
105 12 : W(JSL)=XL(J)-X(J)
106 12 : W(JSU)=XU(J)-X(J)
107 15 : if (W(JSL) >= -RHOBEG) then
108 0 : if (W(JSL) >= ZERO) then
109 0 : X(J)=XL(J)
110 0 : W(JSL)=ZERO
111 0 : W(JSU)=TEMP
112 : else
113 0 : X(J)=XL(J)+RHOBEG
114 0 : W(JSL)=-RHOBEG
115 0 : W(JSU)=DMAX1(XU(J)-X(J),RHOBEG)
116 : end if
117 12 : else if (W(JSU) <= RHOBEG) then
118 0 : if (W(JSU) <= ZERO) then
119 0 : X(J)=XU(J)
120 0 : W(JSL)=-TEMP
121 0 : W(JSU)=ZERO
122 : else
123 0 : X(J)=XU(J)-RHOBEG
124 0 : W(JSL)=DMIN1(XL(J)-X(J),-RHOBEG)
125 0 : W(JSU)=RHOBEG
126 : end if
127 : end if
128 : end do
129 : !
130 : ! Make the call of BOBYQB.
131 : !
132 : CALL BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB),
133 : & W(IXP),W(IFV),W(IXO),W(IGO),W(IHQ),W(IPQ),W(IBMAT),W(IZMAT),
134 : & NDIM,W(ISL),W(ISU),W(IXN),W(IXA),W(ID),W(IVL),W(IW),CALFUN,
135 3 : & max_valid_value)
136 3 : 40 RETURN
137 : end subroutine do_BOBYQA
138 :
139 :
140 3 : subroutine BOBYQB (N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,
141 3 : & MAXFUN,XBASE,XPT,FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,
142 : & SL,SU,XNEW,XALT,D,VLAG,W,CALFUN,max_valid_value)
143 : implicit real(dp) (A-H,O-Z)
144 : integer :: N, NPT, IPRINT, MAXFUN, NDIM
145 : dimension X(:),XL(:),XU(:),XBASE(*),XPT(NPT,*),FVAL(*),
146 : & XOPT(*),GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),
147 : & SL(*),SU(*),XNEW(*),XALT(*),D(*),VLAG(*),W(*)
148 : interface
149 : #include "num_bobyqa_proc.dek"
150 : end interface
151 : real(dp), intent(in) :: max_valid_value
152 :
153 : logical :: do_replace
154 : !
155 : ! The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN
156 : ! are identical to the corresponding arguments in subroutine BOBYQA.
157 : ! XBASE holds a shift of origin that should reduce the contributions
158 : ! from rounding errors to values of the model and Lagrange functions.
159 : ! XPT is a two-dimensional array that holds the coordinates of the
160 : ! interpolation points relative to XBASE.
161 : ! FVAL holds the values of F at the interpolation points.
162 : ! XOPT is set to the displacement from XBASE of the trust region centre.
163 : ! GOPT holds the gradient of the quadratic model at XBASE+XOPT.
164 : ! HQ holds the explicit second derivatives of the quadratic model.
165 : ! PQ contains the parameters of the implicit second derivatives of the
166 : ! quadratic model.
167 : ! BMAT holds the last N columns of H.
168 : ! ZMAT holds the factorization of the leading NPT by NPT submatrix of H,
169 : ! this factorization being ZMAT times ZMAT^T, which provides both the
170 : ! correct rank and positive semi-definiteness.
171 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
172 : ! SL and SU hold the differences XL-XBASE and XU-XBASE, respectively.
173 : ! All the components of every XOPT are going to satisfy the bounds
174 : ! SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when
175 : ! XOPT is on a constraint boundary.
176 : ! XNEW is chosen by subroutine TRSBOX or ALTMOV. Usually XBASE+XNEW is the
177 : ! vector of variables for the next call of CALFUN. XNEW also satisfies
178 : ! the SL and SU constraints in the way that has just been mentioned.
179 : ! XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW
180 : ! in order to increase the denominator in the updating of UPDATE.
181 : ! D is reserved for a trial step from XOPT, which is usually XNEW-XOPT.
182 : ! VLAG contains the values of the Lagrange functions at a new point X.
183 : ! They are part of a product that requires VLAG to be of length NDIM.
184 : ! W is a one-dimensional array that is used for working space. Its length
185 : ! must be at least 3*NDIM = 3*(NPT+N).
186 : !
187 : ! Set some constants.
188 : !
189 3 : HALF=0.5D0
190 3 : ONE=1.0D0
191 3 : TEN=10.0D0
192 3 : TENTH=0.1D0
193 3 : TWO=2.0D0
194 3 : ZERO=0.0D0
195 3 : NP=N+1
196 3 : NPTM=NPT-NP
197 3 : NH=(N*NP)/2
198 : !
199 : ! The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
200 : ! BMAT and ZMAT for the first iteration, with the corresponding values of
201 : ! of NF and KOPT, which are the number of calls of CALFUN so far and the
202 : ! index of the interpolation point at the trust region centre. Then the
203 : ! initial XOPT is set too. The branch to label 720 occurs if MAXFUN is
204 : ! less than NPT. GOPT will be updated if KOPT is different from KBASE.
205 : !
206 : CALL PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,XPT,
207 3 : & FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT,CALFUN)
208 6 : XOPTSQ=ZERO
209 15 : do I=1,N
210 12 : XOPT(I)=XPT(KOPT,I)
211 15 : XOPTSQ=XOPTSQ+XOPT(I)**2
212 : end do
213 6 : FSAVE=FVAL(1)
214 3 : if (NF < NPT) then
215 0 : if (IPRINT > 0) PRINT 390
216 : GOTO 720
217 : end if
218 3 : KBASE=1
219 : !
220 : ! Complete the settings that are required for the iterative procedure.
221 : !
222 6 : RHO=RHOBEG
223 6 : DELTA=RHO
224 3 : NRESC=NF
225 3 : NTRITS=0
226 6 : DIFFA=ZERO
227 6 : DIFFB=ZERO
228 3 : ITEST=0
229 3 : NFSAV=NF
230 : !
231 : ! Update GOPT if necessary before the first iteration and after each
232 : ! call of RESCUE that makes a call of CALFUN.
233 : !
234 3 : 20 if (KOPT /= KBASE) then
235 3 : IH=0
236 15 : do J=1,N
237 49 : do I=1,J
238 34 : IH=IH+1
239 34 : if (I < J) GOPT(J)=GOPT(J)+HQ(IH)*XOPT(I)
240 46 : GOPT(I)=GOPT(I)+HQ(IH)*XOPT(J)
241 : end do
242 : end do
243 3 : if (NF > NPT) then
244 0 : do K=1,NPT
245 3 : TEMP=ZERO
246 0 : do J=1,N
247 0 : TEMP=TEMP+XPT(K,J)*XOPT(J)
248 : end do
249 0 : TEMP=PQ(K)*TEMP
250 0 : do I=1,N
251 0 : GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
252 : end do
253 : end do
254 : end if
255 : end if
256 : !
257 : ! Generate the next point in the trust region that provides a small value
258 : ! of the quadratic model subject to the constraints on the variables.
259 : ! The integer NTRITS is set to the number "trust region" iterations that
260 : ! have occurred since the last "alternative" iteration. If the length
261 : ! of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to
262 : ! label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW.
263 : !
264 : 60 CALL TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,XNEW,D,
265 227 : & W,W(NP),W(NP+N),W(NP+2*N),W(NP+3*N),DSQ,CRVMIN)
266 230 : DNORM=DMIN1(DELTA,DSQRT(DSQ))
267 227 : if (DNORM < HALF*RHO) then
268 46 : NTRITS=-1
269 49 : DISTSQ=(TEN*RHO)**2
270 46 : if (NF <= NFSAV+2) GOTO 650
271 : !
272 : ! The following choice between labels 650 and 680 depends on whether or
273 : ! not our work with the current RHO seems to be complete. Either RHO is
274 : ! decreased or termination occurs if the errors in the quadratic model at
275 : ! the last three interpolation points compare favourably with predictions
276 : ! of likely improvements to the model within distance HALF*RHO of XOPT.
277 : !
278 11 : ERRBIG=DMAX1(DIFFA,DIFFB,DIFFC)
279 11 : FRHOSQ=0.125D0*RHO*RHO
280 8 : if (CRVMIN > ZERO .and. ERRBIG > FRHOSQ*CRVMIN) GOTO 650
281 3 : BDTOL=ERRBIG/RHO
282 0 : do J=1,N
283 3 : BDTEST=BDTOL
284 0 : if (XNEW(J) == SL(J)) BDTEST=W(J)
285 0 : if (XNEW(J) == SU(J)) BDTEST=-W(J)
286 0 : if (BDTEST < BDTOL) then
287 3 : CURV=HQ((J+J*J)/2)
288 0 : do K=1,NPT
289 0 : CURV=CURV+PQ(K)*XPT(K,J)**2
290 : end do
291 0 : BDTEST=BDTEST+HALF*CURV*RHO
292 0 : if (BDTEST < BDTOL) GOTO 650
293 : end if
294 : end do
295 : GOTO 680
296 : end if
297 181 : NTRITS=NTRITS+1
298 : !
299 : ! Severe cancellation is likely to occur if XOPT is too far from XBASE.
300 : ! If the following test holds, then XBASE is shifted so that XOPT becomes
301 : ! zero. The appropriate changes are made to BMAT and to the second
302 : ! derivatives of the current model, beginning with the changes to BMAT
303 : ! that do not depend on ZMAT. VLAG is used temporarily for working space.
304 : !
305 260 : 90 if (DSQ <= 1.0D-3*XOPTSQ) then
306 15 : FRACSQ=0.25D0*XOPTSQ
307 15 : SUMPQ=ZERO
308 136 : do K=1,NPT
309 124 : SUMPQ=SUMPQ+PQ(K)
310 127 : SUM=-HALF*XOPTSQ
311 756 : do I=1,N
312 756 : SUM=SUM+XPT(K,I)*XOPT(I)
313 : end do
314 124 : W(NPT+K)=SUM
315 124 : TEMP=FRACSQ-HALF*SUM
316 768 : do I=1,N
317 632 : W(I)=BMAT(K,I)
318 632 : VLAG(I)=SUM*XPT(K,I)+TEMP*XOPT(I)
319 632 : IP=NPT+I
320 2784 : do J=1,I
321 2660 : BMAT(IP,J)=BMAT(IP,J)+W(I)*VLAG(J)+VLAG(I)*W(J)
322 : end do
323 : end do
324 : end do
325 : !
326 : ! Then the revisions of BMAT that depend on ZMAT are calculated.
327 : !
328 68 : do JJ=1,NPTM
329 3 : SUMZ=ZERO
330 3 : SUMW=ZERO
331 688 : do K=1,NPT
332 632 : SUMZ=SUMZ+ZMAT(K,JJ)
333 632 : VLAG(K)=W(NPT+K)*ZMAT(K,JJ)
334 688 : SUMW=SUMW+VLAG(K)
335 : end do
336 344 : do J=1,N
337 288 : SUM=(FRACSQ*SUMZ-HALF*SUMW)*XOPT(J)
338 3712 : do K=1,NPT
339 3712 : SUM=SUM+VLAG(K)*XPT(K,J)
340 : end do
341 288 : W(J)=SUM
342 3768 : do K=1,NPT
343 3712 : BMAT(K,J)=BMAT(K,J)+SUM*ZMAT(K,JJ)
344 : end do
345 : end do
346 356 : do I=1,N
347 288 : IP=I+NPT
348 288 : TEMP=W(I)
349 1272 : do J=1,I
350 1216 : BMAT(IP,J)=BMAT(IP,J)+TEMP*W(J)
351 : end do
352 : end do
353 : end do
354 : !
355 : ! The following instructions complete the shift, including the changes
356 : ! to the second derivative parameters of the quadratic model.
357 : !
358 12 : IH=0
359 68 : do J=1,N
360 56 : W(J)=-HALF*SUMPQ*XOPT(J)
361 688 : do K=1,NPT
362 632 : W(J)=W(J)+PQ(K)*XPT(K,J)
363 688 : XPT(K,J)=XPT(K,J)-XOPT(J)
364 : end do
365 240 : do I=1,J
366 172 : IH=IH+1
367 172 : HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
368 228 : BMAT(NPT+I,J)=BMAT(NPT+J,I)
369 : end do
370 : end do
371 68 : do I=1,N
372 56 : XBASE(I)=XBASE(I)+XOPT(I)
373 56 : XNEW(I)=XNEW(I)-XOPT(I)
374 56 : SL(I)=SL(I)-XOPT(I)
375 56 : SU(I)=SU(I)-XOPT(I)
376 68 : XOPT(I)=ZERO
377 : end do
378 : XOPTSQ=ZERO
379 : end if
380 260 : if (NTRITS == 0) GOTO 210
381 0 : GOTO 230
382 : !
383 : ! XBASE is also moved to XOPT by a call of RESCUE. This calculation is
384 : ! more expensive than the previous shift, because new matrices BMAT and
385 : ! ZMAT are generated from scratch, which may include the replacement of
386 : ! interpolation points whose positions seem to be causing near linear
387 : ! dependence in the interpolation conditions. Therefore RESCUE is called
388 : ! only if rounding errors have reduced by at least a factor of two the
389 : ! denominator of the formula for updating the H matrix. It provides a
390 : ! useful safeguard, but is not invoked in most applications of BOBYQA.
391 : !
392 0 : 190 NFSAV=NF
393 0 : KBASE=KOPT
394 : CALL RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,FVAL,
395 : & XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,KOPT,
396 0 : & VLAG,W,W(N+NP),W(NDIM+NP),CALFUN)
397 : !
398 : ! XOPT is updated now in case the branch below to label 720 is taken.
399 : ! Any updating of GOPT occurs after the branch below to label 20, which
400 : ! leads to a trust region iteration as does the branch to label 60.
401 : !
402 0 : XOPTSQ=ZERO
403 0 : if (KOPT /= KBASE) then
404 0 : do I=1,N
405 0 : XOPT(I)=XPT(KOPT,I)
406 0 : XOPTSQ=XOPTSQ+XOPT(I)**2
407 : end do
408 : end if
409 0 : if (NF < 0) then
410 0 : NF=MAXFUN
411 0 : if (IPRINT > 0) PRINT 390
412 : GOTO 720
413 : end if
414 0 : NRESC=NF
415 0 : if (NFSAV < NF) then
416 : NFSAV=NF
417 : GOTO 20
418 : end if
419 0 : if (NTRITS > 0) GOTO 60
420 : !
421 : ! Pick two alternative vectors of variables, relative to XBASE, that
422 : ! are suitable as new positions of the KNEW-th interpolation point.
423 : ! Firstly, XNEW is set to the point on a line through XOPT and another
424 : ! interpolation point that minimizes the predicted value of the next
425 : ! denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL
426 : ! and SU bounds. Secondly, XALT is set to the best feasible point on
427 : ! a constrained version of the Cauchy step of the KNEW-th Lagrange
428 : ! function, the corresponding value of the square of this function
429 : ! being returned in CAUCHY. The choice between these alternatives is
430 : ! going to be made when the denominator is calculated.
431 : !
432 : 210 CALL ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
433 79 : & KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,W,W(NP),W(NDIM+1))
434 463 : do I=1,N
435 463 : D(I)=XNEW(I)-XOPT(I)
436 : end do
437 : !
438 : ! Calculate VLAG and BETA for the current choice of D. The scalar
439 : ! product of D with XPT(K,.) is going to be held in W(NPT+K) for
440 : ! use when VQUAD is calculated.
441 : !
442 3596 : 230 do K=1,NPT
443 3301 : SUMA=ZERO
444 3301 : SUMB=ZERO
445 3298 : SUM=ZERO
446 20958 : do J=1,N
447 17660 : SUMA=SUMA+XPT(K,J)*D(J)
448 17660 : SUMB=SUMB+XPT(K,J)*XOPT(J)
449 20958 : SUM=SUM+BMAT(K,J)*D(J)
450 : end do
451 3298 : W(K)=SUMA*(HALF*SUMA+SUMB)
452 3298 : VLAG(K)=SUM
453 3596 : W(NPT+K)=SUMA
454 : end do
455 301 : BETA=ZERO
456 1798 : do JJ=1,NPTM
457 : SUM=ZERO
458 19160 : do K=1,NPT
459 19160 : SUM=SUM+ZMAT(K,JJ)*W(K)
460 : end do
461 1500 : BETA=BETA-SUM*SUM
462 19458 : do K=1,NPT
463 19160 : VLAG(K)=VLAG(K)+SUM*ZMAT(K,JJ)
464 : end do
465 : end do
466 298 : DSQ=ZERO
467 301 : BSUM=ZERO
468 301 : DX=ZERO
469 1798 : do J=1,N
470 1500 : DSQ=DSQ+D(J)**2
471 1500 : SUM=ZERO
472 19160 : do K=1,NPT
473 19160 : SUM=SUM+W(K)*BMAT(K,J)
474 : end do
475 1500 : BSUM=BSUM+SUM*D(J)
476 1500 : JP=NPT+J
477 9580 : do I=1,N
478 9580 : SUM=SUM+BMAT(JP,I)*D(I)
479 : end do
480 1500 : VLAG(JP)=SUM
481 1500 : BSUM=BSUM+SUM*D(J)
482 1798 : DX=DX+D(J)*XOPT(J)
483 : end do
484 298 : BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
485 298 : VLAG(KOPT)=VLAG(KOPT)+ONE
486 : !
487 : ! If NTRITS is zero, the denominator may be increased by replacing
488 : ! the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if
489 : ! rounding errors have damaged the chosen denominator.
490 : !
491 298 : if (NTRITS == 0) then
492 120 : DENOM=VLAG(KNEW)**2+ALPHA*BETA
493 117 : if (DENOM < CAUCHY .and. CAUCHY > ZERO) then
494 236 : do I=1,N
495 198 : XNEW(I)=XALT(I)
496 236 : D(I)=XNEW(I)-XOPT(I)
497 : end do
498 38 : CAUCHY=ZERO
499 38 : GOTO 230
500 : end if
501 79 : if (DENOM <= HALF*VLAG(KNEW)**2) then
502 0 : if (NF > NRESC) GOTO 190
503 0 : if (IPRINT > 0) PRINT 320
504 : 320 FORMAT (/5X,'Return from BOBYQA because of much cancellation in a denominator.')
505 : GOTO 720
506 : end if
507 : !
508 : ! Alternatively, if NTRITS is positive, then set KNEW to the index of
509 : ! the next interpolation point to be deleted to make room for a trust
510 : ! region step. Again RESCUE may be called if rounding errors have damaged
511 : ! the chosen denominator, which is the reason for attempting to select
512 : ! KNEW before calculating the next value of the objective function.
513 : !
514 : else
515 184 : DELSQ=DELTA*DELTA
516 184 : SCADEN=ZERO
517 184 : BIGLSQ=ZERO
518 181 : KNEW=0
519 2198 : do K=1,NPT
520 2017 : if (K == KOPT) CYCLE
521 3 : HDIAG=ZERO
522 11780 : do JJ=1,NPTM
523 11780 : HDIAG=HDIAG+ZMAT(K,JJ)**2
524 : end do
525 1839 : DEN=BETA*HDIAG+VLAG(K)**2
526 1836 : DISTSQ=ZERO
527 11780 : do J=1,N
528 11780 : DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
529 : end do
530 1836 : TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
531 1836 : if (TEMP*DEN > SCADEN) then
532 497 : SCADEN=TEMP*DEN
533 497 : KNEW=K
534 497 : DENOM=DEN
535 : end if
536 2198 : BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
537 : end do
538 181 : if (SCADEN <= HALF*BIGLSQ) then
539 0 : if (NF > NRESC) GOTO 190
540 0 : if (IPRINT > 0) PRINT 320
541 : GOTO 720
542 : end if
543 : end if
544 : !
545 : ! Put the variables for the next calculation of the objective function
546 : ! in XNEW, with any adjustments for the bounds.
547 : !
548 : !
549 : ! Calculate the value of the objective function at XBASE+XNEW, unless
550 : ! the limit on the number of calculations of F has been reached.
551 : !
552 1577 : 360 do I=1,N
553 1314 : X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XNEW(I)),XU(I))
554 1314 : if (XNEW(I) == SL(I)) X(I)=XL(I)
555 1577 : if (XNEW(I) == SU(I)) X(I)=XU(I)
556 : end do
557 263 : if (NF >= MAXFUN) then
558 0 : if (IPRINT > 0) PRINT 390
559 : 390 FORMAT (/4X,'Return from BOBYQA because CALFUN has been called MAXFUN times.')
560 : GOTO 720
561 : end if
562 263 : NF=NF+1
563 263 : CALL CALFUN (N,X,F)
564 263 : if (IPRINT == 3) then
565 0 : PRINT 400, NF,F,(X(I),I=1,N)
566 : 400 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,
567 : & ' The corresponding X is:'/(2X,5D15.6))
568 : end if
569 263 : if (NTRITS == -1) then
570 3 : FSAVE=F
571 3 : GOTO 720
572 : end if
573 : !
574 : ! Use the quadratic model to predict the change in F due to the step D,
575 : ! and set DIFF to the error of this prediction.
576 : !
577 263 : FOPT=FVAL(KOPT)
578 263 : VQUAD=ZERO
579 260 : IH=0
580 1562 : do J=1,N
581 1302 : VQUAD=VQUAD+D(J)*GOPT(J)
582 5711 : do I=1,J
583 4149 : IH=IH+1
584 4149 : TEMP=D(I)*D(J)
585 4149 : if (I == J) TEMP=HALF*TEMP
586 5451 : VQUAD=VQUAD+HQ(IH)*TEMP
587 : end do
588 : end do
589 3124 : do K=1,NPT
590 3124 : VQUAD=VQUAD+HALF*PQ(K)*W(NPT+K)**2
591 : end do
592 263 : DIFF=F-FOPT-VQUAD
593 260 : DIFFC=DIFFB
594 260 : DIFFB=DIFFA
595 260 : DIFFA=DABS(DIFF)
596 260 : if (DNORM > RHO) NFSAV=NF
597 : !
598 : ! Pick the next value of DELTA after a trust region step.
599 : !
600 260 : if (NTRITS > 0) then
601 181 : if (VQUAD >= ZERO) then
602 0 : if (IPRINT > 0) PRINT 430
603 : 430 FORMAT (/4X,'Return from BOBYQA because a trust',
604 : & ' region step has failed to reduce Q.')
605 : GOTO 720
606 : end if
607 184 : RATIO=(F-FOPT)/VQUAD
608 181 : if (RATIO <= TENTH) then
609 54 : DELTA=DMIN1(HALF*DELTA,DNORM)
610 127 : else if (RATIO. LE. 0.7D0) then
611 31 : DELTA=DMAX1(HALF*DELTA,DNORM)
612 : else
613 96 : DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
614 : end if
615 181 : if (DELTA <= 1.5D0*RHO) DELTA=RHO
616 : !
617 : ! Recalculate KNEW and DENOM if the new F is less than FOPT.
618 : !
619 181 : if (F < FOPT) then
620 133 : KSAV=KNEW
621 136 : DENSAV=DENOM
622 133 : DELSQ=DELTA*DELTA
623 133 : SCADEN=ZERO
624 133 : BIGLSQ=ZERO
625 133 : KNEW=0
626 1598 : do K=1,NPT
627 : HDIAG=ZERO
628 9275 : do JJ=1,NPTM
629 9275 : HDIAG=HDIAG+ZMAT(K,JJ)**2
630 : end do
631 1465 : DEN=BETA*HDIAG+VLAG(K)**2
632 1465 : DISTSQ=ZERO
633 9275 : do J=1,N
634 9275 : DISTSQ=DISTSQ+(XPT(K,J)-XNEW(J))**2
635 : end do
636 1465 : TEMP=DMAX1(ONE,(DISTSQ/DELSQ)**2)
637 1465 : if (TEMP*DEN > SCADEN) then
638 363 : SCADEN=TEMP*DEN
639 363 : KNEW=K
640 363 : DENOM=DEN
641 : end if
642 1598 : BIGLSQ=DMAX1(BIGLSQ,TEMP*VLAG(K)**2)
643 : end do
644 133 : if (SCADEN <= HALF*BIGLSQ) then
645 0 : KNEW=KSAV
646 0 : DENOM=DENSAV
647 : end if
648 : end if
649 : end if
650 : !
651 : ! Update BMAT and ZMAT, so that the KNEW-th interpolation point can be
652 : ! moved. Also update the second derivative terms of the model.
653 : !
654 260 : CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
655 260 : IH=0
656 263 : PQOLD=PQ(KNEW)
657 260 : PQ(KNEW)=ZERO
658 1562 : do I=1,N
659 1302 : TEMP=PQOLD*XPT(KNEW,I)
660 5711 : do J=1,I
661 4149 : IH=IH+1
662 5451 : HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
663 : end do
664 : end do
665 1562 : do JJ=1,NPTM
666 1302 : TEMP=DIFF*ZMAT(KNEW,JJ)
667 16856 : do K=1,NPT
668 16596 : PQ(K)=PQ(K)+TEMP*ZMAT(K,JJ)
669 : end do
670 : end do
671 : !
672 : ! Include the new interpolation point, and make the changes to GOPT at
673 : ! the old XOPT that are caused by the updating of the quadratic model.
674 : !
675 260 : FVAL(KNEW)=F
676 1562 : do I=1,N
677 1302 : XPT(KNEW,I)=XNEW(I)
678 1562 : W(I)=BMAT(KNEW,I)
679 : end do
680 3124 : do K=1,NPT
681 : SUMA=ZERO
682 18158 : do JJ=1,NPTM
683 18158 : SUMA=SUMA+ZMAT(KNEW,JJ)*ZMAT(K,JJ)
684 : end do
685 : SUMB=ZERO
686 18158 : do J=1,N
687 18158 : SUMB=SUMB+XPT(K,J)*XOPT(J)
688 : end do
689 2864 : TEMP=SUMA*SUMB
690 18418 : do I=1,N
691 18158 : W(I)=W(I)+TEMP*XPT(K,I)
692 : end do
693 : end do
694 1562 : do I=1,N
695 1562 : GOPT(I)=GOPT(I)+DIFF*W(I)
696 : end do
697 : !
698 : ! Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT.
699 : !
700 260 : if (F < FOPT) then
701 151 : KOPT=KNEW
702 151 : XOPTSQ=ZERO
703 151 : IH=0
704 919 : do J=1,N
705 768 : XOPT(J)=XNEW(J)
706 768 : XOPTSQ=XOPTSQ+XOPT(J)**2
707 3383 : do I=1,J
708 2464 : IH=IH+1
709 2464 : if (I < J) GOPT(J)=GOPT(J)+HQ(IH)*D(I)
710 3232 : GOPT(I)=GOPT(I)+HQ(IH)*D(J)
711 : end do
712 : end do
713 1838 : do K=1,NPT
714 : TEMP=ZERO
715 10775 : do J=1,N
716 10775 : TEMP=TEMP+XPT(K,J)*D(J)
717 : end do
718 1687 : TEMP=PQ(K)*TEMP
719 10926 : do I=1,N
720 10775 : GOPT(I)=GOPT(I)+TEMP*XPT(K,I)
721 : end do
722 : end do
723 : end if
724 : !
725 : ! Calculate the parameters of the least Frobenius norm interpolant to
726 : ! the current data, the gradient of this interpolant at XOPT being put
727 : ! into VLAG(NPT+I), I=1,2,...,N.
728 : !
729 260 : if (NTRITS > 0) then
730 2198 : do K=1,NPT
731 2017 : VLAG(K)=FVAL(K)-FVAL(KOPT)
732 2198 : W(K)=ZERO
733 : end do
734 1099 : do J=1,NPTM
735 : SUM=ZERO
736 11961 : do K=1,NPT
737 11780 : SUM=SUM+ZMAT(K,J)*VLAG(K)
738 : end do
739 : end do
740 2198 : do K=1,NPT
741 2198 : W(K)=W(K)+SUM*ZMAT(K,J)
742 : end do
743 2198 : do K=1,NPT
744 : SUM=ZERO
745 12879 : do J=1,N
746 12879 : SUM=SUM+XPT(K,J)*XOPT(J)
747 : end do
748 2017 : W(K+NPT)=W(K)
749 2198 : W(K)=SUM*W(K)
750 : end do
751 3 : GQSQ=ZERO
752 3 : GISQ=ZERO
753 1099 : do I=1,N
754 : SUM=ZERO
755 11780 : do K=1,NPT
756 11780 : SUM=SUM+BMAT(K,I)*VLAG(K)+XPT(K,I)*W(K)
757 : end do
758 918 : if (XOPT(I) == SL(I)) then
759 0 : GQSQ=GQSQ+DMIN1(ZERO,GOPT(I))**2
760 0 : GISQ=GISQ+DMIN1(ZERO,SUM)**2
761 918 : else if (XOPT(I) == SU(I)) then
762 0 : GQSQ=GQSQ+DMAX1(ZERO,GOPT(I))**2
763 0 : GISQ=GISQ+DMAX1(ZERO,SUM)**2
764 : else
765 918 : GQSQ=GQSQ+GOPT(I)**2
766 918 : GISQ=GISQ+SUM*SUM
767 : end if
768 1099 : VLAG(NPT+I)=SUM
769 : end do
770 : !
771 : ! Test whether to replace the new quadratic model by the least Frobenius
772 : ! norm interpolant, making the replacement if the test is satisfied.
773 : !
774 181 : ITEST=ITEST+1
775 181 : if (GQSQ < TEN*GISQ) ITEST=0
776 181 : do_replace = (ITEST >= 3)
777 0 : if (.not. do_replace) then ! check for "invalid" value
778 2198 : do k=1,npt
779 2198 : if (fval(k) > max_valid_value) then
780 : do_replace = .true.
781 : !write(*,*) 'bobyqa value > max valid', fval(k), max_valid_value
782 : exit
783 : end if
784 : end do
785 : end if
786 181 : if (do_replace) then
787 : !stop 'bobyqa: do_replace'
788 0 : do I=1,MAX0(NPT,NH)
789 0 : if (I <= N) GOPT(I)=VLAG(NPT+I)
790 0 : if (I <= NPT) PQ(I)=W(NPT+I)
791 0 : if (I <= NH) HQ(I)=ZERO
792 0 : ITEST=0
793 : end do
794 : end if
795 : end if
796 : !
797 : ! If a trust region step has provided a sufficient decrease in F, then
798 : ! branch for another trust region calculation. The case NTRITS=0 occurs
799 : ! when the new interpolation point was reached by an alternative step.
800 : !
801 260 : if (NTRITS == 0) GOTO 60
802 181 : if (F <= FOPT+TENTH*VQUAD) GOTO 60
803 : !
804 : ! Alternatively, find out if the interpolation points are close enough
805 : ! to the best point so far.
806 : !
807 100 : DISTSQ=DMAX1((TWO*DELTA)**2,(TEN*RHO)**2)
808 100 : 650 KNEW=0
809 1128 : do K=1,NPT
810 1028 : SUM=ZERO
811 6260 : do J=1,N
812 6260 : SUM=SUM+(XPT(K,J)-XOPT(J))**2
813 : end do
814 1128 : if (SUM > DISTSQ) then
815 162 : KNEW=K
816 162 : DISTSQ=SUM
817 : end if
818 : end do
819 : !
820 : ! If KNEW is positive, then ALTMOV finds alternative new positions for
821 : ! the KNEW-th interpolation point within distance ADELT of XOPT. It is
822 : ! reached via label 90. Otherwise, there is a branch to label 60 for
823 : ! another trust region iteration, unless the calculations with the
824 : ! current RHO are complete.
825 : !
826 100 : if (KNEW > 0) then
827 3 : DIST=DSQRT(DISTSQ)
828 79 : if (NTRITS == -1) then
829 32 : DELTA=DMIN1(TENTH*DELTA,HALF*DIST)
830 32 : if (DELTA <= 1.5D0*RHO) DELTA=RHO
831 : end if
832 79 : NTRITS=0
833 79 : ADELT=DMAX1(DMIN1(TENTH*DIST,DELTA),RHO)
834 79 : DSQ=ADELT*ADELT
835 79 : GOTO 90
836 : end if
837 21 : if (NTRITS == -1) GOTO 680
838 7 : if (RATIO > ZERO) GOTO 60
839 6 : if (DMAX1(DELTA,DNORM) > RHO) GOTO 60
840 : !
841 : ! The calculations with the current value of RHO are complete. Pick the
842 : ! next values of RHO and DELTA.
843 : !
844 18 : 680 if (RHO > RHOEND) then
845 15 : DELTA=HALF*RHO
846 15 : RATIO=RHO/RHOEND
847 15 : if (RATIO <= 16.0D0) then
848 3 : RHO=RHOEND
849 12 : else if (RATIO <= 250.0D0) then
850 3 : RHO=DSQRT(RATIO)*RHOEND
851 : else
852 9 : RHO=TENTH*RHO
853 : end if
854 15 : DELTA=DMAX1(DELTA,RHO)
855 15 : if (IPRINT >= 2) then
856 0 : if (IPRINT >= 3) PRINT 690
857 : 690 FORMAT (5X)
858 0 : PRINT 700, RHO,NF
859 : 700 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of',
860 : & ' function values =',I6)
861 0 : PRINT 710, FVAL(KOPT),(XBASE(I)+XOPT(I),I=1,N)
862 : 710 FORMAT (4X,'Least value of F =',1PD23.15,9X,
863 : & 'The corresponding X is:'/(2X,5D15.6))
864 : end if
865 15 : NTRITS=0
866 15 : NFSAV=NF
867 15 : GOTO 60
868 : end if
869 : !
870 : ! Return from the calculation, after another Newton-Raphson step, if
871 : ! it is too short to have been tried before.
872 : !
873 3 : if (NTRITS == -1) GOTO 360
874 3 : 720 if (FVAL(KOPT) <= FSAVE) then
875 0 : do I=1,N
876 0 : X(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XOPT(I)),XU(I))
877 0 : if (XOPT(I) == SL(I)) X(I)=XL(I)
878 0 : if (XOPT(I) == SU(I)) X(I)=XU(I)
879 : end do
880 0 : F=FVAL(KOPT)
881 : end if
882 3 : if (IPRINT >= 1) then
883 0 : PRINT 740, NF
884 : 740 FORMAT (/4X,'At the return from BOBYQA',5X,'Number of function values =',I6)
885 0 : PRINT 710, F,(X(I),I=1,N)
886 : end if
887 3 : RETURN
888 : end subroutine BOBYQB
889 :
890 :
891 79 : subroutine ALTMOV (N,NPT,XPT,XOPT,BMAT,ZMAT,NDIM,SL,SU,KOPT,
892 : & KNEW,ADELT,XNEW,XALT,ALPHA,CAUCHY,GLAG,HCOL,W)
893 : implicit real(dp) (A-H,O-Z)
894 : integer :: N, NPT, NDIM, KOPT, KNEW
895 : dimension XPT(NPT,*),XOPT(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),
896 : & SU(*),XNEW(*),XALT(*),GLAG(*),HCOL(*),W(*)
897 : !
898 : ! The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have
899 : ! the same meanings as the corresponding arguments of BOBYQB.
900 : ! KOPT is the index of the optimal interpolation point.
901 : ! KNEW is the index of the interpolation point that is going to be moved.
902 : ! ADELT is the current trust region bound.
903 : ! XNEW will be set to a suitable new position for the interpolation point
904 : ! XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region
905 : ! bounds and it should provide a large denominator in the next call of
906 : ! UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the
907 : ! straight lines through XOPT and another interpolation point.
908 : ! XALT also provides a large value of the modulus of the KNEW-th Lagrange
909 : ! function subject to the constraints that have been mentioned, its main
910 : ! difference from XNEW being that XALT-XOPT is a constrained version of
911 : ! the Cauchy step within the trust region. An exception is that XALT is
912 : ! not calculated if all components of GLAG (see below) are zero.
913 : ! ALPHA will be set to the KNEW-th diagonal element of the H matrix.
914 : ! CAUCHY will be set to the square of the KNEW-th Lagrange function at
915 : ! the step XALT-XOPT from XOPT for the vector XALT that is returned,
916 : ! except that CAUCHY is set to zero if XALT is not calculated.
917 : ! GLAG is a working space vector of length N for the gradient of the
918 : ! KNEW-th Lagrange function at XOPT.
919 : ! HCOL is a working space vector of length NPT for the second derivative
920 : ! coefficients of the KNEW-th Lagrange function.
921 : ! W is a working space vector of length 2N that is going to hold the
922 : ! constrained Cauchy step from XOPT of the Lagrange function, followed
923 : ! by the downhill version of XALT when the uphill step is calculated.
924 : !
925 : ! Set the first NPT components of W to the leading elements of the
926 : ! KNEW-th column of the H matrix.
927 : !
928 79 : HALF=0.5D0
929 79 : ONE=1.0D0
930 79 : ZERO=0.0D0
931 79 : CONST=ONE+DSQRT(2.0D0)
932 926 : do K=1,NPT
933 926 : HCOL(K)=ZERO
934 : end do
935 463 : do J=1,NPT-N-1
936 463 : TEMP=ZMAT(KNEW,J)
937 4895 : do K=1,NPT
938 4816 : HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J)
939 : end do
940 : end do
941 79 : ALPHA=HCOL(KNEW)
942 158 : HA=HALF*ALPHA
943 : !
944 : ! Calculate the gradient of the KNEW-th Lagrange function at XOPT.
945 : !
946 463 : do I=1,N
947 463 : GLAG(I)=BMAT(KNEW,I)
948 : end do
949 926 : do K=1,NPT
950 : TEMP=ZERO
951 5279 : do J=1,N
952 5279 : TEMP=TEMP+XPT(K,J)*XOPT(J)
953 : end do
954 847 : TEMP=HCOL(K)*TEMP
955 5358 : do I=1,N
956 5279 : GLAG(I)=GLAG(I)+TEMP*XPT(K,I)
957 : end do
958 : end do
959 : !
960 : ! Search for a large denominator along the straight lines through XOPT
961 : ! and another interpolation point. SLBD and SUBD will be lower and upper
962 : ! bounds on the step along each of these lines in turn. PREDSQ will be
963 : ! set to the square of the predicted denominator for each line. PRESAV
964 : ! will be set to the largest admissible value of PREDSQ that occurs.
965 : !
966 79 : PRESAV=ZERO
967 926 : do K=1,NPT
968 847 : if (K == KOPT) CYCLE
969 79 : DDERIV=ZERO
970 79 : DISTSQ=ZERO
971 4816 : do I=1,N
972 4048 : TEMP=XPT(K,I)-XOPT(I)
973 4048 : DDERIV=DDERIV+GLAG(I)*TEMP
974 4816 : DISTSQ=DISTSQ+TEMP*TEMP
975 : end do
976 847 : SUBD=ADELT/DSQRT(DISTSQ)
977 847 : SLBD=-SUBD
978 768 : ILBD=0
979 768 : IUBD=0
980 847 : SUMIN=DMIN1(ONE,SUBD)
981 : !
982 : ! Revise SLBD and SUBD if necessary because of the bounds in SL and SU.
983 : !
984 4816 : do I=1,N
985 4048 : TEMP=XPT(K,I)-XOPT(I)
986 4816 : if (TEMP > ZERO) then
987 1949 : if (SLBD*TEMP < SL(I)-XOPT(I)) then
988 0 : SLBD=(SL(I)-XOPT(I))/TEMP
989 0 : ILBD=-I
990 : end if
991 1949 : if (SUBD*TEMP > SU(I)-XOPT(I)) then
992 0 : SUBD=DMAX1(SUMIN,(SU(I)-XOPT(I))/TEMP)
993 0 : IUBD=I
994 : end if
995 2099 : else if (TEMP < ZERO) then
996 2099 : if (SLBD*TEMP > SU(I)-XOPT(I)) then
997 0 : SLBD=(SU(I)-XOPT(I))/TEMP
998 0 : ILBD=I
999 : end if
1000 2099 : if (SUBD*TEMP < SL(I)-XOPT(I)) then
1001 0 : SUBD=DMAX1(SUMIN,(SL(I)-XOPT(I))/TEMP)
1002 0 : IUBD=-I
1003 : end if
1004 : end if
1005 : end do
1006 : !
1007 : ! Seek a large modulus of the KNEW-th Lagrange function when the index
1008 : ! of the other interpolation point on the line through XOPT is KNEW.
1009 : !
1010 768 : if (K == KNEW) then
1011 158 : DIFF=DDERIV-ONE
1012 158 : STEP=SLBD
1013 158 : VLAG=SLBD*(DDERIV-SLBD*DIFF)
1014 79 : ISBD=ILBD
1015 79 : TEMP=SUBD*(DDERIV-SUBD*DIFF)
1016 79 : if (DABS(TEMP) > DABS(VLAG)) then
1017 37 : STEP=SUBD
1018 37 : VLAG=TEMP
1019 37 : ISBD=IUBD
1020 : end if
1021 158 : TEMPD=HALF*DDERIV
1022 158 : TEMPA=TEMPD-DIFF*SLBD
1023 158 : TEMPB=TEMPD-DIFF*SUBD
1024 79 : if (TEMPA*TEMPB < ZERO) then
1025 40 : TEMP=TEMPD*TEMPD/DIFF
1026 40 : if (DABS(TEMP) > DABS(VLAG)) then
1027 0 : STEP=TEMPD/DIFF
1028 0 : VLAG=TEMP
1029 0 : ISBD=0
1030 : end if
1031 : end if
1032 : !
1033 : ! Search along each of the other lines through XOPT and another point.
1034 : !
1035 : else
1036 689 : STEP=SLBD
1037 689 : VLAG=SLBD*(ONE-SLBD)
1038 689 : ISBD=ILBD
1039 689 : TEMP=SUBD*(ONE-SUBD)
1040 689 : if (DABS(TEMP) > DABS(VLAG)) then
1041 0 : STEP=SUBD
1042 0 : VLAG=TEMP
1043 0 : ISBD=IUBD
1044 : end if
1045 689 : if (SUBD > HALF) then
1046 126 : if (DABS(VLAG) < 0.25D0) then
1047 689 : STEP=HALF
1048 689 : VLAG=0.25D0
1049 689 : ISBD=0
1050 : end if
1051 : end if
1052 689 : VLAG=VLAG*DDERIV
1053 : end if
1054 : !
1055 : ! Calculate PREDSQ for the current line search and maintain PRESAV.
1056 : !
1057 768 : TEMP=STEP*(ONE-STEP)*DISTSQ
1058 847 : PREDSQ=VLAG*VLAG*(VLAG*VLAG+HA*TEMP*TEMP)
1059 847 : if (PREDSQ > PRESAV) then
1060 221 : PRESAV=PREDSQ
1061 221 : KSAV=K
1062 300 : STPSAV=STEP
1063 221 : IBDSAV=ISBD
1064 : end if
1065 : end do
1066 : !
1067 : ! Construct XNEW in a way that satisfies the bound constraints exactly.
1068 : !
1069 463 : do I=1,N
1070 384 : TEMP=XOPT(I)+STPSAV*(XPT(KSAV,I)-XOPT(I))
1071 463 : XNEW(I)=DMAX1(SL(I),DMIN1(SU(I),TEMP))
1072 : end do
1073 79 : if (IBDSAV < 0) XNEW(-IBDSAV)=SL(-IBDSAV)
1074 79 : if (IBDSAV > 0) XNEW(IBDSAV)=SU(IBDSAV)
1075 : !
1076 : ! Prepare for the iterative method that assembles the constrained Cauchy
1077 : ! step in W. The sum of squares of the fixed components of W is formed in
1078 : ! WFIXSQ, and the free components of W are set to BIGSTP.
1079 : !
1080 158 : BIGSTP=ADELT+ADELT
1081 79 : IFLAG=0
1082 237 : 100 WFIXSQ=ZERO
1083 237 : GGFREE=ZERO
1084 926 : do I=1,N
1085 768 : W(I)=ZERO
1086 768 : TEMPA=DMIN1(XOPT(I)-SL(I),GLAG(I))
1087 768 : TEMPB=DMAX1(XOPT(I)-SU(I),GLAG(I))
1088 926 : if (TEMPA > ZERO .or. TEMPB < ZERO) then
1089 768 : W(I)=BIGSTP
1090 768 : GGFREE=GGFREE+GLAG(I)**2
1091 : end if
1092 : end do
1093 158 : if (GGFREE == ZERO) then
1094 0 : CAUCHY=ZERO
1095 0 : GOTO 200
1096 : end if
1097 : !
1098 : ! Investigate whether more components of W can be fixed.
1099 : !
1100 158 : 120 TEMP=ADELT*ADELT-WFIXSQ
1101 158 : if (TEMP > ZERO) then
1102 237 : WSQSAV=WFIXSQ
1103 158 : STEP=DSQRT(TEMP/GGFREE)
1104 158 : GGFREE=ZERO
1105 926 : do I=1,N
1106 926 : if (W(I) == BIGSTP) then
1107 768 : TEMP=XOPT(I)-STEP*GLAG(I)
1108 768 : if (TEMP <= SL(I)) then
1109 0 : W(I)=SL(I)-XOPT(I)
1110 0 : WFIXSQ=WFIXSQ+W(I)**2
1111 768 : else if (TEMP >= SU(I)) then
1112 0 : W(I)=SU(I)-XOPT(I)
1113 0 : WFIXSQ=WFIXSQ+W(I)**2
1114 : else
1115 768 : GGFREE=GGFREE+GLAG(I)**2
1116 : end if
1117 : end if
1118 : end do
1119 158 : if (WFIXSQ > WSQSAV .and. GGFREE > ZERO) GOTO 120
1120 : end if
1121 : !
1122 : ! Set the remaining free components of W and all components of XALT,
1123 : ! except that W may be scaled later.
1124 : !
1125 237 : GW=ZERO
1126 926 : do I=1,N
1127 768 : if (W(I) == BIGSTP) then
1128 768 : W(I)=-STEP*GLAG(I)
1129 768 : XALT(I)=DMAX1(SL(I),DMIN1(SU(I),XOPT(I)+W(I)))
1130 0 : else if (W(I) == ZERO) then
1131 0 : XALT(I)=XOPT(I)
1132 0 : else if (GLAG(I) > ZERO) then
1133 0 : XALT(I)=SL(I)
1134 : else
1135 0 : XALT(I)=SU(I)
1136 : end if
1137 926 : GW=GW+GLAG(I)*W(I)
1138 : end do
1139 : !
1140 : ! Set CURV to the curvature of the KNEW-th Lagrange function along W.
1141 : ! Scale W by a factor less than one if that can reduce the modulus of
1142 : ! the Lagrange function at XOPT+W. Set CAUCHY to the final value of
1143 : ! the square of this function.
1144 : !
1145 79 : CURV=ZERO
1146 1852 : do K=1,NPT
1147 : TEMP=ZERO
1148 10558 : do J=1,N
1149 10558 : TEMP=TEMP+XPT(K,J)*W(J)
1150 : end do
1151 1852 : CURV=CURV+HCOL(K)*TEMP*TEMP
1152 : end do
1153 158 : if (IFLAG == 1) CURV=-CURV
1154 158 : if (CURV > -GW .and. CURV < -CONST*GW) then
1155 85 : SCALE=-GW/CURV
1156 30 : do I=1,N
1157 24 : TEMP=XOPT(I)+SCALE*W(I)
1158 30 : XALT(I)=DMAX1(SL(I),DMIN1(SU(I),TEMP))
1159 : end do
1160 6 : CAUCHY=(HALF*GW*SCALE)**2
1161 : else
1162 152 : CAUCHY=(GW+HALF*CURV)**2
1163 : end if
1164 : !
1165 : ! If IFLAG is zero, then XALT is calculated as before after reversing
1166 : ! the sign of GLAG. Thus two XALT vectors become available. The one that
1167 : ! is chosen is the one that gives the larger value of CAUCHY.
1168 : !
1169 158 : if (IFLAG == 0) then
1170 463 : do I=1,N
1171 384 : GLAG(I)=-GLAG(I)
1172 463 : W(N+I)=XALT(I)
1173 : end do
1174 79 : CSAVE=CAUCHY
1175 79 : IFLAG=1
1176 79 : GOTO 100
1177 : end if
1178 79 : if (CSAVE > CAUCHY) then
1179 132 : do I=1,N
1180 132 : XALT(I)=W(N+I)
1181 : end do
1182 20 : CAUCHY=CSAVE
1183 : end if
1184 79 : 200 RETURN
1185 : end subroutine ALTMOV
1186 :
1187 :
1188 6 : subroutine PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE,
1189 3 : & XPT,FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT,CALFUN)
1190 : implicit real(dp) (A-H,O-Z)
1191 : dimension X(:),XL(:),XU(:),XBASE(*),XPT(NPT,*),FVAL(*),GOPT(*),
1192 : & HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*)
1193 : interface
1194 : #include "num_bobyqa_proc.dek"
1195 : end interface
1196 : !
1197 : ! The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the
1198 : ! same as the corresponding arguments in subroutine BOBYQA.
1199 : ! The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU
1200 : ! are the same as the corresponding arguments in BOBYQB, the elements
1201 : ! of SL and SU being set in BOBYQA.
1202 : ! GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but
1203 : ! it is set by PRELIM to the gradient of the quadratic model at XBASE.
1204 : ! If XOPT is nonzero, BOBYQB will change it to its usual value later.
1205 : ! NF is maintained as the number of calls of CALFUN so far.
1206 : ! KOPT will be such that the least calculated value of F so far is at
1207 : ! the point XPT(KOPT,.)+XBASE in the space of the variables.
1208 : !
1209 : ! subroutine PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ,
1210 : ! BMAT and ZMAT for the first iteration, and it maintains the values of
1211 : ! NF and KOPT. The vector X is also changed by PRELIM.
1212 : !
1213 : ! Set some constants.
1214 : !
1215 3 : HALF=0.5D0
1216 3 : ONE=1.0D0
1217 3 : TWO=2.0D0
1218 3 : ZERO=0.0D0
1219 3 : RHOSQ=RHOBEG*RHOBEG
1220 3 : RECIP=ONE/RHOSQ
1221 3 : NP=N+1
1222 : !
1223 : ! Set XBASE to the initial vector of variables, and set the initial
1224 : ! elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
1225 : !
1226 15 : do J=1,N
1227 12 : XBASE(J)=X(J)
1228 136 : do K=1,NPT
1229 136 : XPT(K,J)=ZERO
1230 : end do
1231 195 : do I=1,NDIM
1232 192 : BMAT(I,J)=ZERO
1233 : end do
1234 : end do
1235 37 : do IH=1,(N*NP)/2
1236 37 : HQ(IH)=ZERO
1237 : end do
1238 30 : do K=1,NPT
1239 27 : PQ(K)=ZERO
1240 154 : do J=1,NPT-NP
1241 151 : ZMAT(K,J)=ZERO
1242 : end do
1243 : end do
1244 : !
1245 : ! Begin the initialization procedure. NF becomes one more than the number
1246 : ! of function values so far. The coordinates of the displacement of the
1247 : ! next initial interpolation point from XBASE are set in XPT(NF+1,.).
1248 : !
1249 3 : NF=0
1250 27 : 50 NFM=NF
1251 27 : NFX=NF-N
1252 27 : NF=NF+1
1253 27 : if (NFM <= 2*N) then
1254 27 : if (NFM >= 1 .and. NFM <= N) then
1255 15 : STEPA=RHOBEG
1256 12 : if (SU(NFM) == ZERO) STEPA=-STEPA
1257 12 : XPT(NF,NFM)=STEPA
1258 15 : else if (NFM > N) then
1259 12 : STEPA=XPT(NF-N,NFX)
1260 15 : STEPB=-RHOBEG
1261 12 : if (SL(NFX) == ZERO) STEPB=DMIN1(TWO*RHOBEG,SU(NFX))
1262 12 : if (SU(NFX) == ZERO) STEPB=DMAX1(-TWO*RHOBEG,SL(NFX))
1263 12 : XPT(NF,NFX)=STEPB
1264 : end if
1265 : else
1266 0 : ITEMP=(NFM-NP)/N
1267 0 : JPT=NFM-ITEMP*N-N
1268 0 : IPT=JPT+ITEMP
1269 0 : if (IPT > N) then
1270 0 : ITEMP=JPT
1271 0 : JPT=IPT-N
1272 0 : IPT=ITEMP
1273 : end if
1274 0 : XPT(NF,IPT)=XPT(IPT+1,IPT)
1275 0 : XPT(NF,JPT)=XPT(JPT+1,JPT)
1276 : end if
1277 : !
1278 : ! Calculate the next value of F. The least function value so far and
1279 : ! its index are required.
1280 : !
1281 151 : do J=1,N
1282 124 : X(J)=DMIN1(DMAX1(XL(J),XBASE(J)+XPT(NF,J)),XU(J))
1283 124 : if (XPT(NF,J) == SL(J)) X(J)=XL(J)
1284 151 : if (XPT(NF,J) == SU(J)) X(J)=XU(J)
1285 : end do
1286 27 : CALL CALFUN (N,X,F)
1287 27 : if (IPRINT == 3) then
1288 0 : PRINT 70, NF,F,(X(I),I=1,N)
1289 : 70 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,' The corresponding X is:'/(2X,5D15.6))
1290 : end if
1291 27 : FVAL(NF)=F
1292 27 : if (NF == 1) then
1293 6 : FBEG=F
1294 3 : KOPT=1
1295 24 : else if (F < FVAL(KOPT)) then
1296 8 : KOPT=NF
1297 : end if
1298 : !
1299 : ! Set the nonzero initial elements of BMAT and the quadratic model in the
1300 : ! cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions
1301 : ! of the NF-th and (NF-N)-th interpolation points may be switched, in
1302 : ! order that the function value at the first of them contributes to the
1303 : ! off-diagonal second derivative terms of the initial quadratic model.
1304 : !
1305 27 : if (NF <= 2*N+1) then
1306 27 : if (NF >= 2 .and. NF <= N+1) then
1307 12 : GOPT(NFM)=(F-FBEG)/STEPA
1308 12 : if (NPT < NF+N) then
1309 0 : BMAT(1,NFM)=-ONE/STEPA
1310 0 : BMAT(NF,NFM)=ONE/STEPA
1311 0 : BMAT(NPT+NFM,NFM)=-HALF*RHOSQ
1312 : end if
1313 15 : else if (NF >= N+2) then
1314 12 : IH=(NFX*(NFX+1))/2
1315 15 : TEMP=(F-FBEG)/STEPB
1316 3 : DIFF=STEPB-STEPA
1317 12 : HQ(IH)=TWO*(TEMP-GOPT(NFX))/DIFF
1318 12 : GOPT(NFX)=(GOPT(NFX)*STEPB-TEMP*STEPA)/DIFF
1319 12 : if (STEPA*STEPB < ZERO) then
1320 12 : if (F < FVAL(NF-N)) then
1321 6 : FVAL(NF)=FVAL(NF-N)
1322 6 : FVAL(NF-N)=F
1323 6 : if (KOPT == NF) KOPT=NF-N
1324 6 : XPT(NF-N,NFX)=STEPB
1325 6 : XPT(NF,NFX)=STEPA
1326 : end if
1327 : end if
1328 12 : BMAT(1,NFX)=-(STEPA+STEPB)/(STEPA*STEPB)
1329 12 : BMAT(NF,NFX)=-HALF/XPT(NF-N,NFX)
1330 12 : BMAT(NF-N,NFX)=-BMAT(1,NFX)-BMAT(NF,NFX)
1331 12 : ZMAT(1,NFX)=DSQRT(TWO)/(STEPA*STEPB)
1332 12 : ZMAT(NF,NFX)=DSQRT(HALF)/RHOSQ
1333 12 : ZMAT(NF-N,NFX)=-ZMAT(1,NFX)-ZMAT(NF,NFX)
1334 : end if
1335 : !
1336 : ! Set the off-diagonal second derivatives of the Lagrange functions and
1337 : ! the initial quadratic model.
1338 : !
1339 : else
1340 0 : IH=(IPT*(IPT-1))/2+JPT
1341 0 : ZMAT(1,NFX)=RECIP
1342 0 : ZMAT(NF,NFX)=RECIP
1343 0 : ZMAT(IPT+1,NFX)=-RECIP
1344 0 : ZMAT(JPT+1,NFX)=-RECIP
1345 0 : TEMP=XPT(NF,IPT)*XPT(NF,JPT)
1346 0 : HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/TEMP
1347 : end if
1348 27 : if (NF < NPT .and. NF < MAXFUN) GOTO 50
1349 3 : RETURN
1350 : end subroutine PRELIM
1351 :
1352 :
1353 0 : subroutine RESCUE (N,NPT,XL,XU,IPRINT,MAXFUN,XBASE,XPT,
1354 0 : & FVAL,XOPT,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,DELTA,
1355 : & KOPT,VLAG,PTSAUX,PTSID,W,CALFUN)
1356 : implicit real(dp) (A-H,O-Z)
1357 : dimension XL(:),XU(:),XBASE(*),XPT(NPT,*),FVAL(*),XOPT(*),
1358 : & GOPT(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*),
1359 : & VLAG(*),PTSAUX(2,*),PTSID(*),W(*)
1360 : interface
1361 : #include "num_bobyqa_proc.dek"
1362 : end interface
1363 : !
1364 : ! The arguments N, NPT, XL, XU, IPRINT, MAXFUN, XBASE, XPT, FVAL, XOPT,
1365 : ! GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as
1366 : ! the corresponding arguments of BOBYQB on the entry to RESCUE.
1367 : ! NF is maintained as the number of calls of CALFUN so far, except that
1368 : ! NF is set to -1 if the value of MAXFUN prevents further progress.
1369 : ! KOPT is maintained so that FVAL(KOPT) is the least calculated function
1370 : ! value. Its correct value must be given on entry. It is updated if a
1371 : ! new least function value is found, but the corresponding changes to
1372 : ! XOPT and GOPT have to be made later by the calling program.
1373 : ! DELTA is the current trust region radius.
1374 : ! VLAG is a working space vector that will be used for the values of the
1375 : ! provisional Lagrange functions at each of the interpolation points.
1376 : ! They are part of a product that requires VLAG to be of length NDIM.
1377 : ! PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and
1378 : ! PTSAUX(2,J) specify the two positions of provisional interpolation
1379 : ! points when a nonzero step is taken along e_J (the J-th coordinate
1380 : ! direction) through XBASE+XOPT, as specified below. Usually these
1381 : ! steps have length DELTA, but other lengths are chosen if necessary
1382 : ! in order to satisfy the given bounds on the variables.
1383 : ! PTSID is also a working space array. It has NPT components that denote
1384 : ! provisional new positions of the original interpolation points, in
1385 : ! case changes are needed to restore the linear independence of the
1386 : ! interpolation conditions. The K-th point is a candidate for change
1387 : ! if and only if PTSID(K) is nonzero. In this case let p and q be the
1388 : ! integer parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p
1389 : ! and q are both positive, the step from XBASE+XOPT to the new K-th
1390 : ! interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise
1391 : ! the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or
1392 : ! p=0, respectively.
1393 : ! The first NDIM+NPT elements of the array W are used for working space.
1394 : ! The final elements of BMAT and ZMAT are set in a well-conditioned way
1395 : ! to the values that are appropriate for the new interpolation points.
1396 : ! The elements of GOPT, HQ and PQ are also revised to the values that are
1397 : ! appropriate to the final quadratic model.
1398 : !
1399 : ! Set some constants.
1400 : !
1401 0 : HALF=0.5D0
1402 0 : ONE=1.0D0
1403 0 : ZERO=0.0D0
1404 0 : NP=N+1
1405 0 : SFRAC=HALF/DBLE(NP)
1406 0 : NPTM=NPT-NP
1407 : !
1408 : ! Shift the interpolation points so that XOPT becomes the origin, and set
1409 : ! the elements of ZMAT to zero. The value of SUMPQ is required in the
1410 : ! updating of HQ below. The squares of the distances from XOPT to the
1411 : ! other interpolation points are set at the end of W. Increments of WINC
1412 : ! may be added later to these squares to balance the consideration of
1413 : ! the choice of point that is going to become current.
1414 : !
1415 0 : SUMPQ=ZERO
1416 0 : WINC=ZERO
1417 0 : do K=1,NPT
1418 0 : DISTSQ=ZERO
1419 0 : do J=1,N
1420 0 : XPT(K,J)=XPT(K,J)-XOPT(J)
1421 0 : DISTSQ=DISTSQ+XPT(K,J)**2
1422 : end do
1423 0 : SUMPQ=SUMPQ+PQ(K)
1424 0 : W(NDIM+K)=DISTSQ
1425 0 : WINC=DMAX1(WINC,DISTSQ)
1426 0 : do J=1,NPTM
1427 0 : ZMAT(K,J)=ZERO
1428 : end do
1429 : end do
1430 : !
1431 : ! Update HQ so that HQ and PQ define the second derivatives of the model
1432 : ! after XBASE has been shifted to the trust region centre.
1433 : !
1434 : IH=0
1435 0 : do J=1,N
1436 0 : W(J)=HALF*SUMPQ*XOPT(J)
1437 0 : do K=1,NPT
1438 0 : W(J)=W(J)+PQ(K)*XPT(K,J)
1439 : end do
1440 0 : do I=1,J
1441 0 : IH=IH+1
1442 0 : HQ(IH)=HQ(IH)+W(I)*XOPT(J)+W(J)*XOPT(I)
1443 : end do
1444 : end do
1445 : !
1446 : ! Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and
1447 : ! also set the elements of PTSAUX.
1448 : !
1449 0 : do J=1,N
1450 0 : XBASE(J)=XBASE(J)+XOPT(J)
1451 0 : SL(J)=SL(J)-XOPT(J)
1452 0 : SU(J)=SU(J)-XOPT(J)
1453 0 : XOPT(J)=ZERO
1454 0 : PTSAUX(1,J)=DMIN1(DELTA,SU(J))
1455 0 : PTSAUX(2,J)=DMAX1(-DELTA,SL(J))
1456 0 : if (PTSAUX(1,J)+PTSAUX(2,J) < ZERO) then
1457 0 : TEMP=PTSAUX(1,J)
1458 0 : PTSAUX(1,J)=PTSAUX(2,J)
1459 0 : PTSAUX(2,J)=TEMP
1460 : end if
1461 0 : if (DABS(PTSAUX(2,J)) < HALF*DABS(PTSAUX(1,J))) then
1462 0 : PTSAUX(2,J)=HALF*PTSAUX(1,J)
1463 : end if
1464 0 : do I=1,NDIM
1465 0 : BMAT(I,J)=ZERO
1466 : end do
1467 : end do
1468 0 : FBASE=FVAL(KOPT)
1469 : !
1470 : ! Set the identifiers of the artificial interpolation points that are
1471 : ! along a coordinate direction from XOPT, and set the corresponding
1472 : ! nonzero elements of BMAT and ZMAT.
1473 : !
1474 0 : PTSID(1)=SFRAC
1475 0 : do J=1,N
1476 0 : JP=J+1
1477 0 : JPN=JP+N
1478 0 : PTSID(JP)=DBLE(J)+SFRAC
1479 0 : if (JPN <= NPT) then
1480 0 : PTSID(JPN)=DBLE(J)/DBLE(NP)+SFRAC
1481 0 : TEMP=ONE/(PTSAUX(1,J)-PTSAUX(2,J))
1482 0 : BMAT(JP,J)=-TEMP+ONE/PTSAUX(1,J)
1483 0 : BMAT(JPN,J)=TEMP+ONE/PTSAUX(2,J)
1484 0 : BMAT(1,J)=-BMAT(JP,J)-BMAT(JPN,J)
1485 0 : ZMAT(1,J)=DSQRT(2.0D0)/DABS(PTSAUX(1,J)*PTSAUX(2,J))
1486 0 : ZMAT(JP,J)=ZMAT(1,J)*PTSAUX(2,J)*TEMP
1487 0 : ZMAT(JPN,J)=-ZMAT(1,J)*PTSAUX(1,J)*TEMP
1488 : else
1489 0 : BMAT(1,J)=-ONE/PTSAUX(1,J)
1490 0 : BMAT(JP,J)=ONE/PTSAUX(1,J)
1491 0 : BMAT(J+NPT,J)=-HALF*PTSAUX(1,J)**2
1492 : end if
1493 : end do
1494 : !
1495 : ! Set any remaining identifiers with their nonzero elements of ZMAT.
1496 : !
1497 0 : if (NPT >= N+NP) then
1498 0 : do K=2*NP,NPT
1499 0 : IW=(DBLE(K-NP)-HALF)/DBLE(N)
1500 0 : IP=K-NP-IW*N
1501 0 : IQ=IP+IW
1502 0 : if (IQ > N) IQ=IQ-N
1503 0 : PTSID(K)=DBLE(IP)+DBLE(IQ)/DBLE(NP)+SFRAC
1504 0 : TEMP=ONE/(PTSAUX(1,IP)*PTSAUX(1,IQ))
1505 0 : ZMAT(1,K-NP)=TEMP
1506 0 : ZMAT(IP+1,K-NP)=-TEMP
1507 0 : ZMAT(IQ+1,K-NP)=-TEMP
1508 0 : ZMAT(K,K-NP)=TEMP
1509 : end do
1510 : end if
1511 0 : NREM=NPT
1512 0 : KOLD=1
1513 0 : KNEW=KOPT
1514 : !
1515 : ! Reorder the provisional points in the way that exchanges PTSID(KOLD)
1516 : ! with PTSID(KNEW).
1517 : !
1518 0 : 80 do J=1,N
1519 0 : TEMP=BMAT(KOLD,J)
1520 0 : BMAT(KOLD,J)=BMAT(KNEW,J)
1521 0 : BMAT(KNEW,J)=TEMP
1522 : end do
1523 0 : do J=1,NPTM
1524 0 : TEMP=ZMAT(KOLD,J)
1525 0 : ZMAT(KOLD,J)=ZMAT(KNEW,J)
1526 0 : ZMAT(KNEW,J)=TEMP
1527 : end do
1528 0 : PTSID(KOLD)=PTSID(KNEW)
1529 0 : PTSID(KNEW)=ZERO
1530 0 : W(NDIM+KNEW)=ZERO
1531 0 : NREM=NREM-1
1532 0 : if (KNEW /= KOPT) then
1533 0 : TEMP=VLAG(KOLD)
1534 0 : VLAG(KOLD)=VLAG(KNEW)
1535 0 : VLAG(KNEW)=TEMP
1536 : !
1537 : ! Update the BMAT and ZMAT matrices so that the status of the KNEW-th
1538 : ! interpolation point can be changed from provisional to original. The
1539 : ! branch to label 350 occurs if all the original points are reinstated.
1540 : ! The nonnegative values of W(NDIM+K) are required in the search below.
1541 : !
1542 0 : CALL UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
1543 0 : if (NREM == 0) GOTO 350
1544 0 : do K=1,NPT
1545 0 : W(NDIM+K)=DABS(W(NDIM+K))
1546 : end do
1547 : end if
1548 : !
1549 : ! Pick the index KNEW of an original interpolation point that has not
1550 : ! yet replaced one of the provisional interpolation points, giving
1551 : ! attention to the closeness to XOPT and to previous tries with KNEW.
1552 : !
1553 0 : 120 DSQMIN=ZERO
1554 0 : do K=1,NPT
1555 0 : if (W(NDIM+K) > ZERO) then
1556 0 : if (DSQMIN == ZERO .or. W(NDIM+K) < DSQMIN) then
1557 0 : KNEW=K
1558 0 : DSQMIN=W(NDIM+K)
1559 : end if
1560 : end if
1561 : end do
1562 0 : if (DSQMIN == ZERO) GOTO 260
1563 : !
1564 : ! Form the W-vector of the chosen original interpolation point.
1565 : !
1566 0 : do J=1,N
1567 0 : W(NPT+J)=XPT(KNEW,J)
1568 : end do
1569 0 : do K=1,NPT
1570 0 : SUM=ZERO
1571 0 : if (K == KOPT) then
1572 : ! do nothing
1573 0 : else if (PTSID(K) == ZERO) then
1574 0 : do J=1,N
1575 0 : SUM=SUM+W(NPT+J)*XPT(K,J)
1576 : end do
1577 : else
1578 0 : IP=PTSID(K)
1579 0 : if (IP > 0) SUM=W(NPT+IP)*PTSAUX(1,IP)
1580 0 : IQ=DBLE(NP)*PTSID(K)-DBLE(IP*NP)
1581 0 : if (IQ > 0) then
1582 0 : IW=1
1583 0 : if (IP == 0) IW=2
1584 0 : SUM=SUM+W(NPT+IQ)*PTSAUX(IW,IQ)
1585 : end if
1586 : end if
1587 0 : W(K)=HALF*SUM*SUM
1588 : end do
1589 : !
1590 : ! Calculate VLAG and BETA for the required updating of the H matrix if
1591 : ! XPT(KNEW,.) is reinstated in the set of interpolation points.
1592 : !
1593 0 : do K=1,NPT
1594 : SUM=ZERO
1595 0 : do J=1,N
1596 0 : SUM=SUM+BMAT(K,J)*W(NPT+J)
1597 : end do
1598 0 : VLAG(K)=SUM
1599 : end do
1600 0 : BETA=ZERO
1601 0 : do J=1,NPTM
1602 : SUM=ZERO
1603 0 : do K=1,NPT
1604 0 : SUM=SUM+ZMAT(K,J)*W(K)
1605 : end do
1606 0 : BETA=BETA-SUM*SUM
1607 0 : do K=1,NPT
1608 0 : VLAG(K)=VLAG(K)+SUM*ZMAT(K,J)
1609 : end do
1610 : end do
1611 0 : BSUM=ZERO
1612 : DISTSQ=ZERO
1613 0 : do J=1,N
1614 : SUM=ZERO
1615 0 : do K=1,NPT
1616 0 : SUM=SUM+BMAT(K,J)*W(K)
1617 : end do
1618 0 : JP=J+NPT
1619 0 : BSUM=BSUM+SUM*W(JP)
1620 0 : do IP=NPT+1,NDIM
1621 0 : SUM=SUM+BMAT(IP,J)*W(IP)
1622 : end do
1623 0 : BSUM=BSUM+SUM*W(JP)
1624 0 : VLAG(JP)=SUM
1625 0 : DISTSQ=DISTSQ+XPT(KNEW,J)**2
1626 : end do
1627 0 : BETA=HALF*DISTSQ*DISTSQ+BETA-BSUM
1628 0 : VLAG(KOPT)=VLAG(KOPT)+ONE
1629 : !
1630 : ! KOLD is set to the index of the provisional interpolation point that is
1631 : ! going to be deleted to make way for the KNEW-th original interpolation
1632 : ! point. The choice of KOLD is governed by the avoidance of a small value
1633 : ! of the denominator in the updating calculation of UPDATE.
1634 : !
1635 0 : DENOM=ZERO
1636 0 : VLMXSQ=ZERO
1637 0 : do K=1,NPT
1638 0 : if (PTSID(K) /= ZERO) then
1639 0 : HDIAG=ZERO
1640 0 : do J=1,NPTM
1641 0 : HDIAG=HDIAG+ZMAT(K,J)**2
1642 : end do
1643 0 : DEN=BETA*HDIAG+VLAG(K)**2
1644 0 : if (DEN > DENOM) then
1645 0 : KOLD=K
1646 0 : DENOM=DEN
1647 : end if
1648 : end if
1649 0 : VLMXSQ=DMAX1(VLMXSQ,VLAG(K)**2)
1650 : end do
1651 0 : if (DENOM <= 1.0D-2*VLMXSQ) then
1652 0 : W(NDIM+KNEW)=-W(NDIM+KNEW)-WINC
1653 0 : GOTO 120
1654 : end if
1655 0 : GOTO 80
1656 : !
1657 : ! When label 260 is reached, all the final positions of the interpolation
1658 : ! points have been chosen although any changes have not been included yet
1659 : ! in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart
1660 : ! from the shift of XBASE, the updating of the quadratic model remains to
1661 : ! be done. The following cycle through the new interpolation points begins
1662 : ! by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero,
1663 : ! except that a RETURN occurs if MAXFUN prohibits another value of F.
1664 : !
1665 0 : 260 do KPT=1,NPT
1666 0 : if (PTSID(KPT) == ZERO) CYCLE
1667 0 : if (NF >= MAXFUN) then
1668 0 : NF=-1
1669 0 : GOTO 350
1670 : end if
1671 0 : IH=0
1672 0 : do J=1,N
1673 0 : W(J)=XPT(KPT,J)
1674 0 : XPT(KPT,J)=ZERO
1675 0 : TEMP=PQ(KPT)*W(J)
1676 0 : do I=1,J
1677 0 : IH=IH+1
1678 0 : HQ(IH)=HQ(IH)+TEMP*W(I)
1679 : end do
1680 : end do
1681 0 : PQ(KPT)=ZERO
1682 0 : IP=PTSID(KPT)
1683 0 : IQ=DBLE(NP)*PTSID(KPT)-DBLE(IP*NP)
1684 0 : if (IP > 0) then
1685 0 : XP=PTSAUX(1,IP)
1686 0 : XPT(KPT,IP)=XP
1687 : end if
1688 0 : if (IQ > 0) then
1689 0 : XQ=PTSAUX(1,IQ)
1690 0 : if (IP == 0) XQ=PTSAUX(2,IQ)
1691 0 : XPT(KPT,IQ)=XQ
1692 : end if
1693 : !
1694 : ! Set VQUAD to the value of the current model at the new point.
1695 : !
1696 0 : VQUAD=FBASE
1697 0 : if (IP > 0) then
1698 0 : IHP=(IP+IP*IP)/2
1699 0 : VQUAD=VQUAD+XP*(GOPT(IP)+HALF*XP*HQ(IHP))
1700 : end if
1701 0 : if (IQ > 0) then
1702 0 : IHQ=(IQ+IQ*IQ)/2
1703 0 : VQUAD=VQUAD+XQ*(GOPT(IQ)+HALF*XQ*HQ(IHQ))
1704 0 : if (IP > 0) then
1705 0 : IW=MAX0(IHP,IHQ)-IABS(IP-IQ)
1706 0 : VQUAD=VQUAD+XP*XQ*HQ(IW)
1707 : end if
1708 : end if
1709 0 : do K=1,NPT
1710 0 : TEMP=ZERO
1711 0 : if (IP > 0) TEMP=TEMP+XP*XPT(K,IP)
1712 0 : if (IQ > 0) TEMP=TEMP+XQ*XPT(K,IQ)
1713 0 : VQUAD=VQUAD+HALF*PQ(K)*TEMP*TEMP
1714 : end do
1715 : !
1716 : ! Calculate F at the new interpolation point, and set DIFF to the factor
1717 : ! that is going to multiply the KPT-th Lagrange function when the model
1718 : ! is updated to provide interpolation to the new function value.
1719 : !
1720 0 : do I=1,N
1721 0 : W(I)=DMIN1(DMAX1(XL(I),XBASE(I)+XPT(KPT,I)),XU(I))
1722 0 : if (XPT(KPT,I) == SL(I)) W(I)=XL(I)
1723 0 : if (XPT(KPT,I) == SU(I)) W(I)=XU(I)
1724 : end do
1725 0 : NF=NF+1
1726 0 : CALL CALFUN (N,W,F)
1727 0 : if (IPRINT == 3) then
1728 0 : PRINT 300, NF,F,(W(I),I=1,N)
1729 : 300 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,' The corresponding X is:'/(2X,5D15.6))
1730 : end if
1731 0 : FVAL(KPT)=F
1732 0 : if (F < FVAL(KOPT)) KOPT=KPT
1733 0 : DIFF=F-VQUAD
1734 : !
1735 : ! Update the quadratic model. The RETURN from the subroutine occurs when
1736 : ! all the new interpolation points are included in the model.
1737 : !
1738 0 : do I=1,N
1739 0 : GOPT(I)=GOPT(I)+DIFF*BMAT(KPT,I)
1740 : end do
1741 0 : do K=1,NPT
1742 : SUM=ZERO
1743 0 : do J=1,NPTM
1744 0 : SUM=SUM+ZMAT(K,J)*ZMAT(KPT,J)
1745 : end do
1746 0 : TEMP=DIFF*SUM
1747 0 : if (PTSID(K) == ZERO) then
1748 0 : PQ(K)=PQ(K)+TEMP
1749 : else
1750 0 : IP=PTSID(K)
1751 0 : IQ=DBLE(NP)*PTSID(K)-DBLE(IP*NP)
1752 0 : IHQ=(IQ*IQ+IQ)/2
1753 0 : if (IP == 0) then
1754 0 : HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(2,IQ)**2
1755 : else
1756 0 : IHP=(IP*IP+IP)/2
1757 0 : HQ(IHP)=HQ(IHP)+TEMP*PTSAUX(1,IP)**2
1758 0 : if (IQ > 0) then
1759 0 : HQ(IHQ)=HQ(IHQ)+TEMP*PTSAUX(1,IQ)**2
1760 0 : IW=MAX0(IHP,IHQ)-IABS(IQ-IP)
1761 0 : HQ(IW)=HQ(IW)+TEMP*PTSAUX(1,IP)*PTSAUX(1,IQ)
1762 : end if
1763 : end if
1764 : end if
1765 : end do
1766 0 : PTSID(KPT)=ZERO
1767 : end do
1768 0 : 350 RETURN
1769 : end subroutine RESCUE
1770 :
1771 :
1772 227 : subroutine TRSBOX (N,NPT,XPT,XOPT,GOPT,HQ,PQ,SL,SU,DELTA,
1773 : & XNEW,D,GNEW,XBDI,S,HS,HRED,DSQ,CRVMIN)
1774 : implicit real(dp) (A-H,O-Z)
1775 : dimension XPT(NPT,*),XOPT(*),GOPT(*),HQ(*),PQ(*),SL(*),SU(*),
1776 : & XNEW(*),D(*),GNEW(*),XBDI(*),S(*),HS(*),HRED(*)
1777 : !
1778 : ! The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same
1779 : ! meanings as the corresponding arguments of BOBYQB.
1780 : ! DELTA is the trust region radius for the present calculation, which
1781 : ! seeks a small value of the quadratic model within distance DELTA of
1782 : ! XOPT subject to the bounds on the variables.
1783 : ! XNEW will be set to a new vector of variables that is approximately
1784 : ! the one that minimizes the quadratic model within the trust region
1785 : ! subject to the SL and SU constraints on the variables. It satisfies
1786 : ! as equations the bounds that become active during the calculation.
1787 : ! D is the calculated trial step from XOPT, generated iteratively from an
1788 : ! initial value of zero. Thus XNEW is XOPT+D after the final iteration.
1789 : ! GNEW holds the gradient of the quadratic model at XOPT+D. It is updated
1790 : ! when D is updated.
1791 : ! XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is
1792 : ! set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the
1793 : ! I-th variable has become fixed at a bound, the bound being SL(I) or
1794 : ! SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This
1795 : ! information is accumulated during the construction of XNEW.
1796 : ! The arrays S, HS and HRED are also used for working space. They hold the
1797 : ! current search direction, and the changes in the gradient of Q along S
1798 : ! and the reduced D, respectively, where the reduced D is the same as D,
1799 : ! except that the components of the fixed variables are zero.
1800 : ! DSQ will be set to the square of the length of XNEW-XOPT.
1801 : ! CRVMIN is set to zero if D reaches the trust region boundary. Otherwise
1802 : ! it is set to the least curvature of H that occurs in the conjugate
1803 : ! gradient searches that are not restricted by any constraints. The
1804 : ! value CRVMIN=-1.0D0 is set, however, if all of these searches are
1805 : ! constrained.
1806 : !
1807 : ! A version of the truncated conjugate gradient is applied. If a line
1808 : ! search is restricted by a constraint, then the procedure is restarted,
1809 : ! the values of the variables that are at their bounds being fixed. If
1810 : ! the trust region boundary is reached, then further changes may be made
1811 : ! to D, each one being in the two dimensional space that is spanned
1812 : ! by the current D and the gradient of Q at XOPT+D, staying on the trust
1813 : ! region boundary. Termination occurs when the reduction in Q seems to
1814 : ! be close to the greatest reduction that can be achieved.
1815 : !
1816 : ! Set some constants.
1817 : !
1818 227 : HALF=0.5D0
1819 227 : ONE=1.0D0
1820 227 : ONEMIN=-1.0D0
1821 227 : ZERO=0.0D0
1822 : !
1823 : ! The sign of GOPT(I) gives the sign of the change to the I-th variable
1824 : ! that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether
1825 : ! or not to fix the I-th variable at one of its bounds initially, with
1826 : ! NACT being set to the number of fixed variables. D and GNEW are also
1827 : ! set for the first iteration. DELSQ is the upper bound on the sum of
1828 : ! squares of the free variables. QRED is the reduction in Q so far.
1829 : !
1830 227 : ITERC=0
1831 227 : NACT=0
1832 227 : SQSTP=ZERO
1833 1327 : do I=1,N
1834 1100 : XBDI(I)=ZERO
1835 1100 : if (XOPT(I) <= SL(I)) then
1836 0 : if (GOPT(I) >= ZERO) XBDI(I)=ONEMIN
1837 1100 : else if (XOPT(I) >= SU(I)) then
1838 0 : if (GOPT(I) <= ZERO) XBDI(I)=ONE
1839 : end if
1840 1100 : if (XBDI(I) /= ZERO) NACT=NACT+1
1841 1100 : D(I)=ZERO
1842 1327 : GNEW(I)=GOPT(I)
1843 : end do
1844 454 : DELSQ=DELTA*DELTA
1845 454 : QRED=ZERO
1846 227 : CRVMIN=ONEMIN
1847 : !
1848 : ! Set the next search direction of the conjugate gradient method. It is
1849 : ! the steepest descent direction initially and when the iterations are
1850 : ! restarted because a variable has just been fixed by a bound, and of
1851 : ! course the components of the fixed variables are zero. ITERMAX is an
1852 : ! upper bound on the indices of the conjugate gradient iterations.
1853 : !
1854 227 : 20 BETA=ZERO
1855 913 : 30 STEPSQ=ZERO
1856 4180 : do I=1,N
1857 3494 : if (XBDI(I) /= ZERO) then
1858 0 : S(I)=ZERO
1859 3494 : else if (BETA == ZERO) then
1860 1100 : S(I)=-GNEW(I)
1861 : else
1862 2394 : S(I)=BETA*S(I)-GNEW(I)
1863 : end if
1864 4180 : STEPSQ=STEPSQ+S(I)**2
1865 : end do
1866 686 : if (STEPSQ == ZERO) GOTO 190
1867 686 : if (BETA == ZERO) then
1868 454 : GREDSQ=STEPSQ
1869 227 : ITERMAX=ITERC+N-NACT
1870 : end if
1871 686 : if (GREDSQ*DELSQ <= 1.0D-4*QRED*QRED) GOTO 190
1872 : !
1873 : ! Multiply the search direction by the second derivative matrix of Q and
1874 : ! calculate some scalars for the choice of steplength. Then set BLEN to
1875 : ! the length of the the step to the trust region boundary and STPLEN to
1876 : ! the steplength, ignoring the simple bounds.
1877 : !
1878 686 : GOTO 210
1879 227 : 50 RESID=DELSQ
1880 227 : DS=ZERO
1881 227 : SHS=ZERO
1882 4180 : do I=1,N
1883 4180 : if (XBDI(I) == ZERO) then
1884 3494 : RESID=RESID-D(I)**2
1885 3494 : DS=DS+S(I)*D(I)
1886 3494 : SHS=SHS+S(I)*HS(I)
1887 : end if
1888 : end do
1889 686 : if (RESID <= ZERO) GOTO 90
1890 913 : TEMP=DSQRT(STEPSQ*RESID+DS*DS)
1891 686 : if (DS < ZERO) then
1892 0 : BLEN=(TEMP-DS)/STEPSQ
1893 : else
1894 686 : BLEN=RESID/(TEMP+DS)
1895 : end if
1896 913 : STPLEN=BLEN
1897 686 : if (SHS > ZERO) then
1898 680 : STPLEN=DMIN1(BLEN,GREDSQ/SHS)
1899 : end if
1900 :
1901 : !
1902 : ! Reduce STPLEN if necessary in order to preserve the simple bounds,
1903 : ! letting IACT be the index of the new constrained variable.
1904 : !
1905 : IACT=0
1906 4180 : do I=1,N
1907 4180 : if (S(I) /= ZERO) then
1908 3721 : XSUM=XOPT(I)+D(I)
1909 3494 : if (S(I) > ZERO) then
1910 1713 : TEMP=(SU(I)-XSUM)/S(I)
1911 : else
1912 1781 : TEMP=(SL(I)-XSUM)/S(I)
1913 : end if
1914 3494 : if (TEMP < STPLEN) then
1915 3494 : STPLEN=TEMP
1916 3494 : IACT=I
1917 : end if
1918 : end if
1919 : end do
1920 : !
1921 : ! Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q.
1922 : !
1923 913 : SDEC=ZERO
1924 686 : if (STPLEN > ZERO) then
1925 686 : ITERC=ITERC+1
1926 686 : TEMP=SHS/STEPSQ
1927 686 : if (IACT == 0 .and. TEMP > ZERO) then
1928 680 : CRVMIN=DMIN1(CRVMIN,TEMP)
1929 680 : if (CRVMIN == ONEMIN) CRVMIN=TEMP
1930 : end if
1931 4407 : GGSAV=GREDSQ
1932 : GREDSQ=ZERO
1933 4180 : do I=1,N
1934 3494 : GNEW(I)=GNEW(I)+STPLEN*HS(I)
1935 3494 : if (XBDI(I) == ZERO) GREDSQ=GREDSQ+GNEW(I)**2
1936 4180 : D(I)=D(I)+STPLEN*S(I)
1937 : end do
1938 686 : SDEC=DMAX1(STPLEN*(GGSAV-HALF*STPLEN*SHS),ZERO)
1939 686 : QRED=QRED+SDEC
1940 : end if
1941 : !
1942 : ! Restart the conjugate gradient method if it has hit a new bound.
1943 : !
1944 686 : if (IACT > 0) then
1945 0 : NACT=NACT+1
1946 0 : XBDI(IACT)=ONE
1947 0 : if (S(IACT) < ZERO) XBDI(IACT)=ONEMIN
1948 0 : DELSQ=DELSQ-D(IACT)**2
1949 0 : if (DELSQ <= ZERO) GOTO 90
1950 : GOTO 20
1951 : end if
1952 : !
1953 : ! If STPLEN is less than BLEN, then either apply another conjugate
1954 : ! gradient iteration or RETURN.
1955 : !
1956 686 : if (STPLEN < BLEN) then
1957 573 : if (ITERC == ITERMAX) GOTO 190
1958 485 : if (SDEC <= 0.01D0*QRED) GOTO 190
1959 459 : BETA=GREDSQ/GGSAV
1960 459 : GOTO 30
1961 : end if
1962 113 : 90 CRVMIN=ZERO
1963 : !
1964 : ! Prepare for the alternative iteration by calculating some scalars
1965 : ! and by multiplying the reduced D by the second derivative matrix of
1966 : ! Q, where S holds the reduced D in the call of GGMULT.
1967 : !
1968 113 : 100 if (NACT >= N-1) GOTO 190
1969 227 : DREDSQ=ZERO
1970 227 : DREDG=ZERO
1971 : GREDSQ=ZERO
1972 701 : do I=1,N
1973 701 : if (XBDI(I) == ZERO) then
1974 588 : DREDSQ=DREDSQ+D(I)**2
1975 588 : DREDG=DREDG+D(I)*GNEW(I)
1976 588 : GREDSQ=GREDSQ+GNEW(I)**2
1977 588 : S(I)=D(I)
1978 : else
1979 0 : S(I)=ZERO
1980 : end if
1981 : end do
1982 : ITCSAV=ITERC
1983 137 : GOTO 210
1984 : !
1985 : ! Let the search direction S be a linear combination of the reduced D
1986 : ! and the reduced G that is orthogonal to the reduced D.
1987 : !
1988 250 : 120 ITERC=ITERC+1
1989 250 : TEMP=GREDSQ*DREDSQ-DREDG*DREDG
1990 250 : if (TEMP <= 1.0D-4*QRED*QRED) GOTO 190
1991 241 : TEMP=DSQRT(TEMP)
1992 1549 : do I=1,N
1993 1549 : if (XBDI(I) == ZERO) then
1994 1308 : S(I)=(DREDG*D(I)-DREDSQ*GNEW(I))/TEMP
1995 : else
1996 0 : S(I)=ZERO
1997 : end if
1998 : end do
1999 468 : SREDG=-TEMP
2000 : !
2001 : ! By considering the simple bounds on the variables, calculate an upper
2002 : ! bound on the tangent of half the angle of the alternative iteration,
2003 : ! namely ANGBD, except that, if already a free variable has reached a
2004 : ! bound, there is a branch back to label 100 after fixing that variable.
2005 : !
2006 468 : ANGBD=ONE
2007 241 : IACT=0
2008 1549 : do I=1,N
2009 1549 : if (XBDI(I) == ZERO) then
2010 1535 : TEMPA=XOPT(I)+D(I)-SL(I)
2011 1535 : TEMPB=SU(I)-XOPT(I)-D(I)
2012 1308 : if (TEMPA <= ZERO) then
2013 0 : NACT=NACT+1
2014 0 : XBDI(I)=ONEMIN
2015 0 : GOTO 100
2016 1308 : else if (TEMPB <= ZERO) then
2017 0 : NACT=NACT+1
2018 0 : XBDI(I)=ONE
2019 0 : GOTO 100
2020 : end if
2021 1535 : RATIO=ONE
2022 1535 : SSQ=D(I)**2+S(I)**2
2023 1308 : TEMP=SSQ-(XOPT(I)-SL(I))**2
2024 1308 : if (TEMP > ZERO) then
2025 0 : TEMP=DSQRT(TEMP)-S(I)
2026 0 : if (ANGBD*TEMP > TEMPA) then
2027 0 : ANGBD=TEMPA/TEMP
2028 0 : IACT=I
2029 0 : XSAV=ONEMIN
2030 : end if
2031 : end if
2032 1308 : TEMP=SSQ-(SU(I)-XOPT(I))**2
2033 1308 : if (TEMP > ZERO) then
2034 0 : TEMP=DSQRT(TEMP)+S(I)
2035 0 : if (ANGBD*TEMP > TEMPB) then
2036 0 : ANGBD=TEMPB/TEMP
2037 0 : IACT=I
2038 0 : XSAV=ONE
2039 : end if
2040 : end if
2041 : end if
2042 : end do
2043 : !
2044 : ! Calculate HHD and some curvatures for the alternative iteration.
2045 : !
2046 241 : GOTO 210
2047 : 150 SHS=ZERO
2048 227 : DHS=ZERO
2049 227 : DHD=ZERO
2050 1549 : do I=1,N
2051 1549 : if (XBDI(I) == ZERO) then
2052 1308 : SHS=SHS+S(I)*HS(I)
2053 1308 : DHS=DHS+D(I)*HS(I)
2054 1308 : DHD=DHD+D(I)*HRED(I)
2055 : end if
2056 : end do
2057 : !
2058 : ! Seek the greatest reduction in Q for a range of equally spaced values
2059 : ! of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of
2060 : ! the alternative iteration.
2061 : !
2062 468 : REDMAX=ZERO
2063 241 : ISAV=0
2064 468 : REDSAV=ZERO
2065 241 : IU=17.0D0*ANGBD+3.1D0
2066 5061 : do I=1,IU
2067 5047 : ANGT=ANGBD*DBLE(I)/DBLE(IU)
2068 5047 : STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
2069 4820 : TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
2070 5047 : REDNEW=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
2071 4820 : if (REDNEW > REDMAX) then
2072 : REDMAX=REDNEW
2073 : ISAV=I
2074 227 : RDPREV=REDSAV
2075 4512 : else if (I == ISAV+1) then
2076 241 : RDNEXT=REDNEW
2077 : end if
2078 5061 : REDSAV=REDNEW
2079 : end do
2080 : !
2081 : ! Return if the reduction is zero. Otherwise, set the sine and cosine
2082 : ! of the angle of the alternative iteration, and calculate SDEC.
2083 : !
2084 241 : if (ISAV == 0) GOTO 190
2085 175 : if (ISAV < IU) then
2086 175 : TEMP=(RDNEXT-RDPREV)/(REDMAX+REDMAX-RDPREV-RDNEXT)
2087 175 : ANGT=ANGBD*(DBLE(ISAV)+HALF*TEMP)/DBLE(IU)
2088 : end if
2089 227 : CTH=(ONE-ANGT*ANGT)/(ONE+ANGT*ANGT)
2090 175 : STH=(ANGT+ANGT)/(ONE+ANGT*ANGT)
2091 175 : TEMP=SHS+ANGT*(ANGT*DHD-DHS-DHS)
2092 175 : SDEC=STH*(ANGT*DREDG-SREDG-HALF*STH*TEMP)
2093 175 : if (SDEC <= ZERO) GOTO 190
2094 : !
2095 : ! Update GNEW, D and HRED. If the angle of the alternative iteration
2096 : ! is restricted by a bound on a free variable, that variable is fixed
2097 : ! at the bound.
2098 : !
2099 : DREDG=ZERO
2100 : GREDSQ=ZERO
2101 1119 : do I=1,N
2102 944 : GNEW(I)=GNEW(I)+(CTH-ONE)*HRED(I)+STH*HS(I)
2103 944 : if (XBDI(I) == ZERO) then
2104 944 : D(I)=CTH*D(I)+STH*S(I)
2105 944 : DREDG=DREDG+D(I)*GNEW(I)
2106 944 : GREDSQ=GREDSQ+GNEW(I)**2
2107 : end if
2108 1119 : HRED(I)=CTH*HRED(I)+STH*HS(I)
2109 : end do
2110 175 : QRED=QRED+SDEC
2111 175 : if (IACT > 0 .and. ISAV == IU) then
2112 0 : NACT=NACT+1
2113 0 : XBDI(IACT)=XSAV
2114 0 : GOTO 100
2115 : end if
2116 : !
2117 : ! If SDEC is sufficiently small, then RETURN after setting XNEW to
2118 : ! XOPT+D, giving careful attention to the bounds.
2119 : !
2120 175 : if (SDEC > 0.01D0*QRED) GOTO 120
2121 227 : 190 DSQ=ZERO
2122 1327 : do I=1,N
2123 1100 : XNEW(I)=DMAX1(DMIN1(XOPT(I)+D(I),SU(I)),SL(I))
2124 1100 : if (XBDI(I) == ONEMIN) XNEW(I)=SL(I)
2125 1100 : if (XBDI(I) == ONE) XNEW(I)=SU(I)
2126 1100 : D(I)=XNEW(I)-XOPT(I)
2127 1327 : DSQ=DSQ+D(I)**2
2128 : end do
2129 227 : RETURN
2130 :
2131 : ! The following instructions multiply the current S-vector by the second
2132 : ! derivative matrix of the quadratic model, putting the product in HS.
2133 : ! They are reached from three different parts of the software above and
2134 : ! they can be regarded as an external subroutine.
2135 : !
2136 1040 : 210 IH=0
2137 6430 : do J=1,N
2138 5390 : HS(J)=ZERO
2139 23871 : do I=1,J
2140 17441 : IH=IH+1
2141 17441 : if (I < J) HS(J)=HS(J)+HQ(IH)*S(I)
2142 22831 : HS(I)=HS(I)+HQ(IH)*S(J)
2143 : end do
2144 : end do
2145 12860 : do K=1,NPT
2146 12860 : if (PQ(K) /= ZERO) then
2147 : TEMP=ZERO
2148 75756 : do J=1,N
2149 75756 : TEMP=TEMP+XPT(K,J)*S(J)
2150 : end do
2151 11744 : TEMP=TEMP*PQ(K)
2152 75756 : do I=1,N
2153 75756 : HS(I)=HS(I)+TEMP*XPT(K,I)
2154 : end do
2155 : end if
2156 : end do
2157 1040 : if (CRVMIN /= ZERO) GOTO 50
2158 354 : if (ITERC > ITCSAV) GOTO 150
2159 701 : do I=1,N
2160 701 : HRED(I)=HS(I)
2161 : end do
2162 : GOTO 120
2163 : end subroutine TRSBOX
2164 :
2165 :
2166 260 : subroutine UPDATE (N,NPT,BMAT,ZMAT,NDIM,VLAG,BETA,DENOM,KNEW,W)
2167 : implicit real(dp) (A-H,O-Z)
2168 : dimension BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*)
2169 : !
2170 : ! The arrays BMAT and ZMAT are updated, as required by the new position
2171 : ! of the interpolation point that has the index KNEW. The vector VLAG has
2172 : ! N+NPT components, set on entry to the first NPT and last N components
2173 : ! of the product Hw in equation (4.11) of the Powell (2006) paper on
2174 : ! NEWUOA. Further, BETA is set on entry to the value of the parameter
2175 : ! with that name, and DENOM is set to the denominator of the updating
2176 : ! formula. Elements of ZMAT may be treated as zero if their moduli are
2177 : ! at most ZTEST. The first NDIM elements of W are used for working space.
2178 : !
2179 : ! Set some constants.
2180 : !
2181 260 : ONE=1.0D0
2182 260 : ZERO=0.0D0
2183 260 : NPTM=NPT-N-1
2184 260 : ZTEST=ZERO
2185 3124 : do K=1,NPT
2186 18418 : do J=1,NPTM
2187 18158 : ZTEST=DMAX1(ZTEST,DABS(ZMAT(K,J)))
2188 : end do
2189 : end do
2190 260 : ZTEST=1.0D-20*ZTEST
2191 : !
2192 : ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
2193 : !
2194 260 : JL=1
2195 1302 : do J=2,NPTM
2196 1042 : if (DABS(ZMAT(KNEW,J)) > ZTEST) then
2197 1218 : TEMP=DSQRT(ZMAT(KNEW,1)**2+ZMAT(KNEW,J)**2)
2198 1218 : TEMPA=ZMAT(KNEW,1)/TEMP
2199 1218 : TEMPB=ZMAT(KNEW,J)/TEMP
2200 12400 : do I=1,NPT
2201 11442 : TEMP=TEMPA*ZMAT(I,1)+TEMPB*ZMAT(I,J)
2202 11442 : ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,1)
2203 12400 : ZMAT(I,1)=TEMP
2204 : end do
2205 : end if
2206 1302 : ZMAT(KNEW,J)=ZERO
2207 : end do
2208 :
2209 : ! Put the first NPT components of the KNEW-th column of HLAG into W,
2210 : ! and calculate the parameters of the updating formula.
2211 :
2212 3124 : do I=1,NPT
2213 3124 : W(I)=ZMAT(KNEW,1)*ZMAT(I,1)
2214 : end do
2215 520 : ALPHA=W(KNEW)
2216 260 : TAU=VLAG(KNEW)
2217 260 : VLAG(KNEW)=VLAG(KNEW)-ONE
2218 :
2219 : ! Complete the updating of ZMAT.
2220 :
2221 260 : TEMP=DSQRT(DENOM)
2222 260 : TEMPB=ZMAT(KNEW,1)/TEMP
2223 260 : TEMPA=TAU/TEMP
2224 3124 : do I=1,NPT
2225 3124 : ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
2226 : end do
2227 :
2228 : ! Finally, update the matrix BMAT.
2229 :
2230 1562 : do J=1,N
2231 1302 : JP=NPT+J
2232 1302 : W(JP)=BMAT(KNEW,J)
2233 1302 : TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM
2234 1302 : TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM
2235 21005 : do I=1,JP
2236 19443 : BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I)
2237 20745 : if (I > NPT) BMAT(JP,I-NPT)=BMAT(I,J)
2238 : end do
2239 : end do
2240 260 : RETURN
2241 : end subroutine UPDATE
2242 :
2243 : end module mod_bobyqa
|
