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