1!----------------------------------------------------------------------- 2SUBROUTINE WRITE_DOCUMENTATION() 3 write(6,'(A)') & 4" scriptmini : a fortran program to minimise a script ", & 5" or external program. ", & 6" ", & 7" Written by Joost VandeVondele ", & 8" ", & 9" the script is treated as a black box, that given an ", & 10" input vector x returns the function value f. ", & 11" ", & 12" usage : ", & 13" ", & 14" ./scriptmini ", & 15" ", & 16" inputs : ", & 17" ", & 18" -) the script should be called scriptmini_eval ", & 19" -) the script should read the n real variables ", & 20" from scriptmini_eval.in ", & 21" -) the script should write the function value (1 real number) ", & 22" to 'scriptmini_eval.out' ", & 23" -) input of scriptmini is 'scriptmini.in' with the format: ", & 24" N ", & 25" rhobeg rhoend ", & 26" maxfun ", & 27" iprint ", & 28" x[1] x[2] x[3] ... x[N] ", & 29" ", & 30" where: ", & 31" N : integer : is the number of variables ", & 32" rhobeg : real : initial trust region radius, +- 10% of the ", & 33" largest expected change in the variables ", & 34" rhoend : real : final trust region radius, +- the final ", & 35" uncertainty in the variables ", & 36" maxfun : integer : the maximum number of calls to ", & 37" scriptmini_eval [O(10*N**2)] ", & 38" iprint : integer : output level (0-3) ", & 39" 0 : no output at all (!) ", & 40" ... ", & 41" 3 : info at every step (recommended) ", & 42" x[...] : real : the initial values of the variables " 43END SUBROUTINE WRITE_DOCUMENTATION 44!------------------------------------------------------------------------------- 45 46MODULE Powell_Optimize 47 48! Code converted using TO_F90 by Alan Miller 49! Date: 2002-11-09 Time: 16:58:08 50 51IMPLICIT NONE 52INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) 53 54PRIVATE 55PUBLIC :: uobyqa 56 57 58CONTAINS 59 60 61!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uobyqa.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 62 63SUBROUTINE uobyqa(n, x, rhobeg, rhoend, iprint, maxfun) 64 65INTEGER, INTENT(IN) :: n 66REAL (dp), INTENT(IN OUT) :: x(:) 67REAL (dp), INTENT(IN) :: rhobeg 68REAL (dp), INTENT(IN) :: rhoend 69INTEGER, INTENT(IN) :: iprint 70INTEGER, INTENT(IN) :: maxfun 71 72! This subroutine seeks the least value of a function of many variables, 73! by a trust region method that forms quadratic models by interpolation. 74! The algorithm is described in "UOBYQA: unconstrained optimization by 75! quadratic approximation" by M.J.D. Powell, Report DAMTP 2000/NA14, 76! University of Cambridge. The arguments of the subroutine are as follows. 77 78! N must be set to the number of variables and must be at least two. 79! Initial values of the variables must be set in X(1),X(2),...,X(N). They 80! will be changed to the values that give the least calculated F. 81! RHOBEG and RHOEND must be set to the initial and final values of a trust 82! region radius, so both must be positive with RHOEND<=RHOBEG. Typically 83! RHOBEG should be about one tenth of the greatest expected change to a 84! variable, and RHOEND should indicate the accuracy that is required in 85! the final values of the variables. 86! The value of IPRINT should be set to 0, 1, 2 or 3, which controls the 87! amount of printing. Specifically, there is no output if IPRINT=0 and 88! there is output only at the return if IPRINT=1. Otherwise, each new 89! value of RHO is printed, with the best vector of variables so far and 90! the corresponding value of the objective function. Further, each new 91! value of F with its variables are output if IPRINT=3. 92! MAXFUN must be set to an upper bound on the number of calls of CALFUN. 93! The array W will be used for working space. Its length must be at least 94! ( N**4 + 8*N**3 + 23*N**2 + 42*N + max [ 2*N**2 + 4, 18*N ] ) / 4. 95 96! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to 97! the value of the objective function for the variables X(1),X(2),...,X(N). 98 99INTEGER :: npt 100 101! Partition the working space array, so that different parts of it can be 102! treated separately by the subroutine that performs the main calculation. 103 104npt = (n*n + 3*n + 2) / 2 105CALL uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt) 106RETURN 107END SUBROUTINE uobyqa 108 109!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uobyqb.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 110 111SUBROUTINE uobyqb(n, x, rhobeg, rhoend, iprint, maxfun, npt) 112 113INTEGER, INTENT(IN) :: n 114REAL (dp), INTENT(IN OUT) :: x(:) 115REAL (dp), INTENT(IN) :: rhobeg 116REAL (dp), INTENT(IN) :: rhoend 117INTEGER, INTENT(IN) :: iprint 118INTEGER, INTENT(IN) :: maxfun 119INTEGER, INTENT(IN) :: npt 120 121INTERFACE 122 SUBROUTINE calfun(n, x, f) 123 IMPLICIT NONE 124 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) 125 INTEGER, INTENT(IN) :: n 126 REAL (dp), INTENT(IN) :: x(:) 127 REAL (dp), INTENT(OUT) :: f 128 END SUBROUTINE calfun 129END INTERFACE 130 131! The following arrays were previously passed as arguments: 132 133REAL (dp) :: xbase(n), xopt(n), xnew(n), xpt(npt,n), pq(npt-1) 134REAL (dp) :: pl(npt,npt-1), h(n,n), g(n), d(n), vlag(npt), w(npt) 135 136! The arguments N, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical to 137! the corresponding arguments in SUBROUTINE UOBYQA. 138 139! NPT is set by UOBYQA to (N*N+3*N+2)/2 for the above dimension statement. 140! XBASE will contain a shift of origin that reduces the contributions from 141! rounding errors to values of the model and Lagrange functions. 142! XOPT will be set to the displacement from XBASE of the vector of 143! variables that provides the least calculated F so far. 144! XNEW will be set to the displacement from XBASE of the vector of 145! variables for the current calculation of F. 146! XPT will contain the interpolation point coordinates relative to XBASE. 147! PQ will contain the parameters of the quadratic model. 148! PL will contain the parameters of the Lagrange functions. 149! H will provide the second derivatives that TRSTEP and LAGMAX require. 150! G will provide the first derivatives that TRSTEP and LAGMAX require. 151! D is reserved for trial steps from XOPT, except that it will contain 152! diagonal second derivatives during the initialization procedure. 153! VLAG will contain the values of the Lagrange functions at a new point X. 154! The array W will be used for working space. 155 156REAL (dp) :: half = 0.5_dp, one = 1.0_dp, tol = 0.01_dp, two = 2.0_dp 157REAL (dp) :: zero = 0.0_dp 158REAL (dp) :: ddknew, delta, detrat, diff, distest, dnorm, errtol, estim 159REAL (dp) :: evalue, f, fbase, fopt, fsave, ratio, rho, rhosq, sixthm 160REAL (dp) :: sum, sumg, sumh, temp, tempa, tworsq, vmax, vquad, wmult 161INTEGER :: i, ih, ip, iq, iw, j, jswitch, k, knew, kopt, ksave, ktemp 162INTEGER :: nf, nftest, nnp, nptm 163 164! Set some constants. 165 166nnp = n + n + 1 167nptm = npt - 1 168nftest = MAX(maxfun,1) 169 170! Initialization. NF is the number of function calculations so far. 171 172rho = rhobeg 173rhosq = rho * rho 174nf = 0 175DO i = 1, n 176 xbase(i) = x(i) 177 xpt(1:npt,i) = zero 178END DO 179pl(1:npt,1:nptm) = zero 180 181! The branch to label 120 obtains a new value of the objective function 182! and then there is a branch back to label 50, because the new function 183! value is needed to form the initial quadratic model. The least function 184! value so far and its index are noted below. 185 18650 x(1:n) = xbase(1:n) + xpt(nf+1,1:n) 187GO TO 150 188 18970 IF (nf == 1) THEN 190 fopt = f 191 kopt = nf 192 fbase = f 193 j = 0 194 jswitch = -1 195 ih = n 196ELSE 197 IF (f < fopt) THEN 198 fopt = f 199 kopt = nf 200 END IF 201END IF 202 203! Form the gradient and diagonal second derivatives of the initial 204! quadratic model and Lagrange functions. 205 206IF (nf <= nnp) THEN 207 jswitch = -jswitch 208 IF (jswitch > 0) THEN 209 IF (j >= 1) THEN 210 ih = ih + j 211 IF (w(j) < zero) THEN 212 d(j) = (fsave+f-two*fbase) / rhosq 213 pq(j) = (fsave-f) / (two*rho) 214 pl(1,ih) = -two / rhosq 215 pl(nf-1,j) = half / rho 216 pl(nf-1,ih) = one / rhosq 217 ELSE 218 pq(j) = (4.0D0*fsave-3.0D0*fbase-f) / (two*rho) 219 d(j) = (fbase+f-two*fsave) / rhosq 220 pl(1,j) = -1.5D0 / rho 221 pl(1,ih) = one / rhosq 222 pl(nf-1,j) = two / rho 223 pl(nf-1,ih) = -two / rhosq 224 END IF 225 pq(ih) = d(j) 226 pl(nf,j) = -half / rho 227 pl(nf,ih) = one / rhosq 228 END IF 229 230! Pick the shift from XBASE to the next initial interpolation point 231! that provides diagonal second derivatives. 232 233 IF (j < n) THEN 234 j = j + 1 235 xpt(nf+1,j) = rho 236 END IF 237 ELSE 238 fsave = f 239 IF (f < fbase) THEN 240 w(j) = rho 241 xpt(nf+1,j) = two * rho 242 ELSE 243 w(j) = -rho 244 xpt(nf+1,j) = -rho 245 END IF 246 END IF 247 IF (nf < nnp) GO TO 50 248 249! Form the off-diagonal second derivatives of the initial quadratic model. 250 251 ih = n 252 ip = 1 253 iq = 2 254END IF 255ih = ih + 1 256IF (nf > nnp) THEN 257 temp = one / (w(ip)*w(iq)) 258 tempa = f - fbase - w(ip) * pq(ip) - w(iq) * pq(iq) 259 pq(ih) = (tempa - half*rhosq*(d(ip)+d(iq))) * temp 260 pl(1,ih) = temp 261 iw = ip + ip 262 IF (w(ip) < zero) iw = iw + 1 263 pl(iw,ih) = -temp 264 iw = iq + iq 265 IF (w(iq) < zero) iw = iw + 1 266 pl(iw,ih) = -temp 267 pl(nf,ih) = temp 268 269! Pick the shift from XBASE to the next initial interpolation point 270! that provides off-diagonal second derivatives. 271 272 ip = ip + 1 273END IF 274IF (ip == iq) THEN 275 ih = ih + 1 276 ip = 1 277 iq = iq + 1 278END IF 279IF (nf < npt) THEN 280 xpt(nf+1,ip) = w(ip) 281 xpt(nf+1,iq) = w(iq) 282 GO TO 50 283END IF 284 285! Set parameters to begin the iterations for the current RHO. 286 287sixthm = zero 288delta = rho 28980 tworsq = (two*rho) ** 2 290rhosq = rho * rho 291 292! Form the gradient of the quadratic model at the trust region centre. 293 29490 knew = 0 295ih = n 296DO j = 1, n 297 xopt(j) = xpt(kopt,j) 298 g(j) = pq(j) 299 DO i = 1, j 300 ih = ih + 1 301 g(i) = g(i) + pq(ih) * xopt(j) 302 IF (i < j) g(j) = g(j) + pq(ih) * xopt(i) 303 h(i,j) = pq(ih) 304 END DO 305END DO 306 307! Generate the next trust region step and test its length. Set KNEW 308! to -1 if the purpose of the next F will be to improve conditioning, 309! and also calculate a lower bound on the Hessian term of the model Q. 310 311CALL trstep(n, g, h, delta, tol, d, evalue) 312temp = zero 313DO i = 1, n 314 temp = temp + d(i)**2 315END DO 316dnorm = MIN(delta,SQRT(temp)) 317errtol = -one 318IF (dnorm < half*rho) THEN 319 knew = -1 320 errtol = half * evalue * rho * rho 321 IF (nf <= npt+9) errtol = zero 322 GO TO 290 323END IF 324 325! Calculate the next value of the objective function. 326 327130 DO i = 1, n 328 xnew(i) = xopt(i) + d(i) 329 x(i) = xbase(i) + xnew(i) 330END DO 331150 IF (nf >= nftest) THEN 332 IF (iprint > 0) WRITE(*, 5000) 333 GO TO 420 334END IF 335nf = nf + 1 336CALL calfun(n, x, f) 337IF (iprint == 3) THEN 338 WRITE(*, 5100) nf, f, x(1:n) 339END IF 340IF (nf <= npt) GO TO 70 341IF (knew == -1) GO TO 420 342 343! Use the quadratic model to predict the change in F due to the step D, 344! and find the values of the Lagrange functions at the new point. 345 346vquad = zero 347ih = n 348DO j = 1, n 349 w(j) = d(j) 350 vquad = vquad + w(j) * pq(j) 351 DO i = 1, j 352 ih = ih + 1 353 w(ih) = d(i) * xnew(j) + d(j) * xopt(i) 354 IF (i == j) w(ih) = half * w(ih) 355 vquad = vquad + w(ih) * pq(ih) 356 END DO 357END DO 358DO k = 1, npt 359 temp = zero 360 DO j = 1, nptm 361 temp = temp + w(j) * pl(k,j) 362 END DO 363 vlag(k) = temp 364END DO 365vlag(kopt) = vlag(kopt) + one 366 367! Update SIXTHM, which is a lower bound on one sixth of the greatest 368! third derivative of F. 369 370diff = f - fopt - vquad 371sum = zero 372DO k = 1, npt 373 temp = zero 374 DO i = 1, n 375 temp = temp + (xpt(k,i)-xnew(i)) ** 2 376 END DO 377 temp = SQRT(temp) 378 sum = sum + ABS(temp*temp*temp*vlag(k)) 379END DO 380sixthm = MAX(sixthm, ABS(diff)/sum) 381 382! Update FOPT and XOPT if the new F is the least value of the objective 383! function so far. Then branch if D is not a trust region step. 384 385fsave = fopt 386IF (f < fopt) THEN 387 fopt = f 388 xopt(1:n) = xnew(1:n) 389END IF 390ksave = knew 391IF (knew <= 0) THEN 392 393! Pick the next value of DELTA after a trust region step. 394 395 IF (vquad >= zero) THEN 396 IF (iprint > 0) WRITE(*, 5200) 397 GO TO 420 398 END IF 399 ratio = (f-fsave) / vquad 400 IF (ratio <= 0.1D0) THEN 401 delta = half * dnorm 402 ELSE IF (ratio <= 0.7D0) THEN 403 delta = MAX(half*delta,dnorm) 404 ELSE 405 delta = MAX(delta, 1.25D0*dnorm, dnorm+rho) 406 END IF 407 IF (delta <= 1.5D0*rho) delta = rho 408 409! Set KNEW to the index of the next interpolation point to be deleted. 410 411 ktemp = 0 412 detrat = zero 413 IF (f >= fsave) THEN 414 ktemp = kopt 415 detrat = one 416 END IF 417 DO k = 1, npt 418 sum = zero 419 DO i = 1, n 420 sum = sum + (xpt(k,i)-xopt(i)) ** 2 421 END DO 422 temp = ABS(vlag(k)) 423 IF (sum > rhosq) temp = temp * (sum/rhosq) ** 1.5D0 424 IF (temp > detrat .AND. k /= ktemp) THEN 425 detrat = temp 426 ddknew = sum 427 knew = k 428 END IF 429 END DO 430 IF (knew == 0) GO TO 290 431END IF 432 433! Replace the interpolation point that has index KNEW by the point XNEW, 434! and also update the Lagrange functions and the quadratic model. 435 436DO i = 1, n 437 xpt(knew,i) = xnew(i) 438END DO 439temp = one / vlag(knew) 440DO j = 1, nptm 441 pl(knew,j) = temp * pl(knew,j) 442 pq(j) = pq(j) + diff * pl(knew,j) 443END DO 444DO k = 1, npt 445 IF (k /= knew) THEN 446 temp = vlag(k) 447 DO j = 1, nptm 448 pl(k,j) = pl(k,j) - temp * pl(knew,j) 449 END DO 450 END IF 451END DO 452 453! Update KOPT if F is the least calculated value of the objective function. 454! Then branch for another trust region calculation. The case KSAVE > 0 455! indicates that a model step has just been taken. 456 457IF (f < fsave) THEN 458 kopt = knew 459 GO TO 90 460END IF 461IF (ksave > 0) GO TO 90 462IF (dnorm > two*rho) GO TO 90 463IF (ddknew > tworsq) GO TO 90 464 465! Alternatively, find out if the interpolation points are close 466! enough to the best point so far. 467 468290 DO k = 1, npt 469 w(k) = zero 470 DO i = 1, n 471 w(k) = w(k) + (xpt(k,i)-xopt(i)) ** 2 472 END DO 473END DO 474320 knew = -1 475distest = tworsq 476DO k = 1, npt 477 IF (w(k) > distest) THEN 478 knew = k 479 distest = w(k) 480 END IF 481END DO 482 483! If a point is sufficiently far away, then set the gradient and Hessian 484! of its Lagrange function at the centre of the trust region, and find 485! half the sum of squares of components of the Hessian. 486 487IF (knew > 0) THEN 488 ih = n 489 sumh = zero 490 DO j = 1, n 491 g(j) = pl(knew,j) 492 DO i = 1, j 493 ih = ih + 1 494 temp = pl(knew,ih) 495 g(j) = g(j) + temp * xopt(i) 496 IF (i < j) THEN 497 g(i) = g(i) + temp * xopt(j) 498 sumh = sumh + temp * temp 499 END IF 500 h(i,j) = temp 501 END DO 502 sumh = sumh + half * temp * temp 503 END DO 504 505! If ERRTOL is positive, test whether to replace the interpolation point 506! with index KNEW, using a bound on the maximum modulus of its Lagrange 507! function in the trust region. 508 509 IF (errtol > zero) THEN 510 w(knew) = zero 511 sumg = zero 512 DO i = 1, n 513 sumg = sumg + g(i) ** 2 514 END DO 515 estim = rho * (SQRT(sumg)+rho*SQRT(half*sumh)) 516 wmult = sixthm * distest ** 1.5D0 517 IF (wmult*estim <= errtol) GO TO 320 518 END IF 519 520! If the KNEW-th point may be replaced, then pick a D that gives a large 521! value of the modulus of its Lagrange function within the trust region. 522! Here the vector XNEW is used as temporary working space. 523 524 CALL lagmax(n, g, h, rho, d, xnew, vmax) 525 IF (errtol > zero) THEN 526 IF (wmult*vmax <= errtol) GO TO 320 527 END IF 528 GO TO 130 529END IF 530IF (dnorm > rho) GO TO 90 531 532! Prepare to reduce RHO by shifting XBASE to the best point so far, 533! and make the corresponding changes to the gradients of the Lagrange 534! functions and the quadratic model. 535 536IF (rho > rhoend) THEN 537 ih = n 538 DO j = 1, n 539 xbase(j) = xbase(j) + xopt(j) 540 DO k = 1, npt 541 xpt(k,j) = xpt(k,j) - xopt(j) 542 END DO 543 DO i = 1, j 544 ih = ih + 1 545 pq(i) = pq(i) + pq(ih) * xopt(j) 546 IF (i < j) THEN 547 pq(j) = pq(j) + pq(ih) * xopt(i) 548 DO k = 1, npt 549 pl(k,j) = pl(k,j) + pl(k,ih) * xopt(i) 550 END DO 551 END IF 552 DO k = 1, npt 553 pl(k,i) = pl(k,i) + pl(k,ih) * xopt(j) 554 END DO 555 END DO 556 END DO 557 558! Pick the next values of RHO and DELTA. 559 560 delta = half * rho 561 ratio = rho / rhoend 562 IF (ratio <= 16.0D0) THEN 563 rho = rhoend 564 ELSE IF (ratio <= 250.0D0) THEN 565 rho = SQRT(ratio) * rhoend 566 ELSE 567 rho = 0.1D0 * rho 568 END IF 569 delta = MAX(delta,rho) 570 IF (iprint >= 2) THEN 571 IF (iprint >= 3) WRITE(*, 5300) 572 WRITE(*, 5400) rho, nf 573 WRITE(*, 5500) fopt, xbase(1:n) 574 END IF 575 GO TO 80 576END IF 577 578! Return from the calculation, after another Newton-Raphson step, if 579! it is too short to have been tried before. 580 581IF (errtol >= zero) GO TO 130 582420 IF (fopt <= f) THEN 583 DO i = 1, n 584 x(i) = xbase(i) + xopt(i) 585 END DO 586 f = fopt 587END IF 588IF (iprint >= 1) THEN 589 WRITE(*, 5600) nf 590 WRITE(*, 5500) f, x(1:n) 591END IF 592RETURN 593 5945000 FORMAT (/T5, 'Return from UOBYQA because CALFUN has been', & 595 ' called MAXFUN times') 5965100 FORMAT (/T5, 'Function number',i6,' F =', g18.10, & 597 ' The corresponding X is:'/ (t3, 5g15.6)) 5985200 FORMAT (/T5, 'Return from UOBYQA because a trust', & 599 ' region step has failed to reduce Q') 6005300 FORMAT (' ') 6015400 FORMAT (/T5, 'New RHO =', g11.4, ' Number of function values =',i6) 6025500 FORMAT (T5, 'Least value of F =', g23.15, & 603 ' The corresponding X is:'/ (t3, 5g15.6)) 6045600 FORMAT (/T5, 'At the return from UOBYQA', & 605 ' Number of function values =', i6) 606END SUBROUTINE uobyqb 607 608!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% trstep.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 609 610SUBROUTINE trstep(n, g, h, delta, tol, d, evalue) 611 612INTEGER, INTENT(IN) :: n 613REAL (dp), INTENT(IN) :: g(:) 614REAL (dp), INTENT(IN OUT) :: h(:,:) 615REAL (dp), INTENT(IN) :: delta 616REAL (dp), INTENT(IN) :: tol 617REAL (dp), INTENT(OUT) :: d(:) 618REAL (dp), INTENT(OUT) :: evalue 619 620! N is the number of variables of a quadratic objective function, Q say. 621! G is the gradient of Q at the origin. 622! H is the Hessian matrix of Q. Only the upper triangular and diagonal 623! parts need be set. The lower triangular part is used to store the 624! elements of a Householder similarity transformation. 625! DELTA is the trust region radius, and has to be positive. 626! TOL is the value of a tolerance from the open interval (0,1). 627! D will be set to the calculated vector of variables. 628 629! EVALUE will be set to the least eigenvalue of H if and only if D is a 630! Newton-Raphson step. Then EVALUE will be positive, but otherwise it 631! will be set to zero. 632 633! Let MAXRED be the maximum of Q(0)-Q(D) subject to ||D|| <= DELTA, 634! and let ACTRED be the value of Q(0)-Q(D) that is actually calculated. 635! We take the view that any D is acceptable if it has the properties 636 637! ||D|| <= DELTA and ACTRED <= (1-TOL)*MAXRED. 638 639! The calculation of D is done by the method of Section 2 of the paper 640! by MJDP in the 1997 Dundee Numerical Analysis Conference Proceedings, 641! after transforming H to tridiagonal form. 642 643! The arrays GG, TD, TN, W, PIV and Z will be used for working space. 644REAL (dp) :: gg(n), td(n), tn(n), w(n), piv(n), z(n) 645 646REAL (dp) :: delsq, dhd, dnorm, dsq, dtg, dtz, gam, gnorm, gsq, hnorm 647REAL (dp) :: par, parl, parlest, paru, paruest, phi, phil, phiu, pivksv 648REAL (dp) :: pivot, posdef, scale, shfmax, shfmin, shift, slope, sum 649REAL (dp) :: tdmin, temp, tempa, tempb, wsq, wwsq, wz, zsq 650INTEGER :: i, iterc, j, jp, k, kp, kpp, ksav, ksave, nm 651REAL (dp) :: one = 1.0_dp, two = 2.0_dp, zero = 0.0_dp 652 653! Initialization. 654 655delsq = delta * delta 656evalue = zero 657nm = n - 1 658DO i = 1, n 659 d(i) = zero 660 td(i) = h(i,i) 661 DO j = 1, i 662 h(i,j) = h(j,i) 663 END DO 664END DO 665 666! Apply Householder transformations to obtain a tridiagonal matrix that 667! is similar to H, and put the elements of the Householder vectors in 668! the lower triangular part of H. Further, TD and TN will contain the 669! diagonal and other nonzero elements of the tridiagonal matrix. 670 671DO k = 1, nm 672 kp = k + 1 673 sum = zero 674 IF (kp < n) THEN 675 kpp = kp + 1 676 DO i = kpp, n 677 sum = sum + h(i,k) ** 2 678 END DO 679 END IF 680 IF (sum == zero) THEN 681 tn(k) = h(kp,k) 682 h(kp,k) = zero 683 ELSE 684 temp = h(kp,k) 685 tn(k) = SIGN(SQRT(sum+temp*temp),temp) 686 h(kp,k) = -sum / (temp+tn(k)) 687 temp = SQRT(two/(sum+h(kp,k)**2)) 688 DO i = kp, n 689 w(i) = temp * h(i,k) 690 h(i,k) = w(i) 691 z(i) = td(i) * w(i) 692 END DO 693 wz = zero 694 DO j = kp, nm 695 jp = j + 1 696 DO i = jp, n 697 z(i) = z(i) + h(i,j) * w(j) 698 z(j) = z(j) + h(i,j) * w(i) 699 END DO 700 wz = wz + w(j) * z(j) 701 END DO 702 wz = wz + w(n) * z(n) 703 DO j = kp, n 704 td(j) = td(j) + w(j) * (wz*w(j)-two*z(j)) 705 IF (j < n) THEN 706 jp = j + 1 707 DO i = jp, n 708 h(i,j) = h(i,j) - w(i) * z(j) - w(j) * (z(i)-wz*w(i)) 709 END DO 710 END IF 711 END DO 712 END IF 713END DO 714 715! Form GG by applying the similarity transformation to G. 716 717gsq = zero 718DO i = 1, n 719 gg(i) = g(i) 720 gsq = gsq + g(i) ** 2 721END DO 722gnorm = SQRT(gsq) 723DO k = 1, nm 724 kp = k + 1 725 sum = zero 726 DO i = kp, n 727 sum = sum + gg(i) * h(i,k) 728 END DO 729 DO i = kp, n 730 gg(i) = gg(i) - sum * h(i,k) 731 END DO 732END DO 733 734! Begin the trust region calculation with a tridiagonal matrix by 735! calculating the norm of H. Then treat the case when H is zero. 736 737hnorm = ABS(td(1)) + ABS(tn(1)) 738tdmin = td(1) 739tn(n) = zero 740DO i = 2, n 741 temp = ABS(tn(i-1)) + ABS(td(i)) + ABS(tn(i)) 742 hnorm = MAX(hnorm,temp) 743 tdmin = MIN(tdmin,td(i)) 744END DO 745IF (hnorm == zero) THEN 746 IF (gnorm == zero) GO TO 420 747 scale = delta / gnorm 748 DO i = 1, n 749 d(i) = -scale * gg(i) 750 END DO 751 GO TO 380 752END IF 753 754! Set the initial values of PAR and its bounds. 755 756parl = MAX(zero, -tdmin, gnorm/delta-hnorm) 757parlest = parl 758par = parl 759paru = zero 760paruest = zero 761posdef = zero 762iterc = 0 763 764! Calculate the pivots of the Cholesky factorization of (H+PAR*I). 765 766160 iterc = iterc + 1 767ksav = 0 768piv(1) = td(1) + par 769k = 1 770170 IF (piv(k) > zero) THEN 771 piv(k+1) = td(k+1) + par - tn(k) ** 2 / piv(k) 772ELSE 773 IF (piv(k) < zero .OR. tn(k) /= zero) GO TO 180 774 ksav = k 775 piv(k+1) = td(k+1) + par 776END IF 777k = k + 1 778IF (k < n) GO TO 170 779IF (piv(k) >= zero) THEN 780 IF (piv(k) == zero) ksav = k 781 782! Branch if all the pivots are positive, allowing for the case when 783! G is zero. 784 785 IF (ksav == 0 .AND. gsq > zero) GO TO 250 786 IF (gsq == zero) THEN 787 IF (par == zero) GO TO 380 788 paru = par 789 paruest = par 790 IF (ksav == 0) GO TO 210 791 END IF 792 k = ksav 793END IF 794 795! Set D to a direction of nonpositive curvature of the given tridiagonal 796! matrix, and thus revise PARLEST. 797 798180 d(k) = one 799IF (ABS(tn(k)) <= ABS(piv(k))) THEN 800 dsq = one 801 dhd = piv(k) 802ELSE 803 temp = td(k+1) + par 804 IF (temp <= ABS(piv(k))) THEN 805 d(k+1) = SIGN(one,-tn(k)) 806 dhd = piv(k) + temp - two * ABS(tn(k)) 807 ELSE 808 d(k+1) = -tn(k) / temp 809 dhd = piv(k) + tn(k) * d(k+1) 810 END IF 811 dsq = one + d(k+1) ** 2 812END IF 813190 IF (k > 1) THEN 814 k = k - 1 815 IF (tn(k) /= zero) THEN 816 d(k) = -tn(k) * d(k+1) / piv(k) 817 dsq = dsq + d(k) ** 2 818 GO TO 190 819 END IF 820 d(1:k) = zero 821END IF 822parl = par 823parlest = par - dhd / dsq 824 825! Terminate with D set to a multiple of the current D if the following 826! test suggests that it suitable to do so. 827 828210 temp = paruest 829IF (gsq == zero) temp = temp * (one-tol) 830IF (paruest > zero .AND. parlest >= temp) THEN 831 dtg = DOT_PRODUCT( d(1:n), gg(1:n) ) 832 scale = -SIGN(delta/SQRT(dsq),dtg) 833 d(1:n) = scale * d(1:n) 834 GO TO 380 835END IF 836 837! Pick the value of PAR for the next iteration. 838 839240 IF (paru == zero) THEN 840 par = two * parlest + gnorm / delta 841ELSE 842 par = 0.5D0 * (parl+paru) 843 par = MAX(par,parlest) 844END IF 845IF (paruest > zero) par = MIN(par,paruest) 846GO TO 160 847 848! Calculate D for the current PAR in the positive definite case. 849 850250 w(1) = -gg(1) / piv(1) 851DO i = 2, n 852 w(i) = (-gg(i)-tn(i-1)*w(i-1)) / piv(i) 853END DO 854d(n) = w(n) 855DO i = nm, 1, -1 856 d(i) = w(i) - tn(i) * d(i+1) / piv(i) 857END DO 858 859! Branch if a Newton-Raphson step is acceptable. 860 861dsq = zero 862wsq = zero 863DO i = 1, n 864 dsq = dsq + d(i) ** 2 865 wsq = wsq + piv(i) * w(i) ** 2 866END DO 867IF (par /= zero .OR. dsq > delsq) THEN 868 869! Make the usual test for acceptability of a full trust region step. 870 871 dnorm = SQRT(dsq) 872 phi = one / dnorm - one / delta 873 temp = tol * (one+par*dsq/wsq) - dsq * phi * phi 874 IF (temp >= zero) THEN 875 scale = delta / dnorm 876 DO i = 1, n 877 d(i) = scale * d(i) 878 END DO 879 GO TO 380 880 END IF 881 IF (iterc >= 2 .AND. par <= parl) GO TO 380 882 IF (paru > zero .AND. par >= paru) GO TO 380 883 884! Complete the iteration when PHI is negative. 885 886 IF (phi < zero) THEN 887 parlest = par 888 IF (posdef == one) THEN 889 IF (phi <= phil) GO TO 380 890 slope = (phi-phil) / (par-parl) 891 parlest = par - phi / slope 892 END IF 893 slope = one / gnorm 894 IF (paru > zero) slope = (phiu-phi) / (paru-par) 895 temp = par - phi / slope 896 IF (paruest > zero) temp = MIN(temp,paruest) 897 paruest = temp 898 posdef = one 899 parl = par 900 phil = phi 901 GO TO 240 902 END IF 903 904! If required, calculate Z for the alternative test for convergence. 905 906 IF (posdef == zero) THEN 907 w(1) = one / piv(1) 908 DO i = 2, n 909 temp = -tn(i-1) * w(i-1) 910 w(i) = (SIGN(one,temp)+temp) / piv(i) 911 END DO 912 z(n) = w(n) 913 DO i = nm, 1, -1 914 z(i) = w(i) - tn(i) * z(i+1) / piv(i) 915 END DO 916 wwsq = zero 917 zsq = zero 918 dtz = zero 919 DO i = 1, n 920 wwsq = wwsq + piv(i) * w(i) ** 2 921 zsq = zsq + z(i) ** 2 922 dtz = dtz + d(i) * z(i) 923 END DO 924 925! Apply the alternative test for convergence. 926 927 tempa = ABS(delsq-dsq) 928 tempb = SQRT(dtz*dtz+tempa*zsq) 929 gam = tempa / (SIGN(tempb,dtz)+dtz) 930 temp = tol * (wsq+par*delsq) - gam * gam * wwsq 931 IF (temp >= zero) THEN 932 DO i = 1, n 933 d(i) = d(i) + gam * z(i) 934 END DO 935 GO TO 380 936 END IF 937 parlest = MAX(parlest,par-wwsq/zsq) 938 END IF 939 940! Complete the iteration when PHI is positive. 941 942 slope = one / gnorm 943 IF (paru > zero) THEN 944 IF (phi >= phiu) GO TO 380 945 slope = (phiu-phi) / (paru-par) 946 END IF 947 parlest = MAX(parlest,par-phi/slope) 948 paruest = par 949 IF (posdef == one) THEN 950 slope = (phi-phil) / (par-parl) 951 paruest = par - phi / slope 952 END IF 953 paru = par 954 phiu = phi 955 GO TO 240 956END IF 957 958! Set EVALUE to the least eigenvalue of the second derivative matrix if 959! D is a Newton-Raphson step. SHFMAX will be an upper bound on EVALUE. 960 961shfmin = zero 962pivot = td(1) 963shfmax = pivot 964DO k = 2, n 965 pivot = td(k) - tn(k-1) ** 2 / pivot 966 shfmax = MIN(shfmax,pivot) 967END DO 968 969! Find EVALUE by a bisection method, but occasionally SHFMAX may be 970! adjusted by the rule of false position. 971 972ksave = 0 973350 shift = 0.5D0 * (shfmin+shfmax) 974k = 1 975temp = td(1) - shift 976 977360 IF (temp > zero) THEN 978 piv(k) = temp 979 IF (k < n) THEN 980 temp = td(k+1) - shift - tn(k) ** 2 / temp 981 k = k + 1 982 GO TO 360 983 END IF 984 shfmin = shift 985ELSE 986 IF (k < ksave) GO TO 370 987 IF (k == ksave) THEN 988 IF (pivksv == zero) GO TO 370 989 IF (piv(k)-temp < temp-pivksv) THEN 990 pivksv = temp 991 shfmax = shift 992 ELSE 993 pivksv = zero 994 shfmax = (shift*piv(k) - shfmin*temp) / (piv(k)-temp) 995 END IF 996 ELSE 997 ksave = k 998 pivksv = temp 999 shfmax = shift 1000 END IF 1001END IF 1002IF (shfmin <= 0.99D0*shfmax) GO TO 350 1003370 evalue = shfmin 1004 1005! Apply the inverse Householder transformations to D. 1006 1007380 nm = n - 1 1008DO k = nm, 1, -1 1009 kp = k + 1 1010 sum = zero 1011 DO i = kp, n 1012 sum = sum + d(i) * h(i,k) 1013 END DO 1014 DO i = kp, n 1015 d(i) = d(i) - sum * h(i,k) 1016 END DO 1017END DO 1018 1019! Return from the subroutine. 1020 1021420 RETURN 1022END SUBROUTINE trstep 1023 1024!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lagmax.f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1025 1026SUBROUTINE lagmax(n, g, h, rho, d, v, vmax) 1027 1028INTEGER, INTENT(IN) :: n 1029REAL (dp), INTENT(IN) :: g(:) 1030REAL (dp), INTENT(OUT) :: h(:,:) 1031REAL (dp), INTENT(IN) :: rho 1032REAL (dp), INTENT(OUT) :: d(:) 1033REAL (dp), INTENT(OUT) :: v(:) 1034REAL (dp), INTENT(OUT) :: vmax 1035 1036! N is the number of variables of a quadratic objective function, Q say. 1037! G is the gradient of Q at the origin. 1038! H is the symmetric Hessian matrix of Q. Only the upper triangular and 1039! diagonal parts need be set. 1040! RHO is the trust region radius, and has to be positive. 1041! D will be set to the calculated vector of variables. 1042! The array V will be used for working space. 1043! VMAX will be set to |Q(0)-Q(D)|. 1044 1045! Calculating the D that maximizes |Q(0)-Q(D)| subject to ||D|| <= RHO 1046! requires of order N**3 operations, but sometimes it is adequate if 1047! |Q(0)-Q(D)| is within about 0.9 of its greatest possible value. This 1048! subroutine provides such a solution in only of order N**2 operations, 1049! where the claim of accuracy has been tested by numerical experiments. 1050 1051REAL (dp) :: half = 0.5_dp, one = 1.0_dp, zero = 0.0_dp 1052REAL (dp) :: dd, dhd, dlin, dsq, gd, gg, ghg, gnorm, halfrt, hmax, ratio 1053REAL (dp) :: scale, sum, sumv, temp, tempa, tempb, tempc, tempd, tempv 1054REAL (dp) :: vhg, vhv, vhw, vlin, vmu, vnorm, vsq, vv, wcos, whw, wsin, wsq 1055INTEGER :: i, j, k 1056 1057! Preliminary calculations. 1058 1059halfrt = SQRT(half) 1060 1061! Pick V such that ||HV|| / ||V|| is large. 1062 1063hmax = zero 1064DO i = 1, n 1065 sum = zero 1066 DO j = 1, n 1067 h(j,i) = h(i,j) 1068 sum = sum + h(i,j) ** 2 1069 END DO 1070 IF (sum > hmax) THEN 1071 hmax = sum 1072 k = i 1073 END IF 1074END DO 1075DO j = 1, n 1076 v(j) = h(k,j) 1077END DO 1078 1079! Set D to a vector in the subspace spanned by V and HV that maximizes 1080! |(D,HD)|/(D,D), except that we set D=HV if V and HV are nearly parallel. 1081! The vector that has the name D at label 60 used to be the vector W. 1082 1083vsq = zero 1084vhv = zero 1085dsq = zero 1086DO i = 1, n 1087 vsq = vsq + v(i) ** 2 1088 d(i) = DOT_PRODUCT( h(i,1:n), v(1:n) ) 1089 vhv = vhv + v(i) * d(i) 1090 dsq = dsq + d(i) ** 2 1091END DO 1092IF (vhv*vhv <= 0.9999D0*dsq*vsq) THEN 1093 temp = vhv / vsq 1094 wsq = zero 1095 DO i = 1, n 1096 d(i) = d(i) - temp * v(i) 1097 wsq = wsq + d(i) ** 2 1098 END DO 1099 whw = zero 1100 ratio = SQRT(wsq/vsq) 1101 DO i = 1, n 1102 temp = DOT_PRODUCT( h(i,1:n), d(1:n) ) 1103 whw = whw + temp * d(i) 1104 v(i) = ratio * v(i) 1105 END DO 1106 vhv = ratio * ratio * vhv 1107 vhw = ratio * wsq 1108 temp = half * (whw-vhv) 1109 temp = temp + SIGN(SQRT(temp**2+vhw**2),whw+vhv) 1110 DO i = 1, n 1111 d(i) = vhw * v(i) + temp * d(i) 1112 END DO 1113END IF 1114 1115! We now turn our attention to the subspace spanned by G and D. A multiple 1116! of the current D is returned if that choice seems to be adequate. 1117 1118gg = zero 1119gd = zero 1120dd = zero 1121dhd = zero 1122DO i = 1, n 1123 gg = gg + g(i) ** 2 1124 gd = gd + g(i) * d(i) 1125 dd = dd + d(i) ** 2 1126 sum = DOT_PRODUCT( h(i,1:n), d(1:n) ) 1127 dhd = dhd + sum * d(i) 1128END DO 1129temp = gd / gg 1130vv = zero 1131scale = SIGN(rho/SQRT(dd),gd*dhd) 1132DO i = 1, n 1133 v(i) = d(i) - temp * g(i) 1134 vv = vv + v(i) ** 2 1135 d(i) = scale * d(i) 1136END DO 1137gnorm = SQRT(gg) 1138IF (gnorm*dd <= 0.5D-2*rho*ABS(dhd) .OR. vv/dd <= 1.0D-4) THEN 1139 vmax = ABS(scale*(gd + half*scale*dhd)) 1140 GO TO 170 1141END IF 1142 1143! G and V are now orthogonal in the subspace spanned by G and D. Hence 1144! we generate an orthonormal basis of this subspace such that (D,HV) is 1145! negligible or zero, where D and V will be the basis vectors. 1146 1147ghg = zero 1148vhg = zero 1149vhv = zero 1150DO i = 1, n 1151 sum = DOT_PRODUCT( h(i,1:n), g(1:n) ) 1152 sumv = DOT_PRODUCT( h(i,1:n), v(1:n) ) 1153 ghg = ghg + sum * g(i) 1154 vhg = vhg + sumv * g(i) 1155 vhv = vhv + sumv * v(i) 1156END DO 1157vnorm = SQRT(vv) 1158ghg = ghg / gg 1159vhg = vhg / (vnorm*gnorm) 1160vhv = vhv / vv 1161IF (ABS(vhg) <= 0.01D0*MAX(ABS(ghg),ABS(vhv))) THEN 1162 vmu = ghg - vhv 1163 wcos = one 1164 wsin = zero 1165ELSE 1166 temp = half * (ghg-vhv) 1167 vmu = temp + SIGN(SQRT(temp**2+vhg**2),temp) 1168 temp = SQRT(vmu**2+vhg**2) 1169 wcos = vmu / temp 1170 wsin = vhg / temp 1171END IF 1172tempa = wcos / gnorm 1173tempb = wsin / vnorm 1174tempc = wcos / vnorm 1175tempd = wsin / gnorm 1176DO i = 1, n 1177 d(i) = tempa * g(i) + tempb * v(i) 1178 v(i) = tempc * v(i) - tempd * g(i) 1179END DO 1180 1181! The final D is a multiple of the current D, V, D+V or D-V. We make the 1182! choice from these possibilities that is optimal. 1183 1184dlin = wcos * gnorm / rho 1185vlin = -wsin * gnorm / rho 1186tempa = ABS(dlin) + half * ABS(vmu+vhv) 1187tempb = ABS(vlin) + half * ABS(ghg-vmu) 1188tempc = halfrt * (ABS(dlin)+ABS(vlin)) + 0.25D0 * ABS(ghg+vhv) 1189IF (tempa >= tempb .AND. tempa >= tempc) THEN 1190 tempd = SIGN(rho,dlin*(vmu+vhv)) 1191 tempv = zero 1192ELSE IF (tempb >= tempc) THEN 1193 tempd = zero 1194 tempv = SIGN(rho,vlin*(ghg-vmu)) 1195ELSE 1196 tempd = SIGN(halfrt*rho,dlin*(ghg+vhv)) 1197 tempv = SIGN(halfrt*rho,vlin*(ghg+vhv)) 1198END IF 1199DO i = 1, n 1200 d(i) = tempd * d(i) + tempv * v(i) 1201END DO 1202vmax = rho * rho * MAX(tempa,tempb,tempc) 1203170 RETURN 1204END SUBROUTINE lagmax 1205 1206END MODULE Powell_Optimize 1207 1208!------------------------------------------------------------------------------- 1209! 1210! Main program scriptmini 1211! 1212! reads input and starts optimisation 1213! 1214!------------------------------------------------------------------------------- 1215PROGRAM scriptmini 1216 1217USE Powell_Optimize 1218IMPLICIT NONE 1219INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) 1220 1221REAL (dp) :: rhobeg, rhoend 1222REAL(dp), DIMENSION(:), ALLOCATABLE :: x 1223INTEGER :: iprint, maxfun, n,istat 1224 1225OPEN(UNIT=17,FILE="scriptmini.in",FORM="FORMATTED",STATUS="OLD",IOSTAT=ISTAT) 1226IF (ISTAT.NE.0) THEN 1227 CALL write_documentation() 1228 STOP " Unable to open scriptmini.in " 1229END IF 1230READ(17,*,IOSTAT=ISTAT) N 1231IF (ISTAT.NE.0) THEN 1232 CALL write_documentation() 1233 STOP " Unable to read N in scriptmini.in " 1234END IF 1235ALLOCATE(x(N)) 1236READ(17,*,IOSTAT=ISTAT) rhobeg,rhoend 1237IF (ISTAT.NE.0) THEN 1238 CALL write_documentation() 1239 STOP " Unable to read rhobeg,rhoend in scriptmini.in " 1240END IF 1241READ(17,*,IOSTAT=ISTAT) maxfun 1242IF (ISTAT.NE.0) THEN 1243 CALL write_documentation() 1244 STOP " Unable to read maxfun in scriptmini.in " 1245END IF 1246READ(17,*,IOSTAT=ISTAT) iprint 1247IF (ISTAT.NE.0) THEN 1248 CALL write_documentation() 1249 STOP " Unable to read iprint in scriptmini.in " 1250END IF 1251READ(17,*,IOSTAT=ISTAT) x 1252IF (ISTAT.NE.0) THEN 1253 CALL write_documentation() 1254 STOP " Unable to read x in scriptmini.in " 1255END IF 1256 1257IF (.FALSE.) THEN 1258 CALL uobyqa (n, x, rhobeg, rhoend, iprint, maxfun) 1259ELSE 1260 CALL newuoa (n, x, rhobeg, rhoend, iprint, maxfun) 1261ENDIF 1262 1263DEALLOCATE(x) 1264 1265END PROGRAM scriptmini 1266 1267!------------------------------------------------------------------------------- 1268! 1269! calfun: the actual evaluation of the script 1270! 1271! 1272!------------------------------------------------------------------------------- 1273SUBROUTINE calfun(n, x, f) 1274 1275IMPLICIT NONE 1276INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) 1277 1278INTEGER, INTENT(IN) :: n 1279REAL (dp), INTENT(IN) :: x(:) 1280REAL (dp), INTENT(OUT) :: f 1281 1282character(LEN=40) :: format 1283 1284! write variables on a single line in the file 1285WRITE(format,'(A1,I4.4,A12)') '(',n,'(1X,F30.20))' 1286OPEN(UNIT=17,FILE="scriptmini_eval.in") 1287WRITE(UNIT=17,FMT=format) x(1:n) 1288CLOSE(UNIT=17) 1289 1290! execute scriptmini_eval 1291CALL system("./scriptmini_eval") 1292 1293! read value of the energy back 1294OPEN(UNIT=17,FILE="scriptmini_eval.out") 1295READ(UNIT=17,FMT=*) f 1296CLOSE(UNIT=17) 1297 1298END SUBROUTINE calfun 1299! 1300! 1301! 1302! 1303! 1304 SUBROUTINE NEWUOA (N,X,RHOBEG,RHOEND,IPRINT,MAXFUN) 1305 IMPLICIT REAL*8 (A-H,O-Z) 1306 DIMENSION X(*) 1307 REAL*8, DIMENSION(:), ALLOCATABLE :: W 1308 NPT=2*N+1 1309 ALLOCATE(W((NPT+13)*(NPT+N)+3*N*(N+3)/2)) 1310! 1311! This subroutine seeks the least value of a function of many variab 1312! by a trust region method that forms quadratic models by interpolat 1313! There can be some freedom in the interpolation conditions, which i 1314! taken up by minimizing the Frobenius norm of the change to the sec 1315! derivative of the quadratic model, beginning with a zero matrix. T 1316! arguments of the subroutine are as follows. 1317! 1318! N must be set to the number of variables and must be at least two. 1319! NPT is the number of interpolation conditions. Its value must be i 1320! interval [N+2,(N+1)(N+2)/2]. 1321! Initial values of the variables must be set in X(1),X(2),...,X(N). 1322! will be changed to the values that give the least calculated F. 1323! RHOBEG and RHOEND must be set to the initial and final values of a 1324! region radius, so both must be positive with RHOEND<=RHOBEG. Typ 1325! RHOBEG should be about one tenth of the greatest expected change 1326! variable, and RHOEND should indicate the accuracy that is requir 1327! the final values of the variables. 1328! The value of IPRINT should be set to 0, 1, 2 or 3, which controls 1329! amount of printing. Specifically, there is no output if IPRINT=0 1330! there is output only at the return if IPRINT=1. Otherwise, each 1331! value of RHO is printed, with the best vector of variables so fa 1332! the corresponding value of the objective function. Further, each 1333! value of F with its variables are output if IPRINT=3. 1334! MAXFUN must be set to an upper bound on the number of calls of CAL 1335! The array W will be used for working space. Its length must be at 1336! (NPT+13)*(NPT+N)+3*N*(N+3)/2. 1337! 1338! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must se 1339! the value of the objective function for the variables X(1),X(2),.. 1340! 1341! Partition the working space array, so that different parts of it c 1342! treated separately by the subroutine that performs the main calcul 1343! 1344 NP=N+1 1345 NPTM=NPT-NP 1346 IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN 1347 PRINT 10 1348 10 FORMAT (/4X,'Return from NEWUOA because NPT is not in', & 1349 & ' the required interval') 1350 GO TO 20 1351 END IF 1352 NDIM=NPT+N 1353 IXB=1 1354 IXO=IXB+N 1355 IXN=IXO+N 1356 IXP=IXN+N 1357 IFV=IXP+N*NPT 1358 IGQ=IFV+NPT 1359 IHQ=IGQ+N 1360 IPQ=IHQ+(N*NP)/2 1361 IBMAT=IPQ+NPT 1362 IZMAT=IBMAT+NDIM*N 1363 ID=IZMAT+NPT*NPTM 1364 IVL=ID+N 1365 IW=IVL+NDIM 1366! 1367! The above settings provide a partition of W for subroutine NEWUOB. 1368! The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements 1369! W plus the space that is needed by the last array of NEWUOB. 1370! 1371 CALL NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB), & 1372 & W(IXO),W(IXN),W(IXP),W(IFV),W(IGQ),W(IHQ),W(IPQ),W(IBMAT), & 1373 & W(IZMAT),NDIM,W(ID),W(IVL),W(IW)) 1374 20 RETURN 1375 END 1376 1377 SUBROUTINE NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,XBASE, & 1378 & XOPT,XNEW,XPT,FVAL,GQ,HQ,PQ,BMAT,ZMAT,NDIM,D,VLAG,W) 1379 IMPLICIT REAL*8 (A-H,O-Z) 1380 DIMENSION X(1:N),XBASE(*),XOPT(*),XNEW(*),XPT(NPT,*),FVAL(*), & 1381 & GQ(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),VLAG(*),W(*) 1382! 1383! The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are ide 1384! to the corresponding arguments in SUBROUTINE NEWUOA. 1385! XBASE will hold a shift of origin that should reduce the contribut 1386! from rounding errors to values of the model and Lagrange functio 1387! XOPT will be set to the displacement from XBASE of the vector of 1388! variables that provides the least calculated F so far. 1389! XNEW will be set to the displacement from XBASE of the vector of 1390! variables for the current calculation of F. 1391! XPT will contain the interpolation point coordinates relative to X 1392! FVAL will hold the values of F at the interpolation points. 1393! GQ will hold the gradient of the quadratic model at XBASE. 1394! HQ will hold the explicit second derivatives of the quadratic mode 1395! PQ will contain the parameters of the implicit second derivatives 1396! the quadratic model. 1397! BMAT will hold the last N columns of H. 1398! ZMAT will hold the factorization of the leading NPT by NPT submatr 1399! H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, wh 1400! the elements of DZ are plus or minus one, as specified by IDZ. 1401! NDIM is the first dimension of BMAT and has the value NPT+N. 1402! D is reserved for trial steps from XOPT. 1403! VLAG will contain the values of the Lagrange functions at a new po 1404! They are part of a product that requires VLAG to be of length ND 1405! The array W will be used for working space. Its length must be at 1406! 10*NDIM = 10*(NPT+N). 1407! 1408! Set some constants. 1409! 1410 INTERFACE 1411 SUBROUTINE calfun(n, x, f) 1412 IMPLICIT NONE 1413 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) 1414 INTEGER, INTENT(IN) :: n 1415 REAL (dp), INTENT(IN) :: x(:) 1416 REAL (dp), INTENT(OUT) :: f 1417 END SUBROUTINE calfun 1418 END INTERFACE 1419 1420 HALF=0.5D0 1421 ONE=1.0D0 1422 TENTH=0.1D0 1423 ZERO=0.0D0 1424 NP=N+1 1425 NH=(N*NP)/2 1426 NPTM=NPT-NP 1427 NFTEST=MAX0(MAXFUN,1) 1428! 1429! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. 1430! 1431 DO 20 J=1,N 1432 XBASE(J)=X(J) 1433 DO 10 K=1,NPT 1434 10 XPT(K,J)=ZERO 1435 DO 20 I=1,NDIM 1436 20 BMAT(I,J)=ZERO 1437 DO 30 IH=1,NH 1438 30 HQ(IH)=ZERO 1439 DO 40 K=1,NPT 1440 PQ(K)=ZERO 1441 DO 40 J=1,NPTM 1442 40 ZMAT(K,J)=ZERO 1443! 1444! Begin the initialization procedure. NF becomes one more than the n 1445! of function values so far. The coordinates of the displacement of 1446! next initial interpolation point from XBASE are set in XPT(NF,.). 1447! 1448 RHOSQ=RHOBEG*RHOBEG 1449 RECIP=ONE/RHOSQ 1450 RECIQ=DSQRT(HALF)/RHOSQ 1451 NF=0 1452 50 NFM=NF 1453 NFMM=NF-N 1454 NF=NF+1 1455 IF (NFM .LE. 2*N) THEN 1456 IF (NFM .GE. 1 .AND. NFM .LE. N) THEN 1457 XPT(NF,NFM)=RHOBEG 1458 ELSE IF (NFM .GT. N) THEN 1459 XPT(NF,NFMM)=-RHOBEG 1460 END IF 1461 ELSE 1462 ITEMP=(NFMM-1)/N 1463 JPT=NFM-ITEMP*N-N 1464 IPT=JPT+ITEMP 1465 IF (IPT .GT. N) THEN 1466 ITEMP=JPT 1467 JPT=IPT-N 1468 IPT=ITEMP 1469 END IF 1470 XIPT=RHOBEG 1471 IF (FVAL(IPT+NP) .LT. FVAL(IPT+1)) XIPT=-XIPT 1472 XJPT=RHOBEG 1473 IF (FVAL(JPT+NP) .LT. FVAL(JPT+1)) XJPT=-XJPT 1474 XPT(NF,IPT)=XIPT 1475 XPT(NF,JPT)=XJPT 1476 END IF 1477! 1478! Calculate the next value of F, label 70 being reached immediately 1479! after this calculation. The least function value so far and its in 1480! are required. 1481! 1482 DO 60 J=1,N 1483 60 X(J)=XPT(NF,J)+XBASE(J) 1484 GOTO 310 1485 70 FVAL(NF)=F 1486 IF (NF .EQ. 1) THEN 1487 FBEG=F 1488 FOPT=F 1489 KOPT=1 1490 ELSE IF (F .LT. FOPT) THEN 1491 FOPT=F 1492 KOPT=NF 1493 END IF 1494! 1495! Set the nonzero initial elements of BMAT and the quadratic model i 1496! the cases when NF is at most 2*N+1. 1497! 1498 IF (NFM .LE. 2*N) THEN 1499 IF (NFM .GE. 1 .AND. NFM .LE. N) THEN 1500 GQ(NFM)=(F-FBEG)/RHOBEG 1501 IF (NPT .LT. NF+N) THEN 1502 BMAT(1,NFM)=-ONE/RHOBEG 1503 BMAT(NF,NFM)=ONE/RHOBEG 1504 BMAT(NPT+NFM,NFM)=-HALF*RHOSQ 1505 END IF 1506 ELSE IF (NFM .GT. N) THEN 1507 BMAT(NF-N,NFMM)=HALF/RHOBEG 1508 BMAT(NF,NFMM)=-HALF/RHOBEG 1509 ZMAT(1,NFMM)=-RECIQ-RECIQ 1510 ZMAT(NF-N,NFMM)=RECIQ 1511 ZMAT(NF,NFMM)=RECIQ 1512 IH=(NFMM*(NFMM+1))/2 1513 TEMP=(FBEG-F)/RHOBEG 1514 HQ(IH)=(GQ(NFMM)-TEMP)/RHOBEG 1515 GQ(NFMM)=HALF*(GQ(NFMM)+TEMP) 1516 END IF 1517! 1518! Set the off-diagonal second derivatives of the Lagrange functions 1519! the initial quadratic model. 1520! 1521 ELSE 1522 IH=(IPT*(IPT-1))/2+JPT 1523 IF (XIPT .LT. ZERO) IPT=IPT+N 1524 IF (XJPT .LT. ZERO) JPT=JPT+N 1525 ZMAT(1,NFMM)=RECIP 1526 ZMAT(NF,NFMM)=RECIP 1527 ZMAT(IPT+1,NFMM)=-RECIP 1528 ZMAT(JPT+1,NFMM)=-RECIP 1529 HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/(XIPT*XJPT) 1530 END IF 1531 IF (NF .LT. NPT) GOTO 50 1532! 1533! Begin the iterative procedure, because the initial model is comple 1534! 1535 RHO=RHOBEG 1536 DELTA=RHO 1537 IDZ=1 1538 DIFFA=ZERO 1539 DIFFB=ZERO 1540 ITEST=0 1541 XOPTSQ=ZERO 1542 DO 80 I=1,N 1543 XOPT(I)=XPT(KOPT,I) 1544 80 XOPTSQ=XOPTSQ+XOPT(I)**2 1545 90 NFSAV=NF 1546! 1547! Generate the next trust region step and test its length. Set KNEW 1548! to -1 if the purpose of the next F will be to improve the model. 1549! 1550 100 KNEW=0 1551 CALL TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,D,W,W(NP), & 1552 & W(NP+N),W(NP+2*N),CRVMIN) 1553 DSQ=ZERO 1554 DO 110 I=1,N 1555 110 DSQ=DSQ+D(I)**2 1556 DNORM=DMIN1(DELTA,DSQRT(DSQ)) 1557 IF (DNORM .LT. HALF*RHO) THEN 1558 KNEW=-1 1559 DELTA=TENTH*DELTA 1560 RATIO=-1.0D0 1561 IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO 1562 IF (NF .LE. NFSAV+2) GOTO 460 1563 TEMP=0.125D0*CRVMIN*RHO*RHO 1564 IF (TEMP .LE. DMAX1(DIFFA,DIFFB,DIFFC)) GOTO 460 1565 GOTO 490 1566 END IF 1567! 1568! Shift XBASE if XOPT may be too far from XBASE. First make the chan 1569! to BMAT that do not depend on ZMAT. 1570! 1571 120 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN 1572 TEMPQ=0.25D0*XOPTSQ 1573 DO 140 K=1,NPT 1574 SUM=ZERO 1575 DO 130 I=1,N 1576 130 SUM=SUM+XPT(K,I)*XOPT(I) 1577 TEMP=PQ(K)*SUM 1578 SUM=SUM-HALF*XOPTSQ 1579 W(NPT+K)=SUM 1580 DO 140 I=1,N 1581 GQ(I)=GQ(I)+TEMP*XPT(K,I) 1582 XPT(K,I)=XPT(K,I)-HALF*XOPT(I) 1583 VLAG(I)=BMAT(K,I) 1584 W(I)=SUM*XPT(K,I)+TEMPQ*XOPT(I) 1585 IP=NPT+I 1586 DO 140 J=1,I 1587 140 BMAT(IP,J)=BMAT(IP,J)+VLAG(I)*W(J)+W(I)*VLAG(J) 1588! 1589! Then the revisions of BMAT that depend on ZMAT are calculated. 1590! 1591 DO 180 K=1,NPTM 1592 SUMZ=ZERO 1593 DO 150 I=1,NPT 1594 SUMZ=SUMZ+ZMAT(I,K) 1595 150 W(I)=W(NPT+I)*ZMAT(I,K) 1596 DO 170 J=1,N 1597 SUM=TEMPQ*SUMZ*XOPT(J) 1598 DO 160 I=1,NPT 1599 160 SUM=SUM+W(I)*XPT(I,J) 1600 VLAG(J)=SUM 1601 IF (K .LT. IDZ) SUM=-SUM 1602 DO 170 I=1,NPT 1603 170 BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K) 1604 DO 180 I=1,N 1605 IP=I+NPT 1606 TEMP=VLAG(I) 1607 IF (K .LT. IDZ) TEMP=-TEMP 1608 DO 180 J=1,I 1609 180 BMAT(IP,J)=BMAT(IP,J)+TEMP*VLAG(J) 1610! 1611! The following instructions complete the shift of XBASE, including 1612! the changes to the parameters of the quadratic model. 1613! 1614 IH=0 1615 DO 200 J=1,N 1616 W(J)=ZERO 1617 DO 190 K=1,NPT 1618 W(J)=W(J)+PQ(K)*XPT(K,J) 1619 190 XPT(K,J)=XPT(K,J)-HALF*XOPT(J) 1620 DO 200 I=1,J 1621 IH=IH+1 1622 IF (I .LT. J) GQ(J)=GQ(J)+HQ(IH)*XOPT(I) 1623 GQ(I)=GQ(I)+HQ(IH)*XOPT(J) 1624 HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J) 1625 200 BMAT(NPT+I,J)=BMAT(NPT+J,I) 1626 DO 210 J=1,N 1627 XBASE(J)=XBASE(J)+XOPT(J) 1628 210 XOPT(J)=ZERO 1629 XOPTSQ=ZERO 1630 END IF 1631! 1632! Pick the model step if KNEW is positive. A different choice of D 1633! may be made later, if the choice of D by BIGLAG causes substantial 1634! cancellation in DENOM. 1635! 1636 IF (KNEW .GT. 0) THEN 1637 CALL BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DSTEP, & 1638 & D,ALPHA,VLAG,VLAG(NPT+1),W,W(NP),W(NP+N)) 1639 END IF 1640! 1641! Calculate VLAG and BETA for the current choice of D. The first NPT 1642! components of W_check will be held in W. 1643! 1644 DO 230 K=1,NPT 1645 SUMA=ZERO 1646 SUMB=ZERO 1647 SUM=ZERO 1648 DO 220 J=1,N 1649 SUMA=SUMA+XPT(K,J)*D(J) 1650 SUMB=SUMB+XPT(K,J)*XOPT(J) 1651 220 SUM=SUM+BMAT(K,J)*D(J) 1652 W(K)=SUMA*(HALF*SUMA+SUMB) 1653 230 VLAG(K)=SUM 1654 BETA=ZERO 1655 DO 250 K=1,NPTM 1656 SUM=ZERO 1657 DO 240 I=1,NPT 1658 240 SUM=SUM+ZMAT(I,K)*W(I) 1659 IF (K .LT. IDZ) THEN 1660 BETA=BETA+SUM*SUM 1661 SUM=-SUM 1662 ELSE 1663 BETA=BETA-SUM*SUM 1664 END IF 1665 DO 250 I=1,NPT 1666 250 VLAG(I)=VLAG(I)+SUM*ZMAT(I,K) 1667 BSUM=ZERO 1668 DX=ZERO 1669 DO 280 J=1,N 1670 SUM=ZERO 1671 DO 260 I=1,NPT 1672 260 SUM=SUM+W(I)*BMAT(I,J) 1673 BSUM=BSUM+SUM*D(J) 1674 JP=NPT+J 1675 DO 270 K=1,N 1676 270 SUM=SUM+BMAT(JP,K)*D(K) 1677 VLAG(JP)=SUM 1678 BSUM=BSUM+SUM*D(J) 1679 280 DX=DX+D(J)*XOPT(J) 1680 BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM 1681 VLAG(KOPT)=VLAG(KOPT)+ONE 1682! 1683! If KNEW is positive and if the cancellation in DENOM is unacceptab 1684! then BIGDEN calculates an alternative model step, XNEW being used 1685! working space. 1686! 1687 IF (KNEW .GT. 0) THEN 1688 TEMP=ONE+ALPHA*BETA/VLAG(KNEW)**2 1689 IF (DABS(TEMP) .LE. 0.8D0) THEN 1690 CALL BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT, & 1691 & KNEW,D,W,VLAG,BETA,XNEW,W(NDIM+1),W(6*NDIM+1)) 1692 END IF 1693 END IF 1694! 1695! Calculate the next value of the objective function. 1696! 1697 290 DO 300 I=1,N 1698 XNEW(I)=XOPT(I)+D(I) 1699 300 X(I)=XBASE(I)+XNEW(I) 1700 NF=NF+1 1701 310 IF (NF .GT. NFTEST) THEN 1702 NF=NF-1 1703 IF (IPRINT .GT. 0) PRINT 320 1704 320 FORMAT (/4X,'Return from NEWUOA because CALFUN has been', & 1705 & ' called MAXFUN times.') 1706 GOTO 530 1707 END IF 1708 CALL CALFUN (N,X,F) 1709 IF (IPRINT .EQ. 3) THEN 1710 PRINT 330, NF,F,(X(I),I=1,N) 1711 330 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, & 1712 & ' The corresponding X is:'/(2X,5D15.6)) 1713 END IF 1714 IF (NF .LE. NPT) GOTO 70 1715 IF (KNEW .EQ. -1) GOTO 530 1716! 1717! Use the quadratic model to predict the change in F due to the step 1718! and set DIFF to the error of this prediction. 1719! 1720 VQUAD=ZERO 1721 IH=0 1722 DO 340 J=1,N 1723 VQUAD=VQUAD+D(J)*GQ(J) 1724 DO 340 I=1,J 1725 IH=IH+1 1726 TEMP=D(I)*XNEW(J)+D(J)*XOPT(I) 1727 IF (I .EQ. J) TEMP=HALF*TEMP 1728 340 VQUAD=VQUAD+TEMP*HQ(IH) 1729 DO 350 K=1,NPT 1730 350 VQUAD=VQUAD+PQ(K)*W(K) 1731 DIFF=F-FOPT-VQUAD 1732 DIFFC=DIFFB 1733 DIFFB=DIFFA 1734 DIFFA=DABS(DIFF) 1735 IF (DNORM .GT. RHO) NFSAV=NF 1736! 1737! Update FOPT and XOPT if the new F is the least value of the object 1738! function so far. The branch when KNEW is positive occurs if D is n 1739! a trust region step. 1740! 1741 FSAVE=FOPT 1742 IF (F .LT. FOPT) THEN 1743 FOPT=F 1744 XOPTSQ=ZERO 1745 DO 360 I=1,N 1746 XOPT(I)=XNEW(I) 1747 360 XOPTSQ=XOPTSQ+XOPT(I)**2 1748 END IF 1749 KSAVE=KNEW 1750 IF (KNEW .GT. 0) GOTO 410 1751! 1752! Pick the next value of DELTA after a trust region step. 1753! 1754 IF (VQUAD .GE. ZERO) THEN 1755 IF (IPRINT .GT. 0) PRINT 370 1756 370 FORMAT (/4X,'Return from NEWUOA because a trust', & 1757 & ' region step has failed to reduce Q.') 1758 GOTO 530 1759 END IF 1760 RATIO=(F-FSAVE)/VQUAD 1761 IF (RATIO .LE. TENTH) THEN 1762 DELTA=HALF*DNORM 1763 ELSE IF (RATIO .LE. 0.7D0) THEN 1764 DELTA=DMAX1(HALF*DELTA,DNORM) 1765 ELSE 1766 DELTA=DMAX1(HALF*DELTA,DNORM+DNORM) 1767 END IF 1768 IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO 1769! 1770! Set KNEW to the index of the next interpolation point to be delete 1771! 1772 RHOSQ=DMAX1(TENTH*DELTA,RHO)**2 1773 KTEMP=0 1774 DETRAT=ZERO 1775 IF (F .GE. FSAVE) THEN 1776 KTEMP=KOPT 1777 DETRAT=ONE 1778 END IF 1779 DO 400 K=1,NPT 1780 HDIAG=ZERO 1781 DO 380 J=1,NPTM 1782 TEMP=ONE 1783 IF (J .LT. IDZ) TEMP=-ONE 1784 380 HDIAG=HDIAG+TEMP*ZMAT(K,J)**2 1785 TEMP=DABS(BETA*HDIAG+VLAG(K)**2) 1786 DISTSQ=ZERO 1787 DO 390 J=1,N 1788 390 DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2 1789 IF (DISTSQ .GT. RHOSQ) TEMP=TEMP*(DISTSQ/RHOSQ)**3 1790 IF (TEMP .GT. DETRAT .AND. K .NE. KTEMP) THEN 1791 DETRAT=TEMP 1792 KNEW=K 1793 END IF 1794 400 END DO 1795 IF (KNEW .EQ. 0) GOTO 460 1796! 1797! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point 1798! can be moved. Begin the updating of the quadratic model, starting 1799! with the explicit second derivative term. 1800! 1801 410 CALL UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W) 1802 FVAL(KNEW)=F 1803 IH=0 1804 DO 420 I=1,N 1805 TEMP=PQ(KNEW)*XPT(KNEW,I) 1806 DO 420 J=1,I 1807 IH=IH+1 1808 420 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J) 1809 PQ(KNEW)=ZERO 1810! 1811! Update the other second derivative parameters, and then the gradie 1812! vector of the model. Also include the new interpolation point. 1813! 1814 DO 440 J=1,NPTM 1815 TEMP=DIFF*ZMAT(KNEW,J) 1816 IF (J .LT. IDZ) TEMP=-TEMP 1817 DO 440 K=1,NPT 1818 440 PQ(K)=PQ(K)+TEMP*ZMAT(K,J) 1819 GQSQ=ZERO 1820 DO 450 I=1,N 1821 GQ(I)=GQ(I)+DIFF*BMAT(KNEW,I) 1822 GQSQ=GQSQ+GQ(I)**2 1823 450 XPT(KNEW,I)=XNEW(I) 1824! 1825! If a trust region step makes a small change to the objective funct 1826! then calculate the gradient of the least Frobenius norm interpolan 1827! XBASE, and store it in W, using VLAG for a vector of right hand si 1828! 1829 IF (KSAVE .EQ. 0 .AND. DELTA .EQ. RHO) THEN 1830 IF (DABS(RATIO) .GT. 1.0D-2) THEN 1831 ITEST=0 1832 ELSE 1833 DO 700 K=1,NPT 1834 700 VLAG(K)=FVAL(K)-FVAL(KOPT) 1835 GISQ=ZERO 1836 DO 720 I=1,N 1837 SUM=ZERO 1838 DO 710 K=1,NPT 1839 710 SUM=SUM+BMAT(K,I)*VLAG(K) 1840 GISQ=GISQ+SUM*SUM 1841 720 W(I)=SUM 1842! 1843! Test whether to replace the new quadratic model by the least Frobe 1844! norm interpolant, making the replacement if the test is satisfied. 1845! 1846 ITEST=ITEST+1 1847 IF (GQSQ .LT. 1.0D2*GISQ) ITEST=0 1848 IF (ITEST .GE. 3) THEN 1849 DO 730 I=1,N 1850 730 GQ(I)=W(I) 1851 DO 740 IH=1,NH 1852 740 HQ(IH)=ZERO 1853 DO 760 J=1,NPTM 1854 W(J)=ZERO 1855 DO 750 K=1,NPT 1856 750 W(J)=W(J)+VLAG(K)*ZMAT(K,J) 1857 760 IF (J .LT. IDZ) W(J)=-W(J) 1858 DO 770 K=1,NPT 1859 PQ(K)=ZERO 1860 DO 770 J=1,NPTM 1861 770 PQ(K)=PQ(K)+ZMAT(K,J)*W(J) 1862 ITEST=0 1863 END IF 1864 END IF 1865 END IF 1866 IF (F .LT. FSAVE) KOPT=KNEW 1867! 1868! If a trust region step has provided a sufficient decrease in F, th 1869! branch for another trust region calculation. The case KSAVE>0 occu 1870! when the new function value was calculated by a model step. 1871! 1872 IF (F .LE. FSAVE+TENTH*VQUAD) GOTO 100 1873 IF (KSAVE .GT. 0) GOTO 100 1874! 1875! Alternatively, find out if the interpolation points are close enou 1876! to the best point so far. 1877! 1878 KNEW=0 1879 460 DISTSQ=4.0D0*DELTA*DELTA 1880 DO 480 K=1,NPT 1881 SUM=ZERO 1882 DO 470 J=1,N 1883 470 SUM=SUM+(XPT(K,J)-XOPT(J))**2 1884 IF (SUM .GT. DISTSQ) THEN 1885 KNEW=K 1886 DISTSQ=SUM 1887 END IF 1888 480 END DO 1889! 1890! If KNEW is positive, then set DSTEP, and branch back for the next 1891! iteration, which will generate a "model step". 1892! 1893 IF (KNEW .GT. 0) THEN 1894 DSTEP=DMAX1(DMIN1(TENTH*DSQRT(DISTSQ),HALF*DELTA),RHO) 1895 DSQ=DSTEP*DSTEP 1896 GOTO 120 1897 END IF 1898 IF (RATIO .GT. ZERO) GOTO 100 1899 IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 100 1900! 1901! The calculations with the current value of RHO are complete. Pick 1902! next values of RHO and DELTA. 1903! 1904 490 IF (RHO .GT. RHOEND) THEN 1905 DELTA=HALF*RHO 1906 RATIO=RHO/RHOEND 1907 IF (RATIO .LE. 16.0D0) THEN 1908 RHO=RHOEND 1909 ELSE IF (RATIO .LE. 250.0D0) THEN 1910 RHO=DSQRT(RATIO)*RHOEND 1911 ELSE 1912 RHO=TENTH*RHO 1913 END IF 1914 DELTA=DMAX1(DELTA,RHO) 1915 IF (IPRINT .GE. 2) THEN 1916 IF (IPRINT .GE. 3) PRINT 500 1917 500 FORMAT (5X) 1918 PRINT 510, RHO,NF 1919 510 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of', & 1920 & ' function values =',I6) 1921 PRINT 520, FOPT,(XBASE(I)+XOPT(I),I=1,N) 1922 520 FORMAT (4X,'Least value of F =',1PD23.15,9X, & 1923 & 'The corresponding X is:'/(2X,5D15.6)) 1924 END IF 1925 GOTO 90 1926 END IF 1927! 1928! Return from the calculation, after another Newton-Raphson step, if 1929! it is too short to have been tried before. 1930! 1931 IF (KNEW .EQ. -1) GOTO 290 1932 530 IF (FOPT .LE. F) THEN 1933 DO 540 I=1,N 1934 540 X(I)=XBASE(I)+XOPT(I) 1935 F=FOPT 1936 END IF 1937 IF (IPRINT .GE. 1) THEN 1938 PRINT 550, NF 1939 550 FORMAT (/4X,'At the return from NEWUOA',5X, & 1940 & 'Number of function values =',I6) 1941 PRINT 520, F,(X(I),I=1,N) 1942 END IF 1943 RETURN 1944 END 1945 1946 SUBROUTINE BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT, & 1947 & KNEW,D,W,VLAG,BETA,S,WVEC,PROD) 1948 IMPLICIT REAL*8 (A-H,O-Z) 1949 DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*), & 1950 & W(*),VLAG(*),S(*),WVEC(NDIM,*),PROD(NDIM,*) 1951 DIMENSION DEN(9),DENEX(9),PAR(9) 1952! 1953! N is the number of variables. 1954! NPT is the number of interpolation equations. 1955! XOPT is the best interpolation point so far. 1956! XPT contains the coordinates of the current interpolation points. 1957! BMAT provides the last N columns of H. 1958! ZMAT and IDZ give a factorization of the first NPT by NPT submatri 1959! NDIM is the first dimension of BMAT and has the value NPT+N. 1960! KOPT is the index of the optimal interpolation point. 1961! KNEW is the index of the interpolation point that is going to be m 1962! D will be set to the step from XOPT to the new point, and on entry 1963! should be the D that was calculated by the last call of BIGLAG. 1964! length of the initial D provides a trust region bound on the fin 1965! W will be set to Wcheck for the final choice of D. 1966! VLAG will be set to Theta*Wcheck+e_b for the final choice of D. 1967! BETA will be set to the value that will occur in the updating form 1968! when the KNEW-th interpolation point is moved to its new positio 1969! S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be us 1970! for working space. 1971! 1972! D is calculated in a way that should provide a denominator with a 1973! modulus in the updating formula when the KNEW-th interpolation poi 1974! shifted to the new position XOPT+D. 1975! 1976! Set some constants. 1977! 1978 HALF=0.5D0 1979 ONE=1.0D0 1980 QUART=0.25D0 1981 TWO=2.0D0 1982 ZERO=0.0D0 1983 TWOPI=8.0D0*DATAN(ONE) 1984 NPTM=NPT-N-1 1985! 1986! Store the first NPT elements of the KNEW-th column of H in W(N+1) 1987! to W(N+NPT). 1988! 1989 DO 10 K=1,NPT 1990 10 W(N+K)=ZERO 1991 DO 20 J=1,NPTM 1992 TEMP=ZMAT(KNEW,J) 1993 IF (J .LT. IDZ) TEMP=-TEMP 1994 DO 20 K=1,NPT 1995 20 W(N+K)=W(N+K)+TEMP*ZMAT(K,J) 1996 ALPHA=W(N+KNEW) 1997! 1998! The initial search direction D is taken from the last call of BIGL 1999! and the initial S is set below, usually to the direction from X_OP 2000! to X_KNEW, but a different direction to an interpolation point may 2001! be chosen, in order to prevent S from being nearly parallel to D. 2002! 2003 DD=ZERO 2004 DS=ZERO 2005 SS=ZERO 2006 XOPTSQ=ZERO 2007 DO 30 I=1,N 2008 DD=DD+D(I)**2 2009 S(I)=XPT(KNEW,I)-XOPT(I) 2010 DS=DS+D(I)*S(I) 2011 SS=SS+S(I)**2 2012 30 XOPTSQ=XOPTSQ+XOPT(I)**2 2013 IF (DS*DS .GT. 0.99D0*DD*SS) THEN 2014 KSAV=KNEW 2015 DTEST=DS*DS/SS 2016 DO 50 K=1,NPT 2017 IF (K .NE. KOPT) THEN 2018 DSTEMP=ZERO 2019 SSTEMP=ZERO 2020 DO 40 I=1,N 2021 DIFF=XPT(K,I)-XOPT(I) 2022 DSTEMP=DSTEMP+D(I)*DIFF 2023 40 SSTEMP=SSTEMP+DIFF*DIFF 2024 IF (DSTEMP*DSTEMP/SSTEMP .LT. DTEST) THEN 2025 KSAV=K 2026 DTEST=DSTEMP*DSTEMP/SSTEMP 2027 DS=DSTEMP 2028 SS=SSTEMP 2029 END IF 2030 END IF 2031 50 CONTINUE 2032 DO 60 I=1,N 2033 60 S(I)=XPT(KSAV,I)-XOPT(I) 2034 END IF 2035 SSDEN=DD*SS-DS*DS 2036 ITERC=0 2037 DENSAV=ZERO 2038! 2039! Begin the iteration by overwriting S with a vector that has the 2040! required length and direction. 2041! 2042 70 ITERC=ITERC+1 2043 TEMP=ONE/DSQRT(SSDEN) 2044 XOPTD=ZERO 2045 XOPTS=ZERO 2046 DO 80 I=1,N 2047 S(I)=TEMP*(DD*S(I)-DS*D(I)) 2048 XOPTD=XOPTD+XOPT(I)*D(I) 2049 80 XOPTS=XOPTS+XOPT(I)*S(I) 2050! 2051! Set the coefficients of the first two terms of BETA. 2052! 2053 TEMPA=HALF*XOPTD*XOPTD 2054 TEMPB=HALF*XOPTS*XOPTS 2055 DEN(1)=DD*(XOPTSQ+HALF*DD)+TEMPA+TEMPB 2056 DEN(2)=TWO*XOPTD*DD 2057 DEN(3)=TWO*XOPTS*DD 2058 DEN(4)=TEMPA-TEMPB 2059 DEN(5)=XOPTD*XOPTS 2060 DO 90 I=6,9 2061 90 DEN(I)=ZERO 2062! 2063! Put the coefficients of Wcheck in WVEC. 2064! 2065 DO 110 K=1,NPT 2066 TEMPA=ZERO 2067 TEMPB=ZERO 2068 TEMPC=ZERO 2069 DO 100 I=1,N 2070 TEMPA=TEMPA+XPT(K,I)*D(I) 2071 TEMPB=TEMPB+XPT(K,I)*S(I) 2072 100 TEMPC=TEMPC+XPT(K,I)*XOPT(I) 2073 WVEC(K,1)=QUART*(TEMPA*TEMPA+TEMPB*TEMPB) 2074 WVEC(K,2)=TEMPA*TEMPC 2075 WVEC(K,3)=TEMPB*TEMPC 2076 WVEC(K,4)=QUART*(TEMPA*TEMPA-TEMPB*TEMPB) 2077 110 WVEC(K,5)=HALF*TEMPA*TEMPB 2078 DO 120 I=1,N 2079 IP=I+NPT 2080 WVEC(IP,1)=ZERO 2081 WVEC(IP,2)=D(I) 2082 WVEC(IP,3)=S(I) 2083 WVEC(IP,4)=ZERO 2084 120 WVEC(IP,5)=ZERO 2085! 2086! Put the coefficients of THETA*Wcheck in PROD. 2087! 2088 DO 190 JC=1,5 2089 NW=NPT 2090 IF (JC .EQ. 2 .OR. JC .EQ. 3) NW=NDIM 2091 DO 130 K=1,NPT 2092 130 PROD(K,JC)=ZERO 2093 DO 150 J=1,NPTM 2094 SUM=ZERO 2095 DO 140 K=1,NPT 2096 140 SUM=SUM+ZMAT(K,J)*WVEC(K,JC) 2097 IF (J .LT. IDZ) SUM=-SUM 2098 DO 150 K=1,NPT 2099 150 PROD(K,JC)=PROD(K,JC)+SUM*ZMAT(K,J) 2100 IF (NW .EQ. NDIM) THEN 2101 DO 170 K=1,NPT 2102 SUM=ZERO 2103 DO 160 J=1,N 2104 160 SUM=SUM+BMAT(K,J)*WVEC(NPT+J,JC) 2105 170 PROD(K,JC)=PROD(K,JC)+SUM 2106 END IF 2107 DO 190 J=1,N 2108 SUM=ZERO 2109 DO 180 I=1,NW 2110 180 SUM=SUM+BMAT(I,J)*WVEC(I,JC) 2111 190 PROD(NPT+J,JC)=SUM 2112! 2113! Include in DEN the part of BETA that depends on THETA. 2114! 2115 DO 210 K=1,NDIM 2116 SUM=ZERO 2117 DO 200 I=1,5 2118 PAR(I)=HALF*PROD(K,I)*WVEC(K,I) 2119 200 SUM=SUM+PAR(I) 2120 DEN(1)=DEN(1)-PAR(1)-SUM 2121 TEMPA=PROD(K,1)*WVEC(K,2)+PROD(K,2)*WVEC(K,1) 2122 TEMPB=PROD(K,2)*WVEC(K,4)+PROD(K,4)*WVEC(K,2) 2123 TEMPC=PROD(K,3)*WVEC(K,5)+PROD(K,5)*WVEC(K,3) 2124 DEN(2)=DEN(2)-TEMPA-HALF*(TEMPB+TEMPC) 2125 DEN(6)=DEN(6)-HALF*(TEMPB-TEMPC) 2126 TEMPA=PROD(K,1)*WVEC(K,3)+PROD(K,3)*WVEC(K,1) 2127 TEMPB=PROD(K,2)*WVEC(K,5)+PROD(K,5)*WVEC(K,2) 2128 TEMPC=PROD(K,3)*WVEC(K,4)+PROD(K,4)*WVEC(K,3) 2129 DEN(3)=DEN(3)-TEMPA-HALF*(TEMPB-TEMPC) 2130 DEN(7)=DEN(7)-HALF*(TEMPB+TEMPC) 2131 TEMPA=PROD(K,1)*WVEC(K,4)+PROD(K,4)*WVEC(K,1) 2132 DEN(4)=DEN(4)-TEMPA-PAR(2)+PAR(3) 2133 TEMPA=PROD(K,1)*WVEC(K,5)+PROD(K,5)*WVEC(K,1) 2134 TEMPB=PROD(K,2)*WVEC(K,3)+PROD(K,3)*WVEC(K,2) 2135 DEN(5)=DEN(5)-TEMPA-HALF*TEMPB 2136 DEN(8)=DEN(8)-PAR(4)+PAR(5) 2137 TEMPA=PROD(K,4)*WVEC(K,5)+PROD(K,5)*WVEC(K,4) 2138 210 DEN(9)=DEN(9)-HALF*TEMPA 2139! 2140! Extend DEN so that it holds all the coefficients of DENOM. 2141! 2142 SUM=ZERO 2143 DO 220 I=1,5 2144 PAR(I)=HALF*PROD(KNEW,I)**2 2145 220 SUM=SUM+PAR(I) 2146 DENEX(1)=ALPHA*DEN(1)+PAR(1)+SUM 2147 TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,2) 2148 TEMPB=PROD(KNEW,2)*PROD(KNEW,4) 2149 TEMPC=PROD(KNEW,3)*PROD(KNEW,5) 2150 DENEX(2)=ALPHA*DEN(2)+TEMPA+TEMPB+TEMPC 2151 DENEX(6)=ALPHA*DEN(6)+TEMPB-TEMPC 2152 TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,3) 2153 TEMPB=PROD(KNEW,2)*PROD(KNEW,5) 2154 TEMPC=PROD(KNEW,3)*PROD(KNEW,4) 2155 DENEX(3)=ALPHA*DEN(3)+TEMPA+TEMPB-TEMPC 2156 DENEX(7)=ALPHA*DEN(7)+TEMPB+TEMPC 2157 TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,4) 2158 DENEX(4)=ALPHA*DEN(4)+TEMPA+PAR(2)-PAR(3) 2159 TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,5) 2160 DENEX(5)=ALPHA*DEN(5)+TEMPA+PROD(KNEW,2)*PROD(KNEW,3) 2161 DENEX(8)=ALPHA*DEN(8)+PAR(4)-PAR(5) 2162 DENEX(9)=ALPHA*DEN(9)+PROD(KNEW,4)*PROD(KNEW,5) 2163! 2164! Seek the value of the angle that maximizes the modulus of DENOM. 2165! 2166 SUM=DENEX(1)+DENEX(2)+DENEX(4)+DENEX(6)+DENEX(8) 2167 DENOLD=SUM 2168 DENMAX=SUM 2169 ISAVE=0 2170 IU=49 2171 TEMP=TWOPI/DBLE(IU+1) 2172 PAR(1)=ONE 2173 DO 250 I=1,IU 2174 ANGLE=DBLE(I)*TEMP 2175 PAR(2)=DCOS(ANGLE) 2176 PAR(3)=DSIN(ANGLE) 2177 DO 230 J=4,8,2 2178 PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1) 2179 230 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2) 2180 SUMOLD=SUM 2181 SUM=ZERO 2182 DO 240 J=1,9 2183 240 SUM=SUM+DENEX(J)*PAR(J) 2184 IF (DABS(SUM) .GT. DABS(DENMAX)) THEN 2185 DENMAX=SUM 2186 ISAVE=I 2187 TEMPA=SUMOLD 2188 ELSE IF (I .EQ. ISAVE+1) THEN 2189 TEMPB=SUM 2190 END IF 2191 250 END DO 2192 IF (ISAVE .EQ. 0) TEMPA=SUM 2193 IF (ISAVE .EQ. IU) TEMPB=DENOLD 2194 STEP=ZERO 2195 IF (TEMPA .NE. TEMPB) THEN 2196 TEMPA=TEMPA-DENMAX 2197 TEMPB=TEMPB-DENMAX 2198 STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB) 2199 END IF 2200 ANGLE=TEMP*(DBLE(ISAVE)+STEP) 2201! 2202! Calculate the new parameters of the denominator, the new VLAG vect 2203! and the new D. Then test for convergence. 2204! 2205 PAR(2)=DCOS(ANGLE) 2206 PAR(3)=DSIN(ANGLE) 2207 DO 260 J=4,8,2 2208 PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1) 2209 260 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2) 2210 BETA=ZERO 2211 DENMAX=ZERO 2212 DO 270 J=1,9 2213 BETA=BETA+DEN(J)*PAR(J) 2214 270 DENMAX=DENMAX+DENEX(J)*PAR(J) 2215 DO 280 K=1,NDIM 2216 VLAG(K)=ZERO 2217 DO 280 J=1,5 2218 280 VLAG(K)=VLAG(K)+PROD(K,J)*PAR(J) 2219 TAU=VLAG(KNEW) 2220 DD=ZERO 2221 TEMPA=ZERO 2222 TEMPB=ZERO 2223 DO 290 I=1,N 2224 D(I)=PAR(2)*D(I)+PAR(3)*S(I) 2225 W(I)=XOPT(I)+D(I) 2226 DD=DD+D(I)**2 2227 TEMPA=TEMPA+D(I)*W(I) 2228 290 TEMPB=TEMPB+W(I)*W(I) 2229 IF (ITERC .GE. N) GOTO 340 2230 IF (ITERC .GT. 1) DENSAV=DMAX1(DENSAV,DENOLD) 2231 IF (DABS(DENMAX) .LE. 1.1D0*DABS(DENSAV)) GOTO 340 2232 DENSAV=DENMAX 2233! 2234! Set S to half the gradient of the denominator with respect to D. 2235! Then branch for the next iteration. 2236! 2237 DO 300 I=1,N 2238 TEMP=TEMPA*XOPT(I)+TEMPB*D(I)-VLAG(NPT+I) 2239 300 S(I)=TAU*BMAT(KNEW,I)+ALPHA*TEMP 2240 DO 320 K=1,NPT 2241 SUM=ZERO 2242 DO 310 J=1,N 2243 310 SUM=SUM+XPT(K,J)*W(J) 2244 TEMP=(TAU*W(N+K)-ALPHA*VLAG(K))*SUM 2245 DO 320 I=1,N 2246 320 S(I)=S(I)+TEMP*XPT(K,I) 2247 SS=ZERO 2248 DS=ZERO 2249 DO 330 I=1,N 2250 SS=SS+S(I)**2 2251 330 DS=DS+D(I)*S(I) 2252 SSDEN=DD*SS-DS*DS 2253 IF (SSDEN .GE. 1.0D-8*DD*SS) GOTO 70 2254! 2255! Set the vector W before the RETURN from the subroutine. 2256! 2257 340 DO 350 K=1,NDIM 2258 W(K)=ZERO 2259 DO 350 J=1,5 2260 350 W(K)=W(K)+WVEC(K,J)*PAR(J) 2261 VLAG(KOPT)=VLAG(KOPT)+ONE 2262 RETURN 2263 END 2264 2265 SUBROUTINE BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW, & 2266 & DELTA,D,ALPHA,HCOL,GC,GD,S,W) 2267 IMPLICIT REAL*8 (A-H,O-Z) 2268 DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*), & 2269 & HCOL(*),GC(*),GD(*),S(*),W(*) 2270! 2271! N is the number of variables. 2272! NPT is the number of interpolation equations. 2273! XOPT is the best interpolation point so far. 2274! XPT contains the coordinates of the current interpolation points. 2275! BMAT provides the last N columns of H. 2276! ZMAT and IDZ give a factorization of the first NPT by NPT submatri 2277! NDIM is the first dimension of BMAT and has the value NPT+N. 2278! KNEW is the index of the interpolation point that is going to be m 2279! DELTA is the current trust region bound. 2280! D will be set to the step from XOPT to the new point. 2281! ALPHA will be set to the KNEW-th diagonal element of the H matrix. 2282! HCOL, GC, GD, S and W will be used for working space. 2283! 2284! The step D is calculated in a way that attempts to maximize the mo 2285! of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFU 2286! the KNEW-th Lagrange function. 2287! 2288! Set some constants. 2289! 2290 HALF=0.5D0 2291 ONE=1.0D0 2292 ZERO=0.0D0 2293 TWOPI=8.0D0*DATAN(ONE) 2294 DELSQ=DELTA*DELTA 2295 NPTM=NPT-N-1 2296! 2297! Set the first NPT components of HCOL to the leading elements of th 2298! KNEW-th column of H. 2299! 2300 ITERC=0 2301 DO 10 K=1,NPT 2302 10 HCOL(K)=ZERO 2303 DO 20 J=1,NPTM 2304 TEMP=ZMAT(KNEW,J) 2305 IF (J .LT. IDZ) TEMP=-TEMP 2306 DO 20 K=1,NPT 2307 20 HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J) 2308 ALPHA=HCOL(KNEW) 2309! 2310! Set the unscaled initial direction D. Form the gradient of LFUNC a 2311! XOPT, and multiply D by the second derivative matrix of LFUNC. 2312! 2313 DD=ZERO 2314 DO 30 I=1,N 2315 D(I)=XPT(KNEW,I)-XOPT(I) 2316 GC(I)=BMAT(KNEW,I) 2317 GD(I)=ZERO 2318 30 DD=DD+D(I)**2 2319 DO 50 K=1,NPT 2320 TEMP=ZERO 2321 SUM=ZERO 2322 DO 40 J=1,N 2323 TEMP=TEMP+XPT(K,J)*XOPT(J) 2324 40 SUM=SUM+XPT(K,J)*D(J) 2325 TEMP=HCOL(K)*TEMP 2326 SUM=HCOL(K)*SUM 2327 DO 50 I=1,N 2328 GC(I)=GC(I)+TEMP*XPT(K,I) 2329 50 GD(I)=GD(I)+SUM*XPT(K,I) 2330! 2331! Scale D and GD, with a sign change if required. Set S to another 2332! vector in the initial two dimensional subspace. 2333! 2334 GG=ZERO 2335 SP=ZERO 2336 DHD=ZERO 2337 DO 60 I=1,N 2338 GG=GG+GC(I)**2 2339 SP=SP+D(I)*GC(I) 2340 60 DHD=DHD+D(I)*GD(I) 2341 SCALE=DELTA/DSQRT(DD) 2342 IF (SP*DHD .LT. ZERO) SCALE=-SCALE 2343 TEMP=ZERO 2344 IF (SP*SP .GT. 0.99D0*DD*GG) TEMP=ONE 2345 TAU=SCALE*(DABS(SP)+HALF*SCALE*DABS(DHD)) 2346 IF (GG*DELSQ .LT. 0.01D0*TAU*TAU) TEMP=ONE 2347 DO 70 I=1,N 2348 D(I)=SCALE*D(I) 2349 GD(I)=SCALE*GD(I) 2350 70 S(I)=GC(I)+TEMP*GD(I) 2351! 2352! Begin the iteration by overwriting S with a vector that has the 2353! required length and direction, except that termination occurs if 2354! the given D and S are nearly parallel. 2355! 2356 80 ITERC=ITERC+1 2357 DD=ZERO 2358 SP=ZERO 2359 SS=ZERO 2360 DO 90 I=1,N 2361 DD=DD+D(I)**2 2362 SP=SP+D(I)*S(I) 2363 90 SS=SS+S(I)**2 2364 TEMP=DD*SS-SP*SP 2365 IF (TEMP .LE. 1.0D-8*DD*SS) GOTO 160 2366 DENOM=DSQRT(TEMP) 2367 DO 100 I=1,N 2368 S(I)=(DD*S(I)-SP*D(I))/DENOM 2369 100 W(I)=ZERO 2370! 2371! Calculate the coefficients of the objective function on the circle 2372! beginning with the multiplication of S by the second derivative ma 2373! 2374 DO 120 K=1,NPT 2375 SUM=ZERO 2376 DO 110 J=1,N 2377 110 SUM=SUM+XPT(K,J)*S(J) 2378 SUM=HCOL(K)*SUM 2379 DO 120 I=1,N 2380 120 W(I)=W(I)+SUM*XPT(K,I) 2381 CF1=ZERO 2382 CF2=ZERO 2383 CF3=ZERO 2384 CF4=ZERO 2385 CF5=ZERO 2386 DO 130 I=1,N 2387 CF1=CF1+S(I)*W(I) 2388 CF2=CF2+D(I)*GC(I) 2389 CF3=CF3+S(I)*GC(I) 2390 CF4=CF4+D(I)*GD(I) 2391 130 CF5=CF5+S(I)*GD(I) 2392 CF1=HALF*CF1 2393 CF4=HALF*CF4-CF1 2394! 2395! Seek the value of the angle that maximizes the modulus of TAU. 2396! 2397 TAUBEG=CF1+CF2+CF4 2398 TAUMAX=TAUBEG 2399 TAUOLD=TAUBEG 2400 ISAVE=0 2401 IU=49 2402 TEMP=TWOPI/DBLE(IU+1) 2403 DO 140 I=1,IU 2404 ANGLE=DBLE(I)*TEMP 2405 CTH=DCOS(ANGLE) 2406 STH=DSIN(ANGLE) 2407 TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH 2408 IF (DABS(TAU) .GT. DABS(TAUMAX)) THEN 2409 TAUMAX=TAU 2410 ISAVE=I 2411 TEMPA=TAUOLD 2412 ELSE IF (I .EQ. ISAVE+1) THEN 2413 TEMPB=TAU 2414 END IF 2415 140 TAUOLD=TAU 2416 IF (ISAVE .EQ. 0) TEMPA=TAU 2417 IF (ISAVE .EQ. IU) TEMPB=TAUBEG 2418 STEP=ZERO 2419 IF (TEMPA .NE. TEMPB) THEN 2420 TEMPA=TEMPA-TAUMAX 2421 TEMPB=TEMPB-TAUMAX 2422 STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB) 2423 END IF 2424 ANGLE=TEMP*(DBLE(ISAVE)+STEP) 2425! 2426! Calculate the new D and GD. Then test for convergence. 2427! 2428 CTH=DCOS(ANGLE) 2429 STH=DSIN(ANGLE) 2430 TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH 2431 DO 150 I=1,N 2432 D(I)=CTH*D(I)+STH*S(I) 2433 GD(I)=CTH*GD(I)+STH*W(I) 2434 150 S(I)=GC(I)+GD(I) 2435 IF (DABS(TAU) .LE. 1.1D0*DABS(TAUBEG)) GOTO 160 2436 IF (ITERC .LT. N) GOTO 80 2437 160 RETURN 2438 END 2439 2440 SUBROUTINE TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,STEP, & 2441 & D,G,HD,HS,CRVMIN) 2442 IMPLICIT REAL*8 (A-H,O-Z) 2443 DIMENSION XOPT(*),XPT(NPT,*),GQ(*),HQ(*),PQ(*),STEP(*), & 2444 & D(*),G(*),HD(*),HS(*) 2445! 2446! N is the number of variables of a quadratic objective function, Q 2447! The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meani 2448! in order to define the current quadratic model Q. 2449! DELTA is the trust region radius, and has to be positive. 2450! STEP will be set to the calculated trial step. 2451! The arrays D, G, HD and HS will be used for working space. 2452! CRVMIN will be set to the least curvature of H along the conjugate 2453! directions that occur, except that it is set to zero if STEP goe 2454! all the way to the trust region boundary. 2455! 2456! The calculation of STEP begins with the truncated conjugate gradie 2457! method. If the boundary of the trust region is reached, then furth 2458! changes to STEP may be made, each one being in the 2D space spanne 2459! by the current STEP and the corresponding gradient of Q. Thus STEP 2460! should provide a substantial reduction to Q within the trust regio 2461! 2462! Initialization, which includes setting HD to H times XOPT. 2463! 2464 HALF=0.5D0 2465 ZERO=0.0D0 2466 TWOPI=8.0D0*DATAN(1.0D0) 2467 DELSQ=DELTA*DELTA 2468 ITERC=0 2469 ITERMAX=N 2470 ITERSW=ITERMAX 2471 DO 10 I=1,N 2472 10 D(I)=XOPT(I) 2473 GOTO 170 2474! 2475! Prepare for the first line search. 2476! 2477 20 QRED=ZERO 2478 DD=ZERO 2479 DO 30 I=1,N 2480 STEP(I)=ZERO 2481 HS(I)=ZERO 2482 G(I)=GQ(I)+HD(I) 2483 D(I)=-G(I) 2484 30 DD=DD+D(I)**2 2485 CRVMIN=ZERO 2486 IF (DD .EQ. ZERO) GOTO 160 2487 DS=ZERO 2488 SS=ZERO 2489 GG=DD 2490 GGBEG=GG 2491! 2492! Calculate the step to the trust region boundary and the product HD 2493! 2494 40 ITERC=ITERC+1 2495 TEMP=DELSQ-SS 2496 BSTEP=TEMP/(DS+DSQRT(DS*DS+DD*TEMP)) 2497 GOTO 170 2498 50 DHD=ZERO 2499 DO 60 J=1,N 2500 60 DHD=DHD+D(J)*HD(J) 2501! 2502! Update CRVMIN and set the step-length ALPHA. 2503! 2504 ALPHA=BSTEP 2505 IF (DHD .GT. ZERO) THEN 2506 TEMP=DHD/DD 2507 IF (ITERC .EQ. 1) CRVMIN=TEMP 2508 CRVMIN=DMIN1(CRVMIN,TEMP) 2509 ALPHA=DMIN1(ALPHA,GG/DHD) 2510 END IF 2511 QADD=ALPHA*(GG-HALF*ALPHA*DHD) 2512 QRED=QRED+QADD 2513! 2514! Update STEP and HS. 2515! 2516 GGSAV=GG 2517 GG=ZERO 2518 DO 70 I=1,N 2519 STEP(I)=STEP(I)+ALPHA*D(I) 2520 HS(I)=HS(I)+ALPHA*HD(I) 2521 70 GG=GG+(G(I)+HS(I))**2 2522! 2523! Begin another conjugate direction iteration if required. 2524! 2525 IF (ALPHA .LT. BSTEP) THEN 2526 IF (QADD .LE. 0.01D0*QRED) GOTO 160 2527 IF (GG .LE. 1.0D-4*GGBEG) GOTO 160 2528 IF (ITERC .EQ. ITERMAX) GOTO 160 2529 TEMP=GG/GGSAV 2530 DD=ZERO 2531 DS=ZERO 2532 SS=ZERO 2533 DO 80 I=1,N 2534 D(I)=TEMP*D(I)-G(I)-HS(I) 2535 DD=DD+D(I)**2 2536 DS=DS+D(I)*STEP(I) 2537 80 SS=SS+STEP(I)**2 2538 IF (DS .LE. ZERO) GOTO 160 2539 IF (SS .LT. DELSQ) GOTO 40 2540 END IF 2541 CRVMIN=ZERO 2542 ITERSW=ITERC 2543! 2544! Test whether an alternative iteration is required. 2545! 2546 90 IF (GG .LE. 1.0D-4*GGBEG) GOTO 160 2547 SG=ZERO 2548 SHS=ZERO 2549 DO 100 I=1,N 2550 SG=SG+STEP(I)*G(I) 2551 100 SHS=SHS+STEP(I)*HS(I) 2552 SGK=SG+SHS 2553 ANGTEST=SGK/DSQRT(GG*DELSQ) 2554 IF (ANGTEST .LE. -0.99D0) GOTO 160 2555! 2556! Begin the alternative iteration by calculating D and HD and some 2557! scalar products. 2558! 2559 ITERC=ITERC+1 2560 TEMP=DSQRT(DELSQ*GG-SGK*SGK) 2561 TEMPA=DELSQ/TEMP 2562 TEMPB=SGK/TEMP 2563 DO 110 I=1,N 2564 110 D(I)=TEMPA*(G(I)+HS(I))-TEMPB*STEP(I) 2565 GOTO 170 2566 120 DG=ZERO 2567 DHD=ZERO 2568 DHS=ZERO 2569 DO 130 I=1,N 2570 DG=DG+D(I)*G(I) 2571 DHD=DHD+HD(I)*D(I) 2572 130 DHS=DHS+HD(I)*STEP(I) 2573! 2574! Seek the value of the angle that minimizes Q. 2575! 2576 CF=HALF*(SHS-DHD) 2577 QBEG=SG+CF 2578 QSAV=QBEG 2579 QMIN=QBEG 2580 ISAVE=0 2581 IU=49 2582 TEMP=TWOPI/DBLE(IU+1) 2583 DO 140 I=1,IU 2584 ANGLE=DBLE(I)*TEMP 2585 CTH=DCOS(ANGLE) 2586 STH=DSIN(ANGLE) 2587 QNEW=(SG+CF*CTH)*CTH+(DG+DHS*CTH)*STH 2588 IF (QNEW .LT. QMIN) THEN 2589 QMIN=QNEW 2590 ISAVE=I 2591 TEMPA=QSAV 2592 ELSE IF (I .EQ. ISAVE+1) THEN 2593 TEMPB=QNEW 2594 END IF 2595 140 QSAV=QNEW 2596 IF (ISAVE .EQ. ZERO) TEMPA=QNEW 2597 IF (ISAVE .EQ. IU) TEMPB=QBEG 2598 ANGLE=ZERO 2599 IF (TEMPA .NE. TEMPB) THEN 2600 TEMPA=TEMPA-QMIN 2601 TEMPB=TEMPB-QMIN 2602 ANGLE=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB) 2603 END IF 2604 ANGLE=TEMP*(DBLE(ISAVE)+ANGLE) 2605! 2606! Calculate the new STEP and HS. Then test for convergence. 2607! 2608 CTH=DCOS(ANGLE) 2609 STH=DSIN(ANGLE) 2610 REDUC=QBEG-(SG+CF*CTH)*CTH-(DG+DHS*CTH)*STH 2611 GG=ZERO 2612 DO 150 I=1,N 2613 STEP(I)=CTH*STEP(I)+STH*D(I) 2614 HS(I)=CTH*HS(I)+STH*HD(I) 2615 150 GG=GG+(G(I)+HS(I))**2 2616 QRED=QRED+REDUC 2617 RATIO=REDUC/QRED 2618 IF (ITERC .LT. ITERMAX .AND. RATIO .GT. 0.01D0) GOTO 90 2619 160 RETURN 2620! 2621! The following instructions act as a subroutine for setting the vec 2622! HD to the vector D multiplied by the second derivative matrix of Q 2623! They are called from three different places, which are distinguish 2624! by the value of ITERC. 2625! 2626 170 DO 180 I=1,N 2627 180 HD(I)=ZERO 2628 DO 200 K=1,NPT 2629 TEMP=ZERO 2630 DO 190 J=1,N 2631 190 TEMP=TEMP+XPT(K,J)*D(J) 2632 TEMP=TEMP*PQ(K) 2633 DO 200 I=1,N 2634 200 HD(I)=HD(I)+TEMP*XPT(K,I) 2635 IH=0 2636 DO 210 J=1,N 2637 DO 210 I=1,J 2638 IH=IH+1 2639 IF (I .LT. J) HD(J)=HD(J)+HQ(IH)*D(I) 2640 210 HD(I)=HD(I)+HQ(IH)*D(J) 2641 IF (ITERC .EQ. 0) GOTO 20 2642 IF (ITERC .LE. ITERSW) GOTO 50 2643 GOTO 120 2644 END 2645 2646 SUBROUTINE UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W) 2647 IMPLICIT REAL*8 (A-H,O-Z) 2648 DIMENSION BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*) 2649! 2650! The arrays BMAT and ZMAT with IDZ are updated, in order to shift t 2651! interpolation point that has index KNEW. On entry, VLAG contains t 2652! components of the vector Theta*Wcheck+e_b of the updating formula 2653! (6.11), and BETA holds the value of the parameter that has this na 2654! The vector W is used for working space. 2655! 2656! Set some constants. 2657! 2658 ONE=1.0D0 2659 ZERO=0.0D0 2660 NPTM=NPT-N-1 2661! 2662! Apply the rotations that put zeros in the KNEW-th row of ZMAT. 2663! 2664 JL=1 2665 DO 20 J=2,NPTM 2666 IF (J .EQ. IDZ) THEN 2667 JL=IDZ 2668 ELSE IF (ZMAT(KNEW,J) .NE. ZERO) THEN 2669 TEMP=DSQRT(ZMAT(KNEW,JL)**2+ZMAT(KNEW,J)**2) 2670 TEMPA=ZMAT(KNEW,JL)/TEMP 2671 TEMPB=ZMAT(KNEW,J)/TEMP 2672 DO 10 I=1,NPT 2673 TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J) 2674 ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL) 2675 10 ZMAT(I,JL)=TEMP 2676 ZMAT(KNEW,J)=ZERO 2677 END IF 2678 20 END DO 2679! 2680! Put the first NPT components of the KNEW-th column of HLAG into W, 2681! and calculate the parameters of the updating formula. 2682! 2683 TEMPA=ZMAT(KNEW,1) 2684 IF (IDZ .GE. 2) TEMPA=-TEMPA 2685 IF (JL .GT. 1) TEMPB=ZMAT(KNEW,JL) 2686 DO 30 I=1,NPT 2687 W(I)=TEMPA*ZMAT(I,1) 2688 IF (JL .GT. 1) W(I)=W(I)+TEMPB*ZMAT(I,JL) 2689 30 END DO 2690 ALPHA=W(KNEW) 2691 TAU=VLAG(KNEW) 2692 TAUSQ=TAU*TAU 2693 DENOM=ALPHA*BETA+TAUSQ 2694 VLAG(KNEW)=VLAG(KNEW)-ONE 2695! 2696! Complete the updating of ZMAT when there is only one nonzero eleme 2697! in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to 2698! then the first column of ZMAT will be exchanged with another one l 2699! 2700 IFLAG=0 2701 IF (JL .EQ. 1) THEN 2702 TEMP=DSQRT(DABS(DENOM)) 2703 TEMPB=TEMPA/TEMP 2704 TEMPA=TAU/TEMP 2705 DO 40 I=1,NPT 2706 40 ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I) 2707 IF (IDZ .EQ. 1 .AND. TEMP .LT. ZERO) IDZ=2 2708 IF (IDZ .GE. 2 .AND. TEMP .GE. ZERO) IFLAG=1 2709 ELSE 2710! 2711! Complete the updating of ZMAT in the alternative case. 2712! 2713 JA=1 2714 IF (BETA .GE. ZERO) JA=JL 2715 JB=JL+1-JA 2716 TEMP=ZMAT(KNEW,JB)/DENOM 2717 TEMPA=TEMP*BETA 2718 TEMPB=TEMP*TAU 2719 TEMP=ZMAT(KNEW,JA) 2720 SCALA=ONE/DSQRT(DABS(BETA)*TEMP*TEMP+TAUSQ) 2721 SCALB=SCALA*DSQRT(DABS(DENOM)) 2722 DO 50 I=1,NPT 2723 ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I)) 2724 50 ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I)) 2725 IF (DENOM .LE. ZERO) THEN 2726 IF (BETA .LT. ZERO) IDZ=IDZ+1 2727 IF (BETA .GE. ZERO) IFLAG=1 2728 END IF 2729 END IF 2730! 2731! IDZ is reduced in the following case, and usually the first column 2732! of ZMAT is exchanged with a later one. 2733! 2734 IF (IFLAG .EQ. 1) THEN 2735 IDZ=IDZ-1 2736 DO 60 I=1,NPT 2737 TEMP=ZMAT(I,1) 2738 ZMAT(I,1)=ZMAT(I,IDZ) 2739 60 ZMAT(I,IDZ)=TEMP 2740 END IF 2741! 2742! Finally, update the matrix BMAT. 2743! 2744 DO 70 J=1,N 2745 JP=NPT+J 2746 W(JP)=BMAT(KNEW,J) 2747 TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM 2748 TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM 2749 DO 70 I=1,JP 2750 BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I) 2751 IF (I .GT. NPT) BMAT(JP,I-NPT)=BMAT(I,J) 2752 70 CONTINUE 2753 RETURN 2754 END 2755