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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static real c_b9 = 0.f;
   6 static real c_b10 = 1.f;
   7 static integer c__0 = 0;
   8 static integer c__1 = 1;
   9 static integer c__2 = 2;
  10
  11 /* Subroutine */ int ssteqr_(char *compz, integer *n, real *d__, real *e,
  12         real *z__, integer *ldz, real *work, integer *info)
  13 {
  14     /* System generated locals */
  15     integer z_dim1, z_offset, i__1, i__2;
  16     real r__1, r__2;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal), r_sign(real *, real *);
  20
  21     /* Local variables */
  22     real b, c__, f, g;
  23     integer i__, j, k, l, m;
  24     real p, r__, s;
  25     integer l1, ii, mm, lm1, mm1, nm1;
  26     real rt1, rt2, eps;
  27     integer lsv;
  28     real tst, eps2;
  29     integer lend, jtot;
  30     extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)
  31             ;
  32     extern logical lsame_(char *, char *);
  33     real anorm;
  34     extern /* Subroutine */ int slasr_(char *, char *, char *, integer *,
  35             integer *, real *, real *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer *);
  36     integer lendm1, lendp1;
  37     extern /* Subroutine */ int slaev2_(real *, real *, real *, real *, real *
  38 , real *, real *);
  39     extern doublereal slapy2_(real *, real *);
  40     integer iscale;
  41     extern doublereal slamch_(char *);
  42     real safmin;
  43     extern /* Subroutine */ int xerbla_(char *, integer *);
  44     real safmax;
  45     extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
  46             real *, integer *, integer *, real *, integer *, integer *);
  47     integer lendsv;
  48     extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real *
  49 ), slaset_(char *, integer *, integer *, real *, real *, real *,
  50             integer *);
  51     real ssfmin;
  52     integer nmaxit, icompz;
  53     real ssfmax;
  54     extern doublereal slanst_(char *, integer *, real *, real *);
  55     extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);
  56
  57
  58 /*  -- LAPACK routine (version 3.1) -- */
  59 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  60 /*     November 2006 */
  61
  62 /*     .. Scalar Arguments .. */
  63 /*     .. */
  64 /*     .. Array Arguments .. */
  65 /*     .. */
  66
  67 /*  Purpose */
  68 /*  ======= */
  69
  70 /*  SSTEQR computes all eigenvalues and, optionally, eigenvectors of a */
  71 /*  symmetric tridiagonal matrix using the implicit QL or QR method. */
  72 /*  The eigenvectors of a full or band symmetric matrix can also be found */
  73 /*  if SSYTRD or SSPTRD or SSBTRD has been used to reduce this matrix to */
  74 /*  tridiagonal form. */
  75
  76 /*  Arguments */
  77 /*  ========= */
  78
  79 /*  COMPZ   (input) CHARACTER*1 */
  80 /*          = 'N':  Compute eigenvalues only. */
  81 /*          = 'V':  Compute eigenvalues and eigenvectors of the original */
  82 /*                  symmetric matrix.  On entry, Z must contain the */
  83 /*                  orthogonal matrix used to reduce the original matrix */
  84 /*                  to tridiagonal form. */
  85 /*          = 'I':  Compute eigenvalues and eigenvectors of the */
  86 /*                  tridiagonal matrix.  Z is initialized to the identity */
  87 /*                  matrix. */
  88
  89 /*  N       (input) INTEGER */
  90 /*          The order of the matrix.  N >= 0. */
  91
  92 /*  D       (input/output) REAL array, dimension (N) */
  93 /*          On entry, the diagonal elements of the tridiagonal matrix. */
  94 /*          On exit, if INFO = 0, the eigenvalues in ascending order. */
  95
  96 /*  E       (input/output) REAL array, dimension (N-1) */
  97 /*          On entry, the (n-1) subdiagonal elements of the tridiagonal */
  98 /*          matrix. */
  99 /*          On exit, E has been destroyed. */
 100
 101 /*  Z       (input/output) REAL array, dimension (LDZ, N) */
 102 /*          On entry, if  COMPZ = 'V', then Z contains the orthogonal */
 103 /*          matrix used in the reduction to tridiagonal form. */
 104 /*          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the */
 105 /*          orthonormal eigenvectors of the original symmetric matrix, */
 106 /*          and if COMPZ = 'I', Z contains the orthonormal eigenvectors */
 107 /*          of the symmetric tridiagonal matrix. */
 108 /*          If COMPZ = 'N', then Z is not referenced. */
 109
 110 /*  LDZ     (input) INTEGER */
 111 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
 112 /*          eigenvectors are desired, then  LDZ >= max(1,N). */
 113
 114 /*  WORK    (workspace) REAL array, dimension (max(1,2*N-2)) */
 115 /*          If COMPZ = 'N', then WORK is not referenced. */
 116
 117 /*  INFO    (output) INTEGER */
 118 /*          = 0:  successful exit */
 119 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 120 /*          > 0:  the algorithm has failed to find all the eigenvalues in */
 121 /*                a total of 30*N iterations; if INFO = i, then i */
 122 /*                elements of E have not converged to zero; on exit, D */
 123 /*                and E contain the elements of a symmetric tridiagonal */
 124 /*                matrix which is orthogonally similar to the original */
 125 /*                matrix. */
 126
 127 /*  ===================================================================== */
 128
 129 /*     .. Parameters .. */
 130 /*     .. */
 131 /*     .. Local Scalars .. */
 132 /*     .. */
 133 /*     .. External Functions .. */
 134 /*     .. */
 135 /*     .. External Subroutines .. */
 136 /*     .. */
 137 /*     .. Intrinsic Functions .. */
 138 /*     .. */
 139 /*     .. Executable Statements .. */
 140
 141 /*     Test the input parameters. */
 142
 143     /* Parameter adjustments */
 144     --d__;
 145     --e;
 146     z_dim1 = *ldz;
 147     z_offset = 1 + z_dim1;
 148     z__ -= z_offset;
 149     --work;
 150
 151     /* Function Body */
 152     *info = 0;
 153
 154     if (lsame_(compz, "N")) {
 155         icompz = 0;
 156     } else if (lsame_(compz, "V")) {
 157         icompz = 1;
 158     } else if (lsame_(compz, "I")) {
 159         icompz = 2;
 160     } else {
 161         icompz = -1;
 162     }
 163     if (icompz < 0) {
 164         *info = -1;
 165     } else if (*n < 0) {
 166         *info = -2;
 167     } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) {
 168         *info = -6;
 169     }
 170     if (*info != 0) {
 171         i__1 = -(*info);
 172         xerbla_("SSTEQR", &i__1);
 173         return 0;
 174     }
 175
 176 /*     Quick return if possible */
 177
 178     if (*n == 0) {
 179         return 0;
 180     }
 181
 182     if (*n == 1) {
 183         if (icompz == 2) {
 184             z__[z_dim1 + 1] = 1.f;
 185         }
 186         return 0;
 187     }
 188
 189 /*     Determine the unit roundoff and over/underflow thresholds. */
 190
 191     eps = slamch_("E");
 192 /* Computing 2nd power */
 193     r__1 = eps;
 194     eps2 = r__1 * r__1;
 195     safmin = slamch_("S");
 196     safmax = 1.f / safmin;
 197     ssfmax = sqrt(safmax) / 3.f;
 198     ssfmin = sqrt(safmin) / eps2;
 199
 200 /*     Compute the eigenvalues and eigenvectors of the tridiagonal */
 201 /*     matrix. */
 202
 203     if (icompz == 2) {
 204         slaset_("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
 205     }
 206
 207     nmaxit = *n * 30;
 208     jtot = 0;
 209
 210 /*     Determine where the matrix splits and choose QL or QR iteration */
 211 /*     for each block, according to whether top or bottom diagonal */
 212 /*     element is smaller. */
 213
 214     l1 = 1;
 215     nm1 = *n - 1;
 216
 217 L10:
 218     if (l1 > *n) {
 219         goto L160;
 220     }
 221     if (l1 > 1) {
 222         e[l1 - 1] = 0.f;
 223     }
 224     if (l1 <= nm1) {
 225         i__1 = nm1;
 226         for (m = l1; m <= i__1; ++m) {
 227             tst = (r__1 = e[m], dabs(r__1));
 228             if (tst == 0.f) {
 229                 goto L30;
 230             }
 231             if (tst <= sqrt((r__1 = d__[m], dabs(r__1))) * sqrt((r__2 = d__[m
 232                     + 1], dabs(r__2))) * eps) {
 233                 e[m] = 0.f;
 234                 goto L30;
 235             }
 236 /* L20: */
 237         }
 238     }
 239     m = *n;
 240
 241 L30:
 242     l = l1;
 243     lsv = l;
 244     lend = m;
 245     lendsv = lend;
 246     l1 = m + 1;
 247     if (lend == l) {
 248         goto L10;
 249     }
 250
 251 /*     Scale submatrix in rows and columns L to LEND */
 252
 253     i__1 = lend - l + 1;
 254     anorm = slanst_("I", &i__1, &d__[l], &e[l]);
 255     iscale = 0;
 256     if (anorm == 0.f) {
 257         goto L10;
 258     }
 259     if (anorm > ssfmax) {
 260         iscale = 1;
 261         i__1 = lend - l + 1;
 262         slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
 263                 info);
 264         i__1 = lend - l;
 265         slascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
 266                 info);
 267     } else if (anorm < ssfmin) {
 268         iscale = 2;
 269         i__1 = lend - l + 1;
 270         slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
 271                 info);
 272         i__1 = lend - l;
 273         slascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
 274                 info);
 275     }
 276
 277 /*     Choose between QL and QR iteration */
 278
 279     if ((r__1 = d__[lend], dabs(r__1)) < (r__2 = d__[l], dabs(r__2))) {
 280         lend = lsv;
 281         l = lendsv;
 282     }
 283
 284     if (lend > l) {
 285
 286 /*        QL Iteration */
 287
 288 /*        Look for small subdiagonal element. */
 289
 290 L40:
 291         if (l != lend) {
 292             lendm1 = lend - 1;
 293             i__1 = lendm1;
 294             for (m = l; m <= i__1; ++m) {
 295 /* Computing 2nd power */
 296                 r__2 = (r__1 = e[m], dabs(r__1));
 297                 tst = r__2 * r__2;
 298                 if (tst <= eps2 * (r__1 = d__[m], dabs(r__1)) * (r__2 = d__[m
 299                         + 1], dabs(r__2)) + safmin) {
 300                     goto L60;
 301                 }
 302 /* L50: */
 303             }
 304         }
 305
 306         m = lend;
 307
 308 L60:
 309         if (m < lend) {
 310             e[m] = 0.f;
 311         }
 312         p = d__[l];
 313         if (m == l) {
 314             goto L80;
 315         }
 316
 317 /*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 */
 318 /*        to compute its eigensystem. */
 319
 320         if (m == l + 1) {
 321             if (icompz > 0) {
 322                 slaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
 323                 work[l] = c__;
 324                 work[*n - 1 + l] = s;
 325                 slasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
 326                         z__[l * z_dim1 + 1], ldz);
 327             } else {
 328                 slae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
 329             }
 330             d__[l] = rt1;
 331             d__[l + 1] = rt2;
 332             e[l] = 0.f;
 333             l += 2;
 334             if (l <= lend) {
 335                 goto L40;
 336             }
 337             goto L140;
 338         }
 339
 340         if (jtot == nmaxit) {
 341             goto L140;
 342         }
 343         ++jtot;
 344
 345 /*        Form shift. */
 346
 347         g = (d__[l + 1] - p) / (e[l] * 2.f);
 348         r__ = slapy2_(&g, &c_b10);
 349         g = d__[m] - p + e[l] / (g + r_sign(&r__, &g));
 350
 351         s = 1.f;
 352         c__ = 1.f;
 353         p = 0.f;
 354
 355 /*        Inner loop */
 356
 357         mm1 = m - 1;
 358         i__1 = l;
 359         for (i__ = mm1; i__ >= i__1; --i__) {
 360             f = s * e[i__];
 361             b = c__ * e[i__];
 362             slartg_(&g, &f, &c__, &s, &r__);
 363             if (i__ != m - 1) {
 364                 e[i__ + 1] = r__;
 365             }
 366             g = d__[i__ + 1] - p;
 367             r__ = (d__[i__] - g) * s + c__ * 2.f * b;
 368             p = s * r__;
 369             d__[i__ + 1] = g + p;
 370             g = c__ * r__ - b;
 371
 372 /*           If eigenvectors are desired, then save rotations. */
 373
 374             if (icompz > 0) {
 375                 work[i__] = c__;
 376                 work[*n - 1 + i__] = -s;
 377             }
 378
 379 /* L70: */
 380         }
 381
 382 /*        If eigenvectors are desired, then apply saved rotations. */
 383
 384         if (icompz > 0) {
 385             mm = m - l + 1;
 386             slasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l
 387                     * z_dim1 + 1], ldz);
 388         }
 389
 390         d__[l] -= p;
 391         e[l] = g;
 392         goto L40;
 393
 394 /*        Eigenvalue found. */
 395
 396 L80:
 397         d__[l] = p;
 398
 399         ++l;
 400         if (l <= lend) {
 401             goto L40;
 402         }
 403         goto L140;
 404
 405     } else {
 406
 407 /*        QR Iteration */
 408
 409 /*        Look for small superdiagonal element. */
 410
 411 L90:
 412         if (l != lend) {
 413             lendp1 = lend + 1;
 414             i__1 = lendp1;
 415             for (m = l; m >= i__1; --m) {
 416 /* Computing 2nd power */
 417                 r__2 = (r__1 = e[m - 1], dabs(r__1));
 418                 tst = r__2 * r__2;
 419                 if (tst <= eps2 * (r__1 = d__[m], dabs(r__1)) * (r__2 = d__[m
 420                         - 1], dabs(r__2)) + safmin) {
 421                     goto L110;
 422                 }
 423 /* L100: */
 424             }
 425         }
 426
 427         m = lend;
 428
 429 L110:
 430         if (m > lend) {
 431             e[m - 1] = 0.f;
 432         }
 433         p = d__[l];
 434         if (m == l) {
 435             goto L130;
 436         }
 437
 438 /*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 */
 439 /*        to compute its eigensystem. */
 440
 441         if (m == l - 1) {
 442             if (icompz > 0) {
 443                 slaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
 444                         ;
 445                 work[m] = c__;
 446                 work[*n - 1 + m] = s;
 447                 slasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
 448                         z__[(l - 1) * z_dim1 + 1], ldz);
 449             } else {
 450                 slae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
 451             }
 452             d__[l - 1] = rt1;
 453             d__[l] = rt2;
 454             e[l - 1] = 0.f;
 455             l += -2;
 456             if (l >= lend) {
 457                 goto L90;
 458             }
 459             goto L140;
 460         }
 461
 462         if (jtot == nmaxit) {
 463             goto L140;
 464         }
 465         ++jtot;
 466
 467 /*        Form shift. */
 468
 469         g = (d__[l - 1] - p) / (e[l - 1] * 2.f);
 470         r__ = slapy2_(&g, &c_b10);
 471         g = d__[m] - p + e[l - 1] / (g + r_sign(&r__, &g));
 472
 473         s = 1.f;
 474         c__ = 1.f;
 475         p = 0.f;
 476
 477 /*        Inner loop */
 478
 479         lm1 = l - 1;
 480         i__1 = lm1;
 481         for (i__ = m; i__ <= i__1; ++i__) {
 482             f = s * e[i__];
 483             b = c__ * e[i__];
 484             slartg_(&g, &f, &c__, &s, &r__);
 485             if (i__ != m) {
 486                 e[i__ - 1] = r__;
 487             }
 488             g = d__[i__] - p;
 489             r__ = (d__[i__ + 1] - g) * s + c__ * 2.f * b;
 490             p = s * r__;
 491             d__[i__] = g + p;
 492             g = c__ * r__ - b;
 493
 494 /*           If eigenvectors are desired, then save rotations. */
 495
 496             if (icompz > 0) {
 497                 work[i__] = c__;
 498                 work[*n - 1 + i__] = s;
 499             }
 500
 501 /* L120: */
 502         }
 503
 504 /*        If eigenvectors are desired, then apply saved rotations. */
 505
 506         if (icompz > 0) {
 507             mm = l - m + 1;
 508             slasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m
 509                     * z_dim1 + 1], ldz);
 510         }
 511
 512         d__[l] -= p;
 513         e[lm1] = g;
 514         goto L90;
 515
 516 /*        Eigenvalue found. */
 517
 518 L130:
 519         d__[l] = p;
 520
 521         --l;
 522         if (l >= lend) {
 523             goto L90;
 524         }
 525         goto L140;
 526
 527     }
 528
 529 /*     Undo scaling if necessary */
 530
 531 L140:
 532     if (iscale == 1) {
 533         i__1 = lendsv - lsv + 1;
 534         slascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
 535                 n, info);
 536         i__1 = lendsv - lsv;
 537         slascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
 538                 info);
 539     } else if (iscale == 2) {
 540         i__1 = lendsv - lsv + 1;
 541         slascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
 542                 n, info);
 543         i__1 = lendsv - lsv;
 544         slascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
 545                 info);
 546     }
 547
 548 /*     Check for no convergence to an eigenvalue after a total */
 549 /*     of N*MAXIT iterations. */
 550
 551     if (jtot < nmaxit) {
 552         goto L10;
 553     }
 554     i__1 = *n - 1;
 555     for (i__ = 1; i__ <= i__1; ++i__) {
 556         if (e[i__] != 0.f) {
 557             ++(*info);
 558         }
 559 /* L150: */
 560     }
 561     goto L190;
 562
 563 /*     Order eigenvalues and eigenvectors. */
 564
 565 L160:
 566     if (icompz == 0) {
 567
 568 /*        Use Quick Sort */
 569
 570         slasrt_("I", n, &d__[1], info);
 571
 572     } else {
 573
 574 /*        Use Selection Sort to minimize swaps of eigenvectors */
 575
 576         i__1 = *n;
 577         for (ii = 2; ii <= i__1; ++ii) {
 578             i__ = ii - 1;
 579             k = i__;
 580             p = d__[i__];
 581             i__2 = *n;
 582             for (j = ii; j <= i__2; ++j) {
 583                 if (d__[j] < p) {
 584                     k = j;
 585                     p = d__[j];
 586                 }
 587 /* L170: */
 588             }
 589             if (k != i__) {
 590                 d__[k] = d__[i__];
 591                 d__[i__] = p;
 592                 sswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1],
 593                          &c__1);
 594             }
 595 /* L180: */
 596         }
 597     }
 598
 599 L190:
 600     return 0;
 601
 602 /*     End of SSTEQR */
 603
 604 } /* ssteqr_ */