1 #ifndef __POWELL_NEWUOA__
2 #define __POWELL_NEWUOA__
3
4 #undef AFmax
5 #undef AFmin
6 #define AFmax(a,b) (((a)<(b)) ? (b) : (a))
7 #define AFmin(a,b) (((a)>(b)) ? (b) : (a))
8
9 #define DEBUG 0
10
11 /*---------------------------------------------------------------------------*/
12
13 /* powell_newuoa.f -- translated by f2c (version 19961017).
14 You must link the resulting object file with the libraries:
15 -lf2c -lm (in that order)
16 */
17
18 #undef STATIC
19 #define STATIC /*static*/ /*** disable use of static internal variables ***/
20 /*** all local variables instead init to =0 ***/
21
22 #include <stdio.h>
23 #include "f2c.h"
24
25 /* CC SUBROUTINE NEWUOA (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W) */
newuoa_(integer * n,integer * npt,doublereal * x,doublereal * rhobeg,doublereal * rhoend,integer * maxfun,doublereal * w,integer * icode)26 /* Subroutine */ int newuoa_(integer *n, integer *npt, doublereal *x,
27 doublereal *rhobeg, doublereal *rhoend, integer *maxfun, doublereal *
28 w, integer *icode)
29 {
30 STATIC integer ndim=0, nptm=0, ibmat=0, izmat=0, id=0, np=0, iw=0;
31 extern /* Subroutine */ int newuob_(integer *, integer *, doublereal *,
32 doublereal *, doublereal *, integer *, doublereal *, doublereal *,
33 doublereal *, doublereal *, doublereal *, doublereal *,
34 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
35 doublereal *, doublereal *, doublereal *, integer *);
36 STATIC integer igq=0, ihq=0, ixb=0, ifv=0, ipq=0, ivl=0, ixn=0, ixo=0, ixp=0;
37
38
39 /* This subroutine seeks the least value of a function of many variable
40 s,*/
41 /* by a trust region method that forms quadratic models by interpolatio
42 n.*/
43 /* There can be some freedom in the interpolation conditions, which is
44 */
45 /* taken up by minimizing the Frobenius norm of the change to the secon
46 d*/
47 /* derivative of the quadratic model, beginning with a zero matrix. The
48 */
49 /* arguments of the subroutine are as follows. */
50
51 /* N must be set to the number of variables and must be at least two.
52 */
53 /* NPT is the number of interpolation conditions. Its value must be in
54 the*/
55 /* interval [N+2,(N+1)(N+2)/2]. */
56 /* Initial values of the variables must be set in X(1),X(2),...,X(N). T
57 hey*/
58 /* will be changed to the values that give the least calculated F.
59 */
60 /* RHOBEG and RHOEND must be set to the initial and final values of a t
61 rust*/
62 /* region radius, so both must be positive with RHOEND<=RHOBEG. Typic
63 ally*/
64 /* RHOBEG should be about one tenth of the greatest expected change t
65 o a*/
66 /* variable, and RHOEND should indicate the accuracy that is required
67 in*/
68 /* the final values of the variables. */
69 /*CCC The value of IPRINT should be set to 0, 1, 2 or 3, which control
70 s the*/
71 /*CCC amount of printing. Specifically, there is no output if IPRINT
72 =0 and*/
73 /*CCC there is output only at the return if IPRINT=1. Otherwise, eac
74 h new*/
75 /*CCC value of RHO is printed, with the best vector of variables so
76 far and*/
77 /*CCC the corresponding value of the objective function. Further, ea
78 ch new*/
79 /* CCC value of F with its variables are output if IPRINT=3. */
80 /* MAXFUN must be set to an upper bound on the number of calls of CALFU
81 N.*/
82 /* The array W will be used for working space. Its length must be at le
83 ast*/
84 /* (NPT+13)*(NPT+N)+3*N*(N+3)/2. */
85
86 /* SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set
87 F to*/
88 /* the value of the objective function for the variables X(1),X(2),...,
89 X(N).*/
90
91 /* Partition the working space array, so that different parts of it can
92 be*/
93 /* treated separately by the subroutine that performs the main calculat
94 ion.*/
95
96 /* Parameter adjustments */
97 --w;
98 --x;
99
100 /* Function Body */
101 np = *n + 1;
102 nptm = *npt - np;
103 /* CC IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN */
104 /* CC PRINT 10 */
105 /* CC 10 FORMAT (/4X,'Return from NEWUOA because NPT is not in', */
106 /* CC 1 ' the required interval') */
107 /* CC GO TO 20 */
108 /* CC END IF */
109 ndim = *npt + *n;
110 ixb = 1;
111 ixo = ixb + *n;
112 ixn = ixo + *n;
113 ixp = ixn + *n;
114 ifv = ixp + *n * *npt;
115 igq = ifv + *npt;
116 ihq = igq + *n;
117 ipq = ihq + *n * np / 2;
118 ibmat = ipq + *npt;
119 izmat = ibmat + ndim * *n;
120 id = izmat + *npt * nptm;
121 ivl = id + *n;
122 iw = ivl + ndim;
123
124 /* The above settings provide a partition of W for subroutine NEWUOB.
125 */
126 /* The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of
127 */
128 /* W plus the space that is needed by the last array of NEWUOB. */
129
130 /* CC CALL NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB), */
131 newuob_(n, npt, &x[1], rhobeg, rhoend, maxfun, &w[ixb], &w[ixo], &w[ixn],
132 &w[ixp], &w[ifv], &w[igq], &w[ihq], &w[ipq], &w[ibmat], &w[izmat],
133 &ndim, &w[id], &w[ivl], &w[iw], icode);
134 /* L20: */
135 return 0;
136 } /* newuoa_ */
137
138 /* CC SUBROUTINE NEWUOB (N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,XBASE, */
newuob_(integer * n,integer * npt,doublereal * x,doublereal * rhobeg,doublereal * rhoend,integer * maxfun,doublereal * xbase,doublereal * xopt,doublereal * xnew,doublereal * xpt,doublereal * fval,doublereal * gq,doublereal * hq,doublereal * pq,doublereal * bmat,doublereal * zmat,integer * ndim,doublereal * d__,doublereal * vlag,doublereal * w,integer * icode)139 /* Subroutine */ int newuob_(integer *n, integer *npt, doublereal *x,
140 doublereal *rhobeg, doublereal *rhoend, integer *maxfun, doublereal *
141 xbase, doublereal *xopt, doublereal *xnew, doublereal *xpt,
142 doublereal *fval, doublereal *gq, doublereal *hq, doublereal *pq,
143 doublereal *bmat, doublereal *zmat, integer *ndim, doublereal *d__,
144 doublereal *vlag, doublereal *w, integer *icode)
145 {
146 /* System generated locals */
147 integer xpt_dim1=0, xpt_offset=0, bmat_dim1=0, bmat_offset=0, zmat_dim1=0,
148 zmat_offset=0, i__1=0, i__2=0, i__3=0;
149 doublereal d__1=0, d__2=0, d__3=0;
150
151 /* Builtin functions */
152 double sqrt(doublereal);
153
154 /* Local variables */
155 STATIC doublereal fbeg=0, diff=0, half=0, beta=0;
156 STATIC integer nfmm=0;
157 STATIC doublereal gisq=0;
158 STATIC integer knew=0;
159 STATIC doublereal temp=0, suma=0, sumb=0, fopt=0, bsum=0, gqsq=0;
160 STATIC integer kopt=0, nptm=0;
161 STATIC doublereal zero=0, xipt=0, xjpt=0, sumz=0, f=0;
162 STATIC integer i__=0, j=0, k=0;
163 STATIC doublereal diffa=0, diffb=0, diffc=0, hdiag=0, alpha=0, delta=0, recip=0, reciq=0,
164 fsave=0;
165 STATIC integer ksave=0, nfsav=0, itemp=0;
166 STATIC doublereal dnorm=0, ratio=0, dstep=0, tenth=0, vquad=0;
167 STATIC integer ktemp=0;
168 STATIC doublereal tempq=0;
169 STATIC integer itest=0;
170 STATIC doublereal rhosq=0;
171 STATIC integer ih=0, nf=0;
172 extern /* Subroutine */ int biglag_(integer *, integer *, doublereal *,
173 doublereal *, doublereal *, doublereal *, integer *, integer *,
174 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
175 doublereal *, doublereal *, doublereal *, doublereal *);
176 STATIC integer nh=0, ip=0, jp=0;
177 STATIC doublereal dx=0;
178 extern /* Subroutine */ int bigden_(integer *, integer *, doublereal *,
179 doublereal *, doublereal *, doublereal *, integer *, integer *,
180 integer *, integer *, doublereal *, doublereal *, doublereal *,
181 doublereal *, doublereal *, doublereal *, doublereal *);
182 STATIC integer np=0;
183 extern /* Subroutine */ int calfun_(integer *, doublereal *, doublereal *)
184 , update_(integer *, integer *, doublereal *, doublereal *,
185 integer *, integer *, doublereal *, doublereal *, integer *,
186 doublereal *);
187 STATIC doublereal detrat=0, crvmin=0;
188 STATIC integer nftest=0;
189 STATIC doublereal distsq=0;
190 extern /* Subroutine */ int trsapp_(integer *, integer *, doublereal *,
191 doublereal *, doublereal *, doublereal *, doublereal *,
192 doublereal *, doublereal *, doublereal *, doublereal *,
193 doublereal *, doublereal *, doublereal *);
194 STATIC doublereal xoptsq=0;
195 STATIC integer nfm=0;
196 STATIC doublereal one=0;
197 STATIC integer idz=0;
198 STATIC doublereal dsq=0, rho=0;
199 STATIC integer ipt=0, jpt=0;
200 STATIC doublereal sum=0;
201
202
203 /* The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are ident
204 ical*/
205 /* to the corresponding arguments in SUBROUTINE NEWUOA. */
206 /* XBASE will hold a shift of origin that should reduce the contributio
207 ns*/
208 /* from rounding errors to values of the model and Lagrange functions
209 .*/
210 /* XOPT will be set to the displacement from XBASE of the vector of */
211 /* variables that provides the least calculated F so far. */
212 /* XNEW will be set to the displacement from XBASE of the vector of */
213 /* variables for the current calculation of F. */
214 /* XPT will contain the interpolation point coordinates relative to XBA
215 SE.*/
216 /* FVAL will hold the values of F at the interpolation points. */
217 /* GQ will hold the gradient of the quadratic model at XBASE. */
218 /* HQ will hold the explicit second derivatives of the quadratic model.
219 */
220 /* PQ will contain the parameters of the implicit second derivatives of
221 */
222 /* the quadratic model. */
223 /* BMAT will hold the last N columns of H. */
224 /* ZMAT will hold the factorization of the leading NPT by NPT submatrix
225 of*/
226 /* H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, wher
227 e*/
228 /* the elements of DZ are plus or minus one, as specified by IDZ. */
229 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
230 /* D is reserved for trial steps from XOPT. */
231 /* VLAG will contain the values of the Lagrange functions at a new poin
232 t X.*/
233 /* They are part of a product that requires VLAG to be of length NDIM
234 .*/
235 /* The array W will be used for working space. Its length must be at le
236 ast*/
237 /* 10*NDIM = 10*(NPT+N). */
238
239 /* Set some constants. */
240
241 /* Parameter adjustments */
242 zmat_dim1 = *npt;
243 zmat_offset = zmat_dim1 + 1;
244 zmat -= zmat_offset;
245 xpt_dim1 = *npt;
246 xpt_offset = xpt_dim1 + 1;
247 xpt -= xpt_offset;
248 --x;
249 --xbase;
250 --xopt;
251 --xnew;
252 --fval;
253 --gq;
254 --hq;
255 --pq;
256 bmat_dim1 = *ndim;
257 bmat_offset = bmat_dim1 + 1;
258 bmat -= bmat_offset;
259 --d__;
260 --vlag;
261 --w;
262
263 /* Function Body */
264 half = .5;
265 one = 1.;
266 tenth = .1;
267 zero = 0.;
268 np = *n + 1;
269 nh = *n * np / 2;
270 nptm = *npt - np;
271 nftest = AFmax(*maxfun,1);
272
273 /* Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. */
274
275 i__1 = *n;
276 for (j = 1; j <= i__1; ++j) {
277 xbase[j] = x[j];
278 i__2 = *npt;
279 for (k = 1; k <= i__2; ++k) {
280 /* L10: */
281 xpt[k + j * xpt_dim1] = zero;
282 }
283 i__2 = *ndim;
284 for (i__ = 1; i__ <= i__2; ++i__) {
285 /* L20: */
286 bmat[i__ + j * bmat_dim1] = zero;
287 }
288 }
289 i__2 = nh;
290 for (ih = 1; ih <= i__2; ++ih) {
291 /* L30: */
292 hq[ih] = zero;
293 }
294 i__2 = *npt;
295 for (k = 1; k <= i__2; ++k) {
296 pq[k] = zero;
297 i__1 = nptm;
298 for (j = 1; j <= i__1; ++j) {
299 /* L40: */
300 zmat[k + j * zmat_dim1] = zero;
301 }
302 }
303
304 /* Begin the initialization procedure. NF becomes one more than the num
305 ber*/
306 /* of function values so far. The coordinates of the displacement of th
307 e*/
308 /* next initial interpolation point from XBASE are set in XPT(NF,.).
309 */
310
311 rhosq = *rhobeg * *rhobeg;
312 recip = one / rhosq;
313 reciq = sqrt(half) / rhosq;
314 nf = 0;
315 *icode = 0;
316 L50:
317 nfm = nf;
318 nfmm = nf - *n;
319 ++nf;
320 if (nfm <= *n << 1) {
321 if (nfm >= 1 && nfm <= *n) {
322 xpt[nf + nfm * xpt_dim1] = *rhobeg;
323 } else if (nfm > *n) {
324 xpt[nf + nfmm * xpt_dim1] = -(*rhobeg);
325 }
326 } else {
327 itemp = (nfmm - 1) / *n;
328 jpt = nfm - itemp * *n - *n;
329 ipt = jpt + itemp;
330 if (ipt > *n) {
331 itemp = jpt;
332 jpt = ipt - *n;
333 ipt = itemp;
334 }
335 xipt = *rhobeg;
336 if (fval[ipt + np] < fval[ipt + 1]) {
337 xipt = -xipt;
338 }
339 xjpt = *rhobeg;
340 if (fval[jpt + np] < fval[jpt + 1]) {
341 xjpt = -xjpt;
342 }
343 xpt[nf + ipt * xpt_dim1] = xipt;
344 xpt[nf + jpt * xpt_dim1] = xjpt;
345 }
346
347 /* Calculate the next value of F, label 70 being reached immediately
348 */
349 /* after this calculation. The least function value so far and its inde
350 x*/
351 /* are required. */
352
353 i__1 = *n;
354 for (j = 1; j <= i__1; ++j) {
355 /* L60: */
356 x[j] = xpt[nf + j * xpt_dim1] + xbase[j];
357 }
358 goto L310;
359 L70:
360 fval[nf] = f;
361 if (nf == 1) {
362 fbeg = f;
363 fopt = f;
364 if(DEBUG) fprintf(stderr,"\nInitialize fopt=%.14g ;\n",fopt) ;
365 kopt = 1;
366 } else if (f < fopt) {
367 fopt = f;
368 if(DEBUG) fprintf(stderr," fopt=%.14g\n",fopt) ;
369 kopt = nf;
370 }
371
372 /* Set the nonzero initial elements of BMAT and the quadratic model in
373 */
374 /* the cases when NF is at most 2*N+1. */
375
376 if (nfm <= *n << 1) {
377 if (nfm >= 1 && nfm <= *n) {
378 gq[nfm] = (f - fbeg) / *rhobeg;
379 if (*npt < nf + *n) {
380 bmat[nfm * bmat_dim1 + 1] = -one / *rhobeg;
381 bmat[nf + nfm * bmat_dim1] = one / *rhobeg;
382 bmat[*npt + nfm + nfm * bmat_dim1] = -half * rhosq;
383 }
384 } else if (nfm > *n) {
385 bmat[nf - *n + nfmm * bmat_dim1] = half / *rhobeg;
386 bmat[nf + nfmm * bmat_dim1] = -half / *rhobeg;
387 zmat[nfmm * zmat_dim1 + 1] = -reciq - reciq;
388 zmat[nf - *n + nfmm * zmat_dim1] = reciq;
389 zmat[nf + nfmm * zmat_dim1] = reciq;
390 ih = nfmm * (nfmm + 1) / 2;
391 temp = (fbeg - f) / *rhobeg;
392 hq[ih] = (gq[nfmm] - temp) / *rhobeg;
393 gq[nfmm] = half * (gq[nfmm] + temp);
394 }
395
396 /* Set the off-diagonal second derivatives of the Lagrange function
397 s and*/
398 /* the initial quadratic model. */
399
400 } else {
401 ih = ipt * (ipt - 1) / 2 + jpt;
402 if (xipt < zero) {
403 ipt += *n;
404 }
405 if (xjpt < zero) {
406 jpt += *n;
407 }
408 zmat[nfmm * zmat_dim1 + 1] = recip;
409 zmat[nf + nfmm * zmat_dim1] = recip;
410 zmat[ipt + 1 + nfmm * zmat_dim1] = -recip;
411 zmat[jpt + 1 + nfmm * zmat_dim1] = -recip;
412 hq[ih] = (fbeg - fval[ipt + 1] - fval[jpt + 1] + f) / (xipt * xjpt);
413 }
414 if (nf < *npt) {
415 goto L50;
416 }
417
418 /* Begin the iterative procedure, because the initial model is complete
419 .*/
420
421 rho = *rhobeg;
422 delta = rho;
423 if(DEBUG) fprintf(stderr," rho=delta=%.14g\n",delta) ;
424 idz = 1;
425 diffa = zero;
426 diffb = zero;
427 itest = 0;
428 xoptsq = zero;
429 i__1 = *n;
430 for (i__ = 1; i__ <= i__1; ++i__) {
431 xopt[i__] = xpt[kopt + i__ * xpt_dim1];
432 /* L80: */
433 /* Computing 2nd power */
434 d__1 = xopt[i__];
435 xoptsq += d__1 * d__1;
436 }
437 L90:
438 nfsav = nf;
439
440 /* Generate the next trust region step and test its length. Set KNEW
441 */
442 /* to -1 if the purpose of the next F will be to improve the model. */
443
444 L100:
445 knew = 0;
446 trsapp_(n, npt, &xopt[1], &xpt[xpt_offset], &gq[1], &hq[1], &pq[1], &
447 delta, &d__[1], &w[1], &w[np], &w[np + *n], &w[np + (*n << 1)], &
448 crvmin);
449 dsq = zero;
450 i__1 = *n;
451 for (i__ = 1; i__ <= i__1; ++i__) {
452 /* L110: */
453 /* Computing 2nd power */
454 d__1 = d__[i__];
455 dsq += d__1 * d__1;
456 }
457 /* Computing MIN */
458 d__1 = delta, d__2 = sqrt(dsq);
459 if(0&&DEBUG) fprintf(stderr," sqrt(dsq=%.14g)=%.14g\n",dsq,d__2) ;
460 dnorm = AFmin(d__1,d__2);
461 if (dnorm < half * rho) {
462 knew = -1;
463 delta = tenth * delta;
464 if(DEBUG) fprintf(stderr," delta/10") ;
465 ratio = -1.;
466 if (delta <= rho * 1.5) {
467 if(DEBUG) fprintf(stderr," delta=rho") ;
468 delta = rho;
469 }
470 if(DEBUG) fprintf(stderr," delta=%.14g\n",delta) ;
471 if (nf <= nfsav + 2) {
472 goto L460;
473 }
474 temp = crvmin * .125 * rho * rho;
475 /* Computing MAX */
476 d__1 = AFmax(diffa,diffb);
477 if (temp <= AFmax(d__1,diffc)) {
478 goto L460;
479 }
480 goto L490;
481 }
482
483 /* Shift XBASE if XOPT may be too far from XBASE. First make the change
484 s*/
485 /* to BMAT that do not depend on ZMAT. */
486
487 L120:
488 if (dsq <= xoptsq * .001) {
489 tempq = xoptsq * .25;
490 i__1 = *npt;
491 for (k = 1; k <= i__1; ++k) {
492 sum = zero;
493 i__2 = *n;
494 for (i__ = 1; i__ <= i__2; ++i__) {
495 /* L130: */
496 sum += xpt[k + i__ * xpt_dim1] * xopt[i__];
497 }
498 temp = pq[k] * sum;
499 sum -= half * xoptsq;
500 w[*npt + k] = sum;
501 i__2 = *n;
502 for (i__ = 1; i__ <= i__2; ++i__) {
503 gq[i__] += temp * xpt[k + i__ * xpt_dim1];
504 xpt[k + i__ * xpt_dim1] -= half * xopt[i__];
505 vlag[i__] = bmat[k + i__ * bmat_dim1];
506 w[i__] = sum * xpt[k + i__ * xpt_dim1] + tempq * xopt[i__];
507 ip = *npt + i__;
508 i__3 = i__;
509 for (j = 1; j <= i__3; ++j) {
510 /* L140: */
511 bmat[ip + j * bmat_dim1] = bmat[ip + j * bmat_dim1] +
512 vlag[i__] * w[j] + w[i__] * vlag[j];
513 }
514 }
515 }
516
517 /* Then the revisions of BMAT that depend on ZMAT are calculated.
518 */
519
520 i__3 = nptm;
521 for (k = 1; k <= i__3; ++k) {
522 sumz = zero;
523 i__2 = *npt;
524 for (i__ = 1; i__ <= i__2; ++i__) {
525 sumz += zmat[i__ + k * zmat_dim1];
526 /* L150: */
527 w[i__] = w[*npt + i__] * zmat[i__ + k * zmat_dim1];
528 }
529 i__2 = *n;
530 for (j = 1; j <= i__2; ++j) {
531 sum = tempq * sumz * xopt[j];
532 i__1 = *npt;
533 for (i__ = 1; i__ <= i__1; ++i__) {
534 /* L160: */
535 sum += w[i__] * xpt[i__ + j * xpt_dim1];
536 }
537 vlag[j] = sum;
538 if (k < idz) {
539 sum = -sum;
540 }
541 i__1 = *npt;
542 for (i__ = 1; i__ <= i__1; ++i__) {
543 /* L170: */
544 bmat[i__ + j * bmat_dim1] += sum * zmat[i__ + k *
545 zmat_dim1];
546 }
547 }
548 i__1 = *n;
549 for (i__ = 1; i__ <= i__1; ++i__) {
550 ip = i__ + *npt;
551 temp = vlag[i__];
552 if (k < idz) {
553 temp = -temp;
554 }
555 i__2 = i__;
556 for (j = 1; j <= i__2; ++j) {
557 /* L180: */
558 bmat[ip + j * bmat_dim1] += temp * vlag[j];
559 }
560 }
561 }
562
563 /* The following instructions complete the shift of XBASE, includi
564 ng */
565 /* the changes to the parameters of the quadratic model. */
566
567 ih = 0;
568 i__2 = *n;
569 for (j = 1; j <= i__2; ++j) {
570 w[j] = zero;
571 i__1 = *npt;
572 for (k = 1; k <= i__1; ++k) {
573 w[j] += pq[k] * xpt[k + j * xpt_dim1];
574 /* L190: */
575 xpt[k + j * xpt_dim1] -= half * xopt[j];
576 }
577 i__1 = j;
578 for (i__ = 1; i__ <= i__1; ++i__) {
579 ++ih;
580 if (i__ < j) {
581 gq[j] += hq[ih] * xopt[i__];
582 }
583 gq[i__] += hq[ih] * xopt[j];
584 hq[ih] = hq[ih] + w[i__] * xopt[j] + xopt[i__] * w[j];
585 /* L200: */
586 bmat[*npt + i__ + j * bmat_dim1] = bmat[*npt + j + i__ *
587 bmat_dim1];
588 }
589 }
590 i__1 = *n;
591 for (j = 1; j <= i__1; ++j) {
592 xbase[j] += xopt[j];
593 /* L210: */
594 xopt[j] = zero;
595 }
596 xoptsq = zero;
597 }
598
599 /* Pick the model step if KNEW is positive. A different choice of D */
600 /* may be made later, if the choice of D by BIGLAG causes substantial
601 */
602 /* cancellation in DENOM. */
603
604 if (knew > 0) {
605 biglag_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &zmat[
606 zmat_offset], &idz, ndim, &knew, &dstep, &d__[1], &alpha, &
607 vlag[1], &vlag[*npt + 1], &w[1], &w[np], &w[np + *n]);
608 }
609
610 /* Calculate VLAG and BETA for the current choice of D. The first NPT
611 */
612 /* components of W_check will be held in W. */
613
614 i__1 = *npt;
615 for (k = 1; k <= i__1; ++k) {
616 suma = zero;
617 sumb = zero;
618 sum = zero;
619 i__2 = *n;
620 for (j = 1; j <= i__2; ++j) {
621 suma += xpt[k + j * xpt_dim1] * d__[j];
622 sumb += xpt[k + j * xpt_dim1] * xopt[j];
623 /* L220: */
624 sum += bmat[k + j * bmat_dim1] * d__[j];
625 }
626 w[k] = suma * (half * suma + sumb);
627 /* L230: */
628 vlag[k] = sum;
629 }
630 beta = zero;
631 i__1 = nptm;
632 for (k = 1; k <= i__1; ++k) {
633 sum = zero;
634 i__2 = *npt;
635 for (i__ = 1; i__ <= i__2; ++i__) {
636 /* L240: */
637 sum += zmat[i__ + k * zmat_dim1] * w[i__];
638 }
639 if (k < idz) {
640 beta += sum * sum;
641 sum = -sum;
642 } else {
643 beta -= sum * sum;
644 }
645 i__2 = *npt;
646 for (i__ = 1; i__ <= i__2; ++i__) {
647 /* L250: */
648 vlag[i__] += sum * zmat[i__ + k * zmat_dim1];
649 }
650 }
651 bsum = zero;
652 dx = zero;
653 i__2 = *n;
654 for (j = 1; j <= i__2; ++j) {
655 sum = zero;
656 i__1 = *npt;
657 for (i__ = 1; i__ <= i__1; ++i__) {
658 /* L260: */
659 sum += w[i__] * bmat[i__ + j * bmat_dim1];
660 }
661 bsum += sum * d__[j];
662 jp = *npt + j;
663 i__1 = *n;
664 for (k = 1; k <= i__1; ++k) {
665 /* L270: */
666 sum += bmat[jp + k * bmat_dim1] * d__[k];
667 }
668 vlag[jp] = sum;
669 bsum += sum * d__[j];
670 /* L280: */
671 dx += d__[j] * xopt[j];
672 }
673 beta = dx * dx + dsq * (xoptsq + dx + dx + half * dsq) + beta - bsum;
674 vlag[kopt] += one;
675
676 /* If KNEW is positive and if the cancellation in DENOM is unacceptable
677 ,*/
678 /* then BIGDEN calculates an alternative model step, XNEW being used fo
679 r*/
680 /* working space. */
681
682 if (knew > 0) {
683 /* Computing 2nd power */
684 d__1 = vlag[knew];
685 temp = one + alpha * beta / (d__1 * d__1);
686 if (abs(temp) <= .8) {
687 bigden_(n, npt, &xopt[1], &xpt[xpt_offset], &bmat[bmat_offset], &
688 zmat[zmat_offset], &idz, ndim, &kopt, &knew, &d__[1], &w[
689 1], &vlag[1], &beta, &xnew[1], &w[*ndim + 1], &w[*ndim *
690 6 + 1]);
691 }
692 }
693
694 /* Calculate the next value of the objective function. */
695
696 L290:
697 i__2 = *n;
698 for (i__ = 1; i__ <= i__2; ++i__) {
699 xnew[i__] = xopt[i__] + d__[i__];
700 /* L300: */
701 x[i__] = xbase[i__] + xnew[i__];
702 }
703 ++nf;
704 L310:
705 if (nf > nftest) {
706 --nf;
707 /* CC IF (IPRINT .GT. 0) PRINT 320 */
708 /* CC 320 FORMAT (/4X,'Return from NEWUOA because CALFUN has been
709 ', */
710 /* CC 1 ' called MAXFUN times.') */
711 goto L530;
712 }
713 calfun_(n, &x[1], &f);
714 /* CC IF (IPRINT .EQ. 3) THEN */
715 /* CC PRINT 330, NF,F,(X(I),I=1,N) */
716 /* CC 330 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, */
717 /* CC 1 ' The corresponding X is:'/(2X,5D15.6)) */
718 /* CC END IF */
719 if (nf <= *npt) {
720 goto L70;
721 }
722 if (knew == -1) {
723 goto L530;
724 }
725
726 /* Use the quadratic model to predict the change in F due to the step D
727 ,*/
728 /* and set DIFF to the error of this prediction. */
729
730 vquad = zero;
731 ih = 0;
732 i__2 = *n;
733 for (j = 1; j <= i__2; ++j) {
734 vquad += d__[j] * gq[j];
735 i__1 = j;
736 for (i__ = 1; i__ <= i__1; ++i__) {
737 ++ih;
738 temp = d__[i__] * xnew[j] + d__[j] * xopt[i__];
739 if (i__ == j) {
740 temp = half * temp;
741 }
742 /* L340: */
743 vquad += temp * hq[ih];
744 }
745 }
746 i__1 = *npt;
747 for (k = 1; k <= i__1; ++k) {
748 /* L350: */
749 vquad += pq[k] * w[k];
750 }
751 diff = f - fopt - vquad;
752 diffc = diffb;
753 diffb = diffa;
754 diffa = abs(diff);
755 if (dnorm > rho) {
756 nfsav = nf;
757 }
758
759 /* Update FOPT and XOPT if the new F is the least value of the objectiv
760 e*/
761 /* function so far. The branch when KNEW is positive occurs if D is not
762 */
763 /* a trust region step. */
764
765 fsave = fopt;
766 if (f < fopt) {
767 fopt = f;
768 if(DEBUG) fprintf(stderr," fopt=%.14g\n",fopt) ;
769 xoptsq = zero;
770 i__1 = *n;
771 for (i__ = 1; i__ <= i__1; ++i__) {
772 xopt[i__] = xnew[i__];
773 /* L360: */
774 /* Computing 2nd power */
775 d__1 = xopt[i__];
776 xoptsq += d__1 * d__1;
777 }
778 }
779 ksave = knew;
780 if (knew > 0) {
781 goto L410;
782 }
783
784 /* Pick the next value of DELTA after a trust region step. */
785
786 if (vquad >= zero) {
787 /* CC IF (IPRINT .GT. 0) PRINT 370 */
788 /* CC 370 FORMAT (/4X,'Return from NEWUOA because a trust', */
789 /* CC 1 ' region step has failed to reduce Q.') */
790 *icode = -1;
791 goto L530;
792 }
793 ratio = (f - fsave) / vquad;
794 if (ratio <= tenth) {
795 if(DEBUG) fprintf(stderr," delta=dnorm/2") ;
796 delta = half * dnorm;
797 } else if (ratio <= .7) {
798 /* Computing MAX */
799 if(DEBUG) fprintf(stderr," delta=max(delta/2,dnorm)") ;
800 d__1 = half * delta;
801 delta = AFmax(d__1,dnorm);
802 } else {
803 /* Computing MAX */
804 d__1 = half * delta, d__2 = dnorm + dnorm;
805 delta = AFmax(d__1,d__2);
806 if(DEBUG) fprintf(stderr," delta=max(delta/2,dnorm*2)") ;
807 }
808 if (delta <= rho * 1.5) {
809 if(DEBUG) fprintf(stderr," delta=RHO") ;
810 delta = rho;
811 }
812 if(DEBUG) fprintf(stderr," delta=%.14g\n",delta) ;
813
814 /* Set KNEW to the index of the next interpolation point to be deleted.
815 */
816
817 /* Computing MAX */
818 d__2 = tenth * delta;
819 /* Computing 2nd power */
820 d__1 = AFmax(d__2,rho);
821 rhosq = d__1 * d__1;
822 ktemp = 0;
823 detrat = zero;
824 if (f >= fsave) {
825 ktemp = kopt;
826 detrat = one;
827 }
828 i__1 = *npt;
829 for (k = 1; k <= i__1; ++k) {
830 hdiag = zero;
831 i__2 = nptm;
832 for (j = 1; j <= i__2; ++j) {
833 temp = one;
834 if (j < idz) {
835 temp = -one;
836 }
837 /* L380: */
838 /* Computing 2nd power */
839 d__1 = zmat[k + j * zmat_dim1];
840 hdiag += temp * (d__1 * d__1);
841 }
842 /* Computing 2nd power */
843 d__2 = vlag[k];
844 temp = (d__1 = beta * hdiag + d__2 * d__2, abs(d__1));
845 distsq = zero;
846 i__2 = *n;
847 for (j = 1; j <= i__2; ++j) {
848 /* L390: */
849 /* Computing 2nd power */
850 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
851 distsq += d__1 * d__1;
852 }
853 if (distsq > rhosq) {
854 /* Computing 3rd power */
855 d__1 = distsq / rhosq, d__2 = d__1;
856 temp *= d__2 * (d__1 * d__1);
857 }
858 if (temp > detrat && k != ktemp) {
859 detrat = temp;
860 knew = k;
861 }
862 /* L400: */
863 }
864 if (knew == 0) {
865 goto L460;
866 }
867
868 /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
869 */
870 /* can be moved. Begin the updating of the quadratic model, starting
871 */
872 /* with the explicit second derivative term. */
873
874 L410:
875 update_(n, npt, &bmat[bmat_offset], &zmat[zmat_offset], &idz, ndim, &vlag[
876 1], &beta, &knew, &w[1]);
877 fval[knew] = f;
878 ih = 0;
879 i__1 = *n;
880 for (i__ = 1; i__ <= i__1; ++i__) {
881 temp = pq[knew] * xpt[knew + i__ * xpt_dim1];
882 i__2 = i__;
883 for (j = 1; j <= i__2; ++j) {
884 ++ih;
885 /* L420: */
886 hq[ih] += temp * xpt[knew + j * xpt_dim1];
887 }
888 }
889 pq[knew] = zero;
890
891 /* Update the other second derivative parameters, and then the gradient
892 */
893 /* vector of the model. Also include the new interpolation point. */
894
895 i__2 = nptm;
896 for (j = 1; j <= i__2; ++j) {
897 temp = diff * zmat[knew + j * zmat_dim1];
898 if (j < idz) {
899 temp = -temp;
900 }
901 i__1 = *npt;
902 for (k = 1; k <= i__1; ++k) {
903 /* L440: */
904 pq[k] += temp * zmat[k + j * zmat_dim1];
905 }
906 }
907 gqsq = zero;
908 i__1 = *n;
909 for (i__ = 1; i__ <= i__1; ++i__) {
910 gq[i__] += diff * bmat[knew + i__ * bmat_dim1];
911 /* Computing 2nd power */
912 d__1 = gq[i__];
913 gqsq += d__1 * d__1;
914 /* L450: */
915 xpt[knew + i__ * xpt_dim1] = xnew[i__];
916 }
917
918 /* If a trust region step makes a small change to the objective functio
919 n,*/
920 /* then calculate the gradient of the least Frobenius norm interpolant
921 at*/
922 /* XBASE, and store it in W, using VLAG for a vector of right hand side
923 s.*/
924
925 if (ksave == 0 && delta == rho) {
926 if (abs(ratio) > .01) {
927 itest = 0;
928 } else {
929 i__1 = *npt;
930 for (k = 1; k <= i__1; ++k) {
931 /* L700: */
932 vlag[k] = fval[k] - fval[kopt];
933 }
934 gisq = zero;
935 i__1 = *n;
936 for (i__ = 1; i__ <= i__1; ++i__) {
937 sum = zero;
938 i__2 = *npt;
939 for (k = 1; k <= i__2; ++k) {
940 /* L710: */
941 sum += bmat[k + i__ * bmat_dim1] * vlag[k];
942 }
943 gisq += sum * sum;
944 /* L720: */
945 w[i__] = sum;
946 }
947
948 /* Test whether to replace the new quadratic model by the least
949 Frobenius*/
950 /* norm interpolant, making the replacement if the test is sat
951 isfied. */
952
953 ++itest;
954 if (gqsq < gisq * 100.) {
955 itest = 0;
956 }
957 if (itest >= 3) {
958 i__1 = *n;
959 for (i__ = 1; i__ <= i__1; ++i__) {
960 /* L730: */
961 gq[i__] = w[i__];
962 }
963 i__1 = nh;
964 for (ih = 1; ih <= i__1; ++ih) {
965 /* L740: */
966 hq[ih] = zero;
967 }
968 i__1 = nptm;
969 for (j = 1; j <= i__1; ++j) {
970 w[j] = zero;
971 i__2 = *npt;
972 for (k = 1; k <= i__2; ++k) {
973 /* L750: */
974 w[j] += vlag[k] * zmat[k + j * zmat_dim1];
975 }
976 /* L760: */
977 if (j < idz) {
978 w[j] = -w[j];
979 }
980 }
981 i__1 = *npt;
982 for (k = 1; k <= i__1; ++k) {
983 pq[k] = zero;
984 i__2 = nptm;
985 for (j = 1; j <= i__2; ++j) {
986 /* L770: */
987 pq[k] += zmat[k + j * zmat_dim1] * w[j];
988 }
989 }
990 itest = 0;
991 }
992 }
993 }
994 if (f < fsave) {
995 kopt = knew;
996 }
997
998 /* If a trust region step has provided a sufficient decrease in F, then
999 */
1000 /* branch for another trust region calculation. The case KSAVE>0 occurs
1001 */
1002 /* when the new function value was calculated by a model step. */
1003
1004 if (f <= fsave + tenth * vquad) {
1005 goto L100;
1006 }
1007 if (ksave > 0) {
1008 goto L100;
1009 }
1010
1011 /* Alternatively, find out if the interpolation points are close enough
1012 */
1013 /* to the best point so far. */
1014
1015 knew = 0;
1016 L460:
1017 distsq = delta * 4. * delta;
1018 i__2 = *npt;
1019 for (k = 1; k <= i__2; ++k) {
1020 sum = zero;
1021 i__1 = *n;
1022 for (j = 1; j <= i__1; ++j) {
1023 /* L470: */
1024 /* Computing 2nd power */
1025 d__1 = xpt[k + j * xpt_dim1] - xopt[j];
1026 sum += d__1 * d__1;
1027 }
1028 if (sum > distsq) {
1029 knew = k;
1030 distsq = sum;
1031 }
1032 /* L480: */
1033 }
1034
1035 /* If KNEW is positive, then set DSTEP, and branch back for the next
1036 */
1037 /* iteration, which will generate a "model step". */
1038
1039 if (knew > 0) {
1040 /* Computing MAX */
1041 /* Computing MIN */
1042 d__2 = tenth * sqrt(distsq), d__3 = half * delta;
1043 if(0&&DEBUG) fprintf(stderr," sqrt(distsq=%.14g)=%.14g\n",distsq,d__2*10.0) ;
1044 d__1 = AFmin(d__2,d__3);
1045 dstep = AFmax(d__1,rho);
1046 dsq = dstep * dstep;
1047 goto L120;
1048 }
1049 if (ratio > zero) {
1050 goto L100;
1051 }
1052 if (AFmax(delta,dnorm) > rho) {
1053 goto L100;
1054 }
1055
1056 /* The calculations with the current value of RHO are complete. Pick th
1057 e*/
1058 /* next values of RHO and DELTA. */
1059
1060 L490:
1061 if (rho > *rhoend) {
1062 delta = half * rho;
1063 ratio = rho / *rhoend;
1064 if (ratio <= 16.) {
1065 rho = *rhoend;
1066 if(DEBUG) fprintf(stderr," ratio<16") ;
1067 } else if (ratio <= 250.) {
1068 if(DEBUG) fprintf(stderr," ratio<250") ;
1069 rho = sqrt(ratio) * *rhoend;
1070 } else {
1071 if(DEBUG) fprintf(stderr," ratio=other") ;
1072 rho = tenth * rho;
1073 }
1074 delta = AFmax(delta,rho);
1075 if(DEBUG) fprintf(stderr," rho=%.14g delta=%.14g\n",rho,delta) ;
1076 /* CC IF (IPRINT .GE. 2) THEN */
1077 /* CC IF (IPRINT .GE. 3) PRINT 500 */
1078 /* CC 500 FORMAT (5X) */
1079 /* CC PRINT 510, RHO,NF */
1080 /* CC 510 FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of', */
1081 /* CC 1 ' function values =',I6) */
1082 /* CC PRINT 520, FOPT,(XBASE(I)+XOPT(I),I=1,N) */
1083 /* CC 520 FORMAT (4X,'Least value of F =',1PD23.15,9X, */
1084 /* CC 1 'The corresponding X is:'/(2X,5D15.6)) */
1085 /* CC END IF */
1086 goto L90;
1087 }
1088
1089 /* Return from the calculation, after another Newton-Raphson step, if
1090 */
1091 /* it is too short to have been tried before. */
1092
1093 if (knew == -1) {
1094 goto L290;
1095 }
1096 L530:
1097 if (fopt <= f) {
1098 i__2 = *n;
1099 for (i__ = 1; i__ <= i__2; ++i__) {
1100 /* L540: */
1101 x[i__] = xbase[i__] + xopt[i__];
1102 }
1103 f = fopt;
1104 if(DEBUG) fprintf(stderr," ; final set f=%.14g\n",f) ;
1105 }
1106 /* CC IF (IPRINT .GE. 1) THEN */
1107 /* CC PRINT 550, NF */
1108 /* CC 550 FORMAT (/4X,'At the return from NEWUOA',5X, */
1109 /* CC 1 'Number of function values =',I6) */
1110 /* CC PRINT 520, F,(X(I),I=1,N) */
1111 /* CC END IF */
1112 if (*icode == 0) {
1113 *icode = nf;
1114 }
1115 if(DEBUG) fprintf(stderr," ; return f=%.14g with rho=%.14g\n",f,rho) ;
1116 return 0;
1117 } /* newuob_ */
1118
trsapp_(integer * n,integer * npt,doublereal * xopt,doublereal * xpt,doublereal * gq,doublereal * hq,doublereal * pq,doublereal * delta,doublereal * step,doublereal * d__,doublereal * g,doublereal * hd,doublereal * hs,doublereal * crvmin)1119 /* Subroutine */ int trsapp_(integer *n, integer *npt, doublereal *xopt,
1120 doublereal *xpt, doublereal *gq, doublereal *hq, doublereal *pq,
1121 doublereal *delta, doublereal *step, doublereal *d__, doublereal *g,
1122 doublereal *hd, doublereal *hs, doublereal *crvmin)
1123 {
1124 /* System generated locals */
1125 integer xpt_dim1=0, xpt_offset=0, i__1=0, i__2=0;
1126 doublereal d__1=0, d__2=0;
1127
1128 /* Builtin functions */
1129 double atan(doublereal), sqrt(doublereal), cos(doublereal), sin(
1130 doublereal);
1131
1132 /* Local variables */
1133 STATIC doublereal qadd=0, half=0, qbeg=0, qred=0, qmin=0, temp=0, qsav=0, qnew=0, zero=0;
1134 STATIC integer i__=0, j=0;
1135 STATIC doublereal ggbeg=0;
1136 STATIC integer k=0;
1137 STATIC doublereal alpha=0, angle=0, reduc=0;
1138 STATIC integer iterc=0;
1139 STATIC doublereal ggsav=0, delsq=0, tempa=0, tempb=0;
1140 STATIC integer isave=0;
1141 STATIC doublereal bstep=0, ratio=0, twopi=0, dd=0, cf=0, dg=0, gg=0;
1142 STATIC integer ih=0;
1143 STATIC doublereal ds=0, sg=0;
1144 STATIC integer iu=0;
1145 STATIC doublereal ss=0;
1146 STATIC integer itersw=0;
1147 STATIC doublereal dhd=0, dhs=0, cth=0, sgk=0, shs=0, sth=0, angtest=0;
1148 STATIC integer itermax=0;
1149
1150
1151 /* N is the number of variables of a quadratic objective function, Q sa
1152 y.*/
1153 /* The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meaning
1154 s,*/
1155 /* in order to define the current quadratic model Q. */
1156 /* DELTA is the trust region radius, and has to be positive. */
1157 /* STEP will be set to the calculated trial step. */
1158 /* The arrays D, G, HD and HS will be used for working space. */
1159 /* CRVMIN will be set to the least curvature of H along the conjugate
1160 */
1161 /* directions that occur, except that it is set to zero if STEP goes
1162 */
1163 /* all the way to the trust region boundary. */
1164
1165 /* The calculation of STEP begins with the truncated conjugate gradient
1166 */
1167 /* method. If the boundary of the trust region is reached, then further
1168 */
1169 /* changes to STEP may be made, each one being in the 2D space spanned
1170 */
1171 /* by the current STEP and the corresponding gradient of Q. Thus STEP
1172 */
1173 /* should provide a substantial reduction to Q within the trust region.
1174 */
1175
1176 /* Initialization, which includes setting HD to H times XOPT. */
1177
1178 /* Parameter adjustments */
1179 xpt_dim1 = *npt;
1180 xpt_offset = xpt_dim1 + 1;
1181 xpt -= xpt_offset;
1182 --xopt;
1183 --gq;
1184 --hq;
1185 --pq;
1186 --step;
1187 --d__;
1188 --g;
1189 --hd;
1190 --hs;
1191
1192 /* Function Body */
1193 half = .5;
1194 zero = 0.;
1195 twopi = atan(1.) * 8.;
1196 delsq = *delta * *delta;
1197 iterc = 0;
1198 itermax = *n;
1199 itersw = itermax;
1200 i__1 = *n;
1201 for (i__ = 1; i__ <= i__1; ++i__) {
1202 /* L10: */
1203 d__[i__] = xopt[i__];
1204 }
1205 goto L170;
1206
1207 /* Prepare for the first line search. */
1208
1209 L20:
1210 qred = zero;
1211 dd = zero;
1212 i__1 = *n;
1213 for (i__ = 1; i__ <= i__1; ++i__) {
1214 step[i__] = zero;
1215 hs[i__] = zero;
1216 g[i__] = gq[i__] + hd[i__];
1217 d__[i__] = -g[i__];
1218 /* L30: */
1219 /* Computing 2nd power */
1220 d__1 = d__[i__];
1221 dd += d__1 * d__1;
1222 }
1223 *crvmin = zero;
1224 if (dd == zero) {
1225 goto L160;
1226 }
1227 ds = zero;
1228 ss = zero;
1229 gg = dd;
1230 ggbeg = gg;
1231
1232 /* Calculate the step to the trust region boundary and the product HD.
1233 */
1234
1235 L40:
1236 ++iterc;
1237 temp = delsq - ss;
1238 bstep = temp / (ds + sqrt(ds * ds + dd * temp));
1239 goto L170;
1240 L50:
1241 dhd = zero;
1242 i__1 = *n;
1243 for (j = 1; j <= i__1; ++j) {
1244 /* L60: */
1245 dhd += d__[j] * hd[j];
1246 }
1247
1248 /* Update CRVMIN and set the step-length ALPHA. */
1249
1250 alpha = bstep;
1251 if (dhd > zero) {
1252 temp = dhd / dd;
1253 if (iterc == 1) {
1254 *crvmin = temp;
1255 }
1256 *crvmin = AFmin(*crvmin,temp);
1257 /* Computing MIN */
1258 d__1 = alpha, d__2 = gg / dhd;
1259 alpha = AFmin(d__1,d__2);
1260 }
1261 qadd = alpha * (gg - half * alpha * dhd);
1262 qred += qadd;
1263
1264 /* Update STEP and HS. */
1265
1266 ggsav = gg;
1267 gg = zero;
1268 i__1 = *n;
1269 for (i__ = 1; i__ <= i__1; ++i__) {
1270 step[i__] += alpha * d__[i__];
1271 hs[i__] += alpha * hd[i__];
1272 /* L70: */
1273 /* Computing 2nd power */
1274 d__1 = g[i__] + hs[i__];
1275 gg += d__1 * d__1;
1276 }
1277
1278 /* Begin another conjugate direction iteration if required. */
1279
1280 if (alpha < bstep) {
1281 if (qadd <= qred * .01) {
1282 goto L160;
1283 }
1284 if (gg <= ggbeg * 1e-4) {
1285 goto L160;
1286 }
1287 if (iterc == itermax) {
1288 goto L160;
1289 }
1290 temp = gg / ggsav;
1291 dd = zero;
1292 ds = zero;
1293 ss = zero;
1294 i__1 = *n;
1295 for (i__ = 1; i__ <= i__1; ++i__) {
1296 d__[i__] = temp * d__[i__] - g[i__] - hs[i__];
1297 /* Computing 2nd power */
1298 d__1 = d__[i__];
1299 dd += d__1 * d__1;
1300 ds += d__[i__] * step[i__];
1301 /* L80: */
1302 /* Computing 2nd power */
1303 d__1 = step[i__];
1304 ss += d__1 * d__1;
1305 }
1306 if (ds <= zero) {
1307 goto L160;
1308 }
1309 if (ss < delsq) {
1310 goto L40;
1311 }
1312 }
1313 *crvmin = zero;
1314 itersw = iterc;
1315
1316 /* Test whether an alternative iteration is required. */
1317
1318 L90:
1319 if (gg <= ggbeg * 1e-4) {
1320 goto L160;
1321 }
1322 sg = zero;
1323 shs = zero;
1324 i__1 = *n;
1325 for (i__ = 1; i__ <= i__1; ++i__) {
1326 sg += step[i__] * g[i__];
1327 /* L100: */
1328 shs += step[i__] * hs[i__];
1329 }
1330 sgk = sg + shs;
1331 angtest = sgk / sqrt(gg * delsq);
1332 if (angtest <= -.99) {
1333 goto L160;
1334 }
1335
1336 /* Begin the alternative iteration by calculating D and HD and some */
1337 /* scalar products. */
1338
1339 ++iterc;
1340 temp = sqrt(delsq * gg - sgk * sgk);
1341 tempa = delsq / temp;
1342 tempb = sgk / temp;
1343 i__1 = *n;
1344 for (i__ = 1; i__ <= i__1; ++i__) {
1345 /* L110: */
1346 d__[i__] = tempa * (g[i__] + hs[i__]) - tempb * step[i__];
1347 }
1348 goto L170;
1349 L120:
1350 dg = zero;
1351 dhd = zero;
1352 dhs = zero;
1353 i__1 = *n;
1354 for (i__ = 1; i__ <= i__1; ++i__) {
1355 dg += d__[i__] * g[i__];
1356 dhd += hd[i__] * d__[i__];
1357 /* L130: */
1358 dhs += hd[i__] * step[i__];
1359 }
1360
1361 /* Seek the value of the angle that minimizes Q. */
1362
1363 cf = half * (shs - dhd);
1364 qbeg = sg + cf;
1365 qsav = qbeg;
1366 qmin = qbeg;
1367 isave = 0;
1368 iu = 49;
1369 temp = twopi / (doublereal) (iu + 1);
1370 i__1 = iu;
1371 for (i__ = 1; i__ <= i__1; ++i__) {
1372 angle = (doublereal) i__ * temp;
1373 cth = cos(angle);
1374 sth = sin(angle);
1375 qnew = (sg + cf * cth) * cth + (dg + dhs * cth) * sth;
1376 if (qnew < qmin) {
1377 qmin = qnew;
1378 isave = i__;
1379 tempa = qsav;
1380 } else if (i__ == isave + 1) {
1381 tempb = qnew;
1382 }
1383 /* L140: */
1384 qsav = qnew;
1385 }
1386 if ((doublereal) isave == zero) {
1387 tempa = qnew;
1388 }
1389 if (isave == iu) {
1390 tempb = qbeg;
1391 }
1392 angle = zero;
1393 if (tempa != tempb) {
1394 tempa -= qmin;
1395 tempb -= qmin;
1396 angle = half * (tempa - tempb) / (tempa + tempb);
1397 }
1398 angle = temp * ((doublereal) isave + angle);
1399
1400 /* Calculate the new STEP and HS. Then test for convergence. */
1401
1402 cth = cos(angle);
1403 sth = sin(angle);
1404 reduc = qbeg - (sg + cf * cth) * cth - (dg + dhs * cth) * sth;
1405 gg = zero;
1406 i__1 = *n;
1407 for (i__ = 1; i__ <= i__1; ++i__) {
1408 step[i__] = cth * step[i__] + sth * d__[i__];
1409 hs[i__] = cth * hs[i__] + sth * hd[i__];
1410 /* L150: */
1411 /* Computing 2nd power */
1412 d__1 = g[i__] + hs[i__];
1413 gg += d__1 * d__1;
1414 }
1415 qred += reduc;
1416 ratio = reduc / qred;
1417 if (iterc < itermax && ratio > .01) {
1418 goto L90;
1419 }
1420 L160:
1421 return 0;
1422
1423 /* The following instructions act as a subroutine for setting the vecto
1424 r*/
1425 /* HD to the vector D multiplied by the second derivative matrix of Q.
1426 */
1427 /* They are called from three different places, which are distinguished
1428 */
1429 /* by the value of ITERC. */
1430
1431 L170:
1432 i__1 = *n;
1433 for (i__ = 1; i__ <= i__1; ++i__) {
1434 /* L180: */
1435 hd[i__] = zero;
1436 }
1437 i__1 = *npt;
1438 for (k = 1; k <= i__1; ++k) {
1439 temp = zero;
1440 i__2 = *n;
1441 for (j = 1; j <= i__2; ++j) {
1442 /* L190: */
1443 temp += xpt[k + j * xpt_dim1] * d__[j];
1444 }
1445 temp *= pq[k];
1446 i__2 = *n;
1447 for (i__ = 1; i__ <= i__2; ++i__) {
1448 /* L200: */
1449 hd[i__] += temp * xpt[k + i__ * xpt_dim1];
1450 }
1451 }
1452 ih = 0;
1453 i__2 = *n;
1454 for (j = 1; j <= i__2; ++j) {
1455 i__1 = j;
1456 for (i__ = 1; i__ <= i__1; ++i__) {
1457 ++ih;
1458 if (i__ < j) {
1459 hd[j] += hq[ih] * d__[i__];
1460 }
1461 /* L210: */
1462 hd[i__] += hq[ih] * d__[j];
1463 }
1464 }
1465 if (iterc == 0) {
1466 goto L20;
1467 }
1468 if (iterc <= itersw) {
1469 goto L50;
1470 }
1471 goto L120;
1472 } /* trsapp_ */
1473
update_(integer * n,integer * npt,doublereal * bmat,doublereal * zmat,integer * idz,integer * ndim,doublereal * vlag,doublereal * beta,integer * knew,doublereal * w)1474 /* Subroutine */ int update_(integer *n, integer *npt, doublereal *bmat,
1475 doublereal *zmat, integer *idz, integer *ndim, doublereal *vlag,
1476 doublereal *beta, integer *knew, doublereal *w)
1477 {
1478 /* System generated locals */
1479 integer bmat_dim1=0, bmat_offset=0, zmat_dim1=0, zmat_offset=0, i__1=0, i__2=0;
1480 doublereal d__1=0, d__2=0;
1481
1482 /* Builtin functions */
1483 double sqrt(doublereal);
1484
1485 /* Local variables */
1486 STATIC doublereal temp=0;
1487 STATIC integer nptm=0;
1488 STATIC doublereal zero=0;
1489 STATIC integer i__=0, j=0, iflag=0;
1490 STATIC doublereal scala=0, scalb=0, alpha=0, denom=0, tempa=0, tempb=0, tausq=0;
1491 STATIC integer ja=0, jb=0, jl=0, jp=0;
1492 STATIC doublereal one=0, tau=0;
1493
1494
1495 /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift the
1496 */
1497 /* interpolation point that has index KNEW. On entry, VLAG contains the
1498 */
1499 /* components of the vector Theta*Wcheck+e_b of the updating formula
1500 */
1501 /* (6.11), and BETA holds the value of the parameter that has this name
1502 .*/
1503 /* The vector W is used for working space. */
1504
1505 /* Set some constants. */
1506
1507 /* Parameter adjustments */
1508 zmat_dim1 = *npt;
1509 zmat_offset = zmat_dim1 + 1;
1510 zmat -= zmat_offset;
1511 bmat_dim1 = *ndim;
1512 bmat_offset = bmat_dim1 + 1;
1513 bmat -= bmat_offset;
1514 --vlag;
1515 --w;
1516
1517 /* Function Body */
1518 one = 1.;
1519 zero = 0.;
1520 nptm = *npt - *n - 1;
1521
1522 /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */
1523
1524 jl = 1;
1525 i__1 = nptm;
1526 for (j = 2; j <= i__1; ++j) {
1527 if (j == *idz) {
1528 jl = *idz;
1529 } else if (zmat[*knew + j * zmat_dim1] != zero) {
1530 /* Computing 2nd power */
1531 d__1 = zmat[*knew + jl * zmat_dim1];
1532 /* Computing 2nd power */
1533 d__2 = zmat[*knew + j * zmat_dim1];
1534 temp = sqrt(d__1 * d__1 + d__2 * d__2);
1535 tempa = zmat[*knew + jl * zmat_dim1] / temp;
1536 tempb = zmat[*knew + j * zmat_dim1] / temp;
1537 i__2 = *npt;
1538 for (i__ = 1; i__ <= i__2; ++i__) {
1539 temp = tempa * zmat[i__ + jl * zmat_dim1] + tempb * zmat[i__
1540 + j * zmat_dim1];
1541 zmat[i__ + j * zmat_dim1] = tempa * zmat[i__ + j * zmat_dim1]
1542 - tempb * zmat[i__ + jl * zmat_dim1];
1543 /* L10: */
1544 zmat[i__ + jl * zmat_dim1] = temp;
1545 }
1546 zmat[*knew + j * zmat_dim1] = zero;
1547 }
1548 /* L20: */
1549 }
1550
1551 /* Put the first NPT components of the KNEW-th column of HLAG into W,
1552 */
1553 /* and calculate the parameters of the updating formula. */
1554
1555 tempa = zmat[*knew + zmat_dim1];
1556 if (*idz >= 2) {
1557 tempa = -tempa;
1558 }
1559 if (jl > 1) {
1560 tempb = zmat[*knew + jl * zmat_dim1];
1561 }
1562 i__1 = *npt;
1563 for (i__ = 1; i__ <= i__1; ++i__) {
1564 w[i__] = tempa * zmat[i__ + zmat_dim1];
1565 if (jl > 1) {
1566 w[i__] += tempb * zmat[i__ + jl * zmat_dim1];
1567 }
1568 /* L30: */
1569 }
1570 alpha = w[*knew];
1571 tau = vlag[*knew];
1572 tausq = tau * tau;
1573 denom = alpha * *beta + tausq;
1574 vlag[*knew] -= one;
1575
1576 /* Complete the updating of ZMAT when there is only one nonzero element
1577 */
1578 /* in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to o
1579 ne,*/
1580 /* then the first column of ZMAT will be exchanged with another one lat
1581 er.*/
1582
1583 iflag = 0;
1584 if (jl == 1) {
1585 temp = sqrt((abs(denom)));
1586 tempb = tempa / temp;
1587 tempa = tau / temp;
1588 i__1 = *npt;
1589 for (i__ = 1; i__ <= i__1; ++i__) {
1590 /* L40: */
1591 zmat[i__ + zmat_dim1] = tempa * zmat[i__ + zmat_dim1] - tempb *
1592 vlag[i__];
1593 }
1594 if (*idz == 1 && temp < zero) {
1595 *idz = 2;
1596 }
1597 if (*idz >= 2 && temp >= zero) {
1598 iflag = 1;
1599 }
1600 } else {
1601
1602 /* Complete the updating of ZMAT in the alternative case. */
1603
1604 ja = 1;
1605 if (*beta >= zero) {
1606 ja = jl;
1607 }
1608 jb = jl + 1 - ja;
1609 temp = zmat[*knew + jb * zmat_dim1] / denom;
1610 tempa = temp * *beta;
1611 tempb = temp * tau;
1612 temp = zmat[*knew + ja * zmat_dim1];
1613 scala = one / sqrt(abs(*beta) * temp * temp + tausq);
1614 scalb = scala * sqrt((abs(denom)));
1615 i__1 = *npt;
1616 for (i__ = 1; i__ <= i__1; ++i__) {
1617 zmat[i__ + ja * zmat_dim1] = scala * (tau * zmat[i__ + ja *
1618 zmat_dim1] - temp * vlag[i__]);
1619 /* L50: */
1620 zmat[i__ + jb * zmat_dim1] = scalb * (zmat[i__ + jb * zmat_dim1]
1621 - tempa * w[i__] - tempb * vlag[i__]);
1622 }
1623 if (denom <= zero) {
1624 if (*beta < zero) {
1625 ++(*idz);
1626 }
1627 if (*beta >= zero) {
1628 iflag = 1;
1629 }
1630 }
1631 }
1632
1633 /* IDZ is reduced in the following case, and usually the first column
1634 */
1635 /* of ZMAT is exchanged with a later one. */
1636
1637 if (iflag == 1) {
1638 --(*idz);
1639 i__1 = *npt;
1640 for (i__ = 1; i__ <= i__1; ++i__) {
1641 temp = zmat[i__ + zmat_dim1];
1642 zmat[i__ + zmat_dim1] = zmat[i__ + *idz * zmat_dim1];
1643 /* L60: */
1644 zmat[i__ + *idz * zmat_dim1] = temp;
1645 }
1646 }
1647
1648 /* Finally, update the matrix BMAT. */
1649
1650 i__1 = *n;
1651 for (j = 1; j <= i__1; ++j) {
1652 jp = *npt + j;
1653 w[jp] = bmat[*knew + j * bmat_dim1];
1654 tempa = (alpha * vlag[jp] - tau * w[jp]) / denom;
1655 tempb = (-(*beta) * w[jp] - tau * vlag[jp]) / denom;
1656 i__2 = jp;
1657 for (i__ = 1; i__ <= i__2; ++i__) {
1658 bmat[i__ + j * bmat_dim1] = bmat[i__ + j * bmat_dim1] + tempa *
1659 vlag[i__] + tempb * w[i__];
1660 if (i__ > *npt) {
1661 bmat[jp + (i__ - *npt) * bmat_dim1] = bmat[i__ + j *
1662 bmat_dim1];
1663 }
1664 /* L70: */
1665 }
1666 }
1667 return 0;
1668 } /* update_ */
1669
bigden_(integer * n,integer * npt,doublereal * xopt,doublereal * xpt,doublereal * bmat,doublereal * zmat,integer * idz,integer * ndim,integer * kopt,integer * knew,doublereal * d__,doublereal * w,doublereal * vlag,doublereal * beta,doublereal * s,doublereal * wvec,doublereal * prod)1670 /* Subroutine */ int bigden_(integer *n, integer *npt, doublereal *xopt,
1671 doublereal *xpt, doublereal *bmat, doublereal *zmat, integer *idz,
1672 integer *ndim, integer *kopt, integer *knew, doublereal *d__,
1673 doublereal *w, doublereal *vlag, doublereal *beta, doublereal *s,
1674 doublereal *wvec, doublereal *prod)
1675 {
1676 /* System generated locals */
1677 integer xpt_dim1=0, xpt_offset=0, bmat_dim1=0, bmat_offset=0, zmat_dim1=0,
1678 zmat_offset=0, wvec_dim1=0, wvec_offset=0, prod_dim1=0, prod_offset=0, i__1=0,
1679 i__2=0;
1680 doublereal d__1=0;
1681
1682 /* Builtin functions */
1683 double atan(doublereal), sqrt(doublereal), cos(doublereal), sin(
1684 doublereal);
1685
1686 /* Local variables */
1687 STATIC doublereal diff=0, half=0, temp=0;
1688 STATIC integer ksav=0;
1689 STATIC doublereal step=0;
1690 STATIC integer nptm=0;
1691 STATIC doublereal zero=0;
1692 STATIC integer i__=0, j=0, k=0;
1693 STATIC doublereal alpha=0, angle=0, denex[9]={0,0,0,0,0,0,0,0,0};
1694 STATIC integer iterc=0;
1695 STATIC doublereal tempa=0, tempb=0, tempc=0;
1696 STATIC integer isave=0;
1697 STATIC doublereal ssden=0, dtest=0, quart=0, xoptd=0, twopi=0, xopts=0, dd=0;
1698 STATIC integer jc=0;
1699 STATIC doublereal ds=0;
1700 STATIC integer ip=0, iu=0, nw=0;
1701 STATIC doublereal ss=0, denold=0, denmax=0, densav=0, dstemp=0, sumold=0, sstemp=0,
1702 xoptsq=0, den[9]={0,0,0,0,0,0,0,0,0}, one=0,
1703 par[9]={0,0,0,0,0,0,0,0,0}, tau=0, sum=0, two=0;
1704
1705
1706 /* N is the number of variables. */
1707 /* NPT is the number of interpolation equations. */
1708 /* XOPT is the best interpolation point so far. */
1709 /* XPT contains the coordinates of the current interpolation points.
1710 */
1711 /* BMAT provides the last N columns of H. */
1712 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix
1713 of H.*/
1714 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
1715 /* KOPT is the index of the optimal interpolation point. */
1716 /* KNEW is the index of the interpolation point that is going to be mov
1717 ed.*/
1718 /* D will be set to the step from XOPT to the new point, and on entry i
1719 t*/
1720 /* should be the D that was calculated by the last call of BIGLAG. Th
1721 e*/
1722 /* length of the initial D provides a trust region bound on the final
1723 D.*/
1724 /* W will be set to Wcheck for the final choice of D. */
1725 /* VLAG will be set to Theta*Wcheck+e_b for the final choice of D. */
1726 /* BETA will be set to the value that will occur in the updating formul
1727 a*/
1728 /* when the KNEW-th interpolation point is moved to its new position.
1729 */
1730 /* S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used
1731 */
1732 /* for working space. */
1733
1734 /* D is calculated in a way that should provide a denominator with a la
1735 rge*/
1736 /* modulus in the updating formula when the KNEW-th interpolation point
1737 is*/
1738 /* shifted to the new position XOPT+D. */
1739
1740 /* Set some constants. */
1741
1742 /* Parameter adjustments */
1743 zmat_dim1 = *npt;
1744 zmat_offset = zmat_dim1 + 1;
1745 zmat -= zmat_offset;
1746 xpt_dim1 = *npt;
1747 xpt_offset = xpt_dim1 + 1;
1748 xpt -= xpt_offset;
1749 --xopt;
1750 prod_dim1 = *ndim;
1751 prod_offset = prod_dim1 + 1;
1752 prod -= prod_offset;
1753 wvec_dim1 = *ndim;
1754 wvec_offset = wvec_dim1 + 1;
1755 wvec -= wvec_offset;
1756 bmat_dim1 = *ndim;
1757 bmat_offset = bmat_dim1 + 1;
1758 bmat -= bmat_offset;
1759 --d__;
1760 --w;
1761 --vlag;
1762 --s;
1763
1764 /* Function Body */
1765 half = .5;
1766 one = 1.;
1767 quart = .25;
1768 two = 2.;
1769 zero = 0.;
1770 #if 0
1771 twopi = atan(one) * 8.;
1772 #else
1773 twopi = 3.14159265358979323846 * 2.0 ;
1774 #endif
1775 nptm = *npt - *n - 1;
1776
1777 /* Store the first NPT elements of the KNEW-th column of H in W(N+1)
1778 */
1779 /* to W(N+NPT). */
1780
1781 i__1 = *npt;
1782 for (k = 1; k <= i__1; ++k) {
1783 /* L10: */
1784 w[*n + k] = zero;
1785 }
1786 i__1 = nptm;
1787 for (j = 1; j <= i__1; ++j) {
1788 temp = zmat[*knew + j * zmat_dim1];
1789 if (j < *idz) {
1790 temp = -temp;
1791 }
1792 i__2 = *npt;
1793 for (k = 1; k <= i__2; ++k) {
1794 /* L20: */
1795 w[*n + k] += temp * zmat[k + j * zmat_dim1];
1796 }
1797 }
1798 alpha = w[*n + *knew];
1799
1800 /* The initial search direction D is taken from the last call of BIGLAG
1801 ,*/
1802 /* and the initial S is set below, usually to the direction from X_OPT
1803 */
1804 /* to X_KNEW, but a different direction to an interpolation point may
1805 */
1806 /* be chosen, in order to prevent S from being nearly parallel to D.
1807 */
1808
1809 dd = zero;
1810 ds = zero;
1811 ss = zero;
1812 xoptsq = zero;
1813 i__2 = *n;
1814 for (i__ = 1; i__ <= i__2; ++i__) {
1815 /* Computing 2nd power */
1816 d__1 = d__[i__];
1817 dd += d__1 * d__1;
1818 s[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
1819 ds += d__[i__] * s[i__];
1820 /* Computing 2nd power */
1821 d__1 = s[i__];
1822 ss += d__1 * d__1;
1823 /* L30: */
1824 /* Computing 2nd power */
1825 d__1 = xopt[i__];
1826 xoptsq += d__1 * d__1;
1827 }
1828 if (ds * ds > dd * .99 * ss) {
1829 ksav = *knew;
1830 dtest = ds * ds / ss;
1831 i__2 = *npt;
1832 for (k = 1; k <= i__2; ++k) {
1833 if (k != *kopt) {
1834 dstemp = zero;
1835 sstemp = zero;
1836 i__1 = *n;
1837 for (i__ = 1; i__ <= i__1; ++i__) {
1838 diff = xpt[k + i__ * xpt_dim1] - xopt[i__];
1839 dstemp += d__[i__] * diff;
1840 /* L40: */
1841 sstemp += diff * diff;
1842 }
1843 if (dstemp * dstemp / sstemp < dtest) {
1844 ksav = k;
1845 dtest = dstemp * dstemp / sstemp;
1846 ds = dstemp;
1847 ss = sstemp;
1848 }
1849 }
1850 /* L50: */
1851 }
1852 i__2 = *n;
1853 for (i__ = 1; i__ <= i__2; ++i__) {
1854 /* L60: */
1855 s[i__] = xpt[ksav + i__ * xpt_dim1] - xopt[i__];
1856 }
1857 }
1858 ssden = dd * ss - ds * ds;
1859 iterc = 0;
1860 densav = zero;
1861
1862 /* Begin the iteration by overwriting S with a vector that has the */
1863 /* required length and direction. */
1864
1865 L70:
1866 ++iterc;
1867 temp = one / sqrt(ssden);
1868 xoptd = zero;
1869 xopts = zero;
1870 i__2 = *n;
1871 for (i__ = 1; i__ <= i__2; ++i__) {
1872 s[i__] = temp * (dd * s[i__] - ds * d__[i__]);
1873 xoptd += xopt[i__] * d__[i__];
1874 /* L80: */
1875 xopts += xopt[i__] * s[i__];
1876 }
1877
1878 /* Set the coefficients of the first two terms of BETA. */
1879
1880 tempa = half * xoptd * xoptd;
1881 tempb = half * xopts * xopts;
1882 den[0] = dd * (xoptsq + half * dd) + tempa + tempb;
1883 den[1] = two * xoptd * dd;
1884 den[2] = two * xopts * dd;
1885 den[3] = tempa - tempb;
1886 den[4] = xoptd * xopts;
1887 for (i__ = 6; i__ <= 9; ++i__) {
1888 /* L90: */
1889 den[i__ - 1] = zero;
1890 }
1891
1892 /* Put the coefficients of Wcheck in WVEC. */
1893
1894 i__2 = *npt;
1895 for (k = 1; k <= i__2; ++k) {
1896 tempa = zero;
1897 tempb = zero;
1898 tempc = zero;
1899 i__1 = *n;
1900 for (i__ = 1; i__ <= i__1; ++i__) {
1901 tempa += xpt[k + i__ * xpt_dim1] * d__[i__];
1902 tempb += xpt[k + i__ * xpt_dim1] * s[i__];
1903 /* L100: */
1904 tempc += xpt[k + i__ * xpt_dim1] * xopt[i__];
1905 }
1906 wvec[k + wvec_dim1] = quart * (tempa * tempa + tempb * tempb);
1907 wvec[k + (wvec_dim1 << 1)] = tempa * tempc;
1908 wvec[k + wvec_dim1 * 3] = tempb * tempc;
1909 wvec[k + (wvec_dim1 << 2)] = quart * (tempa * tempa - tempb * tempb);
1910 /* L110: */
1911 wvec[k + wvec_dim1 * 5] = half * tempa * tempb;
1912 }
1913 i__2 = *n;
1914 for (i__ = 1; i__ <= i__2; ++i__) {
1915 ip = i__ + *npt;
1916 wvec[ip + wvec_dim1] = zero;
1917 wvec[ip + (wvec_dim1 << 1)] = d__[i__];
1918 wvec[ip + wvec_dim1 * 3] = s[i__];
1919 wvec[ip + (wvec_dim1 << 2)] = zero;
1920 /* L120: */
1921 wvec[ip + wvec_dim1 * 5] = zero;
1922 }
1923
1924 /* Put the coefficents of THETA*Wcheck in PROD. */
1925
1926 for (jc = 1; jc <= 5; ++jc) {
1927 nw = *npt;
1928 if (jc == 2 || jc == 3) {
1929 nw = *ndim;
1930 }
1931 i__2 = *npt;
1932 for (k = 1; k <= i__2; ++k) {
1933 /* L130: */
1934 prod[k + jc * prod_dim1] = zero;
1935 }
1936 i__2 = nptm;
1937 for (j = 1; j <= i__2; ++j) {
1938 sum = zero;
1939 i__1 = *npt;
1940 for (k = 1; k <= i__1; ++k) {
1941 /* L140: */
1942 sum += zmat[k + j * zmat_dim1] * wvec[k + jc * wvec_dim1];
1943 }
1944 if (j < *idz) {
1945 sum = -sum;
1946 }
1947 i__1 = *npt;
1948 for (k = 1; k <= i__1; ++k) {
1949 /* L150: */
1950 prod[k + jc * prod_dim1] += sum * zmat[k + j * zmat_dim1];
1951 }
1952 }
1953 if (nw == *ndim) {
1954 i__1 = *npt;
1955 for (k = 1; k <= i__1; ++k) {
1956 sum = zero;
1957 i__2 = *n;
1958 for (j = 1; j <= i__2; ++j) {
1959 /* L160: */
1960 sum += bmat[k + j * bmat_dim1] * wvec[*npt + j + jc *
1961 wvec_dim1];
1962 }
1963 /* L170: */
1964 prod[k + jc * prod_dim1] += sum;
1965 }
1966 }
1967 i__1 = *n;
1968 for (j = 1; j <= i__1; ++j) {
1969 sum = zero;
1970 i__2 = nw;
1971 for (i__ = 1; i__ <= i__2; ++i__) {
1972 /* L180: */
1973 sum += bmat[i__ + j * bmat_dim1] * wvec[i__ + jc * wvec_dim1];
1974 }
1975 /* L190: */
1976 prod[*npt + j + jc * prod_dim1] = sum;
1977 }
1978 }
1979
1980 /* Include in DEN the part of BETA that depends on THETA. */
1981
1982 i__1 = *ndim;
1983 for (k = 1; k <= i__1; ++k) {
1984 sum = zero;
1985 for (i__ = 1; i__ <= 5; ++i__) {
1986 par[i__ - 1] = half * prod[k + i__ * prod_dim1] * wvec[k + i__ *
1987 wvec_dim1];
1988 /* L200: */
1989 sum += par[i__ - 1];
1990 }
1991 den[0] = den[0] - par[0] - sum;
1992 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 1)] + prod[k + (
1993 prod_dim1 << 1)] * wvec[k + wvec_dim1];
1994 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + (wvec_dim1 << 2)] +
1995 prod[k + (prod_dim1 << 2)] * wvec[k + (wvec_dim1 << 1)];
1996 tempc = prod[k + prod_dim1 * 3] * wvec[k + wvec_dim1 * 5] + prod[k +
1997 prod_dim1 * 5] * wvec[k + wvec_dim1 * 3];
1998 den[1] = den[1] - tempa - half * (tempb + tempc);
1999 den[5] -= half * (tempb - tempc);
2000 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 3] + prod[k +
2001 prod_dim1 * 3] * wvec[k + wvec_dim1];
2002 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 5] + prod[k
2003 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 1)];
2004 tempc = prod[k + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 2)] + prod[k
2005 + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 3];
2006 den[2] = den[2] - tempa - half * (tempb - tempc);
2007 den[6] -= half * (tempb + tempc);
2008 tempa = prod[k + prod_dim1] * wvec[k + (wvec_dim1 << 2)] + prod[k + (
2009 prod_dim1 << 2)] * wvec[k + wvec_dim1];
2010 den[3] = den[3] - tempa - par[1] + par[2];
2011 tempa = prod[k + prod_dim1] * wvec[k + wvec_dim1 * 5] + prod[k +
2012 prod_dim1 * 5] * wvec[k + wvec_dim1];
2013 tempb = prod[k + (prod_dim1 << 1)] * wvec[k + wvec_dim1 * 3] + prod[k
2014 + prod_dim1 * 3] * wvec[k + (wvec_dim1 << 1)];
2015 den[4] = den[4] - tempa - half * tempb;
2016 den[7] = den[7] - par[3] + par[4];
2017 tempa = prod[k + (prod_dim1 << 2)] * wvec[k + wvec_dim1 * 5] + prod[k
2018 + prod_dim1 * 5] * wvec[k + (wvec_dim1 << 2)];
2019 /* L210: */
2020 den[8] -= half * tempa;
2021 }
2022
2023 /* Extend DEN so that it holds all the coefficients of DENOM. */
2024
2025 sum = zero;
2026 for (i__ = 1; i__ <= 5; ++i__) {
2027 /* Computing 2nd power */
2028 d__1 = prod[*knew + i__ * prod_dim1];
2029 par[i__ - 1] = half * (d__1 * d__1);
2030 /* L220: */
2031 sum += par[i__ - 1];
2032 }
2033 denex[0] = alpha * den[0] + par[0] + sum;
2034 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 1)];
2035 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + (prod_dim1 << 2)];
2036 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + prod_dim1 * 5];
2037 denex[1] = alpha * den[1] + tempa + tempb + tempc;
2038 denex[5] = alpha * den[5] + tempb - tempc;
2039 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 3];
2040 tempb = prod[*knew + (prod_dim1 << 1)] * prod[*knew + prod_dim1 * 5];
2041 tempc = prod[*knew + prod_dim1 * 3] * prod[*knew + (prod_dim1 << 2)];
2042 denex[2] = alpha * den[2] + tempa + tempb - tempc;
2043 denex[6] = alpha * den[6] + tempb + tempc;
2044 tempa = two * prod[*knew + prod_dim1] * prod[*knew + (prod_dim1 << 2)];
2045 denex[3] = alpha * den[3] + tempa + par[1] - par[2];
2046 tempa = two * prod[*knew + prod_dim1] * prod[*knew + prod_dim1 * 5];
2047 denex[4] = alpha * den[4] + tempa + prod[*knew + (prod_dim1 << 1)] * prod[
2048 *knew + prod_dim1 * 3];
2049 denex[7] = alpha * den[7] + par[3] - par[4];
2050 denex[8] = alpha * den[8] + prod[*knew + (prod_dim1 << 2)] * prod[*knew +
2051 prod_dim1 * 5];
2052
2053 /* Seek the value of the angle that maximizes the modulus of DENOM. */
2054
2055 sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7];
2056 denold = sum;
2057 denmax = sum;
2058 isave = 0;
2059 iu = 49;
2060 temp = twopi / (doublereal) (iu + 1);
2061 par[0] = one;
2062 i__1 = iu;
2063 for (i__ = 1; i__ <= i__1; ++i__) {
2064 angle = (doublereal) i__ * temp;
2065 par[1] = cos(angle);
2066 par[2] = sin(angle);
2067 if(0&&DEBUG) fprintf(stderr," cos(angle=%.14g)=%.14g sin()=%.14g\n",angle,par[1],par[2]) ;
2068 for (j = 4; j <= 8; j += 2) {
2069 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
2070 /* L230: */
2071 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
2072 }
2073 sumold = sum;
2074 sum = zero;
2075 for (j = 1; j <= 9; ++j) {
2076 /* L240: */
2077 sum += denex[j - 1] * par[j - 1];
2078 }
2079 if (abs(sum) > abs(denmax)) {
2080 denmax = sum;
2081 isave = i__;
2082 tempa = sumold;
2083 } else if (i__ == isave + 1) {
2084 tempb = sum;
2085 }
2086 /* L250: */
2087 }
2088 if (isave == 0) {
2089 tempa = sum;
2090 }
2091 if (isave == iu) {
2092 tempb = denold;
2093 }
2094 step = zero;
2095 if (tempa != tempb) {
2096 tempa -= denmax;
2097 tempb -= denmax;
2098 step = half * (tempa - tempb) / (tempa + tempb);
2099 }
2100 angle = temp * ((doublereal) isave + step);
2101
2102 /* Calculate the new parameters of the denominator, the new VLAG vector
2103 */
2104 /* and the new D. Then test for convergence. */
2105
2106 par[1] = cos(angle);
2107 par[2] = sin(angle);
2108 if(0&&DEBUG) fprintf(stderr," cos(angle=%.14g)=%.14g sin()=%.14g\n",angle,par[1],par[2]) ;
2109 for (j = 4; j <= 8; j += 2) {
2110 par[j - 1] = par[1] * par[j - 3] - par[2] * par[j - 2];
2111 /* L260: */
2112 par[j] = par[1] * par[j - 2] + par[2] * par[j - 3];
2113 }
2114 *beta = zero;
2115 denmax = zero;
2116 for (j = 1; j <= 9; ++j) {
2117 *beta += den[j - 1] * par[j - 1];
2118 /* L270: */
2119 denmax += denex[j - 1] * par[j - 1];
2120 }
2121 i__1 = *ndim;
2122 for (k = 1; k <= i__1; ++k) {
2123 vlag[k] = zero;
2124 for (j = 1; j <= 5; ++j) {
2125 /* L280: */
2126 vlag[k] += prod[k + j * prod_dim1] * par[j - 1];
2127 }
2128 }
2129 tau = vlag[*knew];
2130 dd = zero;
2131 tempa = zero;
2132 tempb = zero;
2133 i__1 = *n;
2134 for (i__ = 1; i__ <= i__1; ++i__) {
2135 d__[i__] = par[1] * d__[i__] + par[2] * s[i__];
2136 w[i__] = xopt[i__] + d__[i__];
2137 /* Computing 2nd power */
2138 d__1 = d__[i__];
2139 dd += d__1 * d__1;
2140 tempa += d__[i__] * w[i__];
2141 /* L290: */
2142 tempb += w[i__] * w[i__];
2143 }
2144 if (iterc >= *n) {
2145 goto L340;
2146 }
2147 if (iterc > 1) {
2148 densav = AFmax(densav,denold);
2149 }
2150 if (abs(denmax) <= abs(densav) * 1.1) {
2151 goto L340;
2152 }
2153 densav = denmax;
2154
2155 /* Set S to half the gradient of the denominator with respect to D. */
2156 /* Then branch for the next iteration. */
2157
2158 i__1 = *n;
2159 for (i__ = 1; i__ <= i__1; ++i__) {
2160 temp = tempa * xopt[i__] + tempb * d__[i__] - vlag[*npt + i__];
2161 /* L300: */
2162 s[i__] = tau * bmat[*knew + i__ * bmat_dim1] + alpha * temp;
2163 }
2164 i__1 = *npt;
2165 for (k = 1; k <= i__1; ++k) {
2166 sum = zero;
2167 i__2 = *n;
2168 for (j = 1; j <= i__2; ++j) {
2169 /* L310: */
2170 sum += xpt[k + j * xpt_dim1] * w[j];
2171 }
2172 temp = (tau * w[*n + k] - alpha * vlag[k]) * sum;
2173 i__2 = *n;
2174 for (i__ = 1; i__ <= i__2; ++i__) {
2175 /* L320: */
2176 s[i__] += temp * xpt[k + i__ * xpt_dim1];
2177 }
2178 }
2179 ss = zero;
2180 ds = zero;
2181 i__2 = *n;
2182 for (i__ = 1; i__ <= i__2; ++i__) {
2183 /* Computing 2nd power */
2184 d__1 = s[i__];
2185 ss += d__1 * d__1;
2186 /* L330: */
2187 ds += d__[i__] * s[i__];
2188 }
2189 ssden = dd * ss - ds * ds;
2190 if (ssden >= dd * 1e-8 * ss) {
2191 goto L70;
2192 }
2193
2194 /* Set the vector W before the RETURN from the subroutine. */
2195
2196 L340:
2197 i__2 = *ndim;
2198 for (k = 1; k <= i__2; ++k) {
2199 w[k] = zero;
2200 for (j = 1; j <= 5; ++j) {
2201 /* L350: */
2202 w[k] += wvec[k + j * wvec_dim1] * par[j - 1];
2203 }
2204 }
2205 vlag[*kopt] += one;
2206 return 0;
2207 } /* bigden_ */
2208
biglag_(integer * n,integer * npt,doublereal * xopt,doublereal * xpt,doublereal * bmat,doublereal * zmat,integer * idz,integer * ndim,integer * knew,doublereal * delta,doublereal * d__,doublereal * alpha,doublereal * hcol,doublereal * gc,doublereal * gd,doublereal * s,doublereal * w)2209 /* Subroutine */ int biglag_(integer *n, integer *npt, doublereal *xopt,
2210 doublereal *xpt, doublereal *bmat, doublereal *zmat, integer *idz,
2211 integer *ndim, integer *knew, doublereal *delta, doublereal *d__,
2212 doublereal *alpha, doublereal *hcol, doublereal *gc, doublereal *gd,
2213 doublereal *s, doublereal *w)
2214 {
2215 /* System generated locals */
2216 integer xpt_dim1=0, xpt_offset=0, bmat_dim1=0, bmat_offset=0, zmat_dim1=0,
2217 zmat_offset=0, i__1=0, i__2=0;
2218 doublereal d__1=0;
2219
2220 /* Builtin functions */
2221 double atan(doublereal), sqrt(doublereal), cos(doublereal), sin(
2222 doublereal);
2223
2224 /* Local variables */
2225 STATIC doublereal half=0, temp=0, step=0;
2226 STATIC integer nptm=0;
2227 STATIC doublereal zero=0;
2228 STATIC integer i__=0, j=0, k=0;
2229 STATIC doublereal angle=0, scale=0, denom=0;
2230 STATIC integer iterc=0, isave=0;
2231 STATIC doublereal delsq=0, tempa=0, tempb=0, twopi=0, dd=0, gg=0;
2232 STATIC integer iu=0;
2233 STATIC doublereal sp=0, ss=0, taubeg=0, tauold=0, cf1=0, cf2=0, cf3=0, cf4=0, cf5=0, taumax=0,
2234 dhd=0, cth=0, one=0, tau=0, sth=0, sum=0;
2235
2236
2237 /* N is the number of variables. */
2238 /* NPT is the number of interpolation equations. */
2239 /* XOPT is the best interpolation point so far. */
2240 /* XPT contains the coordinates of the current interpolation points.
2241 */
2242 /* BMAT provides the last N columns of H. */
2243 /* ZMAT and IDZ give a factorization of the first NPT by NPT submatrix
2244 of H.*/
2245 /* NDIM is the first dimension of BMAT and has the value NPT+N. */
2246 /* KNEW is the index of the interpolation point that is going to be mov
2247 ed.*/
2248 /* DELTA is the current trust region bound. */
2249 /* D will be set to the step from XOPT to the new point. */
2250 /* ALPHA will be set to the KNEW-th diagonal element of the H matrix.
2251 */
2252 /* HCOL, GC, GD, S and W will be used for working space. */
2253
2254 /* The step D is calculated in a way that attempts to maximize the modu
2255 lus*/
2256 /* of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC
2257 is*/
2258 /* the KNEW-th Lagrange function. */
2259
2260 /* Set some constants. */
2261
2262 /* Parameter adjustments */
2263 zmat_dim1 = *npt;
2264 zmat_offset = zmat_dim1 + 1;
2265 zmat -= zmat_offset;
2266 xpt_dim1 = *npt;
2267 xpt_offset = xpt_dim1 + 1;
2268 xpt -= xpt_offset;
2269 --xopt;
2270 bmat_dim1 = *ndim;
2271 bmat_offset = bmat_dim1 + 1;
2272 bmat -= bmat_offset;
2273 --d__;
2274 --hcol;
2275 --gc;
2276 --gd;
2277 --s;
2278 --w;
2279
2280 /* Function Body */
2281 half = .5;
2282 one = 1.;
2283 zero = 0.;
2284 twopi = atan(one) * 8.;
2285 delsq = *delta * *delta;
2286 nptm = *npt - *n - 1;
2287
2288 /* Set the first NPT components of HCOL to the leading elements of the
2289 */
2290 /* KNEW-th column of H. */
2291
2292 iterc = 0;
2293 i__1 = *npt;
2294 for (k = 1; k <= i__1; ++k) {
2295 /* L10: */
2296 hcol[k] = zero;
2297 }
2298 i__1 = nptm;
2299 for (j = 1; j <= i__1; ++j) {
2300 temp = zmat[*knew + j * zmat_dim1];
2301 if (j < *idz) {
2302 temp = -temp;
2303 }
2304 i__2 = *npt;
2305 for (k = 1; k <= i__2; ++k) {
2306 /* L20: */
2307 hcol[k] += temp * zmat[k + j * zmat_dim1];
2308 }
2309 }
2310 *alpha = hcol[*knew];
2311
2312 /* Set the unscaled initial direction D. Form the gradient of LFUNC at
2313 */
2314 /* XOPT, and multiply D by the second derivative matrix of LFUNC. */
2315
2316 dd = zero;
2317 i__2 = *n;
2318 for (i__ = 1; i__ <= i__2; ++i__) {
2319 d__[i__] = xpt[*knew + i__ * xpt_dim1] - xopt[i__];
2320 gc[i__] = bmat[*knew + i__ * bmat_dim1];
2321 gd[i__] = zero;
2322 /* L30: */
2323 /* Computing 2nd power */
2324 d__1 = d__[i__];
2325 dd += d__1 * d__1;
2326 }
2327 i__2 = *npt;
2328 for (k = 1; k <= i__2; ++k) {
2329 temp = zero;
2330 sum = zero;
2331 i__1 = *n;
2332 for (j = 1; j <= i__1; ++j) {
2333 temp += xpt[k + j * xpt_dim1] * xopt[j];
2334 /* L40: */
2335 sum += xpt[k + j * xpt_dim1] * d__[j];
2336 }
2337 temp = hcol[k] * temp;
2338 sum = hcol[k] * sum;
2339 i__1 = *n;
2340 for (i__ = 1; i__ <= i__1; ++i__) {
2341 gc[i__] += temp * xpt[k + i__ * xpt_dim1];
2342 /* L50: */
2343 gd[i__] += sum * xpt[k + i__ * xpt_dim1];
2344 }
2345 }
2346
2347 /* Scale D and GD, with a sign change if required. Set S to another */
2348 /* vector in the initial two dimensional subspace. */
2349
2350 gg = zero;
2351 sp = zero;
2352 dhd = zero;
2353 i__1 = *n;
2354 for (i__ = 1; i__ <= i__1; ++i__) {
2355 /* Computing 2nd power */
2356 d__1 = gc[i__];
2357 gg += d__1 * d__1;
2358 sp += d__[i__] * gc[i__];
2359 /* L60: */
2360 dhd += d__[i__] * gd[i__];
2361 }
2362 scale = *delta / sqrt(dd);
2363 if (sp * dhd < zero) {
2364 scale = -scale;
2365 }
2366 temp = zero;
2367 if (sp * sp > dd * .99 * gg) {
2368 temp = one;
2369 }
2370 tau = scale * (abs(sp) + half * scale * abs(dhd));
2371 if (gg * delsq < tau * .01 * tau) {
2372 temp = one;
2373 }
2374 i__1 = *n;
2375 for (i__ = 1; i__ <= i__1; ++i__) {
2376 d__[i__] = scale * d__[i__];
2377 gd[i__] = scale * gd[i__];
2378 /* L70: */
2379 s[i__] = gc[i__] + temp * gd[i__];
2380 }
2381
2382 /* Begin the iteration by overwriting S with a vector that has the */
2383 /* required length and direction, except that termination occurs if */
2384 /* the given D and S are nearly parallel. */
2385
2386 L80:
2387 ++iterc;
2388 dd = zero;
2389 sp = zero;
2390 ss = zero;
2391 i__1 = *n;
2392 for (i__ = 1; i__ <= i__1; ++i__) {
2393 /* Computing 2nd power */
2394 d__1 = d__[i__];
2395 dd += d__1 * d__1;
2396 sp += d__[i__] * s[i__];
2397 /* L90: */
2398 /* Computing 2nd power */
2399 d__1 = s[i__];
2400 ss += d__1 * d__1;
2401 }
2402 temp = dd * ss - sp * sp;
2403 if (temp <= dd * 1e-8 * ss) {
2404 goto L160;
2405 }
2406 denom = sqrt(temp);
2407 i__1 = *n;
2408 for (i__ = 1; i__ <= i__1; ++i__) {
2409 s[i__] = (dd * s[i__] - sp * d__[i__]) / denom;
2410 /* L100: */
2411 w[i__] = zero;
2412 }
2413
2414 /* Calculate the coefficients of the objective function on the circle,
2415 */
2416 /* beginning with the multiplication of S by the second derivative matr
2417 ix.*/
2418
2419 i__1 = *npt;
2420 for (k = 1; k <= i__1; ++k) {
2421 sum = zero;
2422 i__2 = *n;
2423 for (j = 1; j <= i__2; ++j) {
2424 /* L110: */
2425 sum += xpt[k + j * xpt_dim1] * s[j];
2426 }
2427 sum = hcol[k] * sum;
2428 i__2 = *n;
2429 for (i__ = 1; i__ <= i__2; ++i__) {
2430 /* L120: */
2431 w[i__] += sum * xpt[k + i__ * xpt_dim1];
2432 }
2433 }
2434 cf1 = zero;
2435 cf2 = zero;
2436 cf3 = zero;
2437 cf4 = zero;
2438 cf5 = zero;
2439 i__2 = *n;
2440 for (i__ = 1; i__ <= i__2; ++i__) {
2441 cf1 += s[i__] * w[i__];
2442 cf2 += d__[i__] * gc[i__];
2443 cf3 += s[i__] * gc[i__];
2444 cf4 += d__[i__] * gd[i__];
2445 /* L130: */
2446 cf5 += s[i__] * gd[i__];
2447 }
2448 cf1 = half * cf1;
2449 cf4 = half * cf4 - cf1;
2450
2451 /* Seek the value of the angle that maximizes the modulus of TAU. */
2452
2453 taubeg = cf1 + cf2 + cf4;
2454 taumax = taubeg;
2455 tauold = taubeg;
2456 isave = 0;
2457 iu = 49;
2458 temp = twopi / (doublereal) (iu + 1);
2459 i__2 = iu;
2460 for (i__ = 1; i__ <= i__2; ++i__) {
2461 angle = (doublereal) i__ * temp;
2462 cth = cos(angle);
2463 sth = sin(angle);
2464 if(0&&DEBUG) fprintf(stderr," cos(angle=%.14g)=%.14g sin()=%.14g\n",angle,cth,sth) ;
2465 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
2466 if (abs(tau) > abs(taumax)) {
2467 taumax = tau;
2468 isave = i__;
2469 tempa = tauold;
2470 } else if (i__ == isave + 1) {
2471 tempb = tau;
2472 }
2473 /* L140: */
2474 tauold = tau;
2475 }
2476 if (isave == 0) {
2477 tempa = tau;
2478 }
2479 if (isave == iu) {
2480 tempb = taubeg;
2481 }
2482 step = zero;
2483 if (tempa != tempb) {
2484 tempa -= taumax;
2485 tempb -= taumax;
2486 step = half * (tempa - tempb) / (tempa + tempb);
2487 }
2488 angle = temp * ((doublereal) isave + step);
2489
2490 /* Calculate the new D and GD. Then test for convergence. */
2491
2492 cth = cos(angle);
2493 sth = sin(angle);
2494 if(0&&DEBUG) fprintf(stderr," cos(angle=%.14g)=%.14g sin()=%.14g\n",angle,cth,sth) ;
2495 tau = cf1 + (cf2 + cf4 * cth) * cth + (cf3 + cf5 * cth) * sth;
2496 i__2 = *n;
2497 for (i__ = 1; i__ <= i__2; ++i__) {
2498 d__[i__] = cth * d__[i__] + sth * s[i__];
2499 gd[i__] = cth * gd[i__] + sth * w[i__];
2500 /* L150: */
2501 s[i__] = gc[i__] + gd[i__];
2502 }
2503 if (abs(tau) <= abs(taubeg) * 1.1) {
2504 goto L160;
2505 }
2506 if (iterc < *n) {
2507 goto L80;
2508 }
2509 L160:
2510 return 0;
2511 } /* biglag_ */
2512
2513
2514 /*---------------------------------------------------------------------------*/
2515 #undef AFmax
2516 #undef AFmin
2517 #endif /* __POWELL_NEWUOA__ */
2518