praxis/praxis.c

   1 /* See README */
   2
   3 #include <math.h>
   4 #include "nlopt-util.h"
   5 #include "praxis.h"
   6
   7 /* Common Block Declarations */
   8
   9 struct global_s {
  10     double fx, ldt, dmin__;
  11     int nf, nl;
  12 };
  13
  14 struct q_s {
  15     double v[400]       /* was [20][20] */, q0[20], q1[20], qa, qb, qc, qd0,
  16             qd1, qf1;
  17 };
  18
  19 /* Table of constant values */
  20
  21 static int pow_ii(int x, int n) /* compute x^n, n >= 0 */
  22 {
  23      int p = 1;
  24      while (n > 0) {
  25           if (n & 1) p *= x;
  26           else { n >>= 1; x *= x; }
  27      }
  28      return p;
  29 }
  30
  31 static void minfit_(int *m, int *n, double *machep,
  32              double *tol, double *ab, double *q);
  33 static void min_(int n, int j, int nits, double *d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double *x, double *t, double *machep, double *h__, struct global_s *global_1, struct q_s *q_1);
  34 static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x, int *nf, struct q_s *q_1);
  35 static void sort_(int *m, int *n, double *d__, double *v);
  36 static void quad_(int *n, praxis_func f, void *f_data, double *x, double *t,
  37            double *machep, double *h__, struct global_s *global_1, struct q_s *q_1);
  38
  39 double praxis_(double *t0, double *machep, double *h0,
  40                int *n, double *x, praxis_func f, void *f_data)
  41 {
  42     /* System generated locals */
  43     int i__1, i__2, i__3;
  44     double ret_val, d__1;
  45
  46     /* Global */
  47     struct global_s global_1;
  48     struct q_s q_1;
  49
  50     /* Local variables */
  51     double d__[20], h__;
  52     int i__, j, k;
  53     double s, t, y[20], z__[20], f1;
  54     int k2;
  55     double m2, m4, t2, df, dn;
  56     int kl, ii;
  57     double sf;
  58     int kt;
  59     double sl;
  60     int im1, km1;
  61     double dni, lds;
  62     int ktm;
  63     double scbd;
  64     int idim;
  65     int illc;
  66     int klmk;
  67     double ldfac, large, small, value;
  68     double vlarge;
  69     double vsmall;
  70
  71 /*                             LAST MODIFIED 3/1/73 */
  72
  73 /*     PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(N,X) OF N VARIABLES */
  74 /*     USING THE PRINCIPAL AXIS METHOD.  THE GRADIENT OF THE FUNCTION IS */
  75 /*     NOT REQUIRED. */
  76
  77 /*     FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF */
  78 /*     "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT */
  79 /*     CALCULATING DERIVATIVES" BY RICHARD P BRENT. */
  80
  81 /*     THE PARAMETERS ARE: */
  82 /*     T0       IS A TOLERANCE.  PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X) */
  83 /*              SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN */
  84 /*              NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X). */
  85 /*     MACHEP   IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT */
  86 /*              1 + MACHEP > 1.  MACHEP SHOULD BE 16.**-13 (ABOUT */
  87 /*              2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360. */
  88 /*     H0       IS THE MAXIMUM STEP SIZE.  H0 SHOULD BE SET TO ABOUT THE */
  89 /*              MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM. */
  90 /*              (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF */
  91 /*              CONVERGENCE MAY BE SLOW.) */
  92 /*     N        (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH */
  93 /*              THE FUNCTION DEPENDS. */
  94 /*     X        IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF */
  95 /*              MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM. */
  96 /*     F(N,X)   IS THE FUNCTION TO BE MINIMIZED.  F SHOULD BE A REAL*8 */
  97 /*              FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM. */
  98 /*     THE APPROXIMATING QUADRATIC FORM IS */
  99 /*              Q(X') = F(N,X) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X) */
 100 /*     WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS */
 101 /*              INVERSE(V-TRANSPOSE) * D * INVERSE(V) */
 102 /*     (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY */
 103 /*     OF SECOND DIFFERENCES).  IF F HAS CONTINUOUS SECOND DERIVATIVES */
 104 /*     NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0. */
 105
 106 /*     IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */
 107 /*     TO ZERO. */
 108 /*     THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */
 109 /*     THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */
 110
 111
 112 /* .....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE */
 113 /*     LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND */
 114 /*     IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD. */
 115
 116
 117 /* .....INITIALIZATION..... */
 118 /*     MACHINE DEPENDENT NUMBERS: */
 119
 120     /* Parameter adjustments */
 121     --x;
 122
 123     /* Function Body */
 124     small = *machep * *machep;
 125     vsmall = small * small;
 126     large = 1. / small;
 127     vlarge = 1. / vsmall;
 128     m2 = sqrt(*machep);
 129     m4 = sqrt(m2);
 130
 131 /*     HEURISTIC NUMBERS: */
 132 /*     IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */
 133 /*     POSSIBLE), THEN SET SCBD=10.  OTHERWISE SET SCBD=1. */
 134 /*     IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE. */
 135 /*     OTHERWISE SET ILLC=FALSE. */
 136 /*     KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE */
 137 /*     ALGORITHM TERMINATES.  KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1 */
 138 /*     IS SATISFACTORY. */
 139
 140     scbd = 1.;
 141     illc = 0 /* false */;
 142     ktm = 1;
 143
 144     ldfac = .01;
 145     if (illc) {
 146         ldfac = .1;
 147     }
 148     kt = 0;
 149     global_1.nl = 0;
 150     global_1.nf = 1;
 151     global_1.fx = f(*n, &x[1], f_data);
 152     q_1.qf1 = global_1.fx;
 153     t = small + fabs(*t0);
 154     t2 = t;
 155     global_1.dmin__ = small;
 156     h__ = *h0;
 157     if (h__ < t * 100) {
 158         h__ = t * 100;
 159     }
 160     global_1.ldt = h__;
 161 /* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX..... */
 162     i__1 = *n;
 163     for (i__ = 1; i__ <= i__1; ++i__) {
 164         i__2 = *n;
 165         for (j = 1; j <= i__2; ++j) {
 166 /* L10: */
 167             q_1.v[i__ + j * 20 - 21] = 0.;
 168         }
 169 /* L20: */
 170         q_1.v[i__ + i__ * 20 - 21] = 1.;
 171     }
 172     d__[0] = 0.;
 173     q_1.qd0 = 0.;
 174     i__1 = *n;
 175     for (i__ = 1; i__ <= i__1; ++i__) {
 176         q_1.q0[i__ - 1] = x[i__];
 177 /* L30: */
 178         q_1.q1[i__ - 1] = x[i__];
 179     }
 180
 181 /* .....THE MAIN LOOP STARTS HERE..... */
 182 L40:
 183     sf = d__[0];
 184     d__[0] = 0.;
 185     s = 0.;
 186
 187 /* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */
 188 /*     FX MUST BE PASSED TO MIN BY VALUE. */
 189     value = global_1.fx;
 190     min_(*n, 1, 2, d__, &s, &value, 0, f,f_data, &x[1], &t,
 191             machep, &h__, &global_1, &q_1);
 192     if (s > 0.) {
 193         goto L50;
 194     }
 195     i__1 = *n;
 196     for (i__ = 1; i__ <= i__1; ++i__) {
 197 /* L45: */
 198         q_1.v[i__ - 1] = -q_1.v[i__ - 1];
 199     }
 200 L50:
 201     if (sf > d__[0] * .9 && sf * .9 < d__[0]) {
 202         goto L70;
 203     }
 204     i__1 = *n;
 205     for (i__ = 2; i__ <= i__1; ++i__) {
 206 /* L60: */
 207         d__[i__ - 1] = 0.;
 208     }
 209
 210 /* .....THE INNER LOOP STARTS HERE..... */
 211 L70:
 212     i__1 = *n;
 213     for (k = 2; k <= i__1; ++k) {
 214         i__2 = *n;
 215         for (i__ = 1; i__ <= i__2; ++i__) {
 216 /* L75: */
 217             y[i__ - 1] = x[i__];
 218         }
 219         sf = global_1.fx;
 220         if (kt > 0) {
 221             illc = 1 /* true */;
 222         }
 223 L80:
 224         kl = k;
 225         df = 0.;
 226
 227 /* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */
 228 /*     PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */
 229 /*     DISTRIBUTED IN (0,1). */
 230
 231         if (! illc) {
 232             goto L95;
 233         }
 234         i__2 = *n;
 235         for (i__ = 1; i__ <= i__2; ++i__) {
 236              s = (global_1.ldt * .1 + t2 * pow_ii(10, kt))
 237                   * nlopt_urand(-.5,.5);
 238              /* was: (random_(n) - .5); */
 239             z__[i__ - 1] = s;
 240             i__3 = *n;
 241             for (j = 1; j <= i__3; ++j) {
 242 /* L85: */
 243                 x[j] += s * q_1.v[j + i__ * 20 - 21];
 244             }
 245 /* L90: */
 246         }
 247         global_1.fx = (*f)(*n, &x[1], f_data);
 248         ++global_1.nf;
 249
 250 /* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N) */
 251
 252 L95:
 253         i__2 = *n;
 254         for (k2 = k; k2 <= i__2; ++k2) {
 255             sl = global_1.fx;
 256             s = 0.;
 257             value = global_1.fx;
 258             min_(*n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
 259                     x[1], &t, machep, &h__, &global_1, &q_1);
 260             if (illc) {
 261                 goto L97;
 262             }
 263             s = sl - global_1.fx;
 264             goto L99;
 265 L97:
 266 /* Computing 2nd power */
 267             d__1 = s + z__[k2 - 1];
 268             s = d__[k2 - 1] * (d__1 * d__1);
 269 L99:
 270             if (df > s) {
 271                 goto L105;
 272             }
 273             df = s;
 274             kl = k2;
 275 L105:
 276             ;
 277         }
 278         if (illc || df >= (d__1 = *machep * 100 * global_1.fx, fabs(d__1))) {
 279             goto L110;
 280         }
 281
 282 /* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */
 283 /*     ILLC=TRUE AND START THE INNER LOOP AGAIN..... */
 284
 285         illc = 1 /* true */;
 286         goto L80;
 287 L110:
 288
 289 /* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1) */
 290
 291         km1 = k - 1;
 292         i__2 = km1;
 293         for (k2 = 1; k2 <= i__2; ++k2) {
 294             s = 0.;
 295             value = global_1.fx;
 296             min_(*n, k2, 2, &d__[k2 - 1], &s, &value, 0, f,f_data, &
 297                     x[1], &t, machep, &h__, &global_1, &q_1);
 298 /* L120: */
 299         }
 300         f1 = global_1.fx;
 301         global_1.fx = sf;
 302         lds = 0.;
 303         i__2 = *n;
 304         for (i__ = 1; i__ <= i__2; ++i__) {
 305             sl = x[i__];
 306             x[i__] = y[i__ - 1];
 307             sl -= y[i__ - 1];
 308             y[i__ - 1] = sl;
 309 /* L130: */
 310             lds += sl * sl;
 311         }
 312         lds = sqrt(lds);
 313         if (lds <= small) {
 314             goto L160;
 315         }
 316
 317 /* .....DISCARD DIRECTION V(*,KL). */
 318 /*     IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE" */
 319 /*     DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE..... */
 320
 321         klmk = kl - k;
 322         if (klmk < 1) {
 323             goto L141;
 324         }
 325         i__2 = klmk;
 326         for (ii = 1; ii <= i__2; ++ii) {
 327             i__ = kl - ii;
 328             i__3 = *n;
 329             for (j = 1; j <= i__3; ++j) {
 330 /* L135: */
 331                 q_1.v[j + (i__ + 1) * 20 - 21] = q_1.v[j + i__ * 20 - 21];
 332             }
 333 /* L140: */
 334             d__[i__] = d__[i__ - 1];
 335         }
 336 L141:
 337         d__[k - 1] = 0.;
 338         i__2 = *n;
 339         for (i__ = 1; i__ <= i__2; ++i__) {
 340 /* L145: */
 341             q_1.v[i__ + k * 20 - 21] = y[i__ - 1] / lds;
 342         }
 343
 344 /* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */
 345 /*     THE NORMALIZED VECTOR:  (NEW X) - (0LD X)..... */
 346
 347         value = f1;
 348         min_(*n, k, 4, &d__[k - 1], &lds, &value, 1, f,f_data, &x[1],
 349                  &t, machep, &h__, &global_1, &q_1);
 350         if (lds > 0.) {
 351             goto L160;
 352         }
 353         lds = -lds;
 354         i__2 = *n;
 355         for (i__ = 1; i__ <= i__2; ++i__) {
 356 /* L150: */
 357             q_1.v[i__ + k * 20 - 21] = -q_1.v[i__ + k * 20 - 21];
 358         }
 359 L160:
 360         global_1.ldt = ldfac * global_1.ldt;
 361         if (global_1.ldt < lds) {
 362             global_1.ldt = lds;
 363         }
 364         t2 = 0.;
 365         i__2 = *n;
 366         for (i__ = 1; i__ <= i__2; ++i__) {
 367 /* L165: */
 368 /* Computing 2nd power */
 369             d__1 = x[i__];
 370             t2 += d__1 * d__1;
 371         }
 372         t2 = m2 * sqrt(t2) + t;
 373
 374 /* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
 375 /*     INNER LOOP EXCEEDS HALF THE TOLERANCE..... */
 376
 377         if (global_1.ldt > t2 * .5f) {
 378             kt = -1;
 379         }
 380         ++kt;
 381         if (kt > ktm) {
 382             goto L400;
 383         }
 384 /* L170: */
 385     }
 386 /* .....THE INNER LOOP ENDS HERE. */
 387
 388 /*     TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */
 389
 390 /* L171: */
 391     quad_(n, f,f_data, &x[1], &t, machep, &h__, &global_1, &q_1);
 392     dn = 0.;
 393     i__1 = *n;
 394     for (i__ = 1; i__ <= i__1; ++i__) {
 395         d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
 396         if (dn < d__[i__ - 1]) {
 397             dn = d__[i__ - 1];
 398         }
 399 /* L175: */
 400     }
 401     i__1 = *n;
 402     for (j = 1; j <= i__1; ++j) {
 403         s = d__[j - 1] / dn;
 404         i__2 = *n;
 405         for (i__ = 1; i__ <= i__2; ++i__) {
 406 /* L180: */
 407             q_1.v[i__ + j * 20 - 21] = s * q_1.v[i__ + j * 20 - 21];
 408         }
 409     }
 410
 411 /* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */
 412
 413     if (scbd <= 1.) {
 414         goto L200;
 415     }
 416     s = vlarge;
 417     i__2 = *n;
 418     for (i__ = 1; i__ <= i__2; ++i__) {
 419         sl = 0.;
 420         i__1 = *n;
 421         for (j = 1; j <= i__1; ++j) {
 422 /* L182: */
 423             sl += q_1.v[i__ + j * 20 - 21] * q_1.v[i__ + j * 20 - 21];
 424         }
 425         z__[i__ - 1] = sqrt(sl);
 426         if (z__[i__ - 1] < m4) {
 427             z__[i__ - 1] = m4;
 428         }
 429         if (s > z__[i__ - 1]) {
 430             s = z__[i__ - 1];
 431         }
 432 /* L185: */
 433     }
 434     i__2 = *n;
 435     for (i__ = 1; i__ <= i__2; ++i__) {
 436         sl = s / z__[i__ - 1];
 437         z__[i__ - 1] = 1. / sl;
 438         if (z__[i__ - 1] <= scbd) {
 439             goto L189;
 440         }
 441         sl = 1. / scbd;
 442         z__[i__ - 1] = scbd;
 443 L189:
 444         i__1 = *n;
 445         for (j = 1; j <= i__1; ++j) {
 446 /* L190: */
 447             q_1.v[i__ + j * 20 - 21] = sl * q_1.v[i__ + j * 20 - 21];
 448         }
 449 /* L195: */
 450     }
 451
 452 /* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */
 453 /*     THE MAIN LOOP. */
 454 /*     FIRST TRANSPOSE V FOR MINFIT: */
 455
 456 L200:
 457     i__2 = *n;
 458     for (i__ = 2; i__ <= i__2; ++i__) {
 459         im1 = i__ - 1;
 460         i__1 = im1;
 461         for (j = 1; j <= i__1; ++j) {
 462             s = q_1.v[i__ + j * 20 - 21];
 463             q_1.v[i__ + j * 20 - 21] = q_1.v[j + i__ * 20 - 21];
 464 /* L210: */
 465             q_1.v[j + i__ * 20 - 21] = s;
 466         }
 467 /* L220: */
 468     }
 469
 470 /* .....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V. */
 471 /*     THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE */
 472 /*     APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */
 473 /*     NUMBER..... */
 474
 475     minfit_(&idim, n, machep, &vsmall, q_1.v, d__);
 476
 477 /* .....UNSCALE THE AXES..... */
 478
 479     if (scbd <= 1.) {
 480         goto L250;
 481     }
 482     i__2 = *n;
 483     for (i__ = 1; i__ <= i__2; ++i__) {
 484         s = z__[i__ - 1];
 485         i__1 = *n;
 486         for (j = 1; j <= i__1; ++j) {
 487 /* L225: */
 488             q_1.v[i__ + j * 20 - 21] = s * q_1.v[i__ + j * 20 - 21];
 489         }
 490 /* L230: */
 491     }
 492     i__2 = *n;
 493     for (i__ = 1; i__ <= i__2; ++i__) {
 494         s = 0.;
 495         i__1 = *n;
 496         for (j = 1; j <= i__1; ++j) {
 497 /* L235: */
 498 /* Computing 2nd power */
 499             d__1 = q_1.v[j + i__ * 20 - 21];
 500             s += d__1 * d__1;
 501         }
 502         s = sqrt(s);
 503         d__[i__ - 1] = s * d__[i__ - 1];
 504         s = 1 / s;
 505         i__1 = *n;
 506         for (j = 1; j <= i__1; ++j) {
 507 /* L240: */
 508             q_1.v[j + i__ * 20 - 21] = s * q_1.v[j + i__ * 20 - 21];
 509         }
 510 /* L245: */
 511     }
 512
 513 L250:
 514     i__2 = *n;
 515     for (i__ = 1; i__ <= i__2; ++i__) {
 516         dni = dn * d__[i__ - 1];
 517         if (dni > large) {
 518             goto L265;
 519         }
 520         if (dni < small) {
 521             goto L260;
 522         }
 523         d__[i__ - 1] = 1 / (dni * dni);
 524         goto L270;
 525 L260:
 526         d__[i__ - 1] = vlarge;
 527         goto L270;
 528 L265:
 529         d__[i__ - 1] = vsmall;
 530 L270:
 531         ;
 532     }
 533
 534 /* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */
 535
 536     sort_(&idim, n, d__, q_1.v);
 537     global_1.dmin__ = d__[*n - 1];
 538     if (global_1.dmin__ < small) {
 539         global_1.dmin__ = small;
 540     }
 541     illc = 0 /* false */;
 542     if (m2 * d__[0] > global_1.dmin__) {
 543         illc = 1 /* true */;
 544     }
 545
 546 /* .....THE MAIN LOOP ENDS HERE..... */
 547
 548     goto L40;
 549
 550 /* .....RETURN..... */
 551
 552 L400:
 553     ret_val = global_1.fx;
 554     return ret_val;
 555 } /* praxis_ */
 556
 557 static void minfit_(int *m, int *n, double *machep,
 558         double *tol, double *ab, double *q)
 559 {
 560     /* System generated locals */
 561     int ab_dim1, ab_offset, i__1, i__2, i__3;
 562     double d__1, d__2;
 563
 564     /* Local variables */
 565     double c__, e[20], f, g, h__;
 566     int i__, j, k, l;
 567     double s, x, y, z__;
 568     int l2, ii, kk, kt, ll2, lp1;
 569     double eps, temp;
 570
 571 /* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */
 572 /*   RESTRICTED TO M=N,P=0. */
 573 /*   THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS */
 574 /*   OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V, */
 575 /*   WHERE U IS ANOTHER ORTHOGONAL MATRIX. */
 576 /* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */
 577     /* Parameter adjustments */
 578     --q;
 579     ab_dim1 = *m;
 580     ab_offset = 1 + ab_dim1;
 581     ab -= ab_offset;
 582
 583     /* Function Body */
 584     if (*n == 1) {
 585         goto L200;
 586     }
 587     eps = *machep;
 588     g = 0.;
 589     x = 0.;
 590     i__1 = *n;
 591     for (i__ = 1; i__ <= i__1; ++i__) {
 592         e[i__ - 1] = g;
 593         s = 0.;
 594         l = i__ + 1;
 595         i__2 = *n;
 596         for (j = i__; j <= i__2; ++j) {
 597 /* L1: */
 598 /* Computing 2nd power */
 599             d__1 = ab[j + i__ * ab_dim1];
 600             s += d__1 * d__1;
 601         }
 602         g = 0.;
 603         if (s < *tol) {
 604             goto L4;
 605         }
 606         f = ab[i__ + i__ * ab_dim1];
 607         g = sqrt(s);
 608         if (f >= 0.) {
 609             g = -g;
 610         }
 611         h__ = f * g - s;
 612         ab[i__ + i__ * ab_dim1] = f - g;
 613         if (l > *n) {
 614             goto L4;
 615         }
 616         i__2 = *n;
 617         for (j = l; j <= i__2; ++j) {
 618             f = 0.;
 619             i__3 = *n;
 620             for (k = i__; k <= i__3; ++k) {
 621 /* L2: */
 622                 f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1];
 623             }
 624             f /= h__;
 625             i__3 = *n;
 626             for (k = i__; k <= i__3; ++k) {
 627 /* L3: */
 628                 ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
 629             }
 630         }
 631 L4:
 632         q[i__] = g;
 633         s = 0.;
 634         if (i__ == *n) {
 635             goto L6;
 636         }
 637         i__3 = *n;
 638         for (j = l; j <= i__3; ++j) {
 639 /* L5: */
 640             s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
 641         }
 642 L6:
 643         g = 0.;
 644         if (s < *tol) {
 645             goto L10;
 646         }
 647         if (i__ == *n) {
 648             goto L16;
 649         }
 650         f = ab[i__ + (i__ + 1) * ab_dim1];
 651 L16:
 652         g = sqrt(s);
 653         if (f >= 0.) {
 654             g = -g;
 655         }
 656         h__ = f * g - s;
 657         if (i__ == *n) {
 658             goto L10;
 659         }
 660         ab[i__ + (i__ + 1) * ab_dim1] = f - g;
 661         i__3 = *n;
 662         for (j = l; j <= i__3; ++j) {
 663 /* L7: */
 664             e[j - 1] = ab[i__ + j * ab_dim1] / h__;
 665         }
 666         i__3 = *n;
 667         for (j = l; j <= i__3; ++j) {
 668             s = 0.;
 669             i__2 = *n;
 670             for (k = l; k <= i__2; ++k) {
 671 /* L8: */
 672                 s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1];
 673             }
 674             i__2 = *n;
 675             for (k = l; k <= i__2; ++k) {
 676 /* L9: */
 677                 ab[j + k * ab_dim1] += s * e[k - 1];
 678             }
 679         }
 680 L10:
 681         y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2));
 682 /* L11: */
 683         if (y > x) {
 684             x = y;
 685         }
 686     }
 687 /* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
 688     ab[*n + *n * ab_dim1] = 1.;
 689     g = e[*n - 1];
 690     l = *n;
 691     i__1 = *n;
 692     for (ii = 2; ii <= i__1; ++ii) {
 693         i__ = *n - ii + 1;
 694         if (g == 0.) {
 695             goto L23;
 696         }
 697         h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
 698         i__2 = *n;
 699         for (j = l; j <= i__2; ++j) {
 700 /* L20: */
 701             ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
 702         }
 703         i__2 = *n;
 704         for (j = l; j <= i__2; ++j) {
 705             s = 0.;
 706             i__3 = *n;
 707             for (k = l; k <= i__3; ++k) {
 708 /* L21: */
 709                 s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1];
 710             }
 711             i__3 = *n;
 712             for (k = l; k <= i__3; ++k) {
 713 /* L22: */
 714                 ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
 715             }
 716         }
 717 L23:
 718         i__3 = *n;
 719         for (j = l; j <= i__3; ++j) {
 720             ab[i__ + j * ab_dim1] = 0.;
 721 /* L24: */
 722             ab[j + i__ * ab_dim1] = 0.;
 723         }
 724         ab[i__ + i__ * ab_dim1] = 1.;
 725         g = e[i__ - 1];
 726 /* L25: */
 727         l = i__;
 728     }
 729 /* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
 730 /* L100: */
 731     eps *= x;
 732     i__1 = *n;
 733     for (kk = 1; kk <= i__1; ++kk) {
 734         k = *n - kk + 1;
 735         kt = 0;
 736 L101:
 737         ++kt;
 738         if (kt <= 30) {
 739             goto L102;
 740         }
 741         e[k - 1] = 0.;
 742         /* fprintf(stderr, "QR failed\n"); */
 743 L102:
 744         i__3 = k;
 745         for (ll2 = 1; ll2 <= i__3; ++ll2) {
 746             l2 = k - ll2 + 1;
 747             l = l2;
 748             if ((d__1 = e[l - 1], fabs(d__1)) <= eps) {
 749                 goto L120;
 750             }
 751             if (l == 1) {
 752                 goto L103;
 753             }
 754             if ((d__1 = q[l - 1], fabs(d__1)) <= eps) {
 755                 goto L110;
 756             }
 757 L103:
 758             ;
 759         }
 760 /* ...CANCELLATION OF E(L) IF L>1... */
 761 L110:
 762         c__ = 0.;
 763         s = 1.;
 764         i__3 = k;
 765         for (i__ = l; i__ <= i__3; ++i__) {
 766             f = s * e[i__ - 1];
 767             e[i__ - 1] = c__ * e[i__ - 1];
 768             if (fabs(f) <= eps) {
 769                 goto L120;
 770             }
 771             g = q[i__];
 772 /* ...Q(I) = H = DSQRT(G*G + F*F)... */
 773             if (fabs(f) < fabs(g)) {
 774                 goto L113;
 775             }
 776             if (f != 0.) {
 777                 goto L112;
 778             } else {
 779                 goto L111;
 780             }
 781 L111:
 782             h__ = 0.;
 783             goto L114;
 784 L112:
 785 /* Computing 2nd power */
 786             d__1 = g / f;
 787             h__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 788             goto L114;
 789 L113:
 790 /* Computing 2nd power */
 791             d__1 = f / g;
 792             h__ = fabs(g) * sqrt(d__1 * d__1 + 1);
 793 L114:
 794             q[i__] = h__;
 795             if (h__ != 0.) {
 796                 goto L115;
 797             }
 798             g = 1.;
 799             h__ = 1.;
 800 L115:
 801             c__ = g / h__;
 802 /* L116: */
 803             s = -f / h__;
 804         }
 805 /* ...TEST FOR CONVERGENCE... */
 806 L120:
 807         z__ = q[k];
 808         if (l == k) {
 809             goto L140;
 810         }
 811 /* ...SHIFT FROM BOTTOM 2*2 MINOR... */
 812         x = q[l];
 813         y = q[k - 1];
 814         g = e[k - 2];
 815         h__ = e[k - 1];
 816         f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y);
 817         g = sqrt(f * f + 1.);
 818         temp = f - g;
 819         if (f >= 0.) {
 820             temp = f + g;
 821         }
 822         f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x;
 823 /* ...NEXT QR TRANSFORMATION... */
 824         c__ = 1.;
 825         s = 1.;
 826         lp1 = l + 1;
 827         if (lp1 > k) {
 828             goto L133;
 829         }
 830         i__3 = k;
 831         for (i__ = lp1; i__ <= i__3; ++i__) {
 832             g = e[i__ - 1];
 833             y = q[i__];
 834             h__ = s * g;
 835             g *= c__;
 836             if (fabs(f) < fabs(h__)) {
 837                 goto L123;
 838             }
 839             if (f != 0.) {
 840                 goto L122;
 841             } else {
 842                 goto L121;
 843             }
 844 L121:
 845             z__ = 0.;
 846             goto L124;
 847 L122:
 848 /* Computing 2nd power */
 849             d__1 = h__ / f;
 850             z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 851             goto L124;
 852 L123:
 853 /* Computing 2nd power */
 854             d__1 = f / h__;
 855             z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
 856 L124:
 857             e[i__ - 2] = z__;
 858             if (z__ != 0.) {
 859                 goto L125;
 860             }
 861             f = 1.;
 862             z__ = 1.;
 863 L125:
 864             c__ = f / z__;
 865             s = h__ / z__;
 866             f = x * c__ + g * s;
 867             g = -x * s + g * c__;
 868             h__ = y * s;
 869             y *= c__;
 870             i__2 = *n;
 871             for (j = 1; j <= i__2; ++j) {
 872                 x = ab[j + (i__ - 1) * ab_dim1];
 873                 z__ = ab[j + i__ * ab_dim1];
 874                 ab[j + (i__ - 1) * ab_dim1] = x * c__ + z__ * s;
 875 /* L126: */
 876                 ab[j + i__ * ab_dim1] = -x * s + z__ * c__;
 877             }
 878             if (fabs(f) < fabs(h__)) {
 879                 goto L129;
 880             }
 881             if (f != 0.) {
 882                 goto L128;
 883             } else {
 884                 goto L127;
 885             }
 886 L127:
 887             z__ = 0.;
 888             goto L130;
 889 L128:
 890 /* Computing 2nd power */
 891             d__1 = h__ / f;
 892             z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
 893             goto L130;
 894 L129:
 895 /* Computing 2nd power */
 896             d__1 = f / h__;
 897             z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
 898 L130:
 899             q[i__ - 1] = z__;
 900             if (z__ != 0.) {
 901                 goto L131;
 902             }
 903             f = 1.;
 904             z__ = 1.;
 905 L131:
 906             c__ = f / z__;
 907             s = h__ / z__;
 908             f = c__ * g + s * y;
 909 /* L132: */
 910             x = -s * g + c__ * y;
 911         }
 912 L133:
 913         e[l - 1] = 0.;
 914         e[k - 1] = f;
 915         q[k] = x;
 916         goto L101;
 917 /* ...CONVERGENCE:  Q(K) IS MADE NON-NEGATIVE... */
 918 L140:
 919         if (z__ >= 0.) {
 920             goto L150;
 921         }
 922         q[k] = -z__;
 923         i__3 = *n;
 924         for (j = 1; j <= i__3; ++j) {
 925 /* L141: */
 926             ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
 927         }
 928 L150:
 929         ;
 930     }
 931     return;
 932 L200:
 933     q[1] = ab[ab_dim1 + 1];
 934     ab[ab_dim1 + 1] = 1.;
 935 } /* minfit_ */
 936
 937 static void min_(int n, int j, int nits, double *
 938         d2, double *x1, double *f1, int fk, praxis_func f, void *f_data, double *
 939         x, double *t, double *machep, double *h__, struct global_s *global_1, struct q_s *q_1)
 940 {
 941     /* System generated locals */
 942     int i__1;
 943     double d__1, d__2;
 944
 945     /* Local variables */
 946     int i__, k;
 947     double s, d1, f0, f2, m2, m4, t2, x2, fm;
 948     int dz;
 949     double xm, sf1, sx1;
 950     double temp, small;
 951
 952 /* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */
 953 /*   J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */
 954 /*   DEFINED BY Q0,Q1,X. */
 955 /*   D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */
 956 /*   ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */
 957 /*   ALONG V(*,J) (OR, IF J=0, A CURVE).  ON RETURN, X1 IS THE DISTANCE */
 958 /*   FOUND. */
 959 /*   IF FK=.TRUE., THEN F1 IS FLIN(X1).  OTHERWISE X1 AND F1 ARE IGNORED */
 960 /*   ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */
 961 /*   NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */
 962 /*   THE INTERVAL. */
 963     /* Parameter adjustments */
 964     --x;
 965
 966     /* Function Body */
 967 /* Computing 2nd power */
 968     d__1 = *machep;
 969     small = d__1 * d__1;
 970     m2 = sqrt(*machep);
 971     m4 = sqrt(m2);
 972     sf1 = *f1;
 973     sx1 = *x1;
 974     k = 0;
 975     xm = 0.;
 976     fm = global_1->fx;
 977     f0 = global_1->fx;
 978     dz = *d2 < *machep;
 979 /* ...FIND THE STEP SIZE... */
 980     s = 0.;
 981     i__1 = n;
 982     for (i__ = 1; i__ <= i__1; ++i__) {
 983 /* L1: */
 984 /* Computing 2nd power */
 985         d__1 = x[i__];
 986         s += d__1 * d__1;
 987     }
 988     s = sqrt(s);
 989     temp = *d2;
 990     if (dz) {
 991         temp = global_1->dmin__;
 992     }
 993     t2 = m4 * sqrt(fabs(global_1->fx) / temp + s * global_1->ldt) + m2 *
 994             global_1->ldt;
 995     s = m4 * s + *t;
 996     if (dz && t2 > s) {
 997         t2 = s;
 998     }
 999     t2 = t2 > small ? t2 : small;
1000 /* Computing MIN */
1001     d__1 = t2, d__2 = *h__ * .01;
1002     t2 = d__1 < d__2 ? d__1 : d__2;
1003     if (! (fk) || *f1 > fm) {
1004         goto L2;
1005     }
1006     xm = *x1;
1007     fm = *f1;
1008 L2:
1009     if (fk && fabs(*x1) >= t2) {
1010         goto L3;
1011     }
1012     temp = 1.;
1013     if (*x1 < 0.) {
1014         temp = -1.;
1015     }
1016     *x1 = temp * t2;
1017     *f1 = flin_(n, j, x1, f,f_data, &x[1], &global_1->nf, q_1);
1018 L3:
1019     if (*f1 > fm) {
1020         goto L4;
1021     }
1022     xm = *x1;
1023     fm = *f1;
1024 L4:
1025     if (! dz) {
1026         goto L6;
1027     }
1028 /* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE... */
1029     x2 = -(*x1);
1030     if (f0 >= *f1) {
1031         x2 = *x1 * 2.;
1032     }
1033     f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1);
1034     if (f2 > fm) {
1035         goto L5;
1036     }
1037     xm = x2;
1038     fm = f2;
1039 L5:
1040     *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2));
1041 /* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */
1042 L6:
1043     d1 = (*f1 - f0) / *x1 - *x1 * *d2;
1044     dz = 1 /* true */;
1045 /* ...PREDICT THE MINIMUM... */
1046     if (*d2 > small) {
1047         goto L7;
1048     }
1049     x2 = *h__;
1050     if (d1 >= 0.) {
1051         x2 = -x2;
1052     }
1053     goto L8;
1054 L7:
1055     x2 = d1 * -.5 / *d2;
1056 L8:
1057     if (fabs(x2) <= *h__) {
1058         goto L11;
1059     }
1060     if (x2 <= 0.) {
1061         goto L9;
1062     } else {
1063         goto L10;
1064     }
1065 L9:
1066     x2 = -(*h__);
1067     goto L11;
1068 L10:
1069     x2 = *h__;
1070 /* ...EVALUATE F AT THE PREDICTED MINIMUM... */
1071 L11:
1072     f2 = flin_(n, j, &x2, f,f_data, &x[1], &global_1->nf, q_1);
1073     if (k >= nits || f2 <= f0) {
1074         goto L12;
1075     }
1076 /* ...NO SUCCESS, SO TRY AGAIN... */
1077     ++k;
1078     if (f0 < *f1 && *x1 * x2 > 0.) {
1079         goto L4;
1080     }
1081     x2 *= .5;
1082     goto L11;
1083 /* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */
1084 L12:
1085     ++global_1->nl;
1086     if (f2 <= fm) {
1087         goto L13;
1088     }
1089     x2 = xm;
1090     goto L14;
1091 L13:
1092     fm = f2;
1093 /* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */
1094 L14:
1095     if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small) {
1096         goto L15;
1097     }
1098     *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2));
1099     goto L16;
1100 L15:
1101     if (k > 0) {
1102         *d2 = 0.;
1103     }
1104 L16:
1105     if (*d2 <= small) {
1106         *d2 = small;
1107     }
1108     *x1 = x2;
1109     global_1->fx = fm;
1110     if (sf1 >= global_1->fx) {
1111         goto L17;
1112     }
1113     global_1->fx = sf1;
1114     *x1 = sx1;
1115 /* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
1116 L17:
1117     if (j == 0) {
1118         return;
1119     }
1120     i__1 = n;
1121     for (i__ = 1; i__ <= i__1; ++i__) {
1122 /* L18: */
1123         x[i__] += *x1 * q_1->v[i__ + j * 20 - 21];
1124     }
1125 } /* min_ */
1126
1127 static double flin_(int n, int j, double *l, praxis_func f, void *f_data, double *x,
1128          int *nf, struct q_s *q_1)
1129 {
1130     /* System generated locals */
1131     int i__1;
1132     double ret_val;
1133
1134     /* Local variables */
1135     int i__;
1136     double t[20];
1137
1138 /* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
1139 /*   BY THE SUBROUTINE MIN... */
1140     /* Parameter adjustments */
1141     --x;
1142
1143     /* Function Body */
1144     if (j == 0) {
1145         goto L2;
1146     }
1147 /* ...THE SEARCH IS LINEAR... */
1148     i__1 = n;
1149     for (i__ = 1; i__ <= i__1; ++i__) {
1150 /* L1: */
1151         t[i__ - 1] = x[i__] + *l * q_1->v[i__ + j * 20 - 21];
1152     }
1153     goto L4;
1154 /* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
1155 L2:
1156     q_1->qa = *l * (*l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1));
1157     q_1->qb = (*l + q_1->qd0) * (q_1->qd1 - *l) / (q_1->qd0 * q_1->qd1);
1158     q_1->qc = *l * (*l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1));
1159     i__1 = n;
1160     for (i__ = 1; i__ <= i__1; ++i__) {
1161 /* L3: */
1162         t[i__ - 1] = q_1->qa * q_1->q0[i__ - 1] + q_1->qb * x[i__] + q_1->qc *
1163                 q_1->q1[i__ - 1];
1164     }
1165 /* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */
1166 L4:
1167     ++(*nf);
1168     ret_val = f(n, t, f_data);
1169     return ret_val;
1170 } /* flin_ */
1171
1172 static void sort_(int *m, int *n, double *d__, double *v)
1173 {
1174     /* System generated locals */
1175     int v_dim1, v_offset, i__1, i__2;
1176
1177     /* Local variables */
1178     int i__, j, k;
1179     double s;
1180     int ip1, nm1;
1181
1182 /* ...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE */
1183 /*   CORRESPONDING COLUMNS OF V(N,N). */
1184 /*   M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM. */
1185     /* Parameter adjustments */
1186     v_dim1 = *m;
1187     v_offset = 1 + v_dim1;
1188     v -= v_offset;
1189     --d__;
1190
1191     /* Function Body */
1192     if (*n == 1) {
1193         return;
1194     }
1195     nm1 = *n - 1;
1196     i__1 = nm1;
1197     for (i__ = 1; i__ <= i__1; ++i__) {
1198         k = i__;
1199         s = d__[i__];
1200         ip1 = i__ + 1;
1201         i__2 = *n;
1202         for (j = ip1; j <= i__2; ++j) {
1203             if (d__[j] <= s) {
1204                 goto L1;
1205             }
1206             k = j;
1207             s = d__[j];
1208 L1:
1209             ;
1210         }
1211         if (k <= i__) {
1212             goto L3;
1213         }
1214         d__[k] = d__[i__];
1215         d__[i__] = s;
1216         i__2 = *n;
1217         for (j = 1; j <= i__2; ++j) {
1218             s = v[j + i__ * v_dim1];
1219             v[j + i__ * v_dim1] = v[j + k * v_dim1];
1220 /* L2: */
1221             v[j + k * v_dim1] = s;
1222         }
1223 L3:
1224         ;
1225     }
1226 } /* sort_ */
1227
1228 static void quad_(int *n, praxis_func f, void *f_data, double *x, double *t,
1229                   double *machep, double *h__, struct global_s *global_1, struct q_s *q_1)
1230 {
1231     /* System generated locals */
1232     int i__1;
1233     double d__1;
1234
1235     /* Local variables */
1236     int i__;
1237     double l, s;
1238     double value;
1239
1240 /* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X... */
1241     /* Parameter adjustments */
1242     --x;
1243
1244     /* Function Body */
1245     s = global_1->fx;
1246     global_1->fx = q_1->qf1;
1247     q_1->qf1 = s;
1248     q_1->qd1 = 0.;
1249     i__1 = *n;
1250     for (i__ = 1; i__ <= i__1; ++i__) {
1251         s = x[i__];
1252         l = q_1->q1[i__ - 1];
1253         x[i__] = l;
1254         q_1->q1[i__ - 1] = s;
1255 /* L1: */
1256 /* Computing 2nd power */
1257         d__1 = s - l;
1258         q_1->qd1 += d__1 * d__1;
1259     }
1260     q_1->qd1 = sqrt(q_1->qd1);
1261     l = q_1->qd1;
1262     s = 0.;
1263     if (q_1->qd0 <= 0. || q_1->qd1 <= 0. || global_1->nl < *n * 3 * *n) {
1264         goto L2;
1265     }
1266     value = q_1->qf1;
1267     min_(*n, 0, 2, &s, &l, &value, 1, f,f_data, &x[1], t, machep,
1268             h__, global_1, q_1);
1269     q_1->qa = l * (l - q_1->qd1) / (q_1->qd0 * (q_1->qd0 + q_1->qd1));
1270     q_1->qb = (l + q_1->qd0) * (q_1->qd1 - l) / (q_1->qd0 * q_1->qd1);
1271     q_1->qc = l * (l + q_1->qd0) / (q_1->qd1 * (q_1->qd0 + q_1->qd1));
1272     goto L3;
1273 L2:
1274     global_1->fx = q_1->qf1;
1275     q_1->qa = 0.;
1276     q_1->qb = q_1->qa;
1277     q_1->qc = 1.;
1278 L3:
1279     q_1->qd0 = q_1->qd1;
1280     i__1 = *n;
1281     for (i__ = 1; i__ <= i__1; ++i__) {
1282         s = q_1->q0[i__ - 1];
1283         q_1->q0[i__ - 1] = x[i__];
1284 /* L4: */
1285         x[i__] = q_1->qa * s + q_1->qb * x[i__] + q_1->qc * q_1->q1[i__ - 1];
1286     }
1287 } /* quad_ */