Line data Source code
1 : module mod_newuoa
2 : use const_def, only: dp
3 : use math_lib
4 :
5 : contains
6 :
7 3 : subroutine do_newuoa(N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W,CALFUN,max_valid_value)
8 : implicit real(dp) (A-H,O-Z)
9 : dimension X(*),W(*)
10 : interface
11 : #include "num_newuoa_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 a trust region method that forms quadratic models by interpolation.
17 : ! There can be some freedom in the interpolation conditions, which is
18 : ! taken up by minimizing the Frobenius norm of the change to the second
19 : ! derivative of the quadratic model, beginning with a zero matrix. The
20 : ! arguments of the subroutine are as follows.
21 : !
22 : ! N must be set to the number of variables and must be at least two.
23 : ! NPT is the number of interpolation conditions. Its value must be in the
24 : ! interval [N+2,(N+1)(N+2)/2].
25 : ! Initial values of the variables must be set in X(1),X(2),...,X(N). They
26 : ! will be changed to the values that give the least calculated F.
27 : ! RHOBEG and RHOEND must be set to the initial and final values of a trust
28 : ! region radius, so both must be positive with RHOEND<=RHOBEG. Typically
29 : ! RHOBEG should be about one tenth of the greatest expected change to a
30 : ! variable, and RHOEND should indicate the accuracy that is required in
31 : ! the final values of the variables.
32 : ! The value of IPRINT should be set to 0, 1, 2 or 3, which controls the
33 : ! amount of printing. Specifically, there is no output if IPRINT=0 and
34 : ! there is output only at the return if IPRINT=1. Otherwise, each new
35 : ! value of RHO is printed, with the best vector of variables so far and
36 : ! the corresponding value of the objective function. Further, each new
37 : ! value of F with its variables are output if IPRINT=3.
38 : ! MAXFUN must be set to an upper bound on the number of calls of CALFUN.
39 : ! The array W will be used for working space. Its length must be at least
40 : ! (NPT+13)*(NPT+N)+3*N*(N+3)/2.
41 : !
42 : ! subroutine CALFUN (N,X,F) must be provided by the user. It must set F to
43 : ! the value of the objective function for the variables X(1),X(2),...,X(N).
44 : !
45 : ! Partition the working space array, so that different parts of it can be
46 : ! treated separately by the subroutine that performs the main calculation.
47 : !
48 3 : NP=N+1
49 3 : NPTM=NPT-NP
50 3 : if (NPT < N+2 .or. NPT > ((N+2)*NP)/2) then
51 0 : PRINT 10
52 : 10 FORMAT (/4X,'Return from NEWUOA because NPT is not in the required interval')
53 0 : GOTO 20
54 : end if
55 3 : NDIM=NPT+N
56 3 : IXB=1
57 3 : IXO=IXB+N
58 3 : IXN=IXO+N
59 3 : IXP=IXN+N
60 3 : IFV=IXP+N*NPT
61 3 : IGQ=IFV+NPT
62 3 : IHQ=IGQ+N
63 3 : IPQ=IHQ+(N*NP)/2
64 3 : IBMAT=IPQ+NPT
65 3 : IZMAT=IBMAT+NDIM*N
66 3 : ID=IZMAT+NPT*NPTM
67 3 : IVL=ID+N
68 3 : IW=IVL+NDIM
69 : !
70 : ! The above settings provide a partition of W for subroutine NEWUOB.
71 : ! The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of
72 : ! W plus the space that is needed by the last array of NEWUOB.
73 : !
74 : CALL NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB),
75 : & W(IXO),W(IXN),W(IXP),W(IFV),W(IGQ),W(IHQ),W(IPQ),W(IBMAT),
76 3 : & W(IZMAT),NDIM,W(ID),W(IVL),W(IW),CALFUN,max_valid_value)
77 3 : 20 RETURN
78 : end subroutine do_newuoa
79 :
80 :
81 313 : subroutine NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,XBASE,
82 3 : & XOPT,XNEW,XPT,FVAL,GQ,HQ,PQ,BMAT,ZMAT,NDIM,D,VLAG,W,
83 : & CALFUN,max_valid_value)
84 : implicit real(dp) (A-H,O-Z)
85 : interface
86 : #include "num_newuoa_proc.dek"
87 : end interface
88 : real(dp), intent(in) :: max_valid_value
89 : logical :: do_replace
90 : dimension X(*),XBASE(*),XOPT(*),XNEW(*),XPT(NPT,*),FVAL(*),
91 : & GQ(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),VLAG(*),W(*)
92 : !
93 : ! The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical
94 : ! to the corresponding arguments in subroutine NEWUOA.
95 : ! XBASE will hold a shift of origin that should reduce the contributions
96 : ! from rounding errors to values of the model and Lagrange functions.
97 : ! XOPT will be set to the displacement from XBASE of the vector of
98 : ! variables that provides the least calculated F so far.
99 : ! XNEW will be set to the displacement from XBASE of the vector of
100 : ! variables for the current calculation of F.
101 : ! XPT will contain the interpolation point coordinates relative to XBASE.
102 : ! FVAL will hold the values of F at the interpolation points.
103 : ! GQ will hold the gradient of the quadratic model at XBASE.
104 : ! HQ will hold the explicit second derivatives of the quadratic model.
105 : ! PQ will contain the parameters of the implicit second derivatives of
106 : ! the quadratic model.
107 : ! BMAT will hold the last N columns of H.
108 : ! ZMAT will hold the factorization of the leading NPT by NPT submatrix of
109 : ! H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where
110 : ! the elements of DZ are plus or minus one, as specified by IDZ.
111 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
112 : ! D is reserved for trial steps from XOPT.
113 : ! VLAG will contain the values of the Lagrange functions at a new point X.
114 : ! They are part of a product that requires VLAG to be of length NDIM.
115 : ! The array W will be used for working space. Its length must be at least
116 : ! 10*NDIM = 10*(NPT+N).
117 : !
118 : ! Set some constants.
119 : !
120 3 : HALF=0.5D0
121 3 : ONE=1.0D0
122 3 : TENTH=0.1D0
123 3 : ZERO=0.0D0
124 3 : NP=N+1
125 3 : NH=(N*NP)/2
126 3 : NPTM=NPT-NP
127 3 : NFTEST=MAX0(MAXFUN,1)
128 : !
129 : ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
130 : !
131 15 : do J=1,N
132 12 : XBASE(J)=X(J)
133 136 : do K=1,NPT
134 136 : XPT(K,J)=ZERO
135 : end do
136 195 : do I=1,NDIM
137 192 : BMAT(I,J)=ZERO
138 : end do
139 : end do
140 37 : do IH=1,NH
141 37 : HQ(IH)=ZERO
142 : end do
143 30 : do K=1,NPT
144 27 : PQ(K)=ZERO
145 154 : do J=1,NPTM
146 151 : ZMAT(K,J)=ZERO
147 : end do
148 : end do
149 : !
150 : ! Begin the initialization procedure. NF becomes one more than the number
151 : ! of function values so far. The coordinates of the displacement of the
152 : ! next initial interpolation point from XBASE are set in XPT(NF,.).
153 : !
154 6 : RHOSQ=RHOBEG*RHOBEG
155 6 : RECIP=ONE/RHOSQ
156 6 : RECIQ=DSQRT(HALF)/RHOSQ
157 3 : NF=0
158 27 : 50 NFM=NF
159 27 : NFMM=NF-N
160 27 : NF=NF+1
161 27 : if (NFM <= 2*N) then
162 27 : if (NFM >= 1 .and. NFM <= N) then
163 12 : XPT(NF,NFM)=RHOBEG
164 15 : else if (NFM > N) then
165 12 : XPT(NF,NFMM)=-RHOBEG
166 : end if
167 : else
168 0 : ITEMP=(NFMM-1)/N
169 0 : JPT=NFM-ITEMP*N-N
170 0 : IPT=JPT+ITEMP
171 0 : if (IPT > N) then
172 0 : ITEMP=JPT
173 0 : JPT=IPT-N
174 0 : IPT=ITEMP
175 : end if
176 3 : XIPT=RHOBEG
177 0 : if (FVAL(IPT+NP) < FVAL(IPT+1)) XIPT=-XIPT
178 3 : XJPT=RHOBEG
179 0 : if (FVAL(JPT+NP) < FVAL(JPT+1)) XJPT=-XJPT
180 0 : XPT(NF,IPT)=XIPT
181 0 : XPT(NF,JPT)=XJPT
182 : end if
183 : !
184 : ! Calculate the next value of F, label 70 being reached immediately
185 : ! after this calculation. The least function value so far and its index
186 : ! are required.
187 : !
188 151 : do J=1,N
189 151 : X(J)=XPT(NF,J)+XBASE(J)
190 : end do
191 27 : GOTO 310
192 30 : 70 FVAL(NF)=F
193 27 : if (NF == 1) then
194 6 : FBEG=F
195 6 : FOPT=F
196 3 : KOPT=1
197 24 : else if (F < FOPT) then
198 8 : FOPT=F
199 8 : KOPT=NF
200 : end if
201 : !
202 : ! Set the nonzero initial elements of BMAT and the quadratic model in
203 : ! the cases when NF is at most 2*N+1.
204 : !
205 27 : if (NFM <= 2*N) then
206 27 : if (NFM >= 1 .and. NFM <= N) then
207 12 : GQ(NFM)=(F-FBEG)/RHOBEG
208 12 : if (NPT < NF+N) then
209 0 : BMAT(1,NFM)=-ONE/RHOBEG
210 0 : BMAT(NF,NFM)=ONE/RHOBEG
211 0 : BMAT(NPT+NFM,NFM)=-HALF*RHOSQ
212 : end if
213 15 : else if (NFM > N) then
214 12 : BMAT(NF-N,NFMM)=HALF/RHOBEG
215 12 : BMAT(NF,NFMM)=-HALF/RHOBEG
216 12 : ZMAT(1,NFMM)=-RECIQ-RECIQ
217 12 : ZMAT(NF-N,NFMM)=RECIQ
218 12 : ZMAT(NF,NFMM)=RECIQ
219 12 : IH=(NFMM*(NFMM+1))/2
220 15 : TEMP=(FBEG-F)/RHOBEG
221 12 : HQ(IH)=(GQ(NFMM)-TEMP)/RHOBEG
222 12 : GQ(NFMM)=HALF*(GQ(NFMM)+TEMP)
223 : end if
224 : !
225 : ! Set the off-diagonal second derivatives of the Lagrange functions and
226 : ! the initial quadratic model.
227 : !
228 : else
229 0 : IH=(IPT*(IPT-1))/2+JPT
230 0 : if (XIPT < ZERO) IPT=IPT+N
231 0 : if (XJPT < ZERO) JPT=JPT+N
232 0 : ZMAT(1,NFMM)=RECIP
233 0 : ZMAT(NF,NFMM)=RECIP
234 0 : ZMAT(IPT+1,NFMM)=-RECIP
235 0 : ZMAT(JPT+1,NFMM)=-RECIP
236 0 : HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/(XIPT*XJPT)
237 : end if
238 27 : if (NF < NPT) GOTO 50
239 : !
240 : ! Begin the iterative procedure, because the initial model is complete.
241 : !
242 6 : RHO=RHOBEG
243 6 : DELTA=RHO
244 3 : IDZ=1
245 6 : DIFFA=ZERO
246 6 : DIFFB=ZERO
247 3 : ITEST=0
248 6 : XOPTSQ=ZERO
249 15 : do I=1,N
250 12 : XOPT(I)=XPT(KOPT,I)
251 15 : XOPTSQ=XOPTSQ+XOPT(I)**2
252 : end do
253 18 : 90 NFSAV=NF
254 : !
255 : ! Generate the next trust region step and test its length. Set KNEW
256 : ! to -1 if the purpose of the next F will be to improve the model.
257 : !
258 252 : 100 KNEW=0
259 255 : CALL TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,D,W,W(NP),W(NP+N),W(NP+2*N),CRVMIN)
260 255 : DSQ=ZERO
261 1492 : do I=1,N
262 1492 : DSQ=DSQ+D(I)**2
263 : end do
264 255 : DNORM=DMIN1(DELTA,DSQRT(DSQ))
265 252 : if (DNORM < HALF*RHO) then
266 91 : KNEW=-1
267 91 : DELTA=TENTH*DELTA
268 94 : RATIO=-1.0D0
269 91 : if (DELTA <= 1.5D0*RHO) DELTA=RHO
270 91 : if (NF <= NFSAV+2) GOTO 460
271 62 : TEMP=0.125D0*CRVMIN*RHO*RHO
272 62 : if (TEMP <= DMAX1(DIFFA,DIFFB,DIFFC)) GOTO 460
273 : GOTO 490
274 : end if
275 : !
276 : ! Shift XBASE if XOPT may be too far from XBASE. First make the changes
277 : ! to BMAT that do not depend on ZMAT.
278 : !
279 281 : 120 if (DSQ <= 1.0D-3*XOPTSQ) then
280 12 : TEMPQ=0.25D0*XOPTSQ
281 98 : do K=1,NPT
282 92 : SUM=ZERO
283 529 : do I=1,N
284 529 : SUM=SUM+XPT(K,I)*XOPT(I)
285 : end do
286 89 : TEMP=PQ(K)*SUM
287 89 : SUM=SUM-HALF*XOPTSQ
288 89 : W(NPT+K)=SUM
289 538 : do I=1,N
290 440 : GQ(I)=GQ(I)+TEMP*XPT(K,I)
291 440 : XPT(K,I)=XPT(K,I)-HALF*XOPT(I)
292 440 : VLAG(I)=BMAT(K,I)
293 440 : W(I)=SUM*XPT(K,I)+TEMPQ*XOPT(I)
294 440 : IP=NPT+I
295 1921 : do J=1,I
296 1832 : BMAT(IP,J)=BMAT(IP,J)+VLAG(I)*W(J)+W(I)*VLAG(J)
297 : end do
298 : end do
299 : end do
300 : !
301 : ! Then the revisions of BMAT that depend on ZMAT are calculated.
302 : !
303 49 : do K=1,NPTM
304 3 : SUMZ=ZERO
305 480 : do I=1,NPT
306 440 : SUMZ=SUMZ+ZMAT(I,K)
307 480 : W(I)=W(NPT+I)*ZMAT(I,K)
308 : end do
309 240 : do J=1,N
310 200 : SUM=TEMPQ*SUMZ*XOPT(J)
311 2544 : do I=1,NPT
312 2544 : SUM=SUM+W(I)*XPT(I,J)
313 : end do
314 200 : VLAG(J)=SUM
315 200 : if (K < IDZ) SUM=-SUM
316 2584 : do I=1,NPT
317 2544 : BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K)
318 : end do
319 : end do
320 249 : do I=1,N
321 200 : IP=I+NPT
322 200 : TEMP=VLAG(I)
323 200 : if (K < IDZ) TEMP=-TEMP
324 876 : do J=1,I
325 836 : BMAT(IP,J)=BMAT(IP,J)+TEMP*VLAG(J)
326 : end do
327 : end do
328 : end do
329 : !
330 : ! The following instructions complete the shift of XBASE, including
331 : ! the changes to the parameters of the quadratic model.
332 : !
333 9 : IH=0
334 49 : do J=1,N
335 40 : W(J)=ZERO
336 480 : do K=1,NPT
337 440 : W(J)=W(J)+PQ(K)*XPT(K,J)
338 480 : XPT(K,J)=XPT(K,J)-HALF*XOPT(J)
339 : end do
340 169 : do I=1,J
341 120 : IH=IH+1
342 120 : if (I < J) GQ(J)=GQ(J)+HQ(IH)*XOPT(I)
343 120 : GQ(I)=GQ(I)+HQ(IH)*XOPT(J)
344 120 : HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
345 160 : BMAT(NPT+I,J)=BMAT(NPT+J,I)
346 : end do
347 : end do
348 49 : do J=1,N
349 40 : XBASE(J)=XBASE(J)+XOPT(J)
350 49 : XOPT(J)=ZERO
351 : end do
352 : XOPTSQ=ZERO
353 : end if
354 : !
355 : ! Pick the model step if KNEW is positive. A different choice of D
356 : ! may be made later, if the choice of D by BIGLAG causes substantial
357 : ! cancellation in DENOM.
358 : !
359 281 : if (KNEW > 0) then
360 120 : CALL BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DSTEP,D,ALPHA,VLAG,VLAG(NPT+1),W,W(NP),W(NP+N))
361 : end if
362 : !
363 : ! Calculate VLAG and BETA for the current choice of D. The first NPT
364 : ! components of W_check will be held in W.
365 : !
366 3410 : do K=1,NPT
367 3132 : SUMA=ZERO
368 3132 : SUMB=ZERO
369 3129 : SUM=ZERO
370 19945 : do J=1,N
371 16816 : SUMA=SUMA+XPT(K,J)*D(J)
372 16816 : SUMB=SUMB+XPT(K,J)*XOPT(J)
373 19945 : SUM=SUM+BMAT(K,J)*D(J)
374 : end do
375 3129 : W(K)=SUMA*(HALF*SUMA+SUMB)
376 3410 : VLAG(K)=SUM
377 : end do
378 284 : BETA=ZERO
379 1705 : do K=1,NPTM
380 : SUM=ZERO
381 18240 : do I=1,NPT
382 18240 : SUM=SUM+ZMAT(I,K)*W(I)
383 : end do
384 1424 : if (K < IDZ) then
385 0 : BETA=BETA+SUM*SUM
386 0 : SUM=-SUM
387 : else
388 1424 : BETA=BETA-SUM*SUM
389 : end if
390 18521 : do I=1,NPT
391 18240 : VLAG(I)=VLAG(I)+SUM*ZMAT(I,K)
392 : end do
393 : end do
394 284 : BSUM=ZERO
395 284 : DX=ZERO
396 1705 : do J=1,N
397 : SUM=ZERO
398 18240 : do I=1,NPT
399 18240 : SUM=SUM+W(I)*BMAT(I,J)
400 : end do
401 1424 : BSUM=BSUM+SUM*D(J)
402 1424 : JP=NPT+J
403 9120 : do K=1,N
404 9120 : SUM=SUM+BMAT(JP,K)*D(K)
405 : end do
406 1424 : VLAG(JP)=SUM
407 1424 : BSUM=BSUM+SUM*D(J)
408 1705 : DX=DX+D(J)*XOPT(J)
409 : end do
410 281 : BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
411 281 : VLAG(KOPT)=VLAG(KOPT)+ONE
412 : !
413 : ! If KNEW is positive and if the cancellation in DENOM is unacceptable,
414 : ! then BIGDEN calculates an alternative model step, XNEW being used for
415 : ! working space.
416 : !
417 281 : if (KNEW > 0) then
418 120 : TEMP=ONE+ALPHA*BETA/VLAG(KNEW)**2
419 120 : if (DABS(TEMP) <= 0.8D0) then
420 : CALL BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT,
421 0 : & KNEW,D,W,VLAG,BETA,XNEW,W(NDIM+1),W(6*NDIM+1))
422 : end if
423 : end if
424 : !
425 : ! Calculate the next value of the objective function.
426 : !
427 1713 : 290 do I=1,N
428 1430 : XNEW(I)=XOPT(I)+D(I)
429 1713 : X(I)=XBASE(I)+XNEW(I)
430 : end do
431 283 : NF=NF+1
432 310 : 310 if (NF > NFTEST) then
433 0 : NF=NF-1
434 0 : if (IPRINT > 0) PRINT 320
435 : 320 FORMAT (/4X,'Return from NEWUOA because CALFUN has been called MAXFUN times.')
436 : GOTO 530
437 : end if
438 310 : CALL CALFUN (N,X,F)
439 310 : if (IPRINT == 3) then
440 0 : PRINT 330, NF,F,(X(I),I=1,N)
441 : 330 FORMAT (/4X,'Function number',I6,' F =',1PD18.10,' The corresponding X is:'/(2X,5D15.6))
442 : end if
443 310 : if (NF <= NPT) GOTO 70
444 283 : if (KNEW == -1) GOTO 530
445 : !
446 : ! Use the quadratic model to predict the change in F due to the step D,
447 : ! and set DIFF to the error of this prediction.
448 : !
449 284 : VQUAD=ZERO
450 281 : IH=0
451 1705 : do J=1,N
452 1424 : VQUAD=VQUAD+D(J)*GQ(J)
453 6265 : do I=1,J
454 4560 : IH=IH+1
455 4560 : TEMP=D(I)*XNEW(J)+D(J)*XOPT(I)
456 4560 : if (I == J) TEMP=HALF*TEMP
457 5984 : VQUAD=VQUAD+TEMP*HQ(IH)
458 : end do
459 : end do
460 3410 : do K=1,NPT
461 3410 : VQUAD=VQUAD+PQ(K)*W(K)
462 : end do
463 284 : DIFF=F-FOPT-VQUAD
464 281 : DIFFC=DIFFB
465 281 : DIFFB=DIFFA
466 281 : DIFFA=DABS(DIFF)
467 281 : if (DNORM > RHO) NFSAV=NF
468 : !
469 : ! Update FOPT and XOPT if the new F is the least value of the objective
470 : ! function so far. The branch when KNEW is positive occurs if D is not
471 : ! a trust region step.
472 : !
473 284 : FSAVE=FOPT
474 281 : if (F < FOPT) then
475 131 : FOPT=F
476 131 : XOPTSQ=ZERO
477 803 : do I=1,N
478 672 : XOPT(I)=XNEW(I)
479 803 : XOPTSQ=XOPTSQ+XOPT(I)**2
480 : end do
481 : end if
482 281 : KSAVE=KNEW
483 281 : if (KNEW > 0) GOTO 410
484 : !
485 : ! Pick the next value of DELTA after a trust region step.
486 : !
487 161 : if (VQUAD >= ZERO) then
488 0 : if (IPRINT > 0) PRINT 370
489 : 370 FORMAT (/4X,'Return from NEWUOA because a trust region step has failed to reduce Q.')
490 : GOTO 530
491 : end if
492 161 : RATIO=(F-FSAVE)/VQUAD
493 161 : if (RATIO <= TENTH) then
494 49 : DELTA=HALF*DNORM
495 112 : else if (RATIO. LE. 0.7D0) then
496 40 : DELTA=DMAX1(HALF*DELTA,DNORM)
497 : else
498 72 : DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
499 : end if
500 161 : if (DELTA <= 1.5D0*RHO) DELTA=RHO
501 : !
502 : ! Set KNEW to the index of the next interpolation point to be deleted.
503 : !
504 161 : RHOSQ=DMAX1(TENTH*DELTA,RHO)**2
505 161 : KTEMP=0
506 164 : DETRAT=ZERO
507 161 : if (F >= FSAVE) then
508 47 : KTEMP=KOPT
509 47 : DETRAT=ONE
510 : end if
511 1986 : do K=1,NPT
512 3 : HDIAG=ZERO
513 11761 : do J=1,NPTM
514 9936 : TEMP=ONE
515 9936 : if (J < IDZ) TEMP=-ONE
516 11761 : HDIAG=HDIAG+TEMP*ZMAT(K,J)**2
517 : end do
518 1825 : TEMP=DABS(BETA*HDIAG+VLAG(K)**2)
519 1828 : DISTSQ=ZERO
520 11761 : do J=1,N
521 11761 : DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
522 : end do
523 1825 : if (DISTSQ > RHOSQ) TEMP=TEMP*(DISTSQ/RHOSQ)*(DISTSQ/RHOSQ)*(DISTSQ/RHOSQ)
524 1986 : if (TEMP > DETRAT .and. K /= KTEMP) then
525 478 : DETRAT=TEMP
526 478 : KNEW=K
527 : end if
528 : end do
529 161 : if (KNEW == 0) GOTO 460
530 : !
531 : ! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
532 : ! can be moved. Begin the updating of the quadratic model, starting
533 : ! with the explicit second derivative term.
534 : !
535 281 : 410 CALL UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
536 281 : FVAL(KNEW)=F
537 281 : IH=0
538 1705 : do I=1,N
539 1424 : TEMP=PQ(KNEW)*XPT(KNEW,I)
540 6265 : do J=1,I
541 4560 : IH=IH+1
542 5984 : HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
543 : end do
544 : end do
545 281 : PQ(KNEW)=ZERO
546 : !
547 : ! Update the other second derivative parameters, and then the gradient
548 : ! vector of the model. Also include the new interpolation point.
549 : !
550 1705 : do J=1,NPTM
551 1424 : TEMP=DIFF*ZMAT(KNEW,J)
552 1424 : if (J < IDZ) TEMP=-TEMP
553 18521 : do K=1,NPT
554 18240 : PQ(K)=PQ(K)+TEMP*ZMAT(K,J)
555 : end do
556 : end do
557 3 : GQSQ=ZERO
558 1705 : do I=1,N
559 1424 : GQ(I)=GQ(I)+DIFF*BMAT(KNEW,I)
560 1424 : GQSQ=GQSQ+GQ(I)**2
561 1705 : XPT(KNEW,I)=XNEW(I)
562 : end do
563 : !
564 : ! If a trust region step makes a small change to the objective function,
565 : ! then calculate the gradient of the least Frobenius norm interpolant at
566 : ! XBASE, and store it in W, using VLAG for a vector of right hand sides.
567 : !
568 281 : if (KSAVE == 0 .and. DELTA == RHO) then
569 64 : if (DABS(RATIO) > 1.0D-2) then
570 : ITEST=0
571 : else
572 0 : do K=1,NPT
573 0 : VLAG(K)=FVAL(K)-FVAL(KOPT)
574 : end do
575 3 : GISQ=ZERO
576 0 : do I=1,N
577 : SUM=ZERO
578 0 : do K=1,NPT
579 0 : SUM=SUM+BMAT(K,I)*VLAG(K)
580 : end do
581 0 : GISQ=GISQ+SUM*SUM
582 0 : W(I)=SUM
583 : end do
584 : !
585 : ! Test whether to replace the new quadratic model by the least Frobenius
586 : ! norm interpolant, making the replacement if the test is satisfied.
587 : !
588 0 : ITEST=ITEST+1
589 0 : if (GQSQ < 1.0D2*GISQ) ITEST=0
590 0 : do_replace = (ITEST >= 3)
591 0 : if (.not. do_replace) then ! check for "invalid" value
592 0 : do k=1,npt
593 0 : if (fval(k) > max_valid_value) then
594 : do_replace = .true.
595 : !write(*,*) 'newuoa value > max valid', fval(k), max_valid_value
596 : exit
597 : end if
598 : end do
599 : end if
600 0 : if (do_replace) then
601 0 : do I=1,N
602 0 : GQ(I)=W(I)
603 : end do
604 0 : do IH=1,NH
605 0 : HQ(IH)=ZERO
606 : end do
607 0 : do J=1,NPTM
608 0 : W(J)=ZERO
609 0 : do K=1,NPT
610 0 : W(J)=W(J)+VLAG(K)*ZMAT(K,J)
611 : end do
612 0 : if (J < IDZ) W(J)=-W(J)
613 : end do
614 0 : do K=1,NPT
615 0 : PQ(K)=ZERO
616 0 : do J=1,NPTM
617 0 : PQ(K)=PQ(K)+ZMAT(K,J)*W(J)
618 : end do
619 : end do
620 : ITEST=0
621 : end if
622 : end if
623 : end if
624 281 : if (F < FSAVE) KOPT=KNEW
625 : !
626 : ! If a trust region step has provided a sufficient decrease in F, then
627 : ! branch for another trust region calculation. The case KSAVE>0 occurs
628 : ! when the new function value was calculated by a model step.
629 : !
630 281 : if (F <= FSAVE+TENTH*VQUAD) GOTO 100
631 151 : if (KSAVE > 0) GOTO 100
632 : !
633 : ! Alternatively, find out if the interpolation points are close enough
634 : ! to the best point so far.
635 : !
636 135 : KNEW=0
637 135 : 460 DISTSQ=4.0D0*DELTA*DELTA
638 1586 : do K=1,NPT
639 1451 : SUM=ZERO
640 9045 : do J=1,N
641 9045 : SUM=SUM+(XPT(K,J)-XOPT(J))**2
642 : end do
643 1586 : if (SUM > DISTSQ) then
644 268 : KNEW=K
645 268 : DISTSQ=SUM
646 : end if
647 : end do
648 : !
649 : ! If KNEW is positive, then set DSTEP, and branch back for the next
650 : ! iteration, which will generate a "model step".
651 : !
652 135 : if (KNEW > 0) then
653 120 : DSTEP=DMAX1(DMIN1(TENTH*DSQRT(DISTSQ),HALF*DELTA),RHO)
654 120 : DSQ=DSTEP*DSTEP
655 120 : GOTO 120
656 : end if
657 15 : if (RATIO > ZERO) GOTO 100
658 14 : if (DMAX1(DELTA,DNORM) > RHO) GOTO 100
659 : !
660 : ! The calculations with the current value of RHO are complete. Pick the
661 : ! next values of RHO and DELTA.
662 : !
663 18 : 490 if (RHO > RHOEND) then
664 15 : DELTA=HALF*RHO
665 15 : RATIO=RHO/RHOEND
666 15 : if (RATIO <= 16.0D0) then
667 3 : RHO=RHOEND
668 12 : else if (RATIO <= 250.0D0) then
669 3 : RHO=DSQRT(RATIO)*RHOEND
670 : else
671 9 : RHO=TENTH*RHO
672 : end if
673 15 : DELTA=DMAX1(DELTA,RHO)
674 15 : if (IPRINT >= 2) then
675 0 : if (IPRINT >= 3) PRINT 500
676 : 500 FORMAT (5X)
677 0 : PRINT 510, RHO,NF
678 : 510 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of function values =',I6)
679 0 : PRINT 520, FOPT,(XBASE(I)+XOPT(I),I=1,N)
680 : 520 FORMAT (4X,'Least value of F =',1PD23.15,9X,'The corresponding X is:'/(2X,5D15.6))
681 : end if
682 : GOTO 90
683 : end if
684 : !
685 : ! Return from the calculation, after another Newton-Raphson step, if
686 : ! it is too short to have been tried before.
687 : !
688 3 : if (KNEW == -1) GOTO 290
689 3 : 530 if (FOPT <= F) then
690 7 : do I=1,N
691 7 : X(I)=XBASE(I)+XOPT(I)
692 : end do
693 1 : F=FOPT
694 : end if
695 3 : if (IPRINT >= 1) then
696 0 : PRINT 550, NF
697 : 550 FORMAT (/4X,'At the return from NEWUOA',5X,'Number of function values =',I6)
698 0 : PRINT 520, F,(X(I),I=1,N)
699 : end if
700 3 : RETURN
701 : end subroutine NEWUOB
702 :
703 0 : subroutine BIGDEN(N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT,KNEW,D,W,VLAG,BETA,S,WVEC,PROD)
704 : implicit real(dp) (A-H,O-Z)
705 : dimension XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),
706 : & W(*),VLAG(*),S(*),WVEC(NDIM,*),PROD(NDIM,*)
707 0 : dimension DEN(9),DENEX(9),PAR(9)
708 : !
709 : ! N is the number of variables.
710 : ! NPT is the number of interpolation equations.
711 : ! XOPT is the best interpolation point so far.
712 : ! XPT contains the coordinates of the current interpolation points.
713 : ! BMAT provides the last N columns of H.
714 : ! ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H.
715 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
716 : ! KOPT is the index of the optimal interpolation point.
717 : ! KNEW is the index of the interpolation point that is going to be moved.
718 : ! D will be set to the step from XOPT to the new point, and on entry it
719 : ! should be the D that was calculated by the last call of BIGLAG. The
720 : ! length of the initial D provides a trust region bound on the final D.
721 : ! W will be set to Wcheck for the final choice of D.
722 : ! VLAG will be set to Theta*Wcheck+e_b for the final choice of D.
723 : ! BETA will be set to the value that will occur in the updating formula
724 : ! when the KNEW-th interpolation point is moved to its new position.
725 : ! S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used
726 : ! for working space.
727 : !
728 : ! D is calculated in a way that should provide a denominator with a large
729 : ! modulus in the updating formula when the KNEW-th interpolation point is
730 : ! shifted to the new position XOPT+D.
731 : !
732 : ! Set some constants.
733 : !
734 0 : HALF=0.5D0
735 0 : ONE=1.0D0
736 0 : QUART=0.25D0
737 0 : TWO=2.0D0
738 0 : ZERO=0.0D0
739 0 : TWOPI=8.0D0*atan(ONE)
740 0 : NPTM=NPT-N-1
741 : !
742 : ! Store the first NPT elements of the KNEW-th column of H in W(N+1)
743 : ! to W(N+NPT).
744 : !
745 0 : do K=1,NPT
746 0 : W(N+K)=ZERO
747 : end do
748 0 : do J=1,NPTM
749 0 : TEMP=ZMAT(KNEW,J)
750 0 : if (J < IDZ) TEMP=-TEMP
751 0 : do K=1,NPT
752 0 : W(N+K)=W(N+K)+TEMP*ZMAT(K,J)
753 : end do
754 : end do
755 0 : ALPHA=W(N+KNEW)
756 : !
757 : ! The initial search direction D is taken from the last call of BIGLAG,
758 : ! and the initial S is set below, usually to the direction from X_OPT
759 : ! to X_KNEW, but a different direction to an interpolation point may
760 : ! be chosen, in order to prevent S from being nearly parallel to D.
761 : !
762 0 : DD=ZERO
763 0 : DS=ZERO
764 0 : SS=ZERO
765 0 : XOPTSQ=ZERO
766 0 : do I=1,N
767 0 : DD=DD+D(I)**2
768 0 : S(I)=XPT(KNEW,I)-XOPT(I)
769 0 : DS=DS+D(I)*S(I)
770 0 : SS=SS+S(I)**2
771 0 : XOPTSQ=XOPTSQ+XOPT(I)**2
772 : end do
773 0 : if (DS*DS > 0.99D0*DD*SS) then
774 0 : KSAV=KNEW
775 0 : DTEST=DS*DS/SS
776 0 : do K=1,NPT
777 0 : if (K /= KOPT) then
778 0 : DSTEMP=ZERO
779 0 : SSTEMP=ZERO
780 0 : do I=1,N
781 0 : DIFF=XPT(K,I)-XOPT(I)
782 0 : DSTEMP=DSTEMP+D(I)*DIFF
783 0 : SSTEMP=SSTEMP+DIFF*DIFF
784 : end do
785 0 : if (DSTEMP*DSTEMP/SSTEMP < DTEST) then
786 0 : KSAV=K
787 0 : DTEST=DSTEMP*DSTEMP/SSTEMP
788 0 : DS=DSTEMP
789 0 : SS=SSTEMP
790 : end if
791 : end if
792 : end do
793 0 : do I=1,N
794 0 : S(I)=XPT(KSAV,I)-XOPT(I)
795 : end do
796 : end if
797 0 : SSDEN=DD*SS-DS*DS
798 0 : ITERC=0
799 0 : DENSAV=ZERO
800 : !
801 : ! Begin the iteration by overwriting S with a vector that has the
802 : ! required length and direction.
803 : !
804 0 : 70 ITERC=ITERC+1
805 0 : TEMP=ONE/DSQRT(SSDEN)
806 0 : XOPTD=ZERO
807 0 : XOPTS=ZERO
808 0 : do I=1,N
809 0 : S(I)=TEMP*(DD*S(I)-DS*D(I))
810 0 : XOPTD=XOPTD+XOPT(I)*D(I)
811 0 : XOPTS=XOPTS+XOPT(I)*S(I)
812 : end do
813 : !
814 : ! Set the coefficients of the first two terms of BETA.
815 : !
816 0 : TEMPA=HALF*XOPTD*XOPTD
817 0 : TEMPB=HALF*XOPTS*XOPTS
818 0 : DEN(1)=DD*(XOPTSQ+HALF*DD)+TEMPA+TEMPB
819 0 : DEN(2)=TWO*XOPTD*DD
820 0 : DEN(3)=TWO*XOPTS*DD
821 0 : DEN(4)=TEMPA-TEMPB
822 0 : DEN(5)=XOPTD*XOPTS
823 0 : do I=6,9
824 0 : DEN(I)=ZERO
825 : end do
826 : !
827 : ! Put the coefficients of Wcheck in WVEC.
828 : !
829 0 : do K=1,NPT
830 : TEMPA=ZERO
831 : TEMPB=ZERO
832 0 : TEMPC=ZERO
833 0 : do I=1,N
834 0 : TEMPA=TEMPA+XPT(K,I)*D(I)
835 0 : TEMPB=TEMPB+XPT(K,I)*S(I)
836 0 : TEMPC=TEMPC+XPT(K,I)*XOPT(I)
837 : end do
838 0 : WVEC(K,1)=QUART*(TEMPA*TEMPA+TEMPB*TEMPB)
839 0 : WVEC(K,2)=TEMPA*TEMPC
840 0 : WVEC(K,3)=TEMPB*TEMPC
841 0 : WVEC(K,4)=QUART*(TEMPA*TEMPA-TEMPB*TEMPB)
842 0 : WVEC(K,5)=HALF*TEMPA*TEMPB
843 : end do
844 0 : do I=1,N
845 0 : IP=I+NPT
846 0 : WVEC(IP,1)=ZERO
847 0 : WVEC(IP,2)=D(I)
848 0 : WVEC(IP,3)=S(I)
849 0 : WVEC(IP,4)=ZERO
850 0 : WVEC(IP,5)=ZERO
851 : end do
852 : !
853 : ! Put the coefficients of THETA*Wcheck in PROD.
854 : !
855 0 : do JC=1,5
856 0 : NW=NPT
857 0 : if (JC == 2 .or. JC == 3) NW=NDIM
858 0 : do K=1,NPT
859 0 : PROD(K,JC)=ZERO
860 : end do
861 0 : do J=1,NPTM
862 0 : SUM=ZERO
863 0 : do K=1,NPT
864 0 : SUM=SUM+ZMAT(K,J)*WVEC(K,JC)
865 : end do
866 0 : if (J < IDZ) SUM=-SUM
867 0 : do K=1,NPT
868 0 : PROD(K,JC)=PROD(K,JC)+SUM*ZMAT(K,J)
869 : end do
870 : end do
871 0 : if (NW == NDIM) then
872 0 : do K=1,NPT
873 : SUM=ZERO
874 0 : do J=1,N
875 0 : SUM=SUM+BMAT(K,J)*WVEC(NPT+J,JC)
876 : end do
877 0 : PROD(K,JC)=PROD(K,JC)+SUM
878 : end do
879 : end if
880 0 : do J=1,N
881 : SUM=ZERO
882 0 : do I=1,NW
883 0 : SUM=SUM+BMAT(I,J)*WVEC(I,JC)
884 : end do
885 0 : PROD(NPT+J,JC)=SUM
886 : end do
887 : end do
888 : !
889 : ! Include in DEN the part of BETA that depends on THETA.
890 : !
891 0 : do K=1,NDIM
892 : SUM=ZERO
893 0 : do I=1,5
894 0 : PAR(I)=HALF*PROD(K,I)*WVEC(K,I)
895 0 : SUM=SUM+PAR(I)
896 : end do
897 0 : DEN(1)=DEN(1)-PAR(1)-SUM
898 0 : TEMPA=PROD(K,1)*WVEC(K,2)+PROD(K,2)*WVEC(K,1)
899 0 : TEMPB=PROD(K,2)*WVEC(K,4)+PROD(K,4)*WVEC(K,2)
900 0 : TEMPC=PROD(K,3)*WVEC(K,5)+PROD(K,5)*WVEC(K,3)
901 0 : DEN(2)=DEN(2)-TEMPA-HALF*(TEMPB+TEMPC)
902 0 : DEN(6)=DEN(6)-HALF*(TEMPB-TEMPC)
903 0 : TEMPA=PROD(K,1)*WVEC(K,3)+PROD(K,3)*WVEC(K,1)
904 0 : TEMPB=PROD(K,2)*WVEC(K,5)+PROD(K,5)*WVEC(K,2)
905 0 : TEMPC=PROD(K,3)*WVEC(K,4)+PROD(K,4)*WVEC(K,3)
906 0 : DEN(3)=DEN(3)-TEMPA-HALF*(TEMPB-TEMPC)
907 0 : DEN(7)=DEN(7)-HALF*(TEMPB+TEMPC)
908 0 : TEMPA=PROD(K,1)*WVEC(K,4)+PROD(K,4)*WVEC(K,1)
909 0 : DEN(4)=DEN(4)-TEMPA-PAR(2)+PAR(3)
910 0 : TEMPA=PROD(K,1)*WVEC(K,5)+PROD(K,5)*WVEC(K,1)
911 0 : TEMPB=PROD(K,2)*WVEC(K,3)+PROD(K,3)*WVEC(K,2)
912 0 : DEN(5)=DEN(5)-TEMPA-HALF*TEMPB
913 0 : DEN(8)=DEN(8)-PAR(4)+PAR(5)
914 0 : TEMPA=PROD(K,4)*WVEC(K,5)+PROD(K,5)*WVEC(K,4)
915 0 : DEN(9)=DEN(9)-HALF*TEMPA
916 : end do
917 : !
918 : ! Extend DEN so that it holds all the coefficients of DENOM.
919 : !
920 : SUM=ZERO
921 0 : do I=1,5
922 0 : PAR(I)=HALF*PROD(KNEW,I)**2
923 0 : SUM=SUM+PAR(I)
924 : end do
925 0 : DENEX(1)=ALPHA*DEN(1)+PAR(1)+SUM
926 0 : TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,2)
927 0 : TEMPB=PROD(KNEW,2)*PROD(KNEW,4)
928 0 : TEMPC=PROD(KNEW,3)*PROD(KNEW,5)
929 0 : DENEX(2)=ALPHA*DEN(2)+TEMPA+TEMPB+TEMPC
930 0 : DENEX(6)=ALPHA*DEN(6)+TEMPB-TEMPC
931 0 : TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,3)
932 0 : TEMPB=PROD(KNEW,2)*PROD(KNEW,5)
933 0 : TEMPC=PROD(KNEW,3)*PROD(KNEW,4)
934 0 : DENEX(3)=ALPHA*DEN(3)+TEMPA+TEMPB-TEMPC
935 0 : DENEX(7)=ALPHA*DEN(7)+TEMPB+TEMPC
936 0 : TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,4)
937 0 : DENEX(4)=ALPHA*DEN(4)+TEMPA+PAR(2)-PAR(3)
938 0 : TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,5)
939 0 : DENEX(5)=ALPHA*DEN(5)+TEMPA+PROD(KNEW,2)*PROD(KNEW,3)
940 0 : DENEX(8)=ALPHA*DEN(8)+PAR(4)-PAR(5)
941 0 : DENEX(9)=ALPHA*DEN(9)+PROD(KNEW,4)*PROD(KNEW,5)
942 : !
943 : ! Seek the value of the angle that maximizes the modulus of DENOM.
944 : !
945 0 : SUM=DENEX(1)+DENEX(2)+DENEX(4)+DENEX(6)+DENEX(8)
946 0 : DENOLD=SUM
947 0 : DENMAX=SUM
948 0 : ISAVE=0
949 0 : IU=49
950 0 : TEMP=TWOPI/DBLE(IU+1)
951 0 : PAR(1)=ONE
952 0 : do I=1,IU
953 0 : ANGLE=DBLE(I)*TEMP
954 0 : PAR(2)=cos(ANGLE)
955 0 : PAR(3)=sin(ANGLE)
956 0 : do J=4,8,2
957 0 : PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
958 0 : PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
959 : end do
960 0 : SUMOLD=SUM
961 : SUM=ZERO
962 0 : do J=1,9
963 0 : SUM=SUM+DENEX(J)*PAR(J)
964 : end do
965 0 : if (DABS(SUM) > DABS(DENMAX)) then
966 : DENMAX=SUM
967 : ISAVE=I
968 : TEMPA=SUMOLD
969 0 : else if (I == ISAVE+1) then
970 0 : TEMPB=SUM
971 : end if
972 : end do
973 0 : if (ISAVE == 0) TEMPA=SUM
974 0 : if (ISAVE == IU) TEMPB=DENOLD
975 0 : STEP=ZERO
976 0 : if (TEMPA /= TEMPB) then
977 0 : TEMPA=TEMPA-DENMAX
978 0 : TEMPB=TEMPB-DENMAX
979 0 : STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
980 : end if
981 0 : ANGLE=TEMP*(DBLE(ISAVE)+STEP)
982 : !
983 : ! Calculate the new parameters of the denominator, the new VLAG vector
984 : ! and the new D. Then test for convergence.
985 : !
986 0 : PAR(2)=cos(ANGLE)
987 0 : PAR(3)=sin(ANGLE)
988 0 : do J=4,8,2
989 0 : PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
990 0 : PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
991 : end do
992 0 : BETA=ZERO
993 0 : DENMAX=ZERO
994 0 : do J=1,9
995 0 : BETA=BETA+DEN(J)*PAR(J)
996 0 : DENMAX=DENMAX+DENEX(J)*PAR(J)
997 : end do
998 0 : do K=1,NDIM
999 0 : VLAG(K)=ZERO
1000 0 : do J=1,5
1001 0 : VLAG(K)=VLAG(K)+PROD(K,J)*PAR(J)
1002 : end do
1003 : end do
1004 0 : TAU=VLAG(KNEW)
1005 0 : DD=ZERO
1006 0 : TEMPA=ZERO
1007 0 : TEMPB=ZERO
1008 0 : do I=1,N
1009 0 : D(I)=PAR(2)*D(I)+PAR(3)*S(I)
1010 0 : W(I)=XOPT(I)+D(I)
1011 0 : DD=DD+D(I)**2
1012 0 : TEMPA=TEMPA+D(I)*W(I)
1013 0 : TEMPB=TEMPB+W(I)*W(I)
1014 : end do
1015 0 : if (ITERC >= N) GOTO 340
1016 0 : if (ITERC > 1) DENSAV=DMAX1(DENSAV,DENOLD)
1017 0 : if (DABS(DENMAX) <= 1.1D0*DABS(DENSAV)) GOTO 340
1018 0 : DENSAV=DENMAX
1019 : !
1020 : ! Set S to half the gradient of the denominator with respect to D.
1021 : ! Then branch for the next iteration.
1022 : !
1023 0 : do I=1,N
1024 0 : TEMP=TEMPA*XOPT(I)+TEMPB*D(I)-VLAG(NPT+I)
1025 0 : S(I)=TAU*BMAT(KNEW,I)+ALPHA*TEMP
1026 : end do
1027 0 : do K=1,NPT
1028 : SUM=ZERO
1029 0 : do J=1,N
1030 0 : SUM=SUM+XPT(K,J)*W(J)
1031 : end do
1032 0 : TEMP=(TAU*W(N+K)-ALPHA*VLAG(K))*SUM
1033 0 : do I=1,N
1034 0 : S(I)=S(I)+TEMP*XPT(K,I)
1035 : end do
1036 : end do
1037 : SS=ZERO
1038 : DS=ZERO
1039 0 : do I=1,N
1040 0 : SS=SS+S(I)**2
1041 0 : DS=DS+D(I)*S(I)
1042 : end do
1043 0 : SSDEN=DD*SS-DS*DS
1044 0 : if (SSDEN >= 1.0D-8*DD*SS) GOTO 70
1045 : !
1046 : ! Set the vector W before the RETURN from the subroutine.
1047 : !
1048 0 : 340 do K=1,NDIM
1049 0 : W(K)=ZERO
1050 0 : do J=1,5
1051 0 : W(K)=W(K)+WVEC(K,J)*PAR(J)
1052 : end do
1053 : end do
1054 0 : VLAG(KOPT)=VLAG(KOPT)+ONE
1055 0 : RETURN
1056 : end subroutine BIGDEN
1057 :
1058 :
1059 120 : subroutine BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DELTA,D,ALPHA,HCOL,GC,GD,S,W)
1060 : implicit real(dp) (A-H,O-Z)
1061 : dimension XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),HCOL(*),GC(*),GD(*),S(*),W(*)
1062 : !
1063 : ! N is the number of variables.
1064 : ! NPT is the number of interpolation equations.
1065 : ! XOPT is the best interpolation point so far.
1066 : ! XPT contains the coordinates of the current interpolation points.
1067 : ! BMAT provides the last N columns of H.
1068 : ! ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H.
1069 : ! NDIM is the first dimension of BMAT and has the value NPT+N.
1070 : ! KNEW is the index of the interpolation point that is going to be moved.
1071 : ! DELTA is the current trust region bound.
1072 : ! D will be set to the step from XOPT to the new point.
1073 : ! ALPHA will be set to the KNEW-th diagonal element of the H matrix.
1074 : ! HCOL, GC, GD, S and W will be used for working space.
1075 : !
1076 : ! The step D is calculated in a way that attempts to maximize the modulus
1077 : ! of LFUNC(XOPT+D), subject to the bound ||D|| <= DELTA, where LFUNC is
1078 : ! the KNEW-th Lagrange function.
1079 : !
1080 : ! Set some constants.
1081 : !
1082 120 : HALF=0.5D0
1083 120 : ONE=1.0D0
1084 120 : ZERO=0.0D0
1085 120 : TWOPI=8.0D0*atan(ONE)
1086 120 : DELSQ=DELTA*DELTA
1087 120 : NPTM=NPT-N-1
1088 : !
1089 : ! Set the first NPT components of HCOL to the leading elements of the
1090 : ! KNEW-th column of H.
1091 : !
1092 120 : ITERC=0
1093 1424 : do K=1,NPT
1094 1424 : HCOL(K)=ZERO
1095 : end do
1096 712 : do J=1,NPTM
1097 712 : TEMP=ZMAT(KNEW,J)
1098 592 : if (J < IDZ) TEMP=-TEMP
1099 7592 : do K=1,NPT
1100 7472 : HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J)
1101 : end do
1102 : end do
1103 120 : ALPHA=HCOL(KNEW)
1104 : !
1105 : ! Set the unscaled initial direction D. Form the gradient of LFUNC at
1106 : ! XOPT, and multiply D by the second derivative matrix of LFUNC.
1107 : !
1108 240 : DD=ZERO
1109 712 : do I=1,N
1110 592 : D(I)=XPT(KNEW,I)-XOPT(I)
1111 592 : GC(I)=BMAT(KNEW,I)
1112 592 : GD(I)=ZERO
1113 712 : DD=DD+D(I)**2
1114 : end do
1115 1424 : do K=1,NPT
1116 : TEMP=ZERO
1117 120 : SUM=ZERO
1118 8184 : do J=1,N
1119 6880 : TEMP=TEMP+XPT(K,J)*XOPT(J)
1120 8184 : SUM=SUM+XPT(K,J)*D(J)
1121 : end do
1122 1304 : TEMP=HCOL(K)*TEMP
1123 1304 : SUM=HCOL(K)*SUM
1124 8304 : do I=1,N
1125 6880 : GC(I)=GC(I)+TEMP*XPT(K,I)
1126 8184 : GD(I)=GD(I)+SUM*XPT(K,I)
1127 : end do
1128 : end do
1129 : !
1130 : ! Scale D and GD, with a sign change if required. Set S to another
1131 : ! vector in the initial two dimensional subspace.
1132 : !
1133 120 : GG=ZERO
1134 120 : SP=ZERO
1135 120 : DHD=ZERO
1136 712 : do I=1,N
1137 592 : GG=GG+GC(I)**2
1138 592 : SP=SP+D(I)*GC(I)
1139 712 : DHD=DHD+D(I)*GD(I)
1140 : end do
1141 240 : SCALE=DELTA/DSQRT(DD)
1142 120 : if (SP*DHD < ZERO) SCALE=-SCALE
1143 120 : TEMP=ZERO
1144 120 : if (SP*SP > 0.99D0*DD*GG) TEMP=ONE
1145 240 : TAU=SCALE*(DABS(SP)+HALF*SCALE*DABS(DHD))
1146 120 : if (GG*DELSQ < 0.01D0*TAU*TAU) TEMP=ONE
1147 712 : do I=1,N
1148 592 : D(I)=SCALE*D(I)
1149 592 : GD(I)=SCALE*GD(I)
1150 805 : S(I)=GC(I)+TEMP*GD(I)
1151 : end do
1152 : !
1153 : ! Begin the iteration by overwriting S with a vector that has the
1154 : ! required length and direction, except that termination occurs if
1155 : ! the given D and S are nearly parallel.
1156 : !
1157 213 : 80 ITERC=ITERC+1
1158 213 : DD=ZERO
1159 213 : SP=ZERO
1160 333 : SS=ZERO
1161 1299 : do I=1,N
1162 1086 : DD=DD+D(I)**2
1163 1086 : SP=SP+D(I)*S(I)
1164 1299 : SS=SS+S(I)**2
1165 : end do
1166 213 : TEMP=DD*SS-SP*SP
1167 213 : if (TEMP <= 1.0D-8*DD*SS) GOTO 160
1168 333 : DENOM=DSQRT(TEMP)
1169 1299 : do I=1,N
1170 1086 : S(I)=(DD*S(I)-SP*D(I))/DENOM
1171 1299 : W(I)=ZERO
1172 : end do
1173 : !
1174 : ! Calculate the coefficients of the objective function on the circle,
1175 : ! beginning with the multiplication of S by the second derivative matrix.
1176 : !
1177 2598 : do K=1,NPT
1178 : SUM=ZERO
1179 15239 : do J=1,N
1180 15239 : SUM=SUM+XPT(K,J)*S(J)
1181 : end do
1182 2385 : SUM=HCOL(K)*SUM
1183 15452 : do I=1,N
1184 15239 : W(I)=W(I)+SUM*XPT(K,I)
1185 : end do
1186 : end do
1187 120 : CF1=ZERO
1188 120 : CF2=ZERO
1189 120 : CF3=ZERO
1190 120 : CF4=ZERO
1191 120 : CF5=ZERO
1192 1299 : do I=1,N
1193 1086 : CF1=CF1+S(I)*W(I)
1194 1086 : CF2=CF2+D(I)*GC(I)
1195 1086 : CF3=CF3+S(I)*GC(I)
1196 1086 : CF4=CF4+D(I)*GD(I)
1197 1299 : CF5=CF5+S(I)*GD(I)
1198 : end do
1199 213 : CF1=HALF*CF1
1200 213 : CF4=HALF*CF4-CF1
1201 : !
1202 : ! Seek the value of the angle that maximizes the modulus of TAU.
1203 : !
1204 333 : TAUBEG=CF1+CF2+CF4
1205 333 : TAUMAX=TAUBEG
1206 333 : TAUOLD=TAUBEG
1207 213 : ISAVE=0
1208 213 : IU=49
1209 213 : TEMP=TWOPI/DBLE(IU+1)
1210 10650 : do I=1,IU
1211 10557 : ANGLE=DBLE(I)*TEMP
1212 10557 : CTH=cos(ANGLE)
1213 10557 : STH=sin(ANGLE)
1214 10437 : TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
1215 10437 : if (DABS(TAU) > DABS(TAUMAX)) then
1216 : TAUMAX=TAU
1217 : ISAVE=I
1218 120 : TEMPA=TAUOLD
1219 9621 : else if (I == ISAVE+1) then
1220 235 : TEMPB=TAU
1221 : end if
1222 10650 : TAUOLD=TAU
1223 : end do
1224 213 : if (ISAVE == 0) TEMPA=TAU
1225 184 : if (ISAVE == IU) TEMPB=TAUBEG
1226 120 : STEP=ZERO
1227 213 : if (TEMPA /= TEMPB) then
1228 213 : TEMPA=TEMPA-TAUMAX
1229 213 : TEMPB=TEMPB-TAUMAX
1230 213 : STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
1231 : end if
1232 213 : ANGLE=TEMP*(DBLE(ISAVE)+STEP)
1233 : !
1234 : ! Calculate the new D and GD. Then test for convergence.
1235 : !
1236 213 : CTH=cos(ANGLE)
1237 213 : STH=sin(ANGLE)
1238 213 : TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
1239 1299 : do I=1,N
1240 1086 : D(I)=CTH*D(I)+STH*S(I)
1241 1086 : GD(I)=CTH*GD(I)+STH*W(I)
1242 1299 : S(I)=GC(I)+GD(I)
1243 : end do
1244 213 : if (DABS(TAU) <= 1.1D0*DABS(TAUBEG)) GOTO 160
1245 93 : if (ITERC < N) GOTO 80
1246 120 : 160 RETURN
1247 : end subroutine BIGLAG
1248 :
1249 252 : subroutine TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,STEP,D,G,HD,HS,CRVMIN)
1250 : implicit real(dp) (A-H,O-Z)
1251 : dimension XOPT(*),XPT(NPT,*),GQ(*),HQ(*),PQ(*),STEP(*),D(*),G(*),HD(*),HS(*)
1252 : !
1253 : ! N is the number of variables of a quadratic objective function, Q say.
1254 : ! The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meanings,
1255 : ! in order to define the current quadratic model Q.
1256 : ! DELTA is the trust region radius, and has to be positive.
1257 : ! STEP will be set to the calculated trial step.
1258 : ! The arrays D, G, HD and HS will be used for working space.
1259 : ! CRVMIN will be set to the least curvature of H along the conjugate
1260 : ! directions that occur, except that it is set to zero if STEP goes
1261 : ! all the way to the trust region boundary.
1262 : !
1263 : ! The calculation of STEP begins with the truncated conjugate gradient
1264 : ! method. If the boundary of the trust region is reached, then further
1265 : ! changes to STEP may be made, each one being in the 2D space spanned
1266 : ! by the current STEP and the corresponding gradient of Q. Thus STEP
1267 : ! should provide a substantial reduction to Q within the trust region.
1268 : !
1269 : ! Initialization, which includes setting HD to H times XOPT.
1270 : !
1271 252 : HALF=0.5D0
1272 252 : ZERO=0.0D0
1273 252 : TWOPI=8.0D0*atan(1.0D0)
1274 252 : DELSQ=DELTA*DELTA
1275 252 : ITERC=0
1276 252 : ITERMAX=N
1277 252 : ITERSW=ITERMAX
1278 1492 : do I=1,N
1279 1492 : D(I)=XOPT(I)
1280 : end do
1281 252 : GOTO 170
1282 : !
1283 : ! Prepare for the first line search.
1284 : !
1285 504 : 20 QRED=ZERO
1286 504 : DD=ZERO
1287 1492 : do I=1,N
1288 1240 : STEP(I)=ZERO
1289 1240 : HS(I)=ZERO
1290 1240 : G(I)=GQ(I)+HD(I)
1291 1240 : D(I)=-G(I)
1292 1492 : DD=DD+D(I)**2
1293 : end do
1294 252 : CRVMIN=ZERO
1295 252 : if (DD == ZERO) GOTO 160
1296 252 : DS=ZERO
1297 252 : SS=ZERO
1298 252 : GG=DD
1299 629 : GGBEG=GG
1300 : !
1301 : ! Calculate the step to the trust region boundary and the product HD.
1302 : !
1303 881 : 40 ITERC=ITERC+1
1304 1133 : TEMP=DELSQ-SS
1305 1133 : BSTEP=TEMP/(DS+DSQRT(DS*DS+DD*TEMP))
1306 1762 : GOTO 170
1307 1133 : 50 DHD=ZERO
1308 5445 : do J=1,N
1309 5445 : DHD=DHD+D(J)*HD(J)
1310 : end do
1311 : !
1312 : ! Update CRVMIN and set the step-length ALPHA.
1313 : !
1314 1133 : ALPHA=BSTEP
1315 881 : if (DHD > ZERO) then
1316 879 : TEMP=DHD/DD
1317 879 : if (ITERC == 1) CRVMIN=TEMP
1318 879 : CRVMIN=DMIN1(CRVMIN,TEMP)
1319 879 : ALPHA=DMIN1(ALPHA,GG/DHD)
1320 : end if
1321 1133 : QADD=ALPHA*(GG-HALF*ALPHA*DHD)
1322 881 : QRED=QRED+QADD
1323 : !
1324 : ! Update STEP and HS.
1325 : !
1326 1133 : GGSAV=GG
1327 881 : GG=ZERO
1328 5445 : do I=1,N
1329 4564 : STEP(I)=STEP(I)+ALPHA*D(I)
1330 4564 : HS(I)=HS(I)+ALPHA*HD(I)
1331 5445 : GG=GG+(G(I)+HS(I))**2
1332 : end do
1333 : !
1334 : ! Begin another conjugate direction iteration if required.
1335 : !
1336 881 : if (ALPHA < BSTEP) then
1337 782 : if (QADD <= 0.01D0*QRED) GOTO 160
1338 704 : if (GG <= 1.0D-4*GGBEG) GOTO 160
1339 629 : if (ITERC == ITERMAX) GOTO 160
1340 629 : TEMP=GG/GGSAV
1341 629 : DD=ZERO
1342 629 : DS=ZERO
1343 629 : SS=ZERO
1344 3953 : do I=1,N
1345 3324 : D(I)=TEMP*D(I)-G(I)-HS(I)
1346 3324 : DD=DD+D(I)**2
1347 3324 : DS=DS+D(I)*STEP(I)
1348 3953 : SS=SS+STEP(I)**2
1349 : end do
1350 629 : if (DS <= ZERO) GOTO 160
1351 629 : if (SS < DELSQ) GOTO 40
1352 : end if
1353 99 : CRVMIN=ZERO
1354 99 : ITERSW=ITERC
1355 : !
1356 : ! Test whether an alternative iteration is required.
1357 : !
1358 204 : 90 if (GG <= 1.0D-4*GGBEG) GOTO 160
1359 456 : SG=ZERO
1360 456 : SHS=ZERO
1361 1336 : do I=1,N
1362 1132 : SG=SG+STEP(I)*G(I)
1363 1336 : SHS=SHS+STEP(I)*HS(I)
1364 : end do
1365 456 : SGK=SG+SHS
1366 456 : ANGTEST=SGK/DSQRT(GG*DELSQ)
1367 204 : if (ANGTEST <= -0.99D0) GOTO 160
1368 : !
1369 : ! Begin the alternative iteration by calculating D and HD and some
1370 : ! scalar products.
1371 : !
1372 186 : ITERC=ITERC+1
1373 186 : TEMP=DSQRT(DELSQ*GG-SGK*SGK)
1374 438 : TEMPA=DELSQ/TEMP
1375 438 : TEMPB=SGK/TEMP
1376 1226 : do I=1,N
1377 1226 : D(I)=TEMPA*(G(I)+HS(I))-TEMPB*STEP(I)
1378 : end do
1379 : GOTO 170
1380 252 : 120 DG=ZERO
1381 : DHD=ZERO
1382 252 : DHS=ZERO
1383 1226 : do I=1,N
1384 1040 : DG=DG+D(I)*G(I)
1385 1040 : DHD=DHD+HD(I)*D(I)
1386 1226 : DHS=DHS+HD(I)*STEP(I)
1387 : end do
1388 : !
1389 : ! Seek the value of the angle that minimizes Q.
1390 : !
1391 438 : CF=HALF*(SHS-DHD)
1392 438 : QBEG=SG+CF
1393 438 : QSAV=QBEG
1394 438 : QMIN=QBEG
1395 186 : ISAVE=0
1396 186 : IU=49
1397 186 : TEMP=TWOPI/DBLE(IU+1)
1398 9300 : do I=1,IU
1399 9366 : ANGLE=DBLE(I)*TEMP
1400 9366 : CTH=cos(ANGLE)
1401 9366 : STH=sin(ANGLE)
1402 9366 : QNEW=(SG+CF*CTH)*CTH+(DG+DHS*CTH)*STH
1403 9114 : if (QNEW < QMIN) then
1404 : QMIN=QNEW
1405 : ISAVE=I
1406 : TEMPA=QSAV
1407 8886 : else if (I == ISAVE+1) then
1408 233 : TEMPB=QNEW
1409 : end if
1410 9300 : QSAV=QNEW
1411 : end do
1412 186 : if (ISAVE == ZERO) TEMPA=QNEW
1413 131 : if (ISAVE == IU) TEMPB=QBEG
1414 186 : ANGLE=ZERO
1415 186 : if (TEMPA /= TEMPB) then
1416 186 : TEMPA=TEMPA-QMIN
1417 186 : TEMPB=TEMPB-QMIN
1418 186 : ANGLE=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
1419 : end if
1420 186 : ANGLE=TEMP*(DBLE(ISAVE)+ANGLE)
1421 : !
1422 : ! Calculate the new STEP and HS. Then test for convergence.
1423 : !
1424 186 : CTH=cos(ANGLE)
1425 186 : STH=sin(ANGLE)
1426 438 : REDUC=QBEG-(SG+CF*CTH)*CTH-(DG+DHS*CTH)*STH
1427 186 : GG=ZERO
1428 1226 : do I=1,N
1429 1040 : STEP(I)=CTH*STEP(I)+STH*D(I)
1430 1040 : HS(I)=CTH*HS(I)+STH*HD(I)
1431 1226 : GG=GG+(G(I)+HS(I))**2
1432 : end do
1433 186 : QRED=QRED+REDUC
1434 252 : RATIO=REDUC/QRED
1435 186 : if (ITERC < ITERMAX .and. RATIO > 0.01D0) GOTO 90
1436 252 : 160 RETURN
1437 : !
1438 : ! The following instructions act as a subroutine for setting the vector
1439 : ! HD to the vector D multiplied by the second derivative matrix of Q.
1440 : ! They are called from three different places, which are distinguished
1441 : ! by the value of ITERC.
1442 : !
1443 8163 : 170 do I=1,N
1444 8163 : HD(I)=ZERO
1445 : end do
1446 16326 : do K=1,NPT
1447 : TEMP=ZERO
1448 96859 : do J=1,N
1449 96859 : TEMP=TEMP+XPT(K,J)*D(J)
1450 : end do
1451 15007 : TEMP=TEMP*PQ(K)
1452 98178 : do I=1,N
1453 96859 : HD(I)=HD(I)+TEMP*XPT(K,I)
1454 : end do
1455 : end do
1456 : IH=0
1457 8163 : do J=1,N
1458 30337 : do I=1,J
1459 22174 : IH=IH+1
1460 22174 : if (I < J) HD(J)=HD(J)+HQ(IH)*D(I)
1461 29018 : HD(I)=HD(I)+HQ(IH)*D(J)
1462 : end do
1463 : end do
1464 1319 : if (ITERC == 0) GOTO 20
1465 1067 : if (ITERC <= ITERSW) GOTO 50
1466 : GOTO 120
1467 : end subroutine TRSAPP
1468 :
1469 281 : subroutine UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
1470 : implicit real(dp) (A-H,O-Z)
1471 : dimension BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*)
1472 : !
1473 : ! The arrays BMAT and ZMAT with IDZ are updated, in order to shift the
1474 : ! interpolation point that has index KNEW. On entry, VLAG contains the
1475 : ! components of the vector Theta*Wcheck+e_b of the updating formula
1476 : ! (6.11), and BETA holds the value of the parameter that has this name.
1477 : ! The vector W is used for working space.
1478 : !
1479 : ! Set some constants.
1480 : !
1481 281 : ONE=1.0D0
1482 281 : ZERO=0.0D0
1483 281 : NPTM=NPT-N-1
1484 : !
1485 : ! Apply the rotations that put zeros in the KNEW-th row of ZMAT.
1486 : !
1487 281 : JL=1
1488 1424 : do J=2,NPTM
1489 1424 : if (J == IDZ) then
1490 : JL=IDZ
1491 1143 : else if (ZMAT(KNEW,J) /= ZERO) then
1492 1410 : TEMP=DSQRT(ZMAT(KNEW,JL)**2+ZMAT(KNEW,J)**2)
1493 1410 : TEMPA=ZMAT(KNEW,JL)/TEMP
1494 1410 : TEMPB=ZMAT(KNEW,J)/TEMP
1495 14654 : do I=1,NPT
1496 13525 : TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J)
1497 13525 : ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL)
1498 14654 : ZMAT(I,JL)=TEMP
1499 : end do
1500 1129 : ZMAT(KNEW,J)=ZERO
1501 : end if
1502 : end do
1503 : !
1504 : ! Put the first NPT components of the KNEW-th column of HLAG into W,
1505 : ! and calculate the parameters of the updating formula.
1506 : !
1507 281 : TEMPA=ZMAT(KNEW,1)
1508 281 : if (IDZ >= 2) TEMPA=-TEMPA
1509 281 : if (JL > 1) TEMPB=ZMAT(KNEW,JL)
1510 3410 : do I=1,NPT
1511 3129 : W(I)=TEMPA*ZMAT(I,1)
1512 3410 : if (JL > 1) W(I)=W(I)+TEMPB*ZMAT(I,JL)
1513 : end do
1514 562 : ALPHA=W(KNEW)
1515 562 : TAU=VLAG(KNEW)
1516 562 : TAUSQ=TAU*TAU
1517 562 : DENOM=ALPHA*BETA+TAUSQ
1518 281 : VLAG(KNEW)=VLAG(KNEW)-ONE
1519 : !
1520 : ! Complete the updating of ZMAT when there is only one nonzero element
1521 : ! in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one,
1522 : ! then the first column of ZMAT will be exchanged with another one later.
1523 : !
1524 281 : IFLAG=0
1525 281 : if (JL == 1) then
1526 281 : TEMP=DSQRT(DABS(DENOM))
1527 281 : TEMPB=TEMPA/TEMP
1528 281 : TEMPA=TAU/TEMP
1529 3410 : do I=1,NPT
1530 3410 : ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
1531 : end do
1532 : if (IDZ == 1 .and. TEMP < ZERO) IDZ=2
1533 281 : if (IDZ >= 2 .and. TEMP >= ZERO) IFLAG=1
1534 : else
1535 : !
1536 : ! Complete the updating of ZMAT in the alternative case.
1537 : !
1538 0 : JA=1
1539 0 : if (BETA >= ZERO) JA=JL
1540 0 : JB=JL+1-JA
1541 0 : TEMP=ZMAT(KNEW,JB)/DENOM
1542 0 : TEMPA=TEMP*BETA
1543 0 : TEMPB=TEMP*TAU
1544 0 : TEMP=ZMAT(KNEW,JA)
1545 281 : SCALA=ONE/DSQRT(DABS(BETA)*TEMP*TEMP+TAUSQ)
1546 281 : SCALB=SCALA*DSQRT(DABS(DENOM))
1547 0 : do I=1,NPT
1548 0 : ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I))
1549 0 : ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I))
1550 : end do
1551 0 : if (DENOM <= ZERO) then
1552 0 : if (BETA < ZERO) IDZ=IDZ+1
1553 0 : if (BETA >= ZERO) IFLAG=1
1554 : end if
1555 : end if
1556 : !
1557 : ! IDZ is reduced in the following case, and usually the first column
1558 : ! of ZMAT is exchanged with a later one.
1559 : !
1560 : if (IFLAG == 1) then
1561 0 : IDZ=IDZ-1
1562 0 : do I=1,NPT
1563 0 : TEMP=ZMAT(I,1)
1564 0 : ZMAT(I,1)=ZMAT(I,IDZ)
1565 0 : ZMAT(I,IDZ)=TEMP
1566 : end do
1567 : end if
1568 : !
1569 : ! Finally, update the matrix BMAT.
1570 : !
1571 1705 : do J=1,N
1572 1424 : JP=NPT+J
1573 1424 : W(JP)=BMAT(KNEW,J)
1574 1424 : TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM
1575 1424 : TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM
1576 23081 : do I=1,JP
1577 21376 : BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I)
1578 22800 : if (I > NPT) BMAT(JP,I-NPT)=BMAT(I,J)
1579 : end do
1580 : end do
1581 281 : end subroutine UPDATE
1582 :
1583 : end module mod_newuoa
|
