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