3rdparty/lapack/dsterf.c

   1 /* dsterf.f -- translated by f2c (version 20061008).
   2    You must link the resulting object file with libf2c:
   3         on Microsoft Windows system, link with libf2c.lib;
   4         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
   5         or, if you install libf2c.a in a standard place, with -lf2c -lm
   6         -- in that order, at the end of the command line, as in
   7                 cc *.o -lf2c -lm
   8         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
   9
  10                 http://www.netlib.org/f2c/libf2c.zip
  11 */
  12
  13 #include "clapack.h"
  14
  15
  16 /* Table of constant values */
  17
  18 static integer c__0 = 0;
  19 static integer c__1 = 1;
  20 static doublereal c_b32 = 1.;
  21
  22 /* Subroutine */ int dsterf_(integer *n, doublereal *d__, doublereal *e,
  23         integer *info)
  24 {
  25     /* System generated locals */
  26     integer i__1;
  27     doublereal d__1, d__2, d__3;
  28
  29     /* Builtin functions */
  30     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
  31
  32     /* Local variables */
  33     doublereal c__;
  34     integer i__, l, m;
  35     doublereal p, r__, s;
  36     integer l1;
  37     doublereal bb, rt1, rt2, eps, rte;
  38     integer lsv;
  39     doublereal eps2, oldc;
  40     integer lend, jtot;
  41     extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal
  42             *, doublereal *, doublereal *);
  43     doublereal gamma, alpha, sigma, anorm;
  44     extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *);
  45     integer iscale;
  46     extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
  47             doublereal *, doublereal *, integer *, integer *, doublereal *,
  48             integer *, integer *);
  49     doublereal oldgam, safmin;
  50     extern /* Subroutine */ int xerbla_(char *, integer *);
  51     doublereal safmax;
  52     extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
  53     extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *,
  54             integer *);
  55     integer lendsv;
  56     doublereal ssfmin;
  57     integer nmaxit;
  58     doublereal ssfmax;
  59
  60
  61 /*  -- LAPACK routine (version 3.2) -- */
  62 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  63 /*     November 2006 */
  64
  65 /*     .. Scalar Arguments .. */
  66 /*     .. */
  67 /*     .. Array Arguments .. */
  68 /*     .. */
  69
  70 /*  Purpose */
  71 /*  ======= */
  72
  73 /*  DSTERF computes all eigenvalues of a symmetric tridiagonal matrix */
  74 /*  using the Pal-Walker-Kahan variant of the QL or QR algorithm. */
  75
  76 /*  Arguments */
  77 /*  ========= */
  78
  79 /*  N       (input) INTEGER */
  80 /*          The order of the matrix.  N >= 0. */
  81
  82 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
  83 /*          On entry, the n diagonal elements of the tridiagonal matrix. */
  84 /*          On exit, if INFO = 0, the eigenvalues in ascending order. */
  85
  86 /*  E       (input/output) DOUBLE PRECISION array, dimension (N-1) */
  87 /*          On entry, the (n-1) subdiagonal elements of the tridiagonal */
  88 /*          matrix. */
  89 /*          On exit, E has been destroyed. */
  90
  91 /*  INFO    (output) INTEGER */
  92 /*          = 0:  successful exit */
  93 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
  94 /*          > 0:  the algorithm failed to find all of the eigenvalues in */
  95 /*                a total of 30*N iterations; if INFO = i, then i */
  96 /*                elements of E have not converged to zero. */
  97
  98 /*  ===================================================================== */
  99
 100 /*     .. Parameters .. */
 101 /*     .. */
 102 /*     .. Local Scalars .. */
 103 /*     .. */
 104 /*     .. External Functions .. */
 105 /*     .. */
 106 /*     .. External Subroutines .. */
 107 /*     .. */
 108 /*     .. Intrinsic Functions .. */
 109 /*     .. */
 110 /*     .. Executable Statements .. */
 111
 112 /*     Test the input parameters. */
 113
 114     /* Parameter adjustments */
 115     --e;
 116     --d__;
 117
 118     /* Function Body */
 119     *info = 0;
 120
 121 /*     Quick return if possible */
 122
 123     if (*n < 0) {
 124         *info = -1;
 125         i__1 = -(*info);
 126         xerbla_("DSTERF", &i__1);
 127         return 0;
 128     }
 129     if (*n <= 1) {
 130         return 0;
 131     }
 132
 133 /*     Determine the unit roundoff for this environment. */
 134
 135     eps = dlamch_("E");
 136 /* Computing 2nd power */
 137     d__1 = eps;
 138     eps2 = d__1 * d__1;
 139     safmin = dlamch_("S");
 140     safmax = 1. / safmin;
 141     ssfmax = sqrt(safmax) / 3.;
 142     ssfmin = sqrt(safmin) / eps2;
 143
 144 /*     Compute the eigenvalues of the tridiagonal matrix. */
 145
 146     nmaxit = *n * 30;
 147     sigma = 0.;
 148     jtot = 0;
 149
 150 /*     Determine where the matrix splits and choose QL or QR iteration */
 151 /*     for each block, according to whether top or bottom diagonal */
 152 /*     element is smaller. */
 153
 154     l1 = 1;
 155
 156 L10:
 157     if (l1 > *n) {
 158         goto L170;
 159     }
 160     if (l1 > 1) {
 161         e[l1 - 1] = 0.;
 162     }
 163     i__1 = *n - 1;
 164     for (m = l1; m <= i__1; ++m) {
 165         if ((d__3 = e[m], abs(d__3)) <= sqrt((d__1 = d__[m], abs(d__1))) *
 166                 sqrt((d__2 = d__[m + 1], abs(d__2))) * eps) {
 167             e[m] = 0.;
 168             goto L30;
 169         }
 170 /* L20: */
 171     }
 172     m = *n;
 173
 174 L30:
 175     l = l1;
 176     lsv = l;
 177     lend = m;
 178     lendsv = lend;
 179     l1 = m + 1;
 180     if (lend == l) {
 181         goto L10;
 182     }
 183
 184 /*     Scale submatrix in rows and columns L to LEND */
 185
 186     i__1 = lend - l + 1;
 187     anorm = dlanst_("I", &i__1, &d__[l], &e[l]);
 188     iscale = 0;
 189     if (anorm > ssfmax) {
 190         iscale = 1;
 191         i__1 = lend - l + 1;
 192         dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
 193                 info);
 194         i__1 = lend - l;
 195         dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
 196                 info);
 197     } else if (anorm < ssfmin) {
 198         iscale = 2;
 199         i__1 = lend - l + 1;
 200         dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
 201                 info);
 202         i__1 = lend - l;
 203         dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
 204                 info);
 205     }
 206
 207     i__1 = lend - 1;
 208     for (i__ = l; i__ <= i__1; ++i__) {
 209 /* Computing 2nd power */
 210         d__1 = e[i__];
 211         e[i__] = d__1 * d__1;
 212 /* L40: */
 213     }
 214
 215 /*     Choose between QL and QR iteration */
 216
 217     if ((d__1 = d__[lend], abs(d__1)) < (d__2 = d__[l], abs(d__2))) {
 218         lend = lsv;
 219         l = lendsv;
 220     }
 221
 222     if (lend >= l) {
 223
 224 /*        QL Iteration */
 225
 226 /*        Look for small subdiagonal element. */
 227
 228 L50:
 229         if (l != lend) {
 230             i__1 = lend - 1;
 231             for (m = l; m <= i__1; ++m) {
 232                 if ((d__2 = e[m], abs(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
 233                         + 1], abs(d__1))) {
 234                     goto L70;
 235                 }
 236 /* L60: */
 237             }
 238         }
 239         m = lend;
 240
 241 L70:
 242         if (m < lend) {
 243             e[m] = 0.;
 244         }
 245         p = d__[l];
 246         if (m == l) {
 247             goto L90;
 248         }
 249
 250 /*        If remaining matrix is 2 by 2, use DLAE2 to compute its */
 251 /*        eigenvalues. */
 252
 253         if (m == l + 1) {
 254             rte = sqrt(e[l]);
 255             dlae2_(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
 256             d__[l] = rt1;
 257             d__[l + 1] = rt2;
 258             e[l] = 0.;
 259             l += 2;
 260             if (l <= lend) {
 261                 goto L50;
 262             }
 263             goto L150;
 264         }
 265
 266         if (jtot == nmaxit) {
 267             goto L150;
 268         }
 269         ++jtot;
 270
 271 /*        Form shift. */
 272
 273         rte = sqrt(e[l]);
 274         sigma = (d__[l + 1] - p) / (rte * 2.);
 275         r__ = dlapy2_(&sigma, &c_b32);
 276         sigma = p - rte / (sigma + d_sign(&r__, &sigma));
 277
 278         c__ = 1.;
 279         s = 0.;
 280         gamma = d__[m] - sigma;
 281         p = gamma * gamma;
 282
 283 /*        Inner loop */
 284
 285         i__1 = l;
 286         for (i__ = m - 1; i__ >= i__1; --i__) {
 287             bb = e[i__];
 288             r__ = p + bb;
 289             if (i__ != m - 1) {
 290                 e[i__ + 1] = s * r__;
 291             }
 292             oldc = c__;
 293             c__ = p / r__;
 294             s = bb / r__;
 295             oldgam = gamma;
 296             alpha = d__[i__];
 297             gamma = c__ * (alpha - sigma) - s * oldgam;
 298             d__[i__ + 1] = oldgam + (alpha - gamma);
 299             if (c__ != 0.) {
 300                 p = gamma * gamma / c__;
 301             } else {
 302                 p = oldc * bb;
 303             }
 304 /* L80: */
 305         }
 306
 307         e[l] = s * p;
 308         d__[l] = sigma + gamma;
 309         goto L50;
 310
 311 /*        Eigenvalue found. */
 312
 313 L90:
 314         d__[l] = p;
 315
 316         ++l;
 317         if (l <= lend) {
 318             goto L50;
 319         }
 320         goto L150;
 321
 322     } else {
 323
 324 /*        QR Iteration */
 325
 326 /*        Look for small superdiagonal element. */
 327
 328 L100:
 329         i__1 = lend + 1;
 330         for (m = l; m >= i__1; --m) {
 331             if ((d__2 = e[m - 1], abs(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
 332                     - 1], abs(d__1))) {
 333                 goto L120;
 334             }
 335 /* L110: */
 336         }
 337         m = lend;
 338
 339 L120:
 340         if (m > lend) {
 341             e[m - 1] = 0.;
 342         }
 343         p = d__[l];
 344         if (m == l) {
 345             goto L140;
 346         }
 347
 348 /*        If remaining matrix is 2 by 2, use DLAE2 to compute its */
 349 /*        eigenvalues. */
 350
 351         if (m == l - 1) {
 352             rte = sqrt(e[l - 1]);
 353             dlae2_(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
 354             d__[l] = rt1;
 355             d__[l - 1] = rt2;
 356             e[l - 1] = 0.;
 357             l += -2;
 358             if (l >= lend) {
 359                 goto L100;
 360             }
 361             goto L150;
 362         }
 363
 364         if (jtot == nmaxit) {
 365             goto L150;
 366         }
 367         ++jtot;
 368
 369 /*        Form shift. */
 370
 371         rte = sqrt(e[l - 1]);
 372         sigma = (d__[l - 1] - p) / (rte * 2.);
 373         r__ = dlapy2_(&sigma, &c_b32);
 374         sigma = p - rte / (sigma + d_sign(&r__, &sigma));
 375
 376         c__ = 1.;
 377         s = 0.;
 378         gamma = d__[m] - sigma;
 379         p = gamma * gamma;
 380
 381 /*        Inner loop */
 382
 383         i__1 = l - 1;
 384         for (i__ = m; i__ <= i__1; ++i__) {
 385             bb = e[i__];
 386             r__ = p + bb;
 387             if (i__ != m) {
 388                 e[i__ - 1] = s * r__;
 389             }
 390             oldc = c__;
 391             c__ = p / r__;
 392             s = bb / r__;
 393             oldgam = gamma;
 394             alpha = d__[i__ + 1];
 395             gamma = c__ * (alpha - sigma) - s * oldgam;
 396             d__[i__] = oldgam + (alpha - gamma);
 397             if (c__ != 0.) {
 398                 p = gamma * gamma / c__;
 399             } else {
 400                 p = oldc * bb;
 401             }
 402 /* L130: */
 403         }
 404
 405         e[l - 1] = s * p;
 406         d__[l] = sigma + gamma;
 407         goto L100;
 408
 409 /*        Eigenvalue found. */
 410
 411 L140:
 412         d__[l] = p;
 413
 414         --l;
 415         if (l >= lend) {
 416             goto L100;
 417         }
 418         goto L150;
 419
 420     }
 421
 422 /*     Undo scaling if necessary */
 423
 424 L150:
 425     if (iscale == 1) {
 426         i__1 = lendsv - lsv + 1;
 427         dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
 428                 n, info);
 429     }
 430     if (iscale == 2) {
 431         i__1 = lendsv - lsv + 1;
 432         dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
 433                 n, info);
 434     }
 435
 436 /*     Check for no convergence to an eigenvalue after a total */
 437 /*     of N*MAXIT iterations. */
 438
 439     if (jtot < nmaxit) {
 440         goto L10;
 441     }
 442     i__1 = *n - 1;
 443     for (i__ = 1; i__ <= i__1; ++i__) {
 444         if (e[i__] != 0.) {
 445             ++(*info);
 446         }
 447 /* L160: */
 448     }
 449     goto L180;
 450
 451 /*     Sort eigenvalues in increasing order. */
 452
 453 L170:
 454     dlasrt_("I", n, &d__[1], info);
 455
 456 L180:
 457     return 0;
 458
 459 /*     End of DSTERF */
 460
 461 } /* dsterf_ */