git.maemo.org Git - opencv/blob - 3rdparty/lapack/dsterf.c

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