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