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