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