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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c__2 = 2;
   7 static integer c__10 = 10;
   8 static integer c__3 = 3;
   9 static integer c__4 = 4;
  10 static integer c__11 = 11;
  11
  12 /* Subroutine */ int dlasq2_(integer *n, doublereal *z__, integer *info)
  13 {
  14     /* System generated locals */
  15     integer i__1, i__2, i__3;
  16     doublereal d__1, d__2;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal);
  20
  21     /* Local variables */
  22     doublereal d__, e;
  23     integer k;
  24     doublereal s, t;
  25     integer i0, i4, n0;
  26     doublereal dn;
  27     integer pp;
  28     doublereal dn1, dn2, eps, tau, tol;
  29     integer ipn4;
  30     doublereal tol2;
  31     logical ieee;
  32     integer nbig;
  33     doublereal dmin__, emin, emax;
  34     integer ndiv, iter;
  35     doublereal qmin, temp, qmax, zmax;
  36     integer splt;
  37     doublereal dmin1, dmin2;
  38     integer nfail;
  39     doublereal desig, trace, sigma;
  40     integer iinfo, ttype;
  41     extern /* Subroutine */ int dlazq3_(integer *, integer *, doublereal *,
  42             integer *, doublereal *, doublereal *, doublereal *, doublereal *,
  43              integer *, integer *, integer *, logical *, integer *,
  44             doublereal *, doublereal *, doublereal *, doublereal *,
  45             doublereal *, doublereal *);
  46     extern doublereal dlamch_(char *);
  47     integer iwhila, iwhilb;
  48     doublereal oldemn, safmin;
  49     extern /* Subroutine */ int xerbla_(char *, integer *);
  50     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  51             integer *, integer *);
  52     extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *,
  53             integer *);
  54
  55
  56 /*  -- LAPACK routine (version 3.1) -- */
  57 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  58 /*     November 2006 */
  59
  60 /*     Modified to call DLAZQ3 in place of DLASQ3, 13 Feb 03, SJH. */
  61
  62 /*     .. Scalar Arguments .. */
  63 /*     .. */
  64 /*     .. Array Arguments .. */
  65 /*     .. */
  66
  67 /*  Purpose */
  68 /*  ======= */
  69
  70 /*  DLASQ2 computes all the eigenvalues of the symmetric positive */
  71 /*  definite tridiagonal matrix associated with the qd array Z to high */
  72 /*  relative accuracy are computed to high relative accuracy, in the */
  73 /*  absence of denormalization, underflow and overflow. */
  74
  75 /*  To see the relation of Z to the tridiagonal matrix, let L be a */
  76 /*  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and */
  77 /*  let U be an upper bidiagonal matrix with 1's above and diagonal */
  78 /*  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the */
  79 /*  symmetric tridiagonal to which it is similar. */
  80
  81 /*  Note : DLASQ2 defines a logical variable, IEEE, which is true */
  82 /*  on machines which follow ieee-754 floating-point standard in their */
  83 /*  handling of infinities and NaNs, and false otherwise. This variable */
  84 /*  is passed to DLAZQ3. */
  85
  86 /*  Arguments */
  87 /*  ========= */
  88
  89 /*  N     (input) INTEGER */
  90 /*        The number of rows and columns in the matrix. N >= 0. */
  91
  92 /*  Z     (workspace) DOUBLE PRECISION array, dimension ( 4*N ) */
  93 /*        On entry Z holds the qd array. On exit, entries 1 to N hold */
  94 /*        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the */
  95 /*        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If */
  96 /*        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) */
  97 /*        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of */
  98 /*        shifts that failed. */
  99
 100 /*  INFO  (output) INTEGER */
 101 /*        = 0: successful exit */
 102 /*        < 0: if the i-th argument is a scalar and had an illegal */
 103 /*             value, then INFO = -i, if the i-th argument is an */
 104 /*             array and the j-entry had an illegal value, then */
 105 /*             INFO = -(i*100+j) */
 106 /*        > 0: the algorithm failed */
 107 /*              = 1, a split was marked by a positive value in E */
 108 /*              = 2, current block of Z not diagonalized after 30*N */
 109 /*                   iterations (in inner while loop) */
 110 /*              = 3, termination criterion of outer while loop not met */
 111 /*                   (program created more than N unreduced blocks) */
 112
 113 /*  Further Details */
 114 /*  =============== */
 115 /*  Local Variables: I0:N0 defines a current unreduced segment of Z. */
 116 /*  The shifts are accumulated in SIGMA. Iteration count is in ITER. */
 117 /*  Ping-pong is controlled by PP (alternates between 0 and 1). */
 118
 119 /*  ===================================================================== */
 120
 121 /*     .. Parameters .. */
 122 /*     .. */
 123 /*     .. Local Scalars .. */
 124 /*     .. */
 125 /*     .. External Subroutines .. */
 126 /*     .. */
 127 /*     .. External Functions .. */
 128 /*     .. */
 129 /*     .. Intrinsic Functions .. */
 130 /*     .. */
 131 /*     .. Executable Statements .. */
 132
 133 /*     Test the input arguments. */
 134 /*     (in case DLASQ2 is not called by DLASQ1) */
 135
 136     /* Parameter adjustments */
 137     --z__;
 138
 139     /* Function Body */
 140     *info = 0;
 141     eps = dlamch_("Precision");
 142     safmin = dlamch_("Safe minimum");
 143     tol = eps * 100.;
 144 /* Computing 2nd power */
 145     d__1 = tol;
 146     tol2 = d__1 * d__1;
 147
 148     if (*n < 0) {
 149         *info = -1;
 150         xerbla_("DLASQ2", &c__1);
 151         return 0;
 152     } else if (*n == 0) {
 153         return 0;
 154     } else if (*n == 1) {
 155
 156 /*        1-by-1 case. */
 157
 158         if (z__[1] < 0.) {
 159             *info = -201;
 160             xerbla_("DLASQ2", &c__2);
 161         }
 162         return 0;
 163     } else if (*n == 2) {
 164
 165 /*        2-by-2 case. */
 166
 167         if (z__[2] < 0. || z__[3] < 0.) {
 168             *info = -2;
 169             xerbla_("DLASQ2", &c__2);
 170             return 0;
 171         } else if (z__[3] > z__[1]) {
 172             d__ = z__[3];
 173             z__[3] = z__[1];
 174             z__[1] = d__;
 175         }
 176         z__[5] = z__[1] + z__[2] + z__[3];
 177         if (z__[2] > z__[3] * tol2) {
 178             t = (z__[1] - z__[3] + z__[2]) * .5;
 179             s = z__[3] * (z__[2] / t);
 180             if (s <= t) {
 181                 s = z__[3] * (z__[2] / (t * (sqrt(s / t + 1.) + 1.)));
 182             } else {
 183                 s = z__[3] * (z__[2] / (t + sqrt(t) * sqrt(t + s)));
 184             }
 185             t = z__[1] + (s + z__[2]);
 186             z__[3] *= z__[1] / t;
 187             z__[1] = t;
 188         }
 189         z__[2] = z__[3];
 190         z__[6] = z__[2] + z__[1];
 191         return 0;
 192     }
 193
 194 /*     Check for negative data and compute sums of q's and e's. */
 195
 196     z__[*n * 2] = 0.;
 197     emin = z__[2];
 198     qmax = 0.;
 199     zmax = 0.;
 200     d__ = 0.;
 201     e = 0.;
 202
 203     i__1 = *n - 1 << 1;
 204     for (k = 1; k <= i__1; k += 2) {
 205         if (z__[k] < 0.) {
 206             *info = -(k + 200);
 207             xerbla_("DLASQ2", &c__2);
 208             return 0;
 209         } else if (z__[k + 1] < 0.) {
 210             *info = -(k + 201);
 211             xerbla_("DLASQ2", &c__2);
 212             return 0;
 213         }
 214         d__ += z__[k];
 215         e += z__[k + 1];
 216 /* Computing MAX */
 217         d__1 = qmax, d__2 = z__[k];
 218         qmax = max(d__1,d__2);
 219 /* Computing MIN */
 220         d__1 = emin, d__2 = z__[k + 1];
 221         emin = min(d__1,d__2);
 222 /* Computing MAX */
 223         d__1 = max(qmax,zmax), d__2 = z__[k + 1];
 224         zmax = max(d__1,d__2);
 225 /* L10: */
 226     }
 227     if (z__[(*n << 1) - 1] < 0.) {
 228         *info = -((*n << 1) + 199);
 229         xerbla_("DLASQ2", &c__2);
 230         return 0;
 231     }
 232     d__ += z__[(*n << 1) - 1];
 233 /* Computing MAX */
 234     d__1 = qmax, d__2 = z__[(*n << 1) - 1];
 235     qmax = max(d__1,d__2);
 236     zmax = max(qmax,zmax);
 237
 238 /*     Check for diagonality. */
 239
 240     if (e == 0.) {
 241         i__1 = *n;
 242         for (k = 2; k <= i__1; ++k) {
 243             z__[k] = z__[(k << 1) - 1];
 244 /* L20: */
 245         }
 246         dlasrt_("D", n, &z__[1], &iinfo);
 247         z__[(*n << 1) - 1] = d__;
 248         return 0;
 249     }
 250
 251     trace = d__ + e;
 252
 253 /*     Check for zero data. */
 254
 255     if (trace == 0.) {
 256         z__[(*n << 1) - 1] = 0.;
 257         return 0;
 258     }
 259
 260 /*     Check whether the machine is IEEE conformable. */
 261
 262     ieee = ilaenv_(&c__10, "DLASQ2", "N", &c__1, &c__2, &c__3, &c__4) == 1 && ilaenv_(&c__11, "DLASQ2", "N", &c__1, &c__2,
 263              &c__3, &c__4) == 1;
 264
 265 /*     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). */
 266
 267     for (k = *n << 1; k >= 2; k += -2) {
 268         z__[k * 2] = 0.;
 269         z__[(k << 1) - 1] = z__[k];
 270         z__[(k << 1) - 2] = 0.;
 271         z__[(k << 1) - 3] = z__[k - 1];
 272 /* L30: */
 273     }
 274
 275     i0 = 1;
 276     n0 = *n;
 277
 278 /*     Reverse the qd-array, if warranted. */
 279
 280     if (z__[(i0 << 2) - 3] * 1.5 < z__[(n0 << 2) - 3]) {
 281         ipn4 = i0 + n0 << 2;
 282         i__1 = i0 + n0 - 1 << 1;
 283         for (i4 = i0 << 2; i4 <= i__1; i4 += 4) {
 284             temp = z__[i4 - 3];
 285             z__[i4 - 3] = z__[ipn4 - i4 - 3];
 286             z__[ipn4 - i4 - 3] = temp;
 287             temp = z__[i4 - 1];
 288             z__[i4 - 1] = z__[ipn4 - i4 - 5];
 289             z__[ipn4 - i4 - 5] = temp;
 290 /* L40: */
 291         }
 292     }
 293
 294 /*     Initial split checking via dqd and Li's test. */
 295
 296     pp = 0;
 297
 298     for (k = 1; k <= 2; ++k) {
 299
 300         d__ = z__[(n0 << 2) + pp - 3];
 301         i__1 = (i0 << 2) + pp;
 302         for (i4 = (n0 - 1 << 2) + pp; i4 >= i__1; i4 += -4) {
 303             if (z__[i4 - 1] <= tol2 * d__) {
 304                 z__[i4 - 1] = -0.;
 305                 d__ = z__[i4 - 3];
 306             } else {
 307                 d__ = z__[i4 - 3] * (d__ / (d__ + z__[i4 - 1]));
 308             }
 309 /* L50: */
 310         }
 311
 312 /*        dqd maps Z to ZZ plus Li's test. */
 313
 314         emin = z__[(i0 << 2) + pp + 1];
 315         d__ = z__[(i0 << 2) + pp - 3];
 316         i__1 = (n0 - 1 << 2) + pp;
 317         for (i4 = (i0 << 2) + pp; i4 <= i__1; i4 += 4) {
 318             z__[i4 - (pp << 1) - 2] = d__ + z__[i4 - 1];
 319             if (z__[i4 - 1] <= tol2 * d__) {
 320                 z__[i4 - 1] = -0.;
 321                 z__[i4 - (pp << 1) - 2] = d__;
 322                 z__[i4 - (pp << 1)] = 0.;
 323                 d__ = z__[i4 + 1];
 324             } else if (safmin * z__[i4 + 1] < z__[i4 - (pp << 1) - 2] &&
 325                     safmin * z__[i4 - (pp << 1) - 2] < z__[i4 + 1]) {
 326                 temp = z__[i4 + 1] / z__[i4 - (pp << 1) - 2];
 327                 z__[i4 - (pp << 1)] = z__[i4 - 1] * temp;
 328                 d__ *= temp;
 329             } else {
 330                 z__[i4 - (pp << 1)] = z__[i4 + 1] * (z__[i4 - 1] / z__[i4 - (
 331                         pp << 1) - 2]);
 332                 d__ = z__[i4 + 1] * (d__ / z__[i4 - (pp << 1) - 2]);
 333             }
 334 /* Computing MIN */
 335             d__1 = emin, d__2 = z__[i4 - (pp << 1)];
 336             emin = min(d__1,d__2);
 337 /* L60: */
 338         }
 339         z__[(n0 << 2) - pp - 2] = d__;
 340
 341 /*        Now find qmax. */
 342
 343         qmax = z__[(i0 << 2) - pp - 2];
 344         i__1 = (n0 << 2) - pp - 2;
 345         for (i4 = (i0 << 2) - pp + 2; i4 <= i__1; i4 += 4) {
 346 /* Computing MAX */
 347             d__1 = qmax, d__2 = z__[i4];
 348             qmax = max(d__1,d__2);
 349 /* L70: */
 350         }
 351
 352 /*        Prepare for the next iteration on K. */
 353
 354         pp = 1 - pp;
 355 /* L80: */
 356     }
 357
 358 /*     Initialise variables to pass to DLAZQ3 */
 359
 360     ttype = 0;
 361     dmin1 = 0.;
 362     dmin2 = 0.;
 363     dn = 0.;
 364     dn1 = 0.;
 365     dn2 = 0.;
 366     tau = 0.;
 367
 368     iter = 2;
 369     nfail = 0;
 370     ndiv = n0 - i0 << 1;
 371
 372     i__1 = *n + 1;
 373     for (iwhila = 1; iwhila <= i__1; ++iwhila) {
 374         if (n0 < 1) {
 375             goto L150;
 376         }
 377
 378 /*        While array unfinished do */
 379
 380 /*        E(N0) holds the value of SIGMA when submatrix in I0:N0 */
 381 /*        splits from the rest of the array, but is negated. */
 382
 383         desig = 0.;
 384         if (n0 == *n) {
 385             sigma = 0.;
 386         } else {
 387             sigma = -z__[(n0 << 2) - 1];
 388         }
 389         if (sigma < 0.) {
 390             *info = 1;
 391             return 0;
 392         }
 393
 394 /*        Find last unreduced submatrix's top index I0, find QMAX and */
 395 /*        EMIN. Find Gershgorin-type bound if Q's much greater than E's. */
 396
 397         emax = 0.;
 398         if (n0 > i0) {
 399             emin = (d__1 = z__[(n0 << 2) - 5], abs(d__1));
 400         } else {
 401             emin = 0.;
 402         }
 403         qmin = z__[(n0 << 2) - 3];
 404         qmax = qmin;
 405         for (i4 = n0 << 2; i4 >= 8; i4 += -4) {
 406             if (z__[i4 - 5] <= 0.) {
 407                 goto L100;
 408             }
 409             if (qmin >= emax * 4.) {
 410 /* Computing MIN */
 411                 d__1 = qmin, d__2 = z__[i4 - 3];
 412                 qmin = min(d__1,d__2);
 413 /* Computing MAX */
 414                 d__1 = emax, d__2 = z__[i4 - 5];
 415                 emax = max(d__1,d__2);
 416             }
 417 /* Computing MAX */
 418             d__1 = qmax, d__2 = z__[i4 - 7] + z__[i4 - 5];
 419             qmax = max(d__1,d__2);
 420 /* Computing MIN */
 421             d__1 = emin, d__2 = z__[i4 - 5];
 422             emin = min(d__1,d__2);
 423 /* L90: */
 424         }
 425         i4 = 4;
 426
 427 L100:
 428         i0 = i4 / 4;
 429
 430 /*        Store EMIN for passing to DLAZQ3. */
 431
 432         z__[(n0 << 2) - 1] = emin;
 433
 434 /*        Put -(initial shift) into DMIN. */
 435
 436 /* Computing MAX */
 437         d__1 = 0., d__2 = qmin - sqrt(qmin) * 2. * sqrt(emax);
 438         dmin__ = -max(d__1,d__2);
 439
 440 /*        Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong. */
 441
 442         pp = 0;
 443
 444         nbig = (n0 - i0 + 1) * 30;
 445         i__2 = nbig;
 446         for (iwhilb = 1; iwhilb <= i__2; ++iwhilb) {
 447             if (i0 > n0) {
 448                 goto L130;
 449             }
 450
 451 /*           While submatrix unfinished take a good dqds step. */
 452
 453             dlazq3_(&i0, &n0, &z__[1], &pp, &dmin__, &sigma, &desig, &qmax, &
 454                     nfail, &iter, &ndiv, &ieee, &ttype, &dmin1, &dmin2, &dn, &
 455                     dn1, &dn2, &tau);
 456
 457             pp = 1 - pp;
 458
 459 /*           When EMIN is very small check for splits. */
 460
 461             if (pp == 0 && n0 - i0 >= 3) {
 462                 if (z__[n0 * 4] <= tol2 * qmax || z__[(n0 << 2) - 1] <= tol2 *
 463                          sigma) {
 464                     splt = i0 - 1;
 465                     qmax = z__[(i0 << 2) - 3];
 466                     emin = z__[(i0 << 2) - 1];
 467                     oldemn = z__[i0 * 4];
 468                     i__3 = n0 - 3 << 2;
 469                     for (i4 = i0 << 2; i4 <= i__3; i4 += 4) {
 470                         if (z__[i4] <= tol2 * z__[i4 - 3] || z__[i4 - 1] <=
 471                                 tol2 * sigma) {
 472                             z__[i4 - 1] = -sigma;
 473                             splt = i4 / 4;
 474                             qmax = 0.;
 475                             emin = z__[i4 + 3];
 476                             oldemn = z__[i4 + 4];
 477                         } else {
 478 /* Computing MAX */
 479                             d__1 = qmax, d__2 = z__[i4 + 1];
 480                             qmax = max(d__1,d__2);
 481 /* Computing MIN */
 482                             d__1 = emin, d__2 = z__[i4 - 1];
 483                             emin = min(d__1,d__2);
 484 /* Computing MIN */
 485                             d__1 = oldemn, d__2 = z__[i4];
 486                             oldemn = min(d__1,d__2);
 487                         }
 488 /* L110: */
 489                     }
 490                     z__[(n0 << 2) - 1] = emin;
 491                     z__[n0 * 4] = oldemn;
 492                     i0 = splt + 1;
 493                 }
 494             }
 495
 496 /* L120: */
 497         }
 498
 499         *info = 2;
 500         return 0;
 501
 502 /*        end IWHILB */
 503
 504 L130:
 505
 506 /* L140: */
 507         ;
 508     }
 509
 510     *info = 3;
 511     return 0;
 512
 513 /*     end IWHILA */
 514
 515 L150:
 516
 517 /*     Move q's to the front. */
 518
 519     i__1 = *n;
 520     for (k = 2; k <= i__1; ++k) {
 521         z__[k] = z__[(k << 2) - 3];
 522 /* L160: */
 523     }
 524
 525 /*     Sort and compute sum of eigenvalues. */
 526
 527     dlasrt_("D", n, &z__[1], &iinfo);
 528
 529     e = 0.;
 530     for (k = *n; k >= 1; --k) {
 531         e += z__[k];
 532 /* L170: */
 533     }
 534
 535 /*     Store trace, sum(eigenvalues) and information on performance. */
 536
 537     z__[(*n << 1) + 1] = trace;
 538     z__[(*n << 1) + 2] = e;
 539     z__[(*n << 1) + 3] = (doublereal) iter;
 540 /* Computing 2nd power */
 541     i__1 = *n;
 542     z__[(*n << 1) + 4] = (doublereal) ndiv / (doublereal) (i__1 * i__1);
 543     z__[(*n << 1) + 5] = nfail * 100. / (doublereal) iter;
 544     return 0;
 545
 546 /*     End of DLASQ2 */
 547
 548 } /* dlasq2_ */