1 // Copyright (C) 2009 M.J.D. Powell, Davis E. King (davis@dlib.net) 2 // License: Boost Software License See LICENSE.txt for the full license. 3 #ifndef DLIB_OPTIMIZATIOn_BOBYQA_Hh_ 4 #define DLIB_OPTIMIZATIOn_BOBYQA_Hh_ 5 6 /* 7 The code in this file is derived from Powell's BOBYQA Fortran code. 8 It was created by running f2c on the original Fortran code and then 9 massaging the resulting C code into what you can see below. 10 11 12 The following paper, published in 2009 by Powell, describes the 13 detailed workings of the BOBYQA algorithm. 14 15 The BOBYQA algorithm for bound constrained optimization 16 without derivatives by M.J.D. Powell 17 */ 18 19 #include <algorithm> 20 #include <cmath> 21 #include <memory> 22 23 #include "../matrix.h" 24 #include "optimization_bobyqa_abstract.h" 25 #include "optimization.h" 26 27 // ---------------------------------------------------------------------------------------- 28 29 namespace dlib 30 { 31 32 class bobyqa_failure : public error { bobyqa_failure(const std::string & s)33 public: bobyqa_failure(const std::string& s):error(s){} 34 }; 35 36 // ---------------------------------------------------------------------------------------- 37 38 class bobyqa_implementation 39 { 40 typedef long integer; 41 typedef double doublereal; 42 43 public: 44 45 template < 46 typename funct, 47 typename T, 48 typename U 49 > find_min(const funct & f,T & x,long npt,const U & xl_,const U & xu_,const double rhobeg,const double rhoend,const long max_f_evals)50 double find_min ( 51 const funct& f, 52 T& x, 53 long npt, 54 const U& xl_, 55 const U& xu_, 56 const double rhobeg, 57 const double rhoend, 58 const long max_f_evals 59 ) const 60 { 61 const unsigned long n = x.size(); 62 const unsigned long w_size = (npt+5)*(npt+n)+3*n*(n+5)/2; 63 std::unique_ptr<doublereal[]> w(new doublereal[w_size]); 64 65 // make these temporary matrices becuse U might be some 66 // kind of matrix_exp that doesn't support taking the address 67 // of the first element. 68 matrix<double,0,1> xl(xl_); 69 matrix<double,0,1> xu(xu_); 70 71 72 return bobyqa_ (f, 73 x.size(), 74 npt, 75 &x(0), 76 &xl(0), 77 &xu(0), 78 rhobeg, 79 rhoend, 80 max_f_evals, 81 w.get() ); 82 } 83 84 private: 85 86 87 template <typename funct> bobyqa_(const funct & calfun,const integer n,const integer npt,doublereal * x,const doublereal * xl,const doublereal * xu,const doublereal rhobeg,const doublereal rhoend,const integer maxfun,doublereal * w)88 doublereal bobyqa_( 89 const funct& calfun, 90 const integer n, 91 const integer npt, 92 doublereal *x, 93 const doublereal *xl, 94 const doublereal *xu, 95 const doublereal rhobeg, 96 const doublereal rhoend, 97 const integer maxfun, 98 doublereal *w 99 ) const 100 { 101 102 /* System generated locals */ 103 integer i__1; 104 doublereal d__1, d__2; 105 106 /* Local variables */ 107 integer j, id_, np, iw, igo, ihq, ixb, ixa, ifv, isl, jsl, ipq, ivl, ixn, ixo, ixp, isu, jsu, ndim; 108 doublereal temp, zero; 109 integer ibmat, izmat; 110 111 112 /* This subroutine seeks the least value of a function of many variables, */ 113 /* by applying a trust region method that forms quadratic models by */ 114 /* interpolation. There is usually some freedom in the interpolation */ 115 /* conditions, which is taken up by minimizing the Frobenius norm of */ 116 /* the change to the second derivative of the model, beginning with the */ 117 /* zero matrix. The values of the variables are constrained by upper and */ 118 /* lower bounds. The arguments of the subroutine are as follows. */ 119 120 /* N must be set to the number of variables and must be at least two. */ 121 /* NPT is the number of interpolation conditions. Its value must be in */ 122 /* the interval [N+2,(N+1)(N+2)/2]. Choices that exceed 2*N+1 are not */ 123 /* recommended. */ 124 /* Initial values of the variables must be set in X(1),X(2),...,X(N). They */ 125 /* will be changed to the values that give the least calculated F. */ 126 /* For I=1,2,...,N, XL(I) and XU(I) must provide the lower and upper */ 127 /* bounds, respectively, on X(I). The construction of quadratic models */ 128 /* requires XL(I) to be strictly less than XU(I) for each I. Further, */ 129 /* the contribution to a model from changes to the I-th variable is */ 130 /* damaged severely by rounding errors if XU(I)-XL(I) is too small. */ 131 /* RHOBEG and RHOEND must be set to the initial and final values of a trust */ 132 /* region radius, so both must be positive with RHOEND no greater than */ 133 /* RHOBEG. Typically, RHOBEG should be about one tenth of the greatest */ 134 /* expected change to a variable, while RHOEND should indicate the */ 135 /* accuracy that is required in the final values of the variables. An */ 136 /* error return occurs if any of the differences XU(I)-XL(I), I=1,...,N, */ 137 /* is less than 2*RHOBEG. */ 138 /* MAXFUN must be set to an upper bound on the number of calls of CALFUN. */ 139 /* The array W will be used for working space. Its length must be at least */ 140 /* (NPT+5)*(NPT+N)+3*N*(N+5)/2. */ 141 142 /* Parameter adjustments */ 143 --w; 144 --xu; 145 --xl; 146 --x; 147 148 /* Function Body */ 149 np = n + 1; 150 151 /* Return if the value of NPT is unacceptable. */ 152 if (npt < n + 2 || npt > (n + 2) * np / 2) { 153 throw bobyqa_failure("Return from BOBYQA because NPT is not in the required interval"); 154 //goto L40; 155 } 156 157 /* Partition the working space array, so that different parts of it can */ 158 /* be treated separately during the calculation of BOBYQB. The partition */ 159 /* requires the first (NPT+2)*(NPT+N)+3*N*(N+5)/2 elements of W plus the */ 160 /* space that is taken by the last array in the argument list of BOBYQB. */ 161 162 ndim = npt + n; 163 ixb = 1; 164 ixp = ixb + n; 165 ifv = ixp + n * npt; 166 ixo = ifv + npt; 167 igo = ixo + n; 168 ihq = igo + n; 169 ipq = ihq + n * np / 2; 170 ibmat = ipq + npt; 171 izmat = ibmat + ndim * n; 172 isl = izmat + npt * (npt - np); 173 isu = isl + n; 174 ixn = isu + n; 175 ixa = ixn + n; 176 id_ = ixa + n; 177 ivl = id_ + n; 178 iw = ivl + ndim; 179 180 /* Return if there is insufficient space between the bounds. Modify the */ 181 /* initial X if necessary in order to avoid conflicts between the bounds */ 182 /* and the construction of the first quadratic model. The lower and upper */ 183 /* bounds on moves from the updated X are set now, in the ISL and ISU */ 184 /* partitions of W, in order to provide useful and exact information about */ 185 /* components of X that become within distance RHOBEG from their bounds. */ 186 187 zero = 0.; 188 i__1 = n; 189 for (j = 1; j <= i__1; ++j) { 190 temp = xu[j] - xl[j]; 191 if (temp < rhobeg + rhobeg) { 192 throw bobyqa_failure("Return from BOBYQA because one of the differences in x_lower and x_upper is less than 2*rho_begin"); 193 //goto L40; 194 } 195 jsl = isl + j - 1; 196 jsu = jsl + n; 197 w[jsl] = xl[j] - x[j]; 198 w[jsu] = xu[j] - x[j]; 199 if (w[jsl] >= -(rhobeg)) { 200 if (w[jsl] >= zero) { 201 x[j] = xl[j]; 202 w[jsl] = zero; 203 w[jsu] = temp; 204 } else { 205 x[j] = xl[j] + rhobeg; 206 w[jsl] = -(rhobeg); 207 /* Computing MAX */ 208 d__1 = xu[j] - x[j]; 209 w[jsu] = std::max(d__1,rhobeg); 210 } 211 } else if (w[jsu] <= rhobeg) { 212 if (w[jsu] <= zero) { 213 x[j] = xu[j]; 214 w[jsl] = -temp; 215 w[jsu] = zero; 216 } else { 217 x[j] = xu[j] - rhobeg; 218 /* Computing MIN */ 219 d__1 = xl[j] - x[j], d__2 = -(rhobeg); 220 w[jsl] = std::min(d__1,d__2); 221 w[jsu] = rhobeg; 222 } 223 } 224 /* L30: */ 225 } 226 227 /* Make the call of BOBYQB. */ 228 229 return bobyqb_(calfun, n, npt, &x[1], &xl[1], &xu[1], rhobeg, rhoend, maxfun, &w[ 230 ixb], &w[ixp], &w[ifv], &w[ixo], &w[igo], &w[ihq], &w[ipq], &w[ 231 ibmat], &w[izmat], ndim, &w[isl], &w[isu], &w[ixn], &w[ixa], &w[ 232 id_], &w[ivl], &w[iw]); 233 //L40: 234 ; 235 } /* bobyqa_ */ 236 237 // ---------------------------------------------------------------------------------------- 238 239 template <typename funct> bobyqb_(const funct & calfun,const integer n,const integer npt,doublereal * x,const doublereal * xl,const doublereal * xu,const doublereal rhobeg,const doublereal rhoend,const integer maxfun,doublereal * xbase,doublereal * xpt,doublereal * fval,doublereal * xopt,doublereal * gopt,doublereal * hq,doublereal * pq,doublereal * bmat,doublereal * zmat,const integer ndim,doublereal * sl,doublereal * su,doublereal * xnew,doublereal * xalt,doublereal * d__,doublereal * vlag,doublereal * w)240 doublereal bobyqb_( 241 const funct& calfun, 242 const integer n, 243 const integer npt, 244 doublereal *x, 245 const doublereal *xl, 246 const doublereal *xu, 247 const doublereal rhobeg, 248 const doublereal rhoend, 249 const integer maxfun, 250 doublereal *xbase, 251 doublereal *xpt, 252 doublereal *fval, 253 doublereal *xopt, 254 doublereal *gopt, 255 doublereal *hq, 256 doublereal *pq, 257 doublereal *bmat, 258 doublereal *zmat, 259 const integer ndim, 260 doublereal *sl, 261 doublereal *su, 262 doublereal *xnew, 263 doublereal *xalt, 264 doublereal *d__, 265 doublereal *vlag, 266 doublereal *w 267 ) const 268 { 269 /* System generated locals */ 270 integer xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, 271 zmat_offset, i__1, i__2, i__3; 272 doublereal d__1, d__2, d__3, d__4; 273 274 /* Local variables */ 275 doublereal f = 0; 276 integer i__, j, k, ih, nf, jj, nh, ip, jp; 277 doublereal dx; 278 integer np; 279 doublereal den = 0, one = 0, ten = 0, dsq = 0, rho = 0, sum = 0, two = 0, diff = 0, half = 0, beta = 0, gisq = 0; 280 integer knew = 0; 281 doublereal temp, suma, sumb, bsum, fopt; 282 integer kopt = 0, nptm; 283 doublereal zero, curv; 284 integer ksav; 285 doublereal gqsq = 0, dist = 0, sumw = 0, sumz = 0, diffa = 0, diffb = 0, diffc = 0, hdiag = 0; 286 integer kbase; 287 doublereal alpha = 0, delta = 0, adelt = 0, denom = 0, fsave = 0, bdtol = 0, delsq = 0; 288 integer nresc, nfsav; 289 doublereal ratio = 0, dnorm = 0, vquad = 0, pqold = 0, tenth = 0; 290 integer itest; 291 doublereal sumpq, scaden; 292 doublereal errbig, cauchy=0, fracsq, biglsq, densav; 293 doublereal bdtest; 294 doublereal crvmin, frhosq; 295 doublereal distsq; 296 integer ntrits; 297 doublereal xoptsq; 298 299 300 301 /* The arguments N, NPT, X, XL, XU, RHOBEG, RHOEND, IPRINT and MAXFUN */ 302 /* are identical to the corresponding arguments in SUBROUTINE BOBYQA. */ 303 /* XBASE holds a shift of origin that should reduce the contributions */ 304 /* from rounding errors to values of the model and Lagrange functions. */ 305 /* XPT is a two-dimensional array that holds the coordinates of the */ 306 /* interpolation points relative to XBASE. */ 307 /* FVAL holds the values of F at the interpolation points. */ 308 /* XOPT is set to the displacement from XBASE of the trust region centre. */ 309 /* GOPT holds the gradient of the quadratic model at XBASE+XOPT. */ 310 /* HQ holds the explicit second derivatives of the quadratic model. */ 311 /* PQ contains the parameters of the implicit second derivatives of the */ 312 /* quadratic model. */ 313 /* BMAT holds the last N columns of H. */ 314 /* ZMAT holds the factorization of the leading NPT by NPT submatrix of H, */ 315 /* this factorization being ZMAT times ZMAT^T, which provides both the */ 316 /* correct rank and positive semi-definiteness. */ 317 /* NDIM is the first dimension of BMAT and has the value NPT+N. */ 318 /* SL and SU hold the differences XL-XBASE and XU-XBASE, respectively. */ 319 /* All the components of every XOPT are going to satisfy the bounds */ 320 /* SL(I) .LEQ. XOPT(I) .LEQ. SU(I), with appropriate equalities when */ 321 /* XOPT is on a constraint boundary. */ 322 /* XNEW is chosen by SUBROUTINE TRSBOX or ALTMOV. Usually XBASE+XNEW is the */ 323 /* vector of variables for the next call of CALFUN. XNEW also satisfies */ 324 /* the SL and SU constraints in the way that has just been mentioned. */ 325 /* XALT is an alternative to XNEW, chosen by ALTMOV, that may replace XNEW */ 326 /* in order to increase the denominator in the updating of UPDATE. */ 327 /* D is reserved for a trial step from XOPT, which is usually XNEW-XOPT. */ 328 /* VLAG contains the values of the Lagrange functions at a new point X. */ 329 /* They are part of a product that requires VLAG to be of length NDIM. */ 330 /* W is a one-dimensional array that is used for working space. Its length */ 331 /* must be at least 3*NDIM = 3*(NPT+N). */ 332 333 /* Set some constants. */ 334 335 /* Parameter adjustments */ 336 zmat_dim1 = npt; 337 zmat_offset = 1 + zmat_dim1; 338 zmat -= zmat_offset; 339 xpt_dim1 = npt; 340 xpt_offset = 1 + xpt_dim1; 341 xpt -= xpt_offset; 342 --x; 343 --xl; 344 --xu; 345 --xbase; 346 --fval; 347 --xopt; 348 --gopt; 349 --hq; 350 --pq; 351 bmat_dim1 = ndim; 352 bmat_offset = 1 + bmat_dim1; 353 bmat -= bmat_offset; 354 --sl; 355 --su; 356 --xnew; 357 --xalt; 358 --d__; 359 --vlag; 360 --w; 361 362 /* Function Body */ 363 half = .5; 364 one = 1.; 365 ten = 10.; 366 tenth = .1; 367 two = 2.; 368 zero = 0.; 369 np = n + 1; 370 nptm = npt - np; 371 nh = n * np / 2; 372 373 /* The call of PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */ 374 /* BMAT and ZMAT for the first iteration, with the corresponding values of */ 375 /* of NF and KOPT, which are the number of calls of CALFUN so far and the */ 376 /* index of the interpolation point at the trust region centre. Then the */ 377 /* initial XOPT is set too. The branch to label 720 occurs if MAXFUN is */ 378 /* less than NPT. GOPT will be updated if KOPT is different from KBASE. */ 379 380 prelim_(calfun, n, npt, &x[1], &xl[1], &xu[1], rhobeg, maxfun, &xbase[1], 381 &xpt[xpt_offset], &fval[1], &gopt[1], &hq[1], &pq[1], &bmat[bmat_offset], 382 &zmat[zmat_offset], ndim, &sl[1], &su[1], nf, kopt); 383 xoptsq = zero; 384 i__1 = n; 385 for (i__ = 1; i__ <= i__1; ++i__) { 386 xopt[i__] = xpt[kopt + i__ * xpt_dim1]; 387 /* L10: */ 388 /* Computing 2nd power */ 389 d__1 = xopt[i__]; 390 xoptsq += d__1 * d__1; 391 } 392 fsave = fval[1]; 393 if (nf < npt) { 394 throw bobyqa_failure("Return from BOBYQA because the objective function has been called max_f_evals times."); 395 //goto L720; 396 } 397 kbase = 1; 398 399 /* Complete the settings that are required for the iterative procedure. */ 400 401 rho = rhobeg; 402 delta = rho; 403 nresc = nf; 404 ntrits = 0; 405 diffa = zero; 406 diffb = zero; 407 itest = 0; 408 nfsav = nf; 409 410 /* Update GOPT if necessary before the first iteration and after each */ 411 /* call of RESCUE that makes a call of CALFUN. */ 412 413 L20: 414 if (kopt != kbase) { 415 ih = 0; 416 i__1 = n; 417 for (j = 1; j <= i__1; ++j) { 418 i__2 = j; 419 for (i__ = 1; i__ <= i__2; ++i__) { 420 ++ih; 421 if (i__ < j) { 422 gopt[j] += hq[ih] * xopt[i__]; 423 } 424 /* L30: */ 425 gopt[i__] += hq[ih] * xopt[j]; 426 } 427 } 428 if (nf > npt) { 429 i__2 = npt; 430 for (k = 1; k <= i__2; ++k) { 431 temp = zero; 432 i__1 = n; 433 for (j = 1; j <= i__1; ++j) { 434 /* L40: */ 435 temp += xpt[k + j * xpt_dim1] * xopt[j]; 436 } 437 temp = pq[k] * temp; 438 i__1 = n; 439 for (i__ = 1; i__ <= i__1; ++i__) { 440 /* L50: */ 441 gopt[i__] += temp * xpt[k + i__ * xpt_dim1]; 442 } 443 } 444 } 445 } 446 447 /* Generate the next point in the trust region that provides a small value */ 448 /* of the quadratic model subject to the constraints on the variables. */ 449 /* The integer NTRITS is set to the number "trust region" iterations that */ 450 /* have occurred since the last "alternative" iteration. If the length */ 451 /* of XNEW-XOPT is less than HALF*RHO, however, then there is a branch to */ 452 /* label 650 or 680 with NTRITS=-1, instead of calculating F at XNEW. */ 453 454 L60: 455 trsbox_(n, npt, &xpt[xpt_offset], &xopt[1], &gopt[1], &hq[1], &pq[1], &sl[1], 456 &su[1], delta, &xnew[1], &d__[1], &w[1], &w[np], &w[np + n], 457 &w[np + (n << 1)], &w[np + n * 3], &dsq, &crvmin); 458 /* Computing MIN */ 459 d__1 = delta, d__2 = std::sqrt(dsq); 460 dnorm = std::min(d__1,d__2); 461 if (dnorm < half * rho) { 462 ntrits = -1; 463 /* Computing 2nd power */ 464 d__1 = ten * rho; 465 distsq = d__1 * d__1; 466 if (nf <= nfsav + 2) { 467 goto L650; 468 } 469 470 /* The following choice between labels 650 and 680 depends on whether or */ 471 /* not our work with the current RHO seems to be complete. Either RHO is */ 472 /* decreased or termination occurs if the errors in the quadratic model at */ 473 /* the last three interpolation points compare favourably with predictions */ 474 /* of likely improvements to the model within distance HALF*RHO of XOPT. */ 475 476 /* Computing MAX */ 477 d__1 = std::max(diffa,diffb); 478 errbig = std::max(d__1,diffc); 479 frhosq = rho * .125 * rho; 480 if (crvmin > zero && errbig > frhosq * crvmin) { 481 goto L650; 482 } 483 bdtol = errbig / rho; 484 i__1 = n; 485 for (j = 1; j <= i__1; ++j) { 486 bdtest = bdtol; 487 if (xnew[j] == sl[j]) { 488 bdtest = w[j]; 489 } 490 if (xnew[j] == su[j]) { 491 bdtest = -w[j]; 492 } 493 if (bdtest < bdtol) { 494 curv = hq[(j + j * j) / 2]; 495 i__2 = npt; 496 for (k = 1; k <= i__2; ++k) { 497 /* L70: */ 498 /* Computing 2nd power */ 499 d__1 = xpt[k + j * xpt_dim1]; 500 curv += pq[k] * (d__1 * d__1); 501 } 502 bdtest += half * curv * rho; 503 if (bdtest < bdtol) { 504 goto L650; 505 } 506 } 507 /* L80: */ 508 } 509 goto L680; 510 } 511 ++ntrits; 512 513 /* Severe cancellation is likely to occur if XOPT is too far from XBASE. */ 514 /* If the following test holds, then XBASE is shifted so that XOPT becomes */ 515 /* zero. The appropriate changes are made to BMAT and to the second */ 516 /* derivatives of the current model, beginning with the changes to BMAT */ 517 /* that do not depend on ZMAT. VLAG is used temporarily for working space. */ 518 519 L90: 520 if (dsq <= xoptsq * .001) { 521 fracsq = xoptsq * .25; 522 sumpq = zero; 523 i__1 = npt; 524 for (k = 1; k <= i__1; ++k) { 525 sumpq += pq[k]; 526 sum = -half * xoptsq; 527 i__2 = n; 528 for (i__ = 1; i__ <= i__2; ++i__) { 529 /* L100: */ 530 sum += xpt[k + i__ * xpt_dim1] * xopt[i__]; 531 } 532 w[npt + k] = sum; 533 temp = fracsq - half * sum; 534 i__2 = n; 535 for (i__ = 1; i__ <= i__2; ++i__) { 536 w[i__] = bmat[k + i__ * bmat_dim1]; 537 vlag[i__] = sum * xpt[k + i__ * xpt_dim1] + temp * xopt[i__]; 538 ip = npt + i__; 539 i__3 = i__; 540 for (j = 1; j <= i__3; ++j) { 541 /* L110: */ 542 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] + w[ 543 i__] * vlag[j] + vlag[i__] * w[j]; 544 } 545 } 546 } 547 548 /* Then the revisions of BMAT that depend on ZMAT are calculated. */ 549 550 i__3 = nptm; 551 for (jj = 1; jj <= i__3; ++jj) { 552 sumz = zero; 553 sumw = zero; 554 i__2 = npt; 555 for (k = 1; k <= i__2; ++k) { 556 sumz += zmat[k + jj * zmat_dim1]; 557 vlag[k] = w[npt + k] * zmat[k + jj * zmat_dim1]; 558 /* L120: */ 559 sumw += vlag[k]; 560 } 561 i__2 = n; 562 for (j = 1; j <= i__2; ++j) { 563 sum = (fracsq * sumz - half * sumw) * xopt[j]; 564 i__1 = npt; 565 for (k = 1; k <= i__1; ++k) { 566 /* L130: */ 567 sum += vlag[k] * xpt[k + j * xpt_dim1]; 568 } 569 w[j] = sum; 570 i__1 = npt; 571 for (k = 1; k <= i__1; ++k) { 572 /* L140: */ 573 bmat[k + j * bmat_dim1] += sum * zmat[k + jj * zmat_dim1]; 574 } 575 } 576 i__1 = n; 577 for (i__ = 1; i__ <= i__1; ++i__) { 578 ip = i__ + npt; 579 temp = w[i__]; 580 i__2 = i__; 581 for (j = 1; j <= i__2; ++j) { 582 /* L150: */ 583 bmat[ip + j * bmat_dim1] += temp * w[j]; 584 } 585 } 586 } 587 588 /* The following instructions complete the shift, including the changes */ 589 /* to the second derivative parameters of the quadratic model. */ 590 591 ih = 0; 592 i__2 = n; 593 for (j = 1; j <= i__2; ++j) { 594 w[j] = -half * sumpq * xopt[j]; 595 i__1 = npt; 596 for (k = 1; k <= i__1; ++k) { 597 w[j] += pq[k] * xpt[k + j * xpt_dim1]; 598 /* L160: */ 599 xpt[k + j * xpt_dim1] -= xopt[j]; 600 } 601 i__1 = j; 602 for (i__ = 1; i__ <= i__1; ++i__) { 603 ++ih; 604 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j]; 605 /* L170: */ 606 bmat[npt + i__ + j * bmat_dim1] = bmat[npt + j + i__ * 607 bmat_dim1]; 608 } 609 } 610 i__1 = n; 611 for (i__ = 1; i__ <= i__1; ++i__) { 612 xbase[i__] += xopt[i__]; 613 xnew[i__] -= xopt[i__]; 614 sl[i__] -= xopt[i__]; 615 su[i__] -= xopt[i__]; 616 /* L180: */ 617 xopt[i__] = zero; 618 } 619 xoptsq = zero; 620 } 621 if (ntrits == 0) { 622 goto L210; 623 } 624 goto L230; 625 626 /* XBASE is also moved to XOPT by a call of RESCUE. This calculation is */ 627 /* more expensive than the previous shift, because new matrices BMAT and */ 628 /* ZMAT are generated from scratch, which may include the replacement of */ 629 /* interpolation points whose positions seem to be causing near linear */ 630 /* dependence in the interpolation conditions. Therefore RESCUE is called */ 631 /* only if rounding errors have reduced by at least a factor of two the */ 632 /* denominator of the formula for updating the H matrix. It provides a */ 633 /* useful safeguard, but is not invoked in most applications of BOBYQA. */ 634 635 L190: 636 nfsav = nf; 637 kbase = kopt; 638 rescue_(calfun, n, npt, &xl[1], &xu[1], maxfun, &xbase[1], &xpt[ 639 xpt_offset], &fval[1], &xopt[1], &gopt[1], &hq[1], &pq[1], &bmat[ 640 bmat_offset], &zmat[zmat_offset], ndim, &sl[1], &su[1], nf, delta, 641 kopt, &vlag[1], &w[1], &w[n + np], &w[ndim + np]); 642 643 /* XOPT is updated now in case the branch below to label 720 is taken. */ 644 /* Any updating of GOPT occurs after the branch below to label 20, which */ 645 /* leads to a trust region iteration as does the branch to label 60. */ 646 647 xoptsq = zero; 648 if (kopt != kbase) { 649 i__1 = n; 650 for (i__ = 1; i__ <= i__1; ++i__) { 651 xopt[i__] = xpt[kopt + i__ * xpt_dim1]; 652 /* L200: */ 653 /* Computing 2nd power */ 654 d__1 = xopt[i__]; 655 xoptsq += d__1 * d__1; 656 } 657 } 658 if (nf < 0) { 659 nf = maxfun; 660 throw bobyqa_failure("Return from BOBYQA because the objective function has been called max_f_evals times."); 661 //goto L720; 662 } 663 nresc = nf; 664 if (nfsav < nf) { 665 nfsav = nf; 666 goto L20; 667 } 668 if (ntrits > 0) { 669 goto L60; 670 } 671 672 /* Pick two alternative vectors of variables, relative to XBASE, that */ 673 /* are suitable as new positions of the KNEW-th interpolation point. */ 674 /* Firstly, XNEW is set to the point on a line through XOPT and another */ 675 /* interpolation point that minimizes the predicted value of the next */ 676 /* denominator, subject to ||XNEW - XOPT|| .LEQ. ADELT and to the SL */ 677 /* and SU bounds. Secondly, XALT is set to the best feasible point on */ 678 /* a constrained version of the Cauchy step of the KNEW-th Lagrange */ 679 /* function, the corresponding value of the square of this function */ 680 /* being returned in CAUCHY. The choice between these alternatives is */ 681 /* going to be made when the denominator is calculated. */ 682 683 L210: 684 altmov_(n, npt, &xpt[xpt_offset], &xopt[1], &bmat[bmat_offset], &zmat[zmat_offset], 685 ndim, &sl[1], &su[1], kopt, knew, adelt, &xnew[1], 686 &xalt[1], alpha, cauchy, &w[1], &w[np], &w[ndim + 1]); 687 i__1 = n; 688 for (i__ = 1; i__ <= i__1; ++i__) { 689 /* L220: */ 690 d__[i__] = xnew[i__] - xopt[i__]; 691 } 692 693 /* Calculate VLAG and BETA for the current choice of D. The scalar */ 694 /* product of D with XPT(K,.) is going to be held in W(NPT+K) for */ 695 /* use when VQUAD is calculated. */ 696 697 L230: 698 i__1 = npt; 699 for (k = 1; k <= i__1; ++k) { 700 suma = zero; 701 sumb = zero; 702 sum = zero; 703 i__2 = n; 704 for (j = 1; j <= i__2; ++j) { 705 suma += xpt[k + j * xpt_dim1] * d__[j]; 706 sumb += xpt[k + j * xpt_dim1] * xopt[j]; 707 /* L240: */ 708 sum += bmat[k + j * bmat_dim1] * d__[j]; 709 } 710 w[k] = suma * (half * suma + sumb); 711 vlag[k] = sum; 712 /* L250: */ 713 w[npt + k] = suma; 714 } 715 beta = zero; 716 i__1 = nptm; 717 for (jj = 1; jj <= i__1; ++jj) { 718 sum = zero; 719 i__2 = npt; 720 for (k = 1; k <= i__2; ++k) { 721 /* L260: */ 722 sum += zmat[k + jj * zmat_dim1] * w[k]; 723 } 724 beta -= sum * sum; 725 i__2 = npt; 726 for (k = 1; k <= i__2; ++k) { 727 /* L270: */ 728 vlag[k] += sum * zmat[k + jj * zmat_dim1]; 729 } 730 } 731 dsq = zero; 732 bsum = zero; 733 dx = zero; 734 i__2 = n; 735 for (j = 1; j <= i__2; ++j) { 736 /* Computing 2nd power */ 737 d__1 = d__[j]; 738 dsq += d__1 * d__1; 739 sum = zero; 740 i__1 = npt; 741 for (k = 1; k <= i__1; ++k) { 742 /* L280: */ 743 sum += w[k] * bmat[k + j * bmat_dim1]; 744 } 745 bsum += sum * d__[j]; 746 jp = npt + j; 747 i__1 = n; 748 for (i__ = 1; i__ <= i__1; ++i__) { 749 /* L290: */ 750 sum += bmat[jp + i__ * bmat_dim1] * d__[i__]; 751 } 752 vlag[jp] = sum; 753 bsum += sum * d__[j]; 754 /* L300: */ 755 dx += d__[j] * xopt[j]; 756 } 757 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum; 758 vlag[kopt] += one; 759 760 /* If NTRITS is zero, the denominator may be increased by replacing */ 761 /* the step D of ALTMOV by a Cauchy step. Then RESCUE may be called if */ 762 /* rounding errors have damaged the chosen denominator. */ 763 764 if (ntrits == 0) { 765 /* Computing 2nd power */ 766 d__1 = vlag[knew]; 767 denom = d__1 * d__1 + alpha * beta; 768 if (denom < cauchy && cauchy > zero) { 769 i__2 = n; 770 for (i__ = 1; i__ <= i__2; ++i__) { 771 xnew[i__] = xalt[i__]; 772 /* L310: */ 773 d__[i__] = xnew[i__] - xopt[i__]; 774 } 775 cauchy = zero; 776 goto L230; 777 } 778 /* Computing 2nd power */ 779 d__1 = vlag[knew]; 780 if (denom <= half * (d__1 * d__1)) { 781 if (nf > nresc) { 782 goto L190; 783 } 784 throw bobyqa_failure("Return from BOBYQA because of much cancellation in a denominator."); 785 //goto L720; 786 } 787 788 /* Alternatively, if NTRITS is positive, then set KNEW to the index of */ 789 /* the next interpolation point to be deleted to make room for a trust */ 790 /* region step. Again RESCUE may be called if rounding errors have damaged */ 791 /* the chosen denominator, which is the reason for attempting to select */ 792 /* KNEW before calculating the next value of the objective function. */ 793 794 } else { 795 delsq = delta * delta; 796 scaden = zero; 797 biglsq = zero; 798 knew = 0; 799 i__2 = npt; 800 for (k = 1; k <= i__2; ++k) { 801 if (k == kopt) { 802 goto L350; 803 } 804 hdiag = zero; 805 i__1 = nptm; 806 for (jj = 1; jj <= i__1; ++jj) { 807 /* L330: */ 808 /* Computing 2nd power */ 809 d__1 = zmat[k + jj * zmat_dim1]; 810 hdiag += d__1 * d__1; 811 } 812 /* Computing 2nd power */ 813 d__1 = vlag[k]; 814 den = beta * hdiag + d__1 * d__1; 815 distsq = zero; 816 i__1 = n; 817 for (j = 1; j <= i__1; ++j) { 818 /* L340: */ 819 /* Computing 2nd power */ 820 d__1 = xpt[k + j * xpt_dim1] - xopt[j]; 821 distsq += d__1 * d__1; 822 } 823 /* Computing MAX */ 824 /* Computing 2nd power */ 825 d__3 = distsq / delsq; 826 d__1 = one, d__2 = d__3 * d__3; 827 temp = std::max(d__1,d__2); 828 if (temp * den > scaden) { 829 scaden = temp * den; 830 knew = k; 831 denom = den; 832 } 833 /* Computing MAX */ 834 /* Computing 2nd power */ 835 d__3 = vlag[k]; 836 d__1 = biglsq, d__2 = temp * (d__3 * d__3); 837 biglsq = std::max(d__1,d__2); 838 L350: 839 ; 840 } 841 if (scaden <= half * biglsq) { 842 if (nf > nresc) { 843 goto L190; 844 } 845 throw bobyqa_failure("Return from BOBYQA because of much cancellation in a denominator."); 846 //goto L720; 847 } 848 } 849 850 /* Put the variables for the next calculation of the objective function */ 851 /* in XNEW, with any adjustments for the bounds. */ 852 853 854 /* Calculate the value of the objective function at XBASE+XNEW, unless */ 855 /* the limit on the number of calculations of F has been reached. */ 856 857 L360: 858 i__2 = n; 859 for (i__ = 1; i__ <= i__2; ++i__) { 860 /* Computing MIN */ 861 /* Computing MAX */ 862 d__3 = xl[i__], d__4 = xbase[i__] + xnew[i__]; 863 d__1 = std::max(d__3,d__4), d__2 = xu[i__]; 864 x[i__] = std::min(d__1,d__2); 865 if (xnew[i__] == sl[i__]) { 866 x[i__] = xl[i__]; 867 } 868 if (xnew[i__] == su[i__]) { 869 x[i__] = xu[i__]; 870 } 871 /* L380: */ 872 } 873 if (nf >= maxfun) { 874 throw bobyqa_failure("Return from BOBYQA because the objective function has been called max_f_evals times."); 875 //goto L720; 876 } 877 ++nf; 878 f = calfun(mat(&x[1], n)); 879 if (ntrits == -1) { 880 fsave = f; 881 goto L720; 882 } 883 884 /* Use the quadratic model to predict the change in F due to the step D, */ 885 /* and set DIFF to the error of this prediction. */ 886 887 fopt = fval[kopt]; 888 vquad = zero; 889 ih = 0; 890 i__2 = n; 891 for (j = 1; j <= i__2; ++j) { 892 vquad += d__[j] * gopt[j]; 893 i__1 = j; 894 for (i__ = 1; i__ <= i__1; ++i__) { 895 ++ih; 896 temp = d__[i__] * d__[j]; 897 if (i__ == j) { 898 temp = half * temp; 899 } 900 /* L410: */ 901 vquad += hq[ih] * temp; 902 } 903 } 904 i__1 = npt; 905 for (k = 1; k <= i__1; ++k) { 906 /* L420: */ 907 /* Computing 2nd power */ 908 d__1 = w[npt + k]; 909 vquad += half * pq[k] * (d__1 * d__1); 910 } 911 diff = f - fopt - vquad; 912 diffc = diffb; 913 diffb = diffa; 914 diffa = std::abs(diff); 915 if (dnorm > rho) { 916 nfsav = nf; 917 } 918 919 /* Pick the next value of DELTA after a trust region step. */ 920 921 if (ntrits > 0) { 922 if (vquad >= zero) { 923 throw bobyqa_failure("Return from BOBYQA because a trust region step has failed to reduce Q."); 924 //goto L720; 925 } 926 ratio = (f - fopt) / vquad; 927 if (ratio <= tenth) { 928 /* Computing MIN */ 929 d__1 = half * delta; 930 delta = std::min(d__1,dnorm); 931 } else if (ratio <= .7) { 932 /* Computing MAX */ 933 d__1 = half * delta; 934 delta = std::max(d__1,dnorm); 935 } else { 936 /* Computing MAX */ 937 d__1 = half * delta, d__2 = dnorm + dnorm; 938 delta = std::max(d__1,d__2); 939 } 940 if (delta <= rho * 1.5) { 941 delta = rho; 942 } 943 944 /* Recalculate KNEW and DENOM if the new F is less than FOPT. */ 945 946 if (f < fopt) { 947 ksav = knew; 948 densav = denom; 949 delsq = delta * delta; 950 scaden = zero; 951 biglsq = zero; 952 knew = 0; 953 i__1 = npt; 954 for (k = 1; k <= i__1; ++k) { 955 hdiag = zero; 956 i__2 = nptm; 957 for (jj = 1; jj <= i__2; ++jj) { 958 /* L440: */ 959 /* Computing 2nd power */ 960 d__1 = zmat[k + jj * zmat_dim1]; 961 hdiag += d__1 * d__1; 962 } 963 /* Computing 2nd power */ 964 d__1 = vlag[k]; 965 den = beta * hdiag + d__1 * d__1; 966 distsq = zero; 967 i__2 = n; 968 for (j = 1; j <= i__2; ++j) { 969 /* L450: */ 970 /* Computing 2nd power */ 971 d__1 = xpt[k + j * xpt_dim1] - xnew[j]; 972 distsq += d__1 * d__1; 973 } 974 /* Computing MAX */ 975 /* Computing 2nd power */ 976 d__3 = distsq / delsq; 977 d__1 = one, d__2 = d__3 * d__3; 978 temp = std::max(d__1,d__2); 979 if (temp * den > scaden) { 980 scaden = temp * den; 981 knew = k; 982 denom = den; 983 } 984 /* L460: */ 985 /* Computing MAX */ 986 /* Computing 2nd power */ 987 d__3 = vlag[k]; 988 d__1 = biglsq, d__2 = temp * (d__3 * d__3); 989 biglsq = std::max(d__1,d__2); 990 } 991 if (scaden <= half * biglsq) { 992 knew = ksav; 993 denom = densav; 994 } 995 } 996 } 997 998 /* Update BMAT and ZMAT, so that the KNEW-th interpolation point can be */ 999 /* moved. Also update the second derivative terms of the model. */ 1000 1001 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], 1002 beta, denom, knew, &w[1]); 1003 ih = 0; 1004 pqold = pq[knew]; 1005 pq[knew] = zero; 1006 i__1 = n; 1007 for (i__ = 1; i__ <= i__1; ++i__) { 1008 temp = pqold * xpt[knew + i__ * xpt_dim1]; 1009 i__2 = i__; 1010 for (j = 1; j <= i__2; ++j) { 1011 ++ih; 1012 /* L470: */ 1013 hq[ih] += temp * xpt[knew + j * xpt_dim1]; 1014 } 1015 } 1016 i__2 = nptm; 1017 for (jj = 1; jj <= i__2; ++jj) { 1018 temp = diff * zmat[knew + jj * zmat_dim1]; 1019 i__1 = npt; 1020 for (k = 1; k <= i__1; ++k) { 1021 /* L480: */ 1022 pq[k] += temp * zmat[k + jj * zmat_dim1]; 1023 } 1024 } 1025 1026 /* Include the new interpolation point, and make the changes to GOPT at */ 1027 /* the old XOPT that are caused by the updating of the quadratic model. */ 1028 1029 fval[knew] = f; 1030 i__1 = n; 1031 for (i__ = 1; i__ <= i__1; ++i__) { 1032 xpt[knew + i__ * xpt_dim1] = xnew[i__]; 1033 /* L490: */ 1034 w[i__] = bmat[knew + i__ * bmat_dim1]; 1035 } 1036 i__1 = npt; 1037 for (k = 1; k <= i__1; ++k) { 1038 suma = zero; 1039 i__2 = nptm; 1040 for (jj = 1; jj <= i__2; ++jj) { 1041 /* L500: */ 1042 suma += zmat[knew + jj * zmat_dim1] * zmat[k + jj * zmat_dim1]; 1043 } 1044 sumb = zero; 1045 i__2 = n; 1046 for (j = 1; j <= i__2; ++j) { 1047 /* L510: */ 1048 sumb += xpt[k + j * xpt_dim1] * xopt[j]; 1049 } 1050 temp = suma * sumb; 1051 i__2 = n; 1052 for (i__ = 1; i__ <= i__2; ++i__) { 1053 /* L520: */ 1054 w[i__] += temp * xpt[k + i__ * xpt_dim1]; 1055 } 1056 } 1057 i__2 = n; 1058 for (i__ = 1; i__ <= i__2; ++i__) { 1059 /* L530: */ 1060 gopt[i__] += diff * w[i__]; 1061 } 1062 1063 /* Update XOPT, GOPT and KOPT if the new calculated F is less than FOPT. */ 1064 1065 if (f < fopt) { 1066 kopt = knew; 1067 xoptsq = zero; 1068 ih = 0; 1069 i__2 = n; 1070 for (j = 1; j <= i__2; ++j) { 1071 xopt[j] = xnew[j]; 1072 /* Computing 2nd power */ 1073 d__1 = xopt[j]; 1074 xoptsq += d__1 * d__1; 1075 i__1 = j; 1076 for (i__ = 1; i__ <= i__1; ++i__) { 1077 ++ih; 1078 if (i__ < j) { 1079 gopt[j] += hq[ih] * d__[i__]; 1080 } 1081 /* L540: */ 1082 gopt[i__] += hq[ih] * d__[j]; 1083 } 1084 } 1085 i__1 = npt; 1086 for (k = 1; k <= i__1; ++k) { 1087 temp = zero; 1088 i__2 = n; 1089 for (j = 1; j <= i__2; ++j) { 1090 /* L550: */ 1091 temp += xpt[k + j * xpt_dim1] * d__[j]; 1092 } 1093 temp = pq[k] * temp; 1094 i__2 = n; 1095 for (i__ = 1; i__ <= i__2; ++i__) { 1096 /* L560: */ 1097 gopt[i__] += temp * xpt[k + i__ * xpt_dim1]; 1098 } 1099 } 1100 } 1101 1102 /* Calculate the parameters of the least Frobenius norm interpolant to */ 1103 /* the current data, the gradient of this interpolant at XOPT being put */ 1104 /* into VLAG(NPT+I), I=1,2,...,N. */ 1105 1106 if (ntrits > 0) { 1107 i__2 = npt; 1108 for (k = 1; k <= i__2; ++k) { 1109 vlag[k] = fval[k] - fval[kopt]; 1110 /* L570: */ 1111 w[k] = zero; 1112 } 1113 i__2 = nptm; 1114 for (j = 1; j <= i__2; ++j) { 1115 sum = zero; 1116 i__1 = npt; 1117 for (k = 1; k <= i__1; ++k) { 1118 /* L580: */ 1119 sum += zmat[k + j * zmat_dim1] * vlag[k]; 1120 } 1121 i__1 = npt; 1122 for (k = 1; k <= i__1; ++k) { 1123 /* L590: */ 1124 w[k] += sum * zmat[k + j * zmat_dim1]; 1125 } 1126 } 1127 i__1 = npt; 1128 for (k = 1; k <= i__1; ++k) { 1129 sum = zero; 1130 i__2 = n; 1131 for (j = 1; j <= i__2; ++j) { 1132 /* L600: */ 1133 sum += xpt[k + j * xpt_dim1] * xopt[j]; 1134 } 1135 w[k + npt] = w[k]; 1136 /* L610: */ 1137 w[k] = sum * w[k]; 1138 } 1139 gqsq = zero; 1140 gisq = zero; 1141 i__1 = n; 1142 for (i__ = 1; i__ <= i__1; ++i__) { 1143 sum = zero; 1144 i__2 = npt; 1145 for (k = 1; k <= i__2; ++k) { 1146 /* L620: */ 1147 sum = sum + bmat[k + i__ * bmat_dim1] * vlag[k] + xpt[k + i__ 1148 * xpt_dim1] * w[k]; 1149 } 1150 if (xopt[i__] == sl[i__]) { 1151 /* Computing MIN */ 1152 d__2 = zero, d__3 = gopt[i__]; 1153 /* Computing 2nd power */ 1154 d__1 = std::min(d__2,d__3); 1155 gqsq += d__1 * d__1; 1156 /* Computing 2nd power */ 1157 d__1 = std::min(zero,sum); 1158 gisq += d__1 * d__1; 1159 } else if (xopt[i__] == su[i__]) { 1160 /* Computing MAX */ 1161 d__2 = zero, d__3 = gopt[i__]; 1162 /* Computing 2nd power */ 1163 d__1 = std::max(d__2,d__3); 1164 gqsq += d__1 * d__1; 1165 /* Computing 2nd power */ 1166 d__1 = std::max(zero,sum); 1167 gisq += d__1 * d__1; 1168 } else { 1169 /* Computing 2nd power */ 1170 d__1 = gopt[i__]; 1171 gqsq += d__1 * d__1; 1172 gisq += sum * sum; 1173 } 1174 /* L630: */ 1175 vlag[npt + i__] = sum; 1176 } 1177 1178 /* Test whether to replace the new quadratic model by the least Frobenius */ 1179 /* norm interpolant, making the replacement if the test is satisfied. */ 1180 1181 ++itest; 1182 if (gqsq < ten * gisq) { 1183 itest = 0; 1184 } 1185 if (itest >= 3) { 1186 i__1 = std::max(npt,nh); 1187 for (i__ = 1; i__ <= i__1; ++i__) { 1188 if (i__ <= n) { 1189 gopt[i__] = vlag[npt + i__]; 1190 } 1191 if (i__ <= npt) { 1192 pq[i__] = w[npt + i__]; 1193 } 1194 if (i__ <= nh) { 1195 hq[i__] = zero; 1196 } 1197 itest = 0; 1198 /* L640: */ 1199 } 1200 } 1201 } 1202 1203 /* If a trust region step has provided a sufficient decrease in F, then */ 1204 /* branch for another trust region calculation. The case NTRITS=0 occurs */ 1205 /* when the new interpolation point was reached by an alternative step. */ 1206 1207 if (ntrits == 0) { 1208 goto L60; 1209 } 1210 if (f <= fopt + tenth * vquad) { 1211 goto L60; 1212 } 1213 1214 /* Alternatively, find out if the interpolation points are close enough */ 1215 /* to the best point so far. */ 1216 1217 /* Computing MAX */ 1218 /* Computing 2nd power */ 1219 d__3 = two * delta; 1220 /* Computing 2nd power */ 1221 d__4 = ten * rho; 1222 d__1 = d__3 * d__3, d__2 = d__4 * d__4; 1223 distsq = std::max(d__1,d__2); 1224 L650: 1225 knew = 0; 1226 i__1 = npt; 1227 for (k = 1; k <= i__1; ++k) { 1228 sum = zero; 1229 i__2 = n; 1230 for (j = 1; j <= i__2; ++j) { 1231 /* L660: */ 1232 /* Computing 2nd power */ 1233 d__1 = xpt[k + j * xpt_dim1] - xopt[j]; 1234 sum += d__1 * d__1; 1235 } 1236 if (sum > distsq) { 1237 knew = k; 1238 distsq = sum; 1239 } 1240 /* L670: */ 1241 } 1242 1243 /* If KNEW is positive, then ALTMOV finds alternative new positions for */ 1244 /* the KNEW-th interpolation point within distance ADELT of XOPT. It is */ 1245 /* reached via label 90. Otherwise, there is a branch to label 60 for */ 1246 /* another trust region iteration, unless the calculations with the */ 1247 /* current RHO are complete. */ 1248 1249 if (knew > 0) { 1250 dist = std::sqrt(distsq); 1251 if (ntrits == -1) { 1252 /* Computing MIN */ 1253 d__1 = tenth * delta, d__2 = half * dist; 1254 delta = std::min(d__1,d__2); 1255 if (delta <= rho * 1.5) { 1256 delta = rho; 1257 } 1258 } 1259 ntrits = 0; 1260 /* Computing MAX */ 1261 /* Computing MIN */ 1262 d__2 = tenth * dist; 1263 d__1 = std::min(d__2,delta); 1264 adelt = std::max(d__1,rho); 1265 dsq = adelt * adelt; 1266 goto L90; 1267 } 1268 if (ntrits == -1) { 1269 goto L680; 1270 } 1271 if (ratio > zero) { 1272 goto L60; 1273 } 1274 if (std::max(delta,dnorm) > rho) { 1275 goto L60; 1276 } 1277 1278 /* The calculations with the current value of RHO are complete. Pick the */ 1279 /* next values of RHO and DELTA. */ 1280 1281 L680: 1282 if (rho > rhoend) { 1283 delta = half * rho; 1284 ratio = rho / rhoend; 1285 if (ratio <= 16.) { 1286 rho = rhoend; 1287 } else if (ratio <= 250.) { 1288 rho = std::sqrt(ratio) * rhoend; 1289 } else { 1290 rho = tenth * rho; 1291 } 1292 delta = std::max(delta,rho); 1293 ntrits = 0; 1294 nfsav = nf; 1295 goto L60; 1296 } 1297 1298 /* Return from the calculation, after another Newton-Raphson step, if */ 1299 /* it is too short to have been tried before. */ 1300 1301 if (ntrits == -1) { 1302 goto L360; 1303 } 1304 L720: 1305 if (fval[kopt] <= fsave) { 1306 i__1 = n; 1307 for (i__ = 1; i__ <= i__1; ++i__) { 1308 /* Computing MIN */ 1309 /* Computing MAX */ 1310 d__3 = xl[i__], d__4 = xbase[i__] + xopt[i__]; 1311 d__1 = std::max(d__3,d__4), d__2 = xu[i__]; 1312 x[i__] = std::min(d__1,d__2); 1313 if (xopt[i__] == sl[i__]) { 1314 x[i__] = xl[i__]; 1315 } 1316 if (xopt[i__] == su[i__]) { 1317 x[i__] = xu[i__]; 1318 } 1319 /* L730: */ 1320 } 1321 f = fval[kopt]; 1322 } 1323 1324 return f; 1325 } /* bobyqb_ */ 1326 1327 // ---------------------------------------------------------------------------------------- 1328 altmov_(const integer n,const integer npt,const doublereal * xpt,const doublereal * xopt,const doublereal * bmat,const doublereal * zmat,const integer ndim,const doublereal * sl,const doublereal * su,const integer kopt,const integer knew,const doublereal adelt,doublereal * xnew,doublereal * xalt,doublereal & alpha,doublereal & cauchy,doublereal * glag,doublereal * hcol,doublereal * w)1329 void altmov_( 1330 const integer n, 1331 const integer npt, 1332 const doublereal *xpt, 1333 const doublereal *xopt, 1334 const doublereal *bmat, 1335 const doublereal *zmat, 1336 const integer ndim, 1337 const doublereal *sl, 1338 const doublereal *su, 1339 const integer kopt, 1340 const integer knew, 1341 const doublereal adelt, 1342 doublereal *xnew, 1343 doublereal *xalt, 1344 doublereal& alpha, 1345 doublereal& cauchy, 1346 doublereal *glag, 1347 doublereal *hcol, 1348 doublereal *w 1349 ) const 1350 { 1351 /* System generated locals */ 1352 integer xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, 1353 zmat_offset, i__1, i__2; 1354 doublereal d__1, d__2, d__3, d__4; 1355 1356 1357 /* Local variables */ 1358 integer i__, j, k; 1359 doublereal ha, gw, one, diff, half; 1360 integer ilbd, isbd; 1361 doublereal slbd; 1362 integer iubd; 1363 doublereal vlag, subd, temp; 1364 integer ksav = 0; 1365 doublereal step = 0, zero = 0, curv = 0; 1366 integer iflag; 1367 doublereal scale = 0, csave = 0, tempa = 0, tempb = 0, tempd = 0, const__ = 0, sumin = 0, 1368 ggfree = 0; 1369 integer ibdsav = 0; 1370 doublereal dderiv = 0, bigstp = 0, predsq = 0, presav = 0, distsq = 0, stpsav = 0, wfixsq = 0, wsqsav = 0; 1371 1372 1373 /* The arguments N, NPT, XPT, XOPT, BMAT, ZMAT, NDIM, SL and SU all have */ 1374 /* the same meanings as the corresponding arguments of BOBYQB. */ 1375 /* KOPT is the index of the optimal interpolation point. */ 1376 /* KNEW is the index of the interpolation point that is going to be moved. */ 1377 /* ADELT is the current trust region bound. */ 1378 /* XNEW will be set to a suitable new position for the interpolation point */ 1379 /* XPT(KNEW,.). Specifically, it satisfies the SL, SU and trust region */ 1380 /* bounds and it should provide a large denominator in the next call of */ 1381 /* UPDATE. The step XNEW-XOPT from XOPT is restricted to moves along the */ 1382 /* straight lines through XOPT and another interpolation point. */ 1383 /* XALT also provides a large value of the modulus of the KNEW-th Lagrange */ 1384 /* function subject to the constraints that have been mentioned, its main */ 1385 /* difference from XNEW being that XALT-XOPT is a constrained version of */ 1386 /* the Cauchy step within the trust region. An exception is that XALT is */ 1387 /* not calculated if all components of GLAG (see below) are zero. */ 1388 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix. */ 1389 /* CAUCHY will be set to the square of the KNEW-th Lagrange function at */ 1390 /* the step XALT-XOPT from XOPT for the vector XALT that is returned, */ 1391 /* except that CAUCHY is set to zero if XALT is not calculated. */ 1392 /* GLAG is a working space vector of length N for the gradient of the */ 1393 /* KNEW-th Lagrange function at XOPT. */ 1394 /* HCOL is a working space vector of length NPT for the second derivative */ 1395 /* coefficients of the KNEW-th Lagrange function. */ 1396 /* W is a working space vector of length 2N that is going to hold the */ 1397 /* constrained Cauchy step from XOPT of the Lagrange function, followed */ 1398 /* by the downhill version of XALT when the uphill step is calculated. */ 1399 1400 /* Set the first NPT components of W to the leading elements of the */ 1401 /* KNEW-th column of the H matrix. */ 1402 1403 /* Parameter adjustments */ 1404 zmat_dim1 = npt; 1405 zmat_offset = 1 + zmat_dim1; 1406 zmat -= zmat_offset; 1407 xpt_dim1 = npt; 1408 xpt_offset = 1 + xpt_dim1; 1409 xpt -= xpt_offset; 1410 --xopt; 1411 bmat_dim1 = ndim; 1412 bmat_offset = 1 + bmat_dim1; 1413 bmat -= bmat_offset; 1414 --sl; 1415 --su; 1416 --xnew; 1417 --xalt; 1418 --glag; 1419 --hcol; 1420 --w; 1421 1422 /* Function Body */ 1423 half = .5; 1424 one = 1.; 1425 zero = 0.; 1426 const__ = one + std::sqrt(2.); 1427 i__1 = npt; 1428 for (k = 1; k <= i__1; ++k) { 1429 /* L10: */ 1430 hcol[k] = zero; 1431 } 1432 i__1 = npt - n - 1; 1433 for (j = 1; j <= i__1; ++j) { 1434 temp = zmat[knew + j * zmat_dim1]; 1435 i__2 = npt; 1436 for (k = 1; k <= i__2; ++k) { 1437 /* L20: */ 1438 hcol[k] += temp * zmat[k + j * zmat_dim1]; 1439 } 1440 } 1441 alpha = hcol[knew]; 1442 ha = half * alpha; 1443 1444 /* Calculate the gradient of the KNEW-th Lagrange function at XOPT. */ 1445 1446 i__2 = n; 1447 for (i__ = 1; i__ <= i__2; ++i__) { 1448 /* L30: */ 1449 glag[i__] = bmat[knew + i__ * bmat_dim1]; 1450 } 1451 i__2 = npt; 1452 for (k = 1; k <= i__2; ++k) { 1453 temp = zero; 1454 i__1 = n; 1455 for (j = 1; j <= i__1; ++j) { 1456 /* L40: */ 1457 temp += xpt[k + j * xpt_dim1] * xopt[j]; 1458 } 1459 temp = hcol[k] * temp; 1460 i__1 = n; 1461 for (i__ = 1; i__ <= i__1; ++i__) { 1462 /* L50: */ 1463 glag[i__] += temp * xpt[k + i__ * xpt_dim1]; 1464 } 1465 } 1466 1467 /* Search for a large denominator along the straight lines through XOPT */ 1468 /* and another interpolation point. SLBD and SUBD will be lower and upper */ 1469 /* bounds on the step along each of these lines in turn. PREDSQ will be */ 1470 /* set to the square of the predicted denominator for each line. PRESAV */ 1471 /* will be set to the largest admissible value of PREDSQ that occurs. */ 1472 1473 presav = zero; 1474 i__1 = npt; 1475 for (k = 1; k <= i__1; ++k) { 1476 if (k == kopt) { 1477 goto L80; 1478 } 1479 dderiv = zero; 1480 distsq = zero; 1481 i__2 = n; 1482 for (i__ = 1; i__ <= i__2; ++i__) { 1483 temp = xpt[k + i__ * xpt_dim1] - xopt[i__]; 1484 dderiv += glag[i__] * temp; 1485 /* L60: */ 1486 distsq += temp * temp; 1487 } 1488 subd = adelt / std::sqrt(distsq); 1489 slbd = -subd; 1490 ilbd = 0; 1491 iubd = 0; 1492 sumin = std::min(one,subd); 1493 1494 /* Revise SLBD and SUBD if necessary because of the bounds in SL and SU. */ 1495 1496 i__2 = n; 1497 for (i__ = 1; i__ <= i__2; ++i__) { 1498 temp = xpt[k + i__ * xpt_dim1] - xopt[i__]; 1499 if (temp > zero) { 1500 if (slbd * temp < sl[i__] - xopt[i__]) { 1501 slbd = (sl[i__] - xopt[i__]) / temp; 1502 ilbd = -i__; 1503 } 1504 if (subd * temp > su[i__] - xopt[i__]) { 1505 /* Computing MAX */ 1506 d__1 = sumin, d__2 = (su[i__] - xopt[i__]) / temp; 1507 subd = std::max(d__1,d__2); 1508 iubd = i__; 1509 } 1510 } else if (temp < zero) { 1511 if (slbd * temp > su[i__] - xopt[i__]) { 1512 slbd = (su[i__] - xopt[i__]) / temp; 1513 ilbd = i__; 1514 } 1515 if (subd * temp < sl[i__] - xopt[i__]) { 1516 /* Computing MAX */ 1517 d__1 = sumin, d__2 = (sl[i__] - xopt[i__]) / temp; 1518 subd = std::max(d__1,d__2); 1519 iubd = -i__; 1520 } 1521 } 1522 /* L70: */ 1523 } 1524 1525 /* Seek a large modulus of the KNEW-th Lagrange function when the index */ 1526 /* of the other interpolation point on the line through XOPT is KNEW. */ 1527 1528 if (k == knew) { 1529 diff = dderiv - one; 1530 step = slbd; 1531 vlag = slbd * (dderiv - slbd * diff); 1532 isbd = ilbd; 1533 temp = subd * (dderiv - subd * diff); 1534 if (std::abs(temp) > std::abs(vlag)) { 1535 step = subd; 1536 vlag = temp; 1537 isbd = iubd; 1538 } 1539 tempd = half * dderiv; 1540 tempa = tempd - diff * slbd; 1541 tempb = tempd - diff * subd; 1542 if (tempa * tempb < zero) { 1543 temp = tempd * tempd / diff; 1544 if (std::abs(temp) > std::abs(vlag)) { 1545 step = tempd / diff; 1546 vlag = temp; 1547 isbd = 0; 1548 } 1549 } 1550 1551 /* Search along each of the other lines through XOPT and another point. */ 1552 1553 } else { 1554 step = slbd; 1555 vlag = slbd * (one - slbd); 1556 isbd = ilbd; 1557 temp = subd * (one - subd); 1558 if (std::abs(temp) > std::abs(vlag)) { 1559 step = subd; 1560 vlag = temp; 1561 isbd = iubd; 1562 } 1563 if (subd > half) { 1564 if (std::abs(vlag) < .25) { 1565 step = half; 1566 vlag = .25; 1567 isbd = 0; 1568 } 1569 } 1570 vlag *= dderiv; 1571 } 1572 1573 /* Calculate PREDSQ for the current line search and maintain PRESAV. */ 1574 1575 temp = step * (one - step) * distsq; 1576 predsq = vlag * vlag * (vlag * vlag + ha * temp * temp); 1577 if (predsq > presav) { 1578 presav = predsq; 1579 ksav = k; 1580 stpsav = step; 1581 ibdsav = isbd; 1582 } 1583 L80: 1584 ; 1585 } 1586 1587 /* Construct XNEW in a way that satisfies the bound constraints exactly. */ 1588 1589 i__1 = n; 1590 for (i__ = 1; i__ <= i__1; ++i__) { 1591 temp = xopt[i__] + stpsav * (xpt[ksav + i__ * xpt_dim1] - xopt[i__]); 1592 /* L90: */ 1593 /* Computing MAX */ 1594 /* Computing MIN */ 1595 d__3 = su[i__]; 1596 d__1 = sl[i__], d__2 = std::min(d__3,temp); 1597 xnew[i__] = std::max(d__1,d__2); 1598 } 1599 if (ibdsav < 0) { 1600 xnew[-ibdsav] = sl[-ibdsav]; 1601 } 1602 if (ibdsav > 0) { 1603 xnew[ibdsav] = su[ibdsav]; 1604 } 1605 1606 /* Prepare for the iterative method that assembles the constrained Cauchy */ 1607 /* step in W. The sum of squares of the fixed components of W is formed in */ 1608 /* WFIXSQ, and the free components of W are set to BIGSTP. */ 1609 1610 bigstp = adelt + adelt; 1611 iflag = 0; 1612 L100: 1613 wfixsq = zero; 1614 ggfree = zero; 1615 i__1 = n; 1616 for (i__ = 1; i__ <= i__1; ++i__) { 1617 w[i__] = zero; 1618 /* Computing MIN */ 1619 d__1 = xopt[i__] - sl[i__], d__2 = glag[i__]; 1620 tempa = std::min(d__1,d__2); 1621 /* Computing MAX */ 1622 d__1 = xopt[i__] - su[i__], d__2 = glag[i__]; 1623 tempb = std::max(d__1,d__2); 1624 if (tempa > zero || tempb < zero) { 1625 w[i__] = bigstp; 1626 /* Computing 2nd power */ 1627 d__1 = glag[i__]; 1628 ggfree += d__1 * d__1; 1629 } 1630 /* L110: */ 1631 } 1632 if (ggfree == zero) { 1633 cauchy = zero; 1634 goto L200; 1635 } 1636 1637 /* Investigate whether more components of W can be fixed. */ 1638 1639 L120: 1640 temp = adelt * adelt - wfixsq; 1641 if (temp > zero) { 1642 wsqsav = wfixsq; 1643 step = std::sqrt(temp / ggfree); 1644 ggfree = zero; 1645 i__1 = n; 1646 for (i__ = 1; i__ <= i__1; ++i__) { 1647 if (w[i__] == bigstp) { 1648 temp = xopt[i__] - step * glag[i__]; 1649 if (temp <= sl[i__]) { 1650 w[i__] = sl[i__] - xopt[i__]; 1651 /* Computing 2nd power */ 1652 d__1 = w[i__]; 1653 wfixsq += d__1 * d__1; 1654 } else if (temp >= su[i__]) { 1655 w[i__] = su[i__] - xopt[i__]; 1656 /* Computing 2nd power */ 1657 d__1 = w[i__]; 1658 wfixsq += d__1 * d__1; 1659 } else { 1660 /* Computing 2nd power */ 1661 d__1 = glag[i__]; 1662 ggfree += d__1 * d__1; 1663 } 1664 } 1665 /* L130: */ 1666 } 1667 if (wfixsq > wsqsav && ggfree > zero) { 1668 goto L120; 1669 } 1670 } 1671 1672 /* Set the remaining free components of W and all components of XALT, */ 1673 /* except that W may be scaled later. */ 1674 1675 gw = zero; 1676 i__1 = n; 1677 for (i__ = 1; i__ <= i__1; ++i__) { 1678 if (w[i__] == bigstp) { 1679 w[i__] = -step * glag[i__]; 1680 /* Computing MAX */ 1681 /* Computing MIN */ 1682 d__3 = su[i__], d__4 = xopt[i__] + w[i__]; 1683 d__1 = sl[i__], d__2 = std::min(d__3,d__4); 1684 xalt[i__] = std::max(d__1,d__2); 1685 } else if (w[i__] == zero) { 1686 xalt[i__] = xopt[i__]; 1687 } else if (glag[i__] > zero) { 1688 xalt[i__] = sl[i__]; 1689 } else { 1690 xalt[i__] = su[i__]; 1691 } 1692 /* L140: */ 1693 gw += glag[i__] * w[i__]; 1694 } 1695 1696 /* Set CURV to the curvature of the KNEW-th Lagrange function along W. */ 1697 /* Scale W by a factor less than one if that can reduce the modulus of */ 1698 /* the Lagrange function at XOPT+W. Set CAUCHY to the final value of */ 1699 /* the square of this function. */ 1700 1701 curv = zero; 1702 i__1 = npt; 1703 for (k = 1; k <= i__1; ++k) { 1704 temp = zero; 1705 i__2 = n; 1706 for (j = 1; j <= i__2; ++j) { 1707 /* L150: */ 1708 temp += xpt[k + j * xpt_dim1] * w[j]; 1709 } 1710 /* L160: */ 1711 curv += hcol[k] * temp * temp; 1712 } 1713 if (iflag == 1) { 1714 curv = -curv; 1715 } 1716 if (curv > -gw && curv < -const__ * gw) { 1717 scale = -gw / curv; 1718 i__1 = n; 1719 for (i__ = 1; i__ <= i__1; ++i__) { 1720 temp = xopt[i__] + scale * w[i__]; 1721 /* L170: */ 1722 /* Computing MAX */ 1723 /* Computing MIN */ 1724 d__3 = su[i__]; 1725 d__1 = sl[i__], d__2 = std::min(d__3,temp); 1726 xalt[i__] = std::max(d__1,d__2); 1727 } 1728 /* Computing 2nd power */ 1729 d__1 = half * gw * scale; 1730 cauchy = d__1 * d__1; 1731 } else { 1732 /* Computing 2nd power */ 1733 d__1 = gw + half * curv; 1734 cauchy = d__1 * d__1; 1735 } 1736 1737 /* If IFLAG is zero, then XALT is calculated as before after reversing */ 1738 /* the sign of GLAG. Thus two XALT vectors become available. The one that */ 1739 /* is chosen is the one that gives the larger value of CAUCHY. */ 1740 1741 if (iflag == 0) { 1742 i__1 = n; 1743 for (i__ = 1; i__ <= i__1; ++i__) { 1744 glag[i__] = -glag[i__]; 1745 /* L180: */ 1746 w[n + i__] = xalt[i__]; 1747 } 1748 csave = cauchy; 1749 iflag = 1; 1750 goto L100; 1751 } 1752 if (csave > cauchy) { 1753 i__1 = n; 1754 for (i__ = 1; i__ <= i__1; ++i__) { 1755 /* L190: */ 1756 xalt[i__] = w[n + i__]; 1757 } 1758 cauchy = csave; 1759 } 1760 L200: 1761 ; 1762 } /* altmov_ */ 1763 1764 // ---------------------------------------------------------------------------------------- 1765 1766 template <typename funct> prelim_(const funct & calfun,const integer n,const integer npt,doublereal * x,const doublereal * xl,const doublereal * xu,const doublereal rhobeg,const integer maxfun,doublereal * xbase,doublereal * xpt,doublereal * fval,doublereal * gopt,doublereal * hq,doublereal * pq,doublereal * bmat,doublereal * zmat,const integer ndim,const doublereal * sl,const doublereal * su,integer & nf,integer & kopt)1767 void prelim_( 1768 const funct& calfun, 1769 const integer n, 1770 const integer npt, 1771 doublereal *x, 1772 const doublereal *xl, 1773 const doublereal *xu, 1774 const doublereal rhobeg, 1775 const integer maxfun, 1776 doublereal *xbase, 1777 doublereal *xpt, 1778 doublereal *fval, 1779 doublereal *gopt, 1780 doublereal *hq, 1781 doublereal *pq, 1782 doublereal *bmat, 1783 doublereal *zmat, 1784 const integer ndim, 1785 const doublereal *sl, 1786 const doublereal *su, 1787 integer& nf, 1788 integer& kopt 1789 ) const 1790 { 1791 /* System generated locals */ 1792 integer xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, 1793 zmat_offset, i__1, i__2; 1794 doublereal d__1, d__2, d__3, d__4; 1795 1796 1797 /* Local variables */ 1798 doublereal f; 1799 integer i__, j, k, ih, np, nfm; 1800 doublereal one; 1801 integer nfx = 0, ipt = 0, jpt = 0; 1802 doublereal two = 0, fbeg = 0, diff = 0, half = 0, temp = 0, zero = 0, recip = 0, stepa = 0, stepb = 0; 1803 integer itemp; 1804 doublereal rhosq; 1805 1806 1807 1808 /* The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the */ 1809 /* same as the corresponding arguments in SUBROUTINE BOBYQA. */ 1810 /* The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU */ 1811 /* are the same as the corresponding arguments in BOBYQB, the elements */ 1812 /* of SL and SU being set in BOBYQA. */ 1813 /* GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but */ 1814 /* it is set by PRELIM to the gradient of the quadratic model at XBASE. */ 1815 /* If XOPT is nonzero, BOBYQB will change it to its usual value later. */ 1816 /* NF is maintaned as the number of calls of CALFUN so far. */ 1817 /* KOPT will be such that the least calculated value of F so far is at */ 1818 /* the point XPT(KOPT,.)+XBASE in the space of the variables. */ 1819 1820 /* SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, */ 1821 /* BMAT and ZMAT for the first iteration, and it maintains the values of */ 1822 /* NF and KOPT. The vector X is also changed by PRELIM. */ 1823 1824 /* Set some constants. */ 1825 1826 /* Parameter adjustments */ 1827 zmat_dim1 = npt; 1828 zmat_offset = 1 + zmat_dim1; 1829 zmat -= zmat_offset; 1830 xpt_dim1 = npt; 1831 xpt_offset = 1 + xpt_dim1; 1832 xpt -= xpt_offset; 1833 --x; 1834 --xl; 1835 --xu; 1836 --xbase; 1837 --fval; 1838 --gopt; 1839 --hq; 1840 --pq; 1841 bmat_dim1 = ndim; 1842 bmat_offset = 1 + bmat_dim1; 1843 bmat -= bmat_offset; 1844 --sl; 1845 --su; 1846 1847 /* Function Body */ 1848 half = .5; 1849 one = 1.; 1850 two = 2.; 1851 zero = 0.; 1852 rhosq = rhobeg * rhobeg; 1853 recip = one / rhosq; 1854 np = n + 1; 1855 1856 /* Set XBASE to the initial vector of variables, and set the initial */ 1857 /* elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */ 1858 1859 i__1 = n; 1860 for (j = 1; j <= i__1; ++j) { 1861 xbase[j] = x[j]; 1862 i__2 = npt; 1863 for (k = 1; k <= i__2; ++k) { 1864 /* L10: */ 1865 xpt[k + j * xpt_dim1] = zero; 1866 } 1867 i__2 = ndim; 1868 for (i__ = 1; i__ <= i__2; ++i__) { 1869 /* L20: */ 1870 bmat[i__ + j * bmat_dim1] = zero; 1871 } 1872 } 1873 i__2 = n * np / 2; 1874 for (ih = 1; ih <= i__2; ++ih) { 1875 /* L30: */ 1876 hq[ih] = zero; 1877 } 1878 i__2 = npt; 1879 for (k = 1; k <= i__2; ++k) { 1880 pq[k] = zero; 1881 i__1 = npt - np; 1882 for (j = 1; j <= i__1; ++j) { 1883 /* L40: */ 1884 zmat[k + j * zmat_dim1] = zero; 1885 } 1886 } 1887 1888 /* Begin the initialization procedure. NF becomes one more than the number */ 1889 /* of function values so far. The coordinates of the displacement of the */ 1890 /* next initial interpolation point from XBASE are set in XPT(NF+1,.). */ 1891 1892 nf = 0; 1893 L50: 1894 nfm = nf; 1895 nfx = nf - n; 1896 ++(nf); 1897 if (nfm <= n << 1) { 1898 if (nfm >= 1 && nfm <= n) { 1899 stepa = rhobeg; 1900 if (su[nfm] == zero) { 1901 stepa = -stepa; 1902 } 1903 xpt[nf + nfm * xpt_dim1] = stepa; 1904 } else if (nfm > n) { 1905 stepa = xpt[nf - n + nfx * xpt_dim1]; 1906 stepb = -(rhobeg); 1907 if (sl[nfx] == zero) { 1908 /* Computing MIN */ 1909 d__1 = two * rhobeg, d__2 = su[nfx]; 1910 stepb = std::min(d__1,d__2); 1911 } 1912 if (su[nfx] == zero) { 1913 /* Computing MAX */ 1914 d__1 = -two * rhobeg, d__2 = sl[nfx]; 1915 stepb = std::max(d__1,d__2); 1916 } 1917 xpt[nf + nfx * xpt_dim1] = stepb; 1918 } 1919 } else { 1920 itemp = (nfm - np) / n; 1921 jpt = nfm - itemp * n - n; 1922 ipt = jpt + itemp; 1923 if (ipt > n) { 1924 itemp = jpt; 1925 jpt = ipt - n; 1926 ipt = itemp; 1927 } 1928 xpt[nf + ipt * xpt_dim1] = xpt[ipt + 1 + ipt * xpt_dim1]; 1929 xpt[nf + jpt * xpt_dim1] = xpt[jpt + 1 + jpt * xpt_dim1]; 1930 } 1931 1932 /* Calculate the next value of F. The least function value so far and */ 1933 /* its index are required. */ 1934 1935 i__1 = n; 1936 for (j = 1; j <= i__1; ++j) { 1937 /* Computing MIN */ 1938 /* Computing MAX */ 1939 d__3 = xl[j], d__4 = xbase[j] + xpt[nf + j * xpt_dim1]; 1940 d__1 = std::max(d__3,d__4), d__2 = xu[j]; 1941 x[j] = std::min(d__1,d__2); 1942 if (xpt[nf + j * xpt_dim1] == sl[j]) { 1943 x[j] = xl[j]; 1944 } 1945 if (xpt[nf + j * xpt_dim1] == su[j]) { 1946 x[j] = xu[j]; 1947 } 1948 /* L60: */ 1949 } 1950 f = calfun(mat(&x[1],n)); 1951 fval[nf] = f; 1952 if (nf == 1) { 1953 fbeg = f; 1954 kopt = 1; 1955 } else if (f < fval[kopt]) { 1956 kopt = nf; 1957 } 1958 1959 /* Set the nonzero initial elements of BMAT and the quadratic model in the */ 1960 /* cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions */ 1961 /* of the NF-th and (NF-N)-th interpolation points may be switched, in */ 1962 /* order that the function value at the first of them contributes to the */ 1963 /* off-diagonal second derivative terms of the initial quadratic model. */ 1964 1965 if (nf <= (n << 1) + 1) { 1966 if (nf >= 2 && nf <= n + 1) { 1967 gopt[nfm] = (f - fbeg) / stepa; 1968 if (npt < nf + n) { 1969 bmat[nfm * bmat_dim1 + 1] = -one / stepa; 1970 bmat[nf + nfm * bmat_dim1] = one / stepa; 1971 bmat[npt + nfm + nfm * bmat_dim1] = -half * rhosq; 1972 } 1973 } else if (nf >= n + 2) { 1974 ih = nfx * (nfx + 1) / 2; 1975 temp = (f - fbeg) / stepb; 1976 diff = stepb - stepa; 1977 hq[ih] = two * (temp - gopt[nfx]) / diff; 1978 gopt[nfx] = (gopt[nfx] * stepb - temp * stepa) / diff; 1979 if (stepa * stepb < zero) { 1980 if (f < fval[nf - n]) { 1981 fval[nf] = fval[nf - n]; 1982 fval[nf - n] = f; 1983 if (kopt == nf) { 1984 kopt = nf - n; 1985 } 1986 xpt[nf - n + nfx * xpt_dim1] = stepb; 1987 xpt[nf + nfx * xpt_dim1] = stepa; 1988 } 1989 } 1990 bmat[nfx * bmat_dim1 + 1] = -(stepa + stepb) / (stepa * stepb); 1991 bmat[nf + nfx * bmat_dim1] = -half / xpt[nf - n + nfx * 1992 xpt_dim1]; 1993 bmat[nf - n + nfx * bmat_dim1] = -bmat[nfx * bmat_dim1 + 1] - 1994 bmat[nf + nfx * bmat_dim1]; 1995 zmat[nfx * zmat_dim1 + 1] = std::sqrt(two) / (stepa * stepb); 1996 zmat[nf + nfx * zmat_dim1] = std::sqrt(half) / rhosq; 1997 zmat[nf - n + nfx * zmat_dim1] = -zmat[nfx * zmat_dim1 + 1] - 1998 zmat[nf + nfx * zmat_dim1]; 1999 } 2000 2001 /* Set the off-diagonal second derivatives of the Lagrange functions and */ 2002 /* the initial quadratic model. */ 2003 2004 } else { 2005 ih = ipt * (ipt - 1) / 2 + jpt; 2006 zmat[nfx * zmat_dim1 + 1] = recip; 2007 zmat[nf + nfx * zmat_dim1] = recip; 2008 zmat[ipt + 1 + nfx * zmat_dim1] = -recip; 2009 zmat[jpt + 1 + nfx * zmat_dim1] = -recip; 2010 temp = xpt[nf + ipt * xpt_dim1] * xpt[nf + jpt * xpt_dim1]; 2011 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / temp; 2012 } 2013 if (nf < npt && nf < maxfun) { 2014 goto L50; 2015 } 2016 2017 } /* prelim_ */ 2018 2019 // ---------------------------------------------------------------------------------------- 2020 2021 template <typename funct> rescue_(const funct & calfun,const integer n,const integer npt,const doublereal * xl,const doublereal * xu,const integer maxfun,doublereal * xbase,doublereal * xpt,doublereal * fval,doublereal * xopt,doublereal * gopt,doublereal * hq,doublereal * pq,doublereal * bmat,doublereal * zmat,const integer ndim,doublereal * sl,doublereal * su,integer & nf,const doublereal delta,integer & kopt,doublereal * vlag,doublereal * ptsaux,doublereal * ptsid,doublereal * w)2022 void rescue_ ( 2023 const funct& calfun, 2024 const integer n, 2025 const integer npt, 2026 const doublereal *xl, 2027 const doublereal *xu, 2028 const integer maxfun, 2029 doublereal *xbase, 2030 doublereal *xpt, 2031 doublereal *fval, 2032 doublereal *xopt, 2033 doublereal *gopt, 2034 doublereal *hq, 2035 doublereal *pq, 2036 doublereal *bmat, 2037 doublereal *zmat, 2038 const integer ndim, 2039 doublereal *sl, 2040 doublereal *su, 2041 integer& nf, 2042 const doublereal delta, 2043 integer& kopt, 2044 doublereal *vlag, 2045 doublereal * ptsaux, 2046 doublereal *ptsid, 2047 doublereal *w 2048 ) const 2049 { 2050 /* System generated locals */ 2051 integer xpt_dim1, xpt_offset, bmat_dim1, bmat_offset, zmat_dim1, 2052 zmat_offset, i__1, i__2, i__3; 2053 doublereal d__1, d__2, d__3, d__4; 2054 2055 2056 /* Local variables */ 2057 doublereal f; 2058 integer i__, j, k, ih, jp, ip, iq, np, iw; 2059 doublereal xp = 0, xq = 0, den = 0; 2060 integer ihp = 0; 2061 doublereal one; 2062 integer ihq, jpn, kpt; 2063 doublereal sum = 0, diff = 0, half = 0, beta = 0; 2064 integer kold; 2065 doublereal winc; 2066 integer nrem, knew; 2067 doublereal temp, bsum; 2068 integer nptm; 2069 doublereal zero = 0, hdiag = 0, fbase = 0, sfrac = 0, denom = 0, vquad = 0, sumpq = 0; 2070 doublereal dsqmin, distsq, vlmxsq; 2071 2072 2073 2074 /* The arguments N, NPT, XL, XU, IPRINT, MAXFUN, XBASE, XPT, FVAL, XOPT, */ 2075 /* GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU have the same meanings as */ 2076 /* the corresponding arguments of BOBYQB on the entry to RESCUE. */ 2077 /* NF is maintained as the number of calls of CALFUN so far, except that */ 2078 /* NF is set to -1 if the value of MAXFUN prevents further progress. */ 2079 /* KOPT is maintained so that FVAL(KOPT) is the least calculated function */ 2080 /* value. Its correct value must be given on entry. It is updated if a */ 2081 /* new least function value is found, but the corresponding changes to */ 2082 /* XOPT and GOPT have to be made later by the calling program. */ 2083 /* DELTA is the current trust region radius. */ 2084 /* VLAG is a working space vector that will be used for the values of the */ 2085 /* provisional Lagrange functions at each of the interpolation points. */ 2086 /* They are part of a product that requires VLAG to be of length NDIM. */ 2087 /* PTSAUX is also a working space array. For J=1,2,...,N, PTSAUX(1,J) and */ 2088 /* PTSAUX(2,J) specify the two positions of provisional interpolation */ 2089 /* points when a nonzero step is taken along e_J (the J-th coordinate */ 2090 /* direction) through XBASE+XOPT, as specified below. Usually these */ 2091 /* steps have length DELTA, but other lengths are chosen if necessary */ 2092 /* in order to satisfy the given bounds on the variables. */ 2093 /* PTSID is also a working space array. It has NPT components that denote */ 2094 /* provisional new positions of the original interpolation points, in */ 2095 /* case changes are needed to restore the linear independence of the */ 2096 /* interpolation conditions. The K-th point is a candidate for change */ 2097 /* if and only if PTSID(K) is nonzero. In this case let p and q be the */ 2098 /* integer parts of PTSID(K) and (PTSID(K)-p) multiplied by N+1. If p */ 2099 /* and q are both positive, the step from XBASE+XOPT to the new K-th */ 2100 /* interpolation point is PTSAUX(1,p)*e_p + PTSAUX(1,q)*e_q. Otherwise */ 2101 /* the step is PTSAUX(1,p)*e_p or PTSAUX(2,q)*e_q in the cases q=0 or */ 2102 /* p=0, respectively. */ 2103 /* The first NDIM+NPT elements of the array W are used for working space. */ 2104 /* The final elements of BMAT and ZMAT are set in a well-conditioned way */ 2105 /* to the values that are appropriate for the new interpolation points. */ 2106 /* The elements of GOPT, HQ and PQ are also revised to the values that are */ 2107 /* appropriate to the final quadratic model. */ 2108 2109 /* Set some constants. */ 2110 2111 /* Parameter adjustments */ 2112 zmat_dim1 = npt; 2113 zmat_offset = 1 + zmat_dim1; 2114 zmat -= zmat_offset; 2115 xpt_dim1 = npt; 2116 xpt_offset = 1 + xpt_dim1; 2117 xpt -= xpt_offset; 2118 --xl; 2119 --xu; 2120 --xbase; 2121 --fval; 2122 --xopt; 2123 --gopt; 2124 --hq; 2125 --pq; 2126 bmat_dim1 = ndim; 2127 bmat_offset = 1 + bmat_dim1; 2128 bmat -= bmat_offset; 2129 --sl; 2130 --su; 2131 --vlag; 2132 ptsaux -= 3; 2133 --ptsid; 2134 --w; 2135 2136 /* Function Body */ 2137 half = .5; 2138 one = 1.; 2139 zero = 0.; 2140 np = n + 1; 2141 sfrac = half / (doublereal) np; 2142 nptm = npt - np; 2143 2144 /* Shift the interpolation points so that XOPT becomes the origin, and set */ 2145 /* the elements of ZMAT to zero. The value of SUMPQ is required in the */ 2146 /* updating of HQ below. The squares of the distances from XOPT to the */ 2147 /* other interpolation points are set at the end of W. Increments of WINC */ 2148 /* may be added later to these squares to balance the consideration of */ 2149 /* the choice of point that is going to become current. */ 2150 2151 sumpq = zero; 2152 winc = zero; 2153 i__1 = npt; 2154 for (k = 1; k <= i__1; ++k) { 2155 distsq = zero; 2156 i__2 = n; 2157 for (j = 1; j <= i__2; ++j) { 2158 xpt[k + j * xpt_dim1] -= xopt[j]; 2159 /* L10: */ 2160 /* Computing 2nd power */ 2161 d__1 = xpt[k + j * xpt_dim1]; 2162 distsq += d__1 * d__1; 2163 } 2164 sumpq += pq[k]; 2165 w[ndim + k] = distsq; 2166 winc = std::max(winc,distsq); 2167 i__2 = nptm; 2168 for (j = 1; j <= i__2; ++j) { 2169 /* L20: */ 2170 zmat[k + j * zmat_dim1] = zero; 2171 } 2172 } 2173 2174 /* Update HQ so that HQ and PQ define the second derivatives of the model */ 2175 /* after XBASE has been shifted to the trust region centre. */ 2176 2177 ih = 0; 2178 i__2 = n; 2179 for (j = 1; j <= i__2; ++j) { 2180 w[j] = half * sumpq * xopt[j]; 2181 i__1 = npt; 2182 for (k = 1; k <= i__1; ++k) { 2183 /* L30: */ 2184 w[j] += pq[k] * xpt[k + j * xpt_dim1]; 2185 } 2186 i__1 = j; 2187 for (i__ = 1; i__ <= i__1; ++i__) { 2188 ++ih; 2189 /* L40: */ 2190 hq[ih] = hq[ih] + w[i__] * xopt[j] + w[j] * xopt[i__]; 2191 } 2192 } 2193 2194 /* Shift XBASE, SL, SU and XOPT. Set the elements of BMAT to zero, and */ 2195 /* also set the elements of PTSAUX. */ 2196 2197 i__1 = n; 2198 for (j = 1; j <= i__1; ++j) { 2199 xbase[j] += xopt[j]; 2200 sl[j] -= xopt[j]; 2201 su[j] -= xopt[j]; 2202 xopt[j] = zero; 2203 /* Computing MIN */ 2204 d__1 = delta, d__2 = su[j]; 2205 ptsaux[(j << 1) + 1] = std::min(d__1,d__2); 2206 /* Computing MAX */ 2207 d__1 = -(delta), d__2 = sl[j]; 2208 ptsaux[(j << 1) + 2] = std::max(d__1,d__2); 2209 if (ptsaux[(j << 1) + 1] + ptsaux[(j << 1) + 2] < zero) { 2210 temp = ptsaux[(j << 1) + 1]; 2211 ptsaux[(j << 1) + 1] = ptsaux[(j << 1) + 2]; 2212 ptsaux[(j << 1) + 2] = temp; 2213 } 2214 if ((d__2 = ptsaux[(j << 1) + 2], std::abs(d__2)) < half * (d__1 = ptsaux[( 2215 j << 1) + 1], std::abs(d__1))) { 2216 ptsaux[(j << 1) + 2] = half * ptsaux[(j << 1) + 1]; 2217 } 2218 i__2 = ndim; 2219 for (i__ = 1; i__ <= i__2; ++i__) { 2220 /* L50: */ 2221 bmat[i__ + j * bmat_dim1] = zero; 2222 } 2223 } 2224 fbase = fval[kopt]; 2225 2226 /* Set the identifiers of the artificial interpolation points that are */ 2227 /* along a coordinate direction from XOPT, and set the corresponding */ 2228 /* nonzero elements of BMAT and ZMAT. */ 2229 2230 ptsid[1] = sfrac; 2231 i__2 = n; 2232 for (j = 1; j <= i__2; ++j) { 2233 jp = j + 1; 2234 jpn = jp + n; 2235 ptsid[jp] = (doublereal) j + sfrac; 2236 if (jpn <= npt) { 2237 ptsid[jpn] = (doublereal) j / (doublereal) np + sfrac; 2238 temp = one / (ptsaux[(j << 1) + 1] - ptsaux[(j << 1) + 2]); 2239 bmat[jp + j * bmat_dim1] = -temp + one / ptsaux[(j << 1) + 1]; 2240 bmat[jpn + j * bmat_dim1] = temp + one / ptsaux[(j << 1) + 2]; 2241 bmat[j * bmat_dim1 + 1] = -bmat[jp + j * bmat_dim1] - bmat[jpn + 2242 j * bmat_dim1]; 2243 zmat[j * zmat_dim1 + 1] = std::sqrt(2.) / (d__1 = ptsaux[(j << 1) + 1] 2244 * ptsaux[(j << 1) + 2], std::abs(d__1)); 2245 zmat[jp + j * zmat_dim1] = zmat[j * zmat_dim1 + 1] * ptsaux[(j << 2246 1) + 2] * temp; 2247 zmat[jpn + j * zmat_dim1] = -zmat[j * zmat_dim1 + 1] * ptsaux[(j 2248 << 1) + 1] * temp; 2249 } else { 2250 bmat[j * bmat_dim1 + 1] = -one / ptsaux[(j << 1) + 1]; 2251 bmat[jp + j * bmat_dim1] = one / ptsaux[(j << 1) + 1]; 2252 /* Computing 2nd power */ 2253 d__1 = ptsaux[(j << 1) + 1]; 2254 bmat[j + npt + j * bmat_dim1] = -half * (d__1 * d__1); 2255 } 2256 /* L60: */ 2257 } 2258 2259 /* Set any remaining identifiers with their nonzero elements of ZMAT. */ 2260 2261 if (npt >= n + np) { 2262 i__2 = npt; 2263 for (k = np << 1; k <= i__2; ++k) { 2264 iw = (integer) (((doublereal) (k - np) - half) / (doublereal) (n) 2265 ); 2266 ip = k - np - iw * n; 2267 iq = ip + iw; 2268 if (iq > n) { 2269 iq -= n; 2270 } 2271 ptsid[k] = (doublereal) ip + (doublereal) iq / (doublereal) np + 2272 sfrac; 2273 temp = one / (ptsaux[(ip << 1) + 1] * ptsaux[(iq << 1) + 1]); 2274 zmat[(k - np) * zmat_dim1 + 1] = temp; 2275 zmat[ip + 1 + (k - np) * zmat_dim1] = -temp; 2276 zmat[iq + 1 + (k - np) * zmat_dim1] = -temp; 2277 /* L70: */ 2278 zmat[k + (k - np) * zmat_dim1] = temp; 2279 } 2280 } 2281 nrem = npt; 2282 kold = 1; 2283 knew = kopt; 2284 2285 /* Reorder the provisional points in the way that exchanges PTSID(KOLD) */ 2286 /* with PTSID(KNEW). */ 2287 2288 L80: 2289 i__2 = n; 2290 for (j = 1; j <= i__2; ++j) { 2291 temp = bmat[kold + j * bmat_dim1]; 2292 bmat[kold + j * bmat_dim1] = bmat[knew + j * bmat_dim1]; 2293 /* L90: */ 2294 bmat[knew + j * bmat_dim1] = temp; 2295 } 2296 i__2 = nptm; 2297 for (j = 1; j <= i__2; ++j) { 2298 temp = zmat[kold + j * zmat_dim1]; 2299 zmat[kold + j * zmat_dim1] = zmat[knew + j * zmat_dim1]; 2300 /* L100: */ 2301 zmat[knew + j * zmat_dim1] = temp; 2302 } 2303 ptsid[kold] = ptsid[knew]; 2304 ptsid[knew] = zero; 2305 w[ndim + knew] = zero; 2306 --nrem; 2307 if (knew != kopt) { 2308 temp = vlag[kold]; 2309 vlag[kold] = vlag[knew]; 2310 vlag[knew] = temp; 2311 2312 /* Update the BMAT and ZMAT matrices so that the status of the KNEW-th */ 2313 /* interpolation point can be changed from provisional to original. The */ 2314 /* branch to label 350 occurs if all the original points are reinstated. */ 2315 /* The nonnegative values of W(NDIM+K) are required in the search below. */ 2316 2317 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], ndim, &vlag[1], 2318 beta, denom, knew, &w[1]); 2319 if (nrem == 0) { 2320 goto L350; 2321 } 2322 i__2 = npt; 2323 for (k = 1; k <= i__2; ++k) { 2324 /* L110: */ 2325 w[ndim + k] = (d__1 = w[ndim + k], std::abs(d__1)); 2326 } 2327 } 2328 2329 /* Pick the index KNEW of an original interpolation point that has not */ 2330 /* yet replaced one of the provisional interpolation points, giving */ 2331 /* attention to the closeness to XOPT and to previous tries with KNEW. */ 2332 2333 L120: 2334 dsqmin = zero; 2335 i__2 = npt; 2336 for (k = 1; k <= i__2; ++k) { 2337 if (w[ndim + k] > zero) { 2338 if (dsqmin == zero || w[ndim + k] < dsqmin) { 2339 knew = k; 2340 dsqmin = w[ndim + k]; 2341 } 2342 } 2343 /* L130: */ 2344 } 2345 if (dsqmin == zero) { 2346 goto L260; 2347 } 2348 2349 /* Form the W-vector of the chosen original interpolation point. */ 2350 2351 i__2 = n; 2352 for (j = 1; j <= i__2; ++j) { 2353 /* L140: */ 2354 w[npt + j] = xpt[knew + j * xpt_dim1]; 2355 } 2356 i__2 = npt; 2357 for (k = 1; k <= i__2; ++k) { 2358 sum = zero; 2359 if (k == kopt) { 2360 } else if (ptsid[k] == zero) { 2361 i__1 = n; 2362 for (j = 1; j <= i__1; ++j) { 2363 /* L150: */ 2364 sum += w[npt + j] * xpt[k + j * xpt_dim1]; 2365 } 2366 } else { 2367 ip = (integer) ptsid[k]; 2368 if (ip > 0) { 2369 sum = w[npt + ip] * ptsaux[(ip << 1) + 1]; 2370 } 2371 iq = (integer) ((doublereal) np * ptsid[k] - (doublereal) (ip * 2372 np)); 2373 if (iq > 0) { 2374 iw = 1; 2375 if (ip == 0) { 2376 iw = 2; 2377 } 2378 sum += w[npt + iq] * ptsaux[iw + (iq << 1)]; 2379 } 2380 } 2381 /* L160: */ 2382 w[k] = half * sum * sum; 2383 } 2384 2385 /* Calculate VLAG and BETA for the required updating of the H matrix if */ 2386 /* XPT(KNEW,.) is reinstated in the set of interpolation points. */ 2387 2388 i__2 = npt; 2389 for (k = 1; k <= i__2; ++k) { 2390 sum = zero; 2391 i__1 = n; 2392 for (j = 1; j <= i__1; ++j) { 2393 /* L170: */ 2394 sum += bmat[k + j * bmat_dim1] * w[npt + j]; 2395 } 2396 /* L180: */ 2397 vlag[k] = sum; 2398 } 2399 beta = zero; 2400 i__2 = nptm; 2401 for (j = 1; j <= i__2; ++j) { 2402 sum = zero; 2403 i__1 = npt; 2404 for (k = 1; k <= i__1; ++k) { 2405 /* L190: */ 2406 sum += zmat[k + j * zmat_dim1] * w[k]; 2407 } 2408 beta -= sum * sum; 2409 i__1 = npt; 2410 for (k = 1; k <= i__1; ++k) { 2411 /* L200: */ 2412 vlag[k] += sum * zmat[k + j * zmat_dim1]; 2413 } 2414 } 2415 bsum = zero; 2416 distsq = zero; 2417 i__1 = n; 2418 for (j = 1; j <= i__1; ++j) { 2419 sum = zero; 2420 i__2 = npt; 2421 for (k = 1; k <= i__2; ++k) { 2422 /* L210: */ 2423 sum += bmat[k + j * bmat_dim1] * w[k]; 2424 } 2425 jp = j + npt; 2426 bsum += sum * w[jp]; 2427 i__2 = ndim; 2428 for (ip = npt + 1; ip <= i__2; ++ip) { 2429 /* L220: */ 2430 sum += bmat[ip + j * bmat_dim1] * w[ip]; 2431 } 2432 bsum += sum * w[jp]; 2433 vlag[jp] = sum; 2434 /* L230: */ 2435 /* Computing 2nd power */ 2436 d__1 = xpt[knew + j * xpt_dim1]; 2437 distsq += d__1 * d__1; 2438 } 2439 beta = half * distsq * distsq + beta - bsum; 2440 vlag[kopt] += one; 2441 2442 /* KOLD is set to the index of the provisional interpolation point that is */ 2443 /* going to be deleted to make way for the KNEW-th original interpolation */ 2444 /* point. The choice of KOLD is governed by the avoidance of a small value */ 2445 /* of the denominator in the updating calculation of UPDATE. */ 2446 2447 denom = zero; 2448 vlmxsq = zero; 2449 i__1 = npt; 2450 for (k = 1; k <= i__1; ++k) { 2451 if (ptsid[k] != zero) { 2452 hdiag = zero; 2453 i__2 = nptm; 2454 for (j = 1; j <= i__2; ++j) { 2455 /* L240: */ 2456 /* Computing 2nd power */ 2457 d__1 = zmat[k + j * zmat_dim1]; 2458 hdiag += d__1 * d__1; 2459 } 2460 /* Computing 2nd power */ 2461 d__1 = vlag[k]; 2462 den = beta * hdiag + d__1 * d__1; 2463 if (den > denom) { 2464 kold = k; 2465 denom = den; 2466 } 2467 } 2468 /* L250: */ 2469 /* Computing MAX */ 2470 /* Computing 2nd power */ 2471 d__3 = vlag[k]; 2472 d__1 = vlmxsq, d__2 = d__3 * d__3; 2473 vlmxsq = std::max(d__1,d__2); 2474 } 2475 if (denom <= vlmxsq * .01) { 2476 w[ndim + knew] = -w[ndim + knew] - winc; 2477 goto L120; 2478 } 2479 goto L80; 2480 2481 /* When label 260 is reached, all the final positions of the interpolation */ 2482 /* points have been chosen although any changes have not been included yet */ 2483 /* in XPT. Also the final BMAT and ZMAT matrices are complete, but, apart */ 2484 /* from the shift of XBASE, the updating of the quadratic model remains to */ 2485 /* be done. The following cycle through the new interpolation points begins */ 2486 /* by putting the new point in XPT(KPT,.) and by setting PQ(KPT) to zero, */ 2487 /* except that a RETURN occurs if MAXFUN prohibits another value of F. */ 2488 2489 L260: 2490 i__1 = npt; 2491 for (kpt = 1; kpt <= i__1; ++kpt) { 2492 if (ptsid[kpt] == zero) { 2493 goto L340; 2494 } 2495 if (nf >= maxfun) { 2496 nf = -1; 2497 goto L350; 2498 } 2499 ih = 0; 2500 i__2 = n; 2501 for (j = 1; j <= i__2; ++j) { 2502 w[j] = xpt[kpt + j * xpt_dim1]; 2503 xpt[kpt + j * xpt_dim1] = zero; 2504 temp = pq[kpt] * w[j]; 2505 i__3 = j; 2506 for (i__ = 1; i__ <= i__3; ++i__) { 2507 ++ih; 2508 /* L270: */ 2509 hq[ih] += temp * w[i__]; 2510 } 2511 } 2512 pq[kpt] = zero; 2513 ip = (integer) ptsid[kpt]; 2514 iq = (integer) ((doublereal) np * ptsid[kpt] - (doublereal) (ip * np)) 2515 ; 2516 if (ip > 0) { 2517 xp = ptsaux[(ip << 1) + 1]; 2518 xpt[kpt + ip * xpt_dim1] = xp; 2519 } 2520 if (iq > 0) { 2521 xq = ptsaux[(iq << 1) + 1]; 2522 if (ip == 0) { 2523 xq = ptsaux[(iq << 1) + 2]; 2524 } 2525 xpt[kpt + iq * xpt_dim1] = xq; 2526 } 2527 2528 /* Set VQUAD to the value of the current model at the new point. */ 2529 2530 vquad = fbase; 2531 if (ip > 0) { 2532 ihp = (ip + ip * ip) / 2; 2533 vquad += xp * (gopt[ip] + half * xp * hq[ihp]); 2534 } 2535 if (iq > 0) { 2536 ihq = (iq + iq * iq) / 2; 2537 vquad += xq * (gopt[iq] + half * xq * hq[ihq]); 2538 if (ip > 0) { 2539 iw = std::max(ihp,ihq) - (i__3 = ip - iq, std::abs(i__3)); 2540 vquad += xp * xq * hq[iw]; 2541 } 2542 } 2543 i__3 = npt; 2544 for (k = 1; k <= i__3; ++k) { 2545 temp = zero; 2546 if (ip > 0) { 2547 temp += xp * xpt[k + ip * xpt_dim1]; 2548 } 2549 if (iq > 0) { 2550 temp += xq * xpt[k + iq * xpt_dim1]; 2551 } 2552 /* L280: */ 2553 vquad += half * pq[k] * temp * temp; 2554 } 2555 2556 /* Calculate F at the new interpolation point, and set DIFF to the factor */ 2557 /* that is going to multiply the KPT-th Lagrange function when the model */ 2558 /* is updated to provide interpolation to the new function value. */ 2559 2560 i__3 = n; 2561 for (i__ = 1; i__ <= i__3; ++i__) { 2562 /* Computing MIN */ 2563 /* Computing MAX */ 2564 d__3 = xl[i__], d__4 = xbase[i__] + xpt[kpt + i__ * xpt_dim1]; 2565 d__1 = std::max(d__3,d__4), d__2 = xu[i__]; 2566 w[i__] = std::min(d__1,d__2); 2567 if (xpt[kpt + i__ * xpt_dim1] == sl[i__]) { 2568 w[i__] = xl[i__]; 2569 } 2570 if (xpt[kpt + i__ * xpt_dim1] == su[i__]) { 2571 w[i__] = xu[i__]; 2572 } 2573 /* L290: */ 2574 } 2575 ++(nf); 2576 f = calfun(mat(&w[1],n)); 2577 fval[kpt] = f; 2578 if (f < fval[kopt]) { 2579 kopt = kpt; 2580 } 2581 diff = f - vquad; 2582 2583 /* Update the quadratic model. The RETURN from the subroutine occurs when */ 2584 /* all the new interpolation points are included in the model. */ 2585 2586 i__3 = n; 2587 for (i__ = 1; i__ <= i__3; ++i__) { 2588 /* L310: */ 2589 gopt[i__] += diff * bmat[kpt + i__ * bmat_dim1]; 2590 } 2591 i__3 = npt; 2592 for (k = 1; k <= i__3; ++k) { 2593 sum = zero; 2594 i__2 = nptm; 2595 for (j = 1; j <= i__2; ++j) { 2596 /* L320: */ 2597 sum += zmat[k + j * zmat_dim1] * zmat[kpt + j * zmat_dim1]; 2598 } 2599 temp = diff * sum; 2600 if (ptsid[k] == zero) { 2601 pq[k] += temp; 2602 } else { 2603 ip = (integer) ptsid[k]; 2604 iq = (integer) ((doublereal) np * ptsid[k] - (doublereal) (ip 2605 * np)); 2606 ihq = (iq * iq + iq) / 2; 2607 if (ip == 0) { 2608 /* Computing 2nd power */ 2609 d__1 = ptsaux[(iq << 1) + 2]; 2610 hq[ihq] += temp * (d__1 * d__1); 2611 } else { 2612 ihp = (ip * ip + ip) / 2; 2613 /* Computing 2nd power */ 2614 d__1 = ptsaux[(ip << 1) + 1]; 2615 hq[ihp] += temp * (d__1 * d__1); 2616 if (iq > 0) { 2617 /* Computing 2nd power */ 2618 d__1 = ptsaux[(iq << 1) + 1]; 2619 hq[ihq] += temp * (d__1 * d__1); 2620 iw = std::max(ihp,ihq) - (i__2 = iq - ip, std::abs(i__2)); 2621 hq[iw] += temp * ptsaux[(ip << 1) + 1] * ptsaux[(iq << 2622 1) + 1]; 2623 } 2624 } 2625 } 2626 /* L330: */ 2627 } 2628 ptsid[kpt] = zero; 2629 L340: 2630 ; 2631 } 2632 L350: 2633 ; 2634 } /* rescue_ */ 2635 2636 // ---------------------------------------------------------------------------------------- 2637 trsbox_(const integer n,const integer npt,const doublereal * xpt,const doublereal * xopt,const doublereal * gopt,const doublereal * hq,const doublereal * pq,const doublereal * sl,const doublereal * su,const doublereal delta,doublereal * xnew,doublereal * d__,doublereal * gnew,doublereal * xbdi,doublereal * s,doublereal * hs,doublereal * hred,doublereal * dsq,doublereal * crvmin)2638 void trsbox_( 2639 const integer n, 2640 const integer npt, 2641 const doublereal *xpt, 2642 const doublereal *xopt, 2643 const doublereal *gopt, 2644 const doublereal *hq, 2645 const doublereal *pq, 2646 const doublereal *sl, 2647 const doublereal *su, 2648 const doublereal delta, 2649 doublereal *xnew, 2650 doublereal *d__, 2651 doublereal *gnew, 2652 doublereal *xbdi, 2653 doublereal *s, 2654 doublereal *hs, 2655 doublereal *hred, 2656 doublereal *dsq, 2657 doublereal *crvmin 2658 ) const 2659 { 2660 /* System generated locals */ 2661 integer xpt_dim1, xpt_offset, i__1, i__2; 2662 doublereal d__1, d__2, d__3, d__4; 2663 2664 /* Local variables */ 2665 integer i__, j, k, ih; 2666 doublereal ds; 2667 integer iu; 2668 doublereal dhd, dhs, cth, one, shs, sth, ssq, half, beta, sdec, blen; 2669 integer iact = 0, nact = 0; 2670 doublereal angt, qred; 2671 integer isav; 2672 doublereal temp = 0, zero = 0, xsav = 0, xsum = 0, angbd = 0, dredg = 0, sredg = 0; 2673 integer iterc; 2674 doublereal resid = 0, delsq = 0, ggsav = 0, tempa = 0, tempb = 0, 2675 redmax = 0, dredsq = 0, redsav = 0, onemin = 0, gredsq = 0, rednew = 0; 2676 integer itcsav = 0; 2677 doublereal rdprev = 0, rdnext = 0, stplen = 0, stepsq = 0; 2678 integer itermax = 0; 2679 2680 2681 /* The arguments N, NPT, XPT, XOPT, GOPT, HQ, PQ, SL and SU have the same */ 2682 /* meanings as the corresponding arguments of BOBYQB. */ 2683 /* DELTA is the trust region radius for the present calculation, which */ 2684 /* seeks a small value of the quadratic model within distance DELTA of */ 2685 /* XOPT subject to the bounds on the variables. */ 2686 /* XNEW will be set to a new vector of variables that is approximately */ 2687 /* the one that minimizes the quadratic model within the trust region */ 2688 /* subject to the SL and SU constraints on the variables. It satisfies */ 2689 /* as equations the bounds that become active during the calculation. */ 2690 /* D is the calculated trial step from XOPT, generated iteratively from an */ 2691 /* initial value of zero. Thus XNEW is XOPT+D after the final iteration. */ 2692 /* GNEW holds the gradient of the quadratic model at XOPT+D. It is updated */ 2693 /* when D is updated. */ 2694 /* XBDI is a working space vector. For I=1,2,...,N, the element XBDI(I) is */ 2695 /* set to -1.0, 0.0, or 1.0, the value being nonzero if and only if the */ 2696 /* I-th variable has become fixed at a bound, the bound being SL(I) or */ 2697 /* SU(I) in the case XBDI(I)=-1.0 or XBDI(I)=1.0, respectively. This */ 2698 /* information is accumulated during the construction of XNEW. */ 2699 /* The arrays S, HS and HRED are also used for working space. They hold the */ 2700 /* current search direction, and the changes in the gradient of Q along S */ 2701 /* and the reduced D, respectively, where the reduced D is the same as D, */ 2702 /* except that the components of the fixed variables are zero. */ 2703 /* DSQ will be set to the square of the length of XNEW-XOPT. */ 2704 /* CRVMIN is set to zero if D reaches the trust region boundary. Otherwise */ 2705 /* it is set to the least curvature of H that occurs in the conjugate */ 2706 /* gradient searches that are not restricted by any constraints. The */ 2707 /* value CRVMIN=-1.0D0 is set, however, if all of these searches are */ 2708 /* constrained. */ 2709 2710 /* A version of the truncated conjugate gradient is applied. If a line */ 2711 /* search is restricted by a constraint, then the procedure is restarted, */ 2712 /* the values of the variables that are at their bounds being fixed. If */ 2713 /* the trust region boundary is reached, then further changes may be made */ 2714 /* to D, each one being in the two dimensional space that is spanned */ 2715 /* by the current D and the gradient of Q at XOPT+D, staying on the trust */ 2716 /* region boundary. Termination occurs when the reduction in Q seems to */ 2717 /* be close to the greatest reduction that can be achieved. */ 2718 2719 /* Set some constants. */ 2720 2721 /* Parameter adjustments */ 2722 xpt_dim1 = npt; 2723 xpt_offset = 1 + xpt_dim1; 2724 xpt -= xpt_offset; 2725 --xopt; 2726 --gopt; 2727 --hq; 2728 --pq; 2729 --sl; 2730 --su; 2731 --xnew; 2732 --d__; 2733 --gnew; 2734 --xbdi; 2735 --s; 2736 --hs; 2737 --hred; 2738 2739 /* Function Body */ 2740 half = .5; 2741 one = 1.; 2742 onemin = -1.; 2743 zero = 0.; 2744 2745 /* The sign of GOPT(I) gives the sign of the change to the I-th variable */ 2746 /* that will reduce Q from its value at XOPT. Thus XBDI(I) shows whether */ 2747 /* or not to fix the I-th variable at one of its bounds initially, with */ 2748 /* NACT being set to the number of fixed variables. D and GNEW are also */ 2749 /* set for the first iteration. DELSQ is the upper bound on the sum of */ 2750 /* squares of the free variables. QRED is the reduction in Q so far. */ 2751 2752 iterc = 0; 2753 nact = 0; 2754 i__1 = n; 2755 for (i__ = 1; i__ <= i__1; ++i__) { 2756 xbdi[i__] = zero; 2757 if (xopt[i__] <= sl[i__]) { 2758 if (gopt[i__] >= zero) { 2759 xbdi[i__] = onemin; 2760 } 2761 } else if (xopt[i__] >= su[i__]) { 2762 if (gopt[i__] <= zero) { 2763 xbdi[i__] = one; 2764 } 2765 } 2766 if (xbdi[i__] != zero) { 2767 ++nact; 2768 } 2769 d__[i__] = zero; 2770 /* L10: */ 2771 gnew[i__] = gopt[i__]; 2772 } 2773 delsq = delta * delta; 2774 qred = zero; 2775 *crvmin = onemin; 2776 2777 /* Set the next search direction of the conjugate gradient method. It is */ 2778 /* the steepest descent direction initially and when the iterations are */ 2779 /* restarted because a variable has just been fixed by a bound, and of */ 2780 /* course the components of the fixed variables are zero. ITERMAX is an */ 2781 /* upper bound on the indices of the conjugate gradient iterations. */ 2782 2783 L20: 2784 beta = zero; 2785 L30: 2786 stepsq = zero; 2787 i__1 = n; 2788 for (i__ = 1; i__ <= i__1; ++i__) { 2789 if (xbdi[i__] != zero) { 2790 s[i__] = zero; 2791 } else if (beta == zero) { 2792 s[i__] = -gnew[i__]; 2793 } else { 2794 s[i__] = beta * s[i__] - gnew[i__]; 2795 } 2796 /* L40: */ 2797 /* Computing 2nd power */ 2798 d__1 = s[i__]; 2799 stepsq += d__1 * d__1; 2800 } 2801 if (stepsq == zero) { 2802 goto L190; 2803 } 2804 if (beta == zero) { 2805 gredsq = stepsq; 2806 itermax = iterc + n - nact; 2807 } 2808 if (gredsq * delsq <= qred * 1e-4 * qred) { 2809 goto L190; 2810 } 2811 2812 /* Multiply the search direction by the second derivative matrix of Q and */ 2813 /* calculate some scalars for the choice of steplength. Then set BLEN to */ 2814 /* the length of the the step to the trust region boundary and STPLEN to */ 2815 /* the steplength, ignoring the simple bounds. */ 2816 2817 goto L210; 2818 L50: 2819 resid = delsq; 2820 ds = zero; 2821 shs = zero; 2822 i__1 = n; 2823 for (i__ = 1; i__ <= i__1; ++i__) { 2824 if (xbdi[i__] == zero) { 2825 /* Computing 2nd power */ 2826 d__1 = d__[i__]; 2827 resid -= d__1 * d__1; 2828 ds += s[i__] * d__[i__]; 2829 shs += s[i__] * hs[i__]; 2830 } 2831 /* L60: */ 2832 } 2833 if (resid <= zero) { 2834 goto L90; 2835 } 2836 temp = std::sqrt(stepsq * resid + ds * ds); 2837 if (ds < zero) { 2838 blen = (temp - ds) / stepsq; 2839 } else { 2840 blen = resid / (temp + ds); 2841 } 2842 stplen = blen; 2843 if (shs > zero) { 2844 /* Computing MIN */ 2845 d__1 = blen, d__2 = gredsq / shs; 2846 stplen = std::min(d__1,d__2); 2847 } 2848 2849 /* Reduce STPLEN if necessary in order to preserve the simple bounds, */ 2850 /* letting IACT be the index of the new constrained variable. */ 2851 2852 iact = 0; 2853 i__1 = n; 2854 for (i__ = 1; i__ <= i__1; ++i__) { 2855 if (s[i__] != zero) { 2856 xsum = xopt[i__] + d__[i__]; 2857 if (s[i__] > zero) { 2858 temp = (su[i__] - xsum) / s[i__]; 2859 } else { 2860 temp = (sl[i__] - xsum) / s[i__]; 2861 } 2862 if (temp < stplen) { 2863 stplen = temp; 2864 iact = i__; 2865 } 2866 } 2867 /* L70: */ 2868 } 2869 2870 /* Update CRVMIN, GNEW and D. Set SDEC to the decrease that occurs in Q. */ 2871 2872 sdec = zero; 2873 if (stplen > zero) { 2874 ++iterc; 2875 temp = shs / stepsq; 2876 if (iact == 0 && temp > zero) { 2877 *crvmin = std::min(*crvmin,temp); 2878 if (*crvmin == onemin) { 2879 *crvmin = temp; 2880 } 2881 } 2882 ggsav = gredsq; 2883 gredsq = zero; 2884 i__1 = n; 2885 for (i__ = 1; i__ <= i__1; ++i__) { 2886 gnew[i__] += stplen * hs[i__]; 2887 if (xbdi[i__] == zero) { 2888 /* Computing 2nd power */ 2889 d__1 = gnew[i__]; 2890 gredsq += d__1 * d__1; 2891 } 2892 /* L80: */ 2893 d__[i__] += stplen * s[i__]; 2894 } 2895 /* Computing MAX */ 2896 d__1 = stplen * (ggsav - half * stplen * shs); 2897 sdec = std::max(d__1,zero); 2898 qred += sdec; 2899 } 2900 2901 /* Restart the conjugate gradient method if it has hit a new bound. */ 2902 2903 if (iact > 0) { 2904 ++nact; 2905 xbdi[iact] = one; 2906 if (s[iact] < zero) { 2907 xbdi[iact] = onemin; 2908 } 2909 /* Computing 2nd power */ 2910 d__1 = d__[iact]; 2911 delsq -= d__1 * d__1; 2912 if (delsq <= zero) { 2913 goto L90; 2914 } 2915 goto L20; 2916 } 2917 2918 /* If STPLEN is less than BLEN, then either apply another conjugate */ 2919 /* gradient iteration or RETURN. */ 2920 2921 if (stplen < blen) { 2922 if (iterc == itermax) { 2923 goto L190; 2924 } 2925 if (sdec <= qred * .01) { 2926 goto L190; 2927 } 2928 beta = gredsq / ggsav; 2929 goto L30; 2930 } 2931 L90: 2932 *crvmin = zero; 2933 2934 /* Prepare for the alternative iteration by calculating some scalars */ 2935 /* and by multiplying the reduced D by the second derivative matrix of */ 2936 /* Q, where S holds the reduced D in the call of GGMULT. */ 2937 2938 L100: 2939 if (nact >= n - 1) { 2940 goto L190; 2941 } 2942 dredsq = zero; 2943 dredg = zero; 2944 gredsq = zero; 2945 i__1 = n; 2946 for (i__ = 1; i__ <= i__1; ++i__) { 2947 if (xbdi[i__] == zero) { 2948 /* Computing 2nd power */ 2949 d__1 = d__[i__]; 2950 dredsq += d__1 * d__1; 2951 dredg += d__[i__] * gnew[i__]; 2952 /* Computing 2nd power */ 2953 d__1 = gnew[i__]; 2954 gredsq += d__1 * d__1; 2955 s[i__] = d__[i__]; 2956 } else { 2957 s[i__] = zero; 2958 } 2959 /* L110: */ 2960 } 2961 itcsav = iterc; 2962 goto L210; 2963 2964 /* Let the search direction S be a linear combination of the reduced D */ 2965 /* and the reduced G that is orthogonal to the reduced D. */ 2966 2967 L120: 2968 ++iterc; 2969 temp = gredsq * dredsq - dredg * dredg; 2970 if (temp <= qred * 1e-4 * qred) { 2971 goto L190; 2972 } 2973 temp = std::sqrt(temp); 2974 i__1 = n; 2975 for (i__ = 1; i__ <= i__1; ++i__) { 2976 if (xbdi[i__] == zero) { 2977 s[i__] = (dredg * d__[i__] - dredsq * gnew[i__]) / temp; 2978 } else { 2979 s[i__] = zero; 2980 } 2981 /* L130: */ 2982 } 2983 sredg = -temp; 2984 2985 /* By considering the simple bounds on the variables, calculate an upper */ 2986 /* bound on the tangent of half the angle of the alternative iteration, */ 2987 /* namely ANGBD, except that, if already a free variable has reached a */ 2988 /* bound, there is a branch back to label 100 after fixing that variable. */ 2989 2990 angbd = one; 2991 iact = 0; 2992 i__1 = n; 2993 for (i__ = 1; i__ <= i__1; ++i__) { 2994 if (xbdi[i__] == zero) { 2995 tempa = xopt[i__] + d__[i__] - sl[i__]; 2996 tempb = su[i__] - xopt[i__] - d__[i__]; 2997 if (tempa <= zero) { 2998 ++nact; 2999 xbdi[i__] = onemin; 3000 goto L100; 3001 } else if (tempb <= zero) { 3002 ++nact; 3003 xbdi[i__] = one; 3004 goto L100; 3005 } 3006 /* Computing 2nd power */ 3007 d__1 = d__[i__]; 3008 /* Computing 2nd power */ 3009 d__2 = s[i__]; 3010 ssq = d__1 * d__1 + d__2 * d__2; 3011 /* Computing 2nd power */ 3012 d__1 = xopt[i__] - sl[i__]; 3013 temp = ssq - d__1 * d__1; 3014 if (temp > zero) { 3015 temp = std::sqrt(temp) - s[i__]; 3016 if (angbd * temp > tempa) { 3017 angbd = tempa / temp; 3018 iact = i__; 3019 xsav = onemin; 3020 } 3021 } 3022 /* Computing 2nd power */ 3023 d__1 = su[i__] - xopt[i__]; 3024 temp = ssq - d__1 * d__1; 3025 if (temp > zero) { 3026 temp = std::sqrt(temp) + s[i__]; 3027 if (angbd * temp > tempb) { 3028 angbd = tempb / temp; 3029 iact = i__; 3030 xsav = one; 3031 } 3032 } 3033 } 3034 /* L140: */ 3035 } 3036 3037 /* Calculate HHD and some curvatures for the alternative iteration. */ 3038 3039 goto L210; 3040 L150: 3041 shs = zero; 3042 dhs = zero; 3043 dhd = zero; 3044 i__1 = n; 3045 for (i__ = 1; i__ <= i__1; ++i__) { 3046 if (xbdi[i__] == zero) { 3047 shs += s[i__] * hs[i__]; 3048 dhs += d__[i__] * hs[i__]; 3049 dhd += d__[i__] * hred[i__]; 3050 } 3051 /* L160: */ 3052 } 3053 3054 /* Seek the greatest reduction in Q for a range of equally spaced values */ 3055 /* of ANGT in [0,ANGBD], where ANGT is the tangent of half the angle of */ 3056 /* the alternative iteration. */ 3057 3058 redmax = zero; 3059 isav = 0; 3060 redsav = zero; 3061 iu = (integer) (angbd * 17. + 3.1); 3062 i__1 = iu; 3063 for (i__ = 1; i__ <= i__1; ++i__) { 3064 angt = angbd * (doublereal) i__ / (doublereal) iu; 3065 sth = (angt + angt) / (one + angt * angt); 3066 temp = shs + angt * (angt * dhd - dhs - dhs); 3067 rednew = sth * (angt * dredg - sredg - half * sth * temp); 3068 if (rednew > redmax) { 3069 redmax = rednew; 3070 isav = i__; 3071 rdprev = redsav; 3072 } else if (i__ == isav + 1) { 3073 rdnext = rednew; 3074 } 3075 /* L170: */ 3076 redsav = rednew; 3077 } 3078 3079 /* Return if the reduction is zero. Otherwise, set the sine and cosine */ 3080 /* of the angle of the alternative iteration, and calculate SDEC. */ 3081 3082 if (isav == 0) { 3083 goto L190; 3084 } 3085 if (isav < iu) { 3086 temp = (rdnext - rdprev) / (redmax + redmax - rdprev - rdnext); 3087 angt = angbd * ((doublereal) isav + half * temp) / (doublereal) iu; 3088 } 3089 cth = (one - angt * angt) / (one + angt * angt); 3090 sth = (angt + angt) / (one + angt * angt); 3091 temp = shs + angt * (angt * dhd - dhs - dhs); 3092 sdec = sth * (angt * dredg - sredg - half * sth * temp); 3093 if (sdec <= zero) { 3094 goto L190; 3095 } 3096 3097 /* Update GNEW, D and HRED. If the angle of the alternative iteration */ 3098 /* is restricted by a bound on a free variable, that variable is fixed */ 3099 /* at the bound. */ 3100 3101 dredg = zero; 3102 gredsq = zero; 3103 i__1 = n; 3104 for (i__ = 1; i__ <= i__1; ++i__) { 3105 gnew[i__] = gnew[i__] + (cth - one) * hred[i__] + sth * hs[i__]; 3106 if (xbdi[i__] == zero) { 3107 d__[i__] = cth * d__[i__] + sth * s[i__]; 3108 dredg += d__[i__] * gnew[i__]; 3109 /* Computing 2nd power */ 3110 d__1 = gnew[i__]; 3111 gredsq += d__1 * d__1; 3112 } 3113 /* L180: */ 3114 hred[i__] = cth * hred[i__] + sth * hs[i__]; 3115 } 3116 qred += sdec; 3117 if (iact > 0 && isav == iu) { 3118 ++nact; 3119 xbdi[iact] = xsav; 3120 goto L100; 3121 } 3122 3123 /* If SDEC is sufficiently small, then RETURN after setting XNEW to */ 3124 /* XOPT+D, giving careful attention to the bounds. */ 3125 3126 if (sdec > qred * .01) { 3127 goto L120; 3128 } 3129 L190: 3130 *dsq = zero; 3131 i__1 = n; 3132 for (i__ = 1; i__ <= i__1; ++i__) { 3133 /* Computing MAX */ 3134 /* Computing MIN */ 3135 d__3 = xopt[i__] + d__[i__], d__4 = su[i__]; 3136 d__1 = std::min(d__3,d__4), d__2 = sl[i__]; 3137 xnew[i__] = std::max(d__1,d__2); 3138 if (xbdi[i__] == onemin) { 3139 xnew[i__] = sl[i__]; 3140 } 3141 if (xbdi[i__] == one) { 3142 xnew[i__] = su[i__]; 3143 } 3144 d__[i__] = xnew[i__] - xopt[i__]; 3145 /* L200: */ 3146 /* Computing 2nd power */ 3147 d__1 = d__[i__]; 3148 *dsq += d__1 * d__1; 3149 } 3150 return; 3151 /* The following instructions multiply the current S-vector by the second */ 3152 /* derivative matrix of the quadratic model, putting the product in HS. */ 3153 /* They are reached from three different parts of the software above and */ 3154 /* they can be regarded as an external subroutine. */ 3155 3156 L210: 3157 ih = 0; 3158 i__1 = n; 3159 for (j = 1; j <= i__1; ++j) { 3160 hs[j] = zero; 3161 i__2 = j; 3162 for (i__ = 1; i__ <= i__2; ++i__) { 3163 ++ih; 3164 if (i__ < j) { 3165 hs[j] += hq[ih] * s[i__]; 3166 } 3167 /* L220: */ 3168 hs[i__] += hq[ih] * s[j]; 3169 } 3170 } 3171 i__2 = npt; 3172 for (k = 1; k <= i__2; ++k) { 3173 if (pq[k] != zero) { 3174 temp = zero; 3175 i__1 = n; 3176 for (j = 1; j <= i__1; ++j) { 3177 /* L230: */ 3178 temp += xpt[k + j * xpt_dim1] * s[j]; 3179 } 3180 temp *= pq[k]; 3181 i__1 = n; 3182 for (i__ = 1; i__ <= i__1; ++i__) { 3183 /* L240: */ 3184 hs[i__] += temp * xpt[k + i__ * xpt_dim1]; 3185 } 3186 } 3187 /* L250: */ 3188 } 3189 if (*crvmin != zero) { 3190 goto L50; 3191 } 3192 if (iterc > itcsav) { 3193 goto L150; 3194 } 3195 i__2 = n; 3196 for (i__ = 1; i__ <= i__2; ++i__) { 3197 /* L260: */ 3198 hred[i__] = hs[i__]; 3199 } 3200 goto L120; 3201 } /* trsbox_ */ 3202 3203 // ---------------------------------------------------------------------------------------- 3204 update_(const integer n,const integer npt,doublereal * bmat,doublereal * zmat,const integer ndim,doublereal * vlag,const doublereal beta,const doublereal denom,const integer knew,doublereal * w)3205 void update_( 3206 const integer n, 3207 const integer npt, 3208 doublereal *bmat, 3209 doublereal *zmat, 3210 const integer ndim, 3211 doublereal *vlag, 3212 const doublereal beta, 3213 const doublereal denom, 3214 const integer knew, 3215 doublereal *w 3216 ) const 3217 { 3218 /* System generated locals */ 3219 integer bmat_dim1, bmat_offset, zmat_dim1, zmat_offset, i__1, i__2; 3220 doublereal d__1, d__2, d__3; 3221 3222 /* Local variables */ 3223 integer i__, j, k, jp; 3224 doublereal one, tau, temp; 3225 integer nptm; 3226 doublereal zero, alpha, tempa, tempb, ztest; 3227 3228 3229 /* The arrays BMAT and ZMAT are updated, as required by the new position */ 3230 /* of the interpolation point that has the index KNEW. The vector VLAG has */ 3231 /* N+NPT components, set on entry to the first NPT and last N components */ 3232 /* of the product Hw in equation (4.11) of the Powell (2006) paper on */ 3233 /* NEWUOA. Further, BETA is set on entry to the value of the parameter */ 3234 /* with that name, and DENOM is set to the denominator of the updating */ 3235 /* formula. Elements of ZMAT may be treated as zero if their moduli are */ 3236 /* at most ZTEST. The first NDIM elements of W are used for working space. */ 3237 3238 /* Set some constants. */ 3239 3240 /* Parameter adjustments */ 3241 zmat_dim1 = npt; 3242 zmat_offset = 1 + zmat_dim1; 3243 zmat -= zmat_offset; 3244 bmat_dim1 = ndim; 3245 bmat_offset = 1 + bmat_dim1; 3246 bmat -= bmat_offset; 3247 --vlag; 3248 --w; 3249 3250 /* Function Body */ 3251 one = 1.; 3252 zero = 0.; 3253 nptm = npt - n - 1; 3254 ztest = zero; 3255 i__1 = npt; 3256 for (k = 1; k <= i__1; ++k) { 3257 i__2 = nptm; 3258 for (j = 1; j <= i__2; ++j) { 3259 /* L10: */ 3260 /* Computing MAX */ 3261 d__2 = ztest, d__3 = (d__1 = zmat[k + j * zmat_dim1], std::abs(d__1)); 3262 ztest = std::max(d__2,d__3); 3263 } 3264 } 3265 ztest *= 1e-20; 3266 3267 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */ 3268 3269 i__2 = nptm; 3270 for (j = 2; j <= i__2; ++j) { 3271 if ((d__1 = zmat[knew + j * zmat_dim1], std::abs(d__1)) > ztest) { 3272 /* Computing 2nd power */ 3273 d__1 = zmat[knew + zmat_dim1]; 3274 /* Computing 2nd power */ 3275 d__2 = zmat[knew + j * zmat_dim1]; 3276 temp = std::sqrt(d__1 * d__1 + d__2 * d__2); 3277 tempa = zmat[knew + zmat_dim1] / temp; 3278 tempb = zmat[knew + j * zmat_dim1] / temp; 3279 i__1 = npt; 3280 for (i__ = 1; i__ <= i__1; ++i__) { 3281 temp = tempa * zmat[i__ + zmat_dim1] + tempb * zmat[i__ + j * 3282 zmat_dim1]; 3283 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1] 3284 - tempb * zmat[i__ + zmat_dim1]; 3285 /* L20: */ 3286 zmat[i__ + zmat_dim1] = temp; 3287 } 3288 } 3289 zmat[knew + j * zmat_dim1] = zero; 3290 /* L30: */ 3291 } 3292 3293 /* Put the first NPT components of the KNEW-th column of HLAG into W, */ 3294 /* and calculate the parameters of the updating formula. */ 3295 3296 i__2 = npt; 3297 for (i__ = 1; i__ <= i__2; ++i__) { 3298 w[i__] = zmat[knew + zmat_dim1] * zmat[i__ + zmat_dim1]; 3299 /* L40: */ 3300 } 3301 alpha = w[knew]; 3302 tau = vlag[knew]; 3303 vlag[knew] -= one; 3304 3305 /* Complete the updating of ZMAT. */ 3306 3307 temp = std::sqrt(denom); 3308 tempb = zmat[knew + zmat_dim1] / temp; 3309 tempa = tau / temp; 3310 i__2 = npt; 3311 for (i__ = 1; i__ <= i__2; ++i__) { 3312 /* L50: */ 3313 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb * vlag[ 3314 i__]; 3315 } 3316 3317 /* Finally, update the matrix BMAT. */ 3318 3319 i__2 = n; 3320 for (j = 1; j <= i__2; ++j) { 3321 jp = npt + j; 3322 w[jp] = bmat[knew + j * bmat_dim1]; 3323 tempa = (alpha * vlag[jp] - tau * w[jp]) / denom; 3324 tempb = (-(beta) * w[jp] - tau * vlag[jp]) / denom; 3325 i__1 = jp; 3326 for (i__ = 1; i__ <= i__1; ++i__) { 3327 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa * 3328 vlag[i__] + tempb * w[i__]; 3329 if (i__ > npt) { 3330 bmat[jp + (i__ - npt) * bmat_dim1] = bmat[i__ + j * 3331 bmat_dim1]; 3332 } 3333 /* L60: */ 3334 } 3335 } 3336 } /* update_ */ 3337 }; 3338 3339 // ---------------------------------------------------------------------------------------- 3340 3341 template < 3342 typename funct, 3343 typename T, 3344 typename U 3345 > find_min_bobyqa(const funct & f,T & x,long npt,const U & x_lower,const U & x_upper,const double rho_begin,const double rho_end,const long max_f_evals)3346 double find_min_bobyqa ( 3347 const funct& f, 3348 T& x, 3349 long npt, 3350 const U& x_lower, 3351 const U& x_upper, 3352 const double rho_begin, 3353 const double rho_end, 3354 const long max_f_evals 3355 ) 3356 { 3357 // The starting point (i.e. x) must be a column vector. 3358 COMPILE_TIME_ASSERT(T::NC <= 1); 3359 3360 // check the requirements. Also split the assert up so that the error message isn't huge. 3361 DLIB_CASSERT(is_col_vector(x) && is_col_vector(x_lower) && is_col_vector(x_upper) && 3362 x.size() == x_lower.size() && x_lower.size() == x_upper.size() && 3363 x.size() > 1 && max_f_evals > 1, 3364 "\tdouble find_min_bobyqa()" 3365 << "\n\t Invalid arguments have been given to this function" 3366 << "\n\t is_col_vector(x): " << is_col_vector(x) 3367 << "\n\t is_col_vector(x_lower): " << is_col_vector(x_lower) 3368 << "\n\t is_col_vector(x_upper): " << is_col_vector(x_upper) 3369 << "\n\t x.size(): " << x.size() 3370 << "\n\t x_lower.size(): " << x_lower.size() 3371 << "\n\t x_upper.size(): " << x_upper.size() 3372 << "\n\t max_f_evals: " << max_f_evals 3373 ); 3374 3375 DLIB_CASSERT(x.size() + 2 <= npt && npt <= (x.size()+1)*(x.size()+2)/2 && 3376 0 < rho_end && rho_end < rho_begin && 3377 min(x_upper - x_lower) > 2*rho_begin && 3378 min(x - x_lower) >= 0 && min(x_upper - x) >= 0, 3379 "\tdouble find_min_bobyqa()" 3380 << "\n\t Invalid arguments have been given to this function" 3381 << "\n\t ntp in valid range: " << (x.size() + 2 <= npt && npt <= (x.size()+1)*(x.size()+2)/2) 3382 << "\n\t npt: " << npt 3383 << "\n\t rho_begin: " << rho_begin 3384 << "\n\t rho_end: " << rho_end 3385 << "\n\t min(x_upper - x_lower) > 2*rho_begin: " << (min(x_upper - x_lower) > 2*rho_begin) 3386 << "\n\t min(x - x_lower) >= 0 && min(x_upper - x) >= 0: " << (min(x - x_lower) >= 0 && min(x_upper - x) >= 0) 3387 ); 3388 3389 3390 bobyqa_implementation impl; 3391 return impl.find_min(f, x, npt, x_lower, x_upper, rho_begin, rho_end, max_f_evals); 3392 } 3393 3394 // ---------------------------------------------------------------------------------------- 3395 3396 template < 3397 typename funct, 3398 typename T, 3399 typename U 3400 > find_max_bobyqa(const funct & f,T & x,long npt,const U & x_lower,const U & x_upper,const double rho_begin,const double rho_end,const long max_f_evals)3401 double find_max_bobyqa ( 3402 const funct& f, 3403 T& x, 3404 long npt, 3405 const U& x_lower, 3406 const U& x_upper, 3407 const double rho_begin, 3408 const double rho_end, 3409 const long max_f_evals 3410 ) 3411 { 3412 // The starting point (i.e. x) must be a column vector. 3413 COMPILE_TIME_ASSERT(T::NC <= 1); 3414 3415 return -find_min_bobyqa(negate_function(f), x, npt, x_lower, x_upper, rho_begin, rho_end, max_f_evals); 3416 } 3417 3418 // ---------------------------------------------------------------------------------------- 3419 3420 } 3421 3422 #endif // DLIB_OPTIMIZATIOn_BOBYQA_Hh_ 3423 3424