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

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