1      SUBROUTINE BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT,
2     1  KNEW,D,W,VLAG,BETA,S,WVEC,PROD)
3      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4      DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),
5     1  W(*),VLAG(*),S(*),WVEC(NDIM,*),PROD(NDIM,*)
6      DIMENSION DEN(9),DENEX(9),PAR(9)
7C
8C     N is the number of variables.
9C     NPT is the number of interpolation equations.
10C     XOPT is the best interpolation point so far.
11C     XPT contains the coordinates of the current interpolation points.
12C     BMAT provides the last N columns of H.
13C     ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H.
14C     NDIM is the first dimension of BMAT and has the value NPT+N.
15C     KOPT is the index of the optimal interpolation point.
16C     KNEW is the index of the interpolation point that is going to be moved.
17C     D will be set to the step from XOPT to the new point, and on entry it
18C       should be the D that was calculated by the last call of BIGLAG. The
19C       length of the initial D provides a trust region bound on the final D.
20C     W will be set to Wcheck for the final choice of D.
21C     VLAG will be set to Theta*Wcheck+e_b for the final choice of D.
22C     BETA will be set to the value that will occur in the updating formula
23C       when the KNEW-th interpolation point is moved to its new position.
24C     S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used
25C       for working space.
26C
27C     D is calculated in a way that should provide a denominator with a large
28C     modulus in the updating formula when the KNEW-th interpolation point is
29C     shifted to the new position XOPT+D.
30C
31C     Set some constants.
32C
33      HALF=0.5D0
34      ONE=1.0D0
35      QUART=0.25D0
36      TWO=2.0D0
37      ZERO=0.0D0
38      TWOPI=8.0D0*DATAN(ONE)
39      NPTM=NPT-N-1
40C
41C     Store the first NPT elements of the KNEW-th column of H in W(N+1)
42C     to W(N+NPT).
43C
44      DO 10 K=1,NPT
45   10 W(N+K)=ZERO
46      DO 20 J=1,NPTM
47      TEMP=ZMAT(KNEW,J)
48      IF (J .LT. IDZ) TEMP=-TEMP
49      DO 20 K=1,NPT
50   20 W(N+K)=W(N+K)+TEMP*ZMAT(K,J)
51      ALPHA=W(N+KNEW)
52C
53C     The initial search direction D is taken from the last call of BIGLAG,
54C     and the initial S is set below, usually to the direction from X_OPT
55C     to X_KNEW, but a different direction to an interpolation point may
56C     be chosen, in order to prevent S from being nearly parallel to D.
57C
58      DD=ZERO
59      DS=ZERO
60      SS=ZERO
61      XOPTSQ=ZERO
62      DO 30 I=1,N
63      DD=DD+D(I)**2
64      S(I)=XPT(KNEW,I)-XOPT(I)
65      DS=DS+D(I)*S(I)
66      SS=SS+S(I)**2
67   30 XOPTSQ=XOPTSQ+XOPT(I)**2
68      IF (DS*DS .GT. 0.99D0*DD*SS) THEN
69          KSAV=KNEW
70          DTEST=DS*DS/SS
71          DO 50 K=1,NPT
72          IF (K .NE. KOPT) THEN
73              DSTEMP=ZERO
74              SSTEMP=ZERO
75              DO 40 I=1,N
76              DIFF=XPT(K,I)-XOPT(I)
77              DSTEMP=DSTEMP+D(I)*DIFF
78   40         SSTEMP=SSTEMP+DIFF*DIFF
79              IF (DSTEMP*DSTEMP/SSTEMP .LT. DTEST) THEN
80                  KSAV=K
81                  DTEST=DSTEMP*DSTEMP/SSTEMP
82                  DS=DSTEMP
83                  SS=SSTEMP
84              END IF
85          END IF
86   50     CONTINUE
87          DO 60 I=1,N
88   60     S(I)=XPT(KSAV,I)-XOPT(I)
89      END IF
90      SSDEN=DD*SS-DS*DS
91      ITERC=0
92      DENSAV=ZERO
93C
94C     Begin the iteration by overwriting S with a vector that has the
95C     required length and direction.
96C
97   70 ITERC=ITERC+1
98      TEMP=ONE/DSQRT(SSDEN)
99      XOPTD=ZERO
100      XOPTS=ZERO
101      DO 80 I=1,N
102      S(I)=TEMP*(DD*S(I)-DS*D(I))
103      XOPTD=XOPTD+XOPT(I)*D(I)
104   80 XOPTS=XOPTS+XOPT(I)*S(I)
105C
106C     Set the coefficients of the first two terms of BETA.
107C
108      TEMPA=HALF*XOPTD*XOPTD
109      TEMPB=HALF*XOPTS*XOPTS
110      DEN(1)=DD*(XOPTSQ+HALF*DD)+TEMPA+TEMPB
111      DEN(2)=TWO*XOPTD*DD
112      DEN(3)=TWO*XOPTS*DD
113      DEN(4)=TEMPA-TEMPB
114      DEN(5)=XOPTD*XOPTS
115      DO 90 I=6,9
116   90 DEN(I)=ZERO
117C
118C     Put the coefficients of Wcheck in WVEC.
119C
120      DO 110 K=1,NPT
121      TEMPA=ZERO
122      TEMPB=ZERO
123      TEMPC=ZERO
124      DO 100 I=1,N
125      TEMPA=TEMPA+XPT(K,I)*D(I)
126      TEMPB=TEMPB+XPT(K,I)*S(I)
127  100 TEMPC=TEMPC+XPT(K,I)*XOPT(I)
128      WVEC(K,1)=QUART*(TEMPA*TEMPA+TEMPB*TEMPB)
129      WVEC(K,2)=TEMPA*TEMPC
130      WVEC(K,3)=TEMPB*TEMPC
131      WVEC(K,4)=QUART*(TEMPA*TEMPA-TEMPB*TEMPB)
132  110 WVEC(K,5)=HALF*TEMPA*TEMPB
133      DO 120 I=1,N
134      IP=I+NPT
135      WVEC(IP,1)=ZERO
136      WVEC(IP,2)=D(I)
137      WVEC(IP,3)=S(I)
138      WVEC(IP,4)=ZERO
139  120 WVEC(IP,5)=ZERO
140C
141C     Put the coefficents of THETA*Wcheck in PROD.
142C
143      DO 190 JC=1,5
144      NW=NPT
145      IF (JC .EQ. 2 .OR. JC .EQ. 3) NW=NDIM
146      DO 130 K=1,NPT
147  130 PROD(K,JC)=ZERO
148      DO 150 J=1,NPTM
149      SUM=ZERO
150      DO 140 K=1,NPT
151  140 SUM=SUM+ZMAT(K,J)*WVEC(K,JC)
152      IF (J .LT. IDZ) SUM=-SUM
153      DO 150 K=1,NPT
154  150 PROD(K,JC)=PROD(K,JC)+SUM*ZMAT(K,J)
155      IF (NW .EQ. NDIM) THEN
156          DO 170 K=1,NPT
157          SUM=ZERO
158          DO 160 J=1,N
159  160     SUM=SUM+BMAT(K,J)*WVEC(NPT+J,JC)
160  170     PROD(K,JC)=PROD(K,JC)+SUM
161      END IF
162      DO 190 J=1,N
163      SUM=ZERO
164      DO 180 I=1,NW
165  180 SUM=SUM+BMAT(I,J)*WVEC(I,JC)
166  190 PROD(NPT+J,JC)=SUM
167C
168C     Include in DEN the part of BETA that depends on THETA.
169C
170      DO 210 K=1,NDIM
171      SUM=ZERO
172      DO 200 I=1,5
173      PAR(I)=HALF*PROD(K,I)*WVEC(K,I)
174  200 SUM=SUM+PAR(I)
175      DEN(1)=DEN(1)-PAR(1)-SUM
176      TEMPA=PROD(K,1)*WVEC(K,2)+PROD(K,2)*WVEC(K,1)
177      TEMPB=PROD(K,2)*WVEC(K,4)+PROD(K,4)*WVEC(K,2)
178      TEMPC=PROD(K,3)*WVEC(K,5)+PROD(K,5)*WVEC(K,3)
179      DEN(2)=DEN(2)-TEMPA-HALF*(TEMPB+TEMPC)
180      DEN(6)=DEN(6)-HALF*(TEMPB-TEMPC)
181      TEMPA=PROD(K,1)*WVEC(K,3)+PROD(K,3)*WVEC(K,1)
182      TEMPB=PROD(K,2)*WVEC(K,5)+PROD(K,5)*WVEC(K,2)
183      TEMPC=PROD(K,3)*WVEC(K,4)+PROD(K,4)*WVEC(K,3)
184      DEN(3)=DEN(3)-TEMPA-HALF*(TEMPB-TEMPC)
185      DEN(7)=DEN(7)-HALF*(TEMPB+TEMPC)
186      TEMPA=PROD(K,1)*WVEC(K,4)+PROD(K,4)*WVEC(K,1)
187      DEN(4)=DEN(4)-TEMPA-PAR(2)+PAR(3)
188      TEMPA=PROD(K,1)*WVEC(K,5)+PROD(K,5)*WVEC(K,1)
189      TEMPB=PROD(K,2)*WVEC(K,3)+PROD(K,3)*WVEC(K,2)
190      DEN(5)=DEN(5)-TEMPA-HALF*TEMPB
191      DEN(8)=DEN(8)-PAR(4)+PAR(5)
192      TEMPA=PROD(K,4)*WVEC(K,5)+PROD(K,5)*WVEC(K,4)
193  210 DEN(9)=DEN(9)-HALF*TEMPA
194C
195C     Extend DEN so that it holds all the coefficients of DENOM.
196C
197      SUM=ZERO
198      DO 220 I=1,5
199      PAR(I)=HALF*PROD(KNEW,I)**2
200  220 SUM=SUM+PAR(I)
201      DENEX(1)=ALPHA*DEN(1)+PAR(1)+SUM
202      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,2)
203      TEMPB=PROD(KNEW,2)*PROD(KNEW,4)
204      TEMPC=PROD(KNEW,3)*PROD(KNEW,5)
205      DENEX(2)=ALPHA*DEN(2)+TEMPA+TEMPB+TEMPC
206      DENEX(6)=ALPHA*DEN(6)+TEMPB-TEMPC
207      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,3)
208      TEMPB=PROD(KNEW,2)*PROD(KNEW,5)
209      TEMPC=PROD(KNEW,3)*PROD(KNEW,4)
210      DENEX(3)=ALPHA*DEN(3)+TEMPA+TEMPB-TEMPC
211      DENEX(7)=ALPHA*DEN(7)+TEMPB+TEMPC
212      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,4)
213      DENEX(4)=ALPHA*DEN(4)+TEMPA+PAR(2)-PAR(3)
214      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,5)
215      DENEX(5)=ALPHA*DEN(5)+TEMPA+PROD(KNEW,2)*PROD(KNEW,3)
216      DENEX(8)=ALPHA*DEN(8)+PAR(4)-PAR(5)
217      DENEX(9)=ALPHA*DEN(9)+PROD(KNEW,4)*PROD(KNEW,5)
218C
219C     Seek the value of the angle that maximizes the modulus of DENOM.
220C
221      SUM=DENEX(1)+DENEX(2)+DENEX(4)+DENEX(6)+DENEX(8)
222      DENOLD=SUM
223      DENMAX=SUM
224      ISAVE=0
225      IU=49
226      TEMP=TWOPI/DFLOAT(IU+1)
227      PAR(1)=ONE
228      DO 250 I=1,IU
229      ANGLE=DFLOAT(I)*TEMP
230      PAR(2)=DCOS(ANGLE)
231      PAR(3)=DSIN(ANGLE)
232      DO 230 J=4,8,2
233      PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
234  230 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
235      SUMOLD=SUM
236      SUM=ZERO
237      DO 240 J=1,9
238  240 SUM=SUM+DENEX(J)*PAR(J)
239      IF (DABS(SUM) .GT. DABS(DENMAX)) THEN
240          DENMAX=SUM
241          ISAVE=I
242          TEMPA=SUMOLD
243      ELSE IF (I .EQ. ISAVE+1) THEN
244          TEMPB=SUM
245      END IF
246  250 CONTINUE
247      IF (ISAVE .EQ. 0) TEMPA=SUM
248      IF (ISAVE .EQ. IU) TEMPB=DENOLD
249      STEP=ZERO
250      IF (TEMPA .NE. TEMPB) THEN
251          TEMPA=TEMPA-DENMAX
252          TEMPB=TEMPB-DENMAX
253          STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
254      END IF
255      ANGLE=TEMP*(DFLOAT(ISAVE)+STEP)
256C
257C     Calculate the new parameters of the denominator, the new VLAG vector
258C     and the new D. Then test for convergence.
259C
260      PAR(2)=DCOS(ANGLE)
261      PAR(3)=DSIN(ANGLE)
262      DO 260 J=4,8,2
263      PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
264  260 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
265      BETA=ZERO
266      DENMAX=ZERO
267      DO 270 J=1,9
268      BETA=BETA+DEN(J)*PAR(J)
269  270 DENMAX=DENMAX+DENEX(J)*PAR(J)
270      DO 280 K=1,NDIM
271      VLAG(K)=ZERO
272      DO 280 J=1,5
273  280 VLAG(K)=VLAG(K)+PROD(K,J)*PAR(J)
274      TAU=VLAG(KNEW)
275      DD=ZERO
276      TEMPA=ZERO
277      TEMPB=ZERO
278      DO 290 I=1,N
279      D(I)=PAR(2)*D(I)+PAR(3)*S(I)
280      W(I)=XOPT(I)+D(I)
281      DD=DD+D(I)**2
282      TEMPA=TEMPA+D(I)*W(I)
283  290 TEMPB=TEMPB+W(I)*W(I)
284      IF (ITERC .GE. N) GOTO 340
285      IF (ITERC .GT. 1) DENSAV=DMAX1(DENSAV,DENOLD)
286      IF (DABS(DENMAX) .LE. 1.1D0*DABS(DENSAV)) GOTO 340
287      DENSAV=DENMAX
288C
289C     Set S to half the gradient of the denominator with respect to D.
290C     Then branch for the next iteration.
291C
292      DO 300 I=1,N
293      TEMP=TEMPA*XOPT(I)+TEMPB*D(I)-VLAG(NPT+I)
294  300 S(I)=TAU*BMAT(KNEW,I)+ALPHA*TEMP
295      DO 320 K=1,NPT
296      SUM=ZERO
297      DO 310 J=1,N
298  310 SUM=SUM+XPT(K,J)*W(J)
299      TEMP=(TAU*W(N+K)-ALPHA*VLAG(K))*SUM
300      DO 320 I=1,N
301  320 S(I)=S(I)+TEMP*XPT(K,I)
302      SS=ZERO
303      DS=ZERO
304      DO 330 I=1,N
305      SS=SS+S(I)**2
306  330 DS=DS+D(I)*S(I)
307      SSDEN=DD*SS-DS*DS
308      IF (SSDEN .GE. 1.0D-8*DD*SS) GOTO 70
309C
310C     Set the vector W before the RETURN from the subroutine.
311C
312  340 DO 350 K=1,NDIM
313      W(K)=ZERO
314      DO 350 J=1,5
315  350 W(K)=W(K)+WVEC(K,J)*PAR(J)
316      VLAG(KOPT)=VLAG(KOPT)+ONE
317      RETURN
318      END
319
320