vcs.maemo.org Git - opencv/blob - 3rdparty/lapack/dlar1v.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlar1v_(integer *n, integer *b1, integer *bn, doublereal
   4         *lambda, doublereal *d__, doublereal *l, doublereal *ld, doublereal *
   5         lld, doublereal *pivmin, doublereal *gaptol, doublereal *z__, logical
   6         *wantnc, integer *negcnt, doublereal *ztz, doublereal *mingma,
   7         integer *r__, integer *isuppz, doublereal *nrminv, doublereal *resid,
   8         doublereal *rqcorr, doublereal *work)
   9 {
  10     /* System generated locals */
  11     integer i__1;
  12     doublereal d__1, d__2, d__3;
  13
  14     /* Builtin functions */
  15     double sqrt(doublereal);
  16
  17     /* Local variables */
  18     integer i__;
  19     doublereal s;
  20     integer r1, r2;
  21     doublereal eps, tmp;
  22     integer neg1, neg2, indp, inds;
  23     doublereal dplus;
  24     extern doublereal dlamch_(char *);
  25     extern logical disnan_(doublereal *);
  26     integer indlpl, indumn;
  27     doublereal dminus;
  28     logical sawnan1, sawnan2;
  29
  30
  31 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  32 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  33 /*     November 2006 */
  34
  35 /*     .. Scalar Arguments .. */
  36 /*     .. */
  37 /*     .. Array Arguments .. */
  38 /*     .. */
  39
  40 /*  Purpose */
  41 /*  ======= */
  42
  43 /*  DLAR1V computes the (scaled) r-th column of the inverse of */
  44 /*  the sumbmatrix in rows B1 through BN of the tridiagonal matrix */
  45 /*  L D L^T - sigma I. When sigma is close to an eigenvalue, the */
  46 /*  computed vector is an accurate eigenvector. Usually, r corresponds */
  47 /*  to the index where the eigenvector is largest in magnitude. */
  48 /*  The following steps accomplish this computation : */
  49 /*  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T, */
  50 /*  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T, */
  51 /*  (c) Computation of the diagonal elements of the inverse of */
  52 /*      L D L^T - sigma I by combining the above transforms, and choosing */
  53 /*      r as the index where the diagonal of the inverse is (one of the) */
  54 /*      largest in magnitude. */
  55 /*  (d) Computation of the (scaled) r-th column of the inverse using the */
  56 /*      twisted factorization obtained by combining the top part of the */
  57 /*      the stationary and the bottom part of the progressive transform. */
  58
  59 /*  Arguments */
  60 /*  ========= */
  61
  62 /*  N        (input) INTEGER */
  63 /*           The order of the matrix L D L^T. */
  64
  65 /*  B1       (input) INTEGER */
  66 /*           First index of the submatrix of L D L^T. */
  67
  68 /*  BN       (input) INTEGER */
  69 /*           Last index of the submatrix of L D L^T. */
  70
  71 /*  LAMBDA    (input) DOUBLE PRECISION */
  72 /*           The shift. In order to compute an accurate eigenvector, */
  73 /*           LAMBDA should be a good approximation to an eigenvalue */
  74 /*           of L D L^T. */
  75
  76 /*  L        (input) DOUBLE PRECISION array, dimension (N-1) */
  77 /*           The (n-1) subdiagonal elements of the unit bidiagonal matrix */
  78 /*           L, in elements 1 to N-1. */
  79
  80 /*  D        (input) DOUBLE PRECISION array, dimension (N) */
  81 /*           The n diagonal elements of the diagonal matrix D. */
  82
  83 /*  LD       (input) DOUBLE PRECISION array, dimension (N-1) */
  84 /*           The n-1 elements L(i)*D(i). */
  85
  86 /*  LLD      (input) DOUBLE PRECISION array, dimension (N-1) */
  87 /*           The n-1 elements L(i)*L(i)*D(i). */
  88
  89 /*  PIVMIN   (input) DOUBLE PRECISION */
  90 /*           The minimum pivot in the Sturm sequence. */
  91
  92 /*  GAPTOL   (input) DOUBLE PRECISION */
  93 /*           Tolerance that indicates when eigenvector entries are negligible */
  94 /*           w.r.t. their contribution to the residual. */
  95
  96 /*  Z        (input/output) DOUBLE PRECISION array, dimension (N) */
  97 /*           On input, all entries of Z must be set to 0. */
  98 /*           On output, Z contains the (scaled) r-th column of the */
  99 /*           inverse. The scaling is such that Z(R) equals 1. */
 100
 101 /*  WANTNC   (input) LOGICAL */
 102 /*           Specifies whether NEGCNT has to be computed. */
 103
 104 /*  NEGCNT   (output) INTEGER */
 105 /*           If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin */
 106 /*           in the  matrix factorization L D L^T, and NEGCNT = -1 otherwise. */
 107
 108 /*  ZTZ      (output) DOUBLE PRECISION */
 109 /*           The square of the 2-norm of Z. */
 110
 111 /*  MINGMA   (output) DOUBLE PRECISION */
 112 /*           The reciprocal of the largest (in magnitude) diagonal */
 113 /*           element of the inverse of L D L^T - sigma I. */
 114
 115 /*  R        (input/output) INTEGER */
 116 /*           The twist index for the twisted factorization used to */
 117 /*           compute Z. */
 118 /*           On input, 0 <= R <= N. If R is input as 0, R is set to */
 119 /*           the index where (L D L^T - sigma I)^{-1} is largest */
 120 /*           in magnitude. If 1 <= R <= N, R is unchanged. */
 121 /*           On output, R contains the twist index used to compute Z. */
 122 /*           Ideally, R designates the position of the maximum entry in the */
 123 /*           eigenvector. */
 124
 125 /*  ISUPPZ   (output) INTEGER array, dimension (2) */
 126 /*           The support of the vector in Z, i.e., the vector Z is */
 127 /*           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ). */
 128
 129 /*  NRMINV   (output) DOUBLE PRECISION */
 130 /*           NRMINV = 1/SQRT( ZTZ ) */
 131
 132 /*  RESID    (output) DOUBLE PRECISION */
 133 /*           The residual of the FP vector. */
 134 /*           RESID = ABS( MINGMA )/SQRT( ZTZ ) */
 135
 136 /*  RQCORR   (output) DOUBLE PRECISION */
 137 /*           The Rayleigh Quotient correction to LAMBDA. */
 138 /*           RQCORR = MINGMA*TMP */
 139
 140 /*  WORK     (workspace) DOUBLE PRECISION array, dimension (4*N) */
 141
 142 /*  Further Details */
 143 /*  =============== */
 144
 145 /*  Based on contributions by */
 146 /*     Beresford Parlett, University of California, Berkeley, USA */
 147 /*     Jim Demmel, University of California, Berkeley, USA */
 148 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 149 /*     Osni Marques, LBNL/NERSC, USA */
 150 /*     Christof Voemel, University of California, Berkeley, USA */
 151
 152 /*  ===================================================================== */
 153
 154 /*     .. Parameters .. */
 155 /*     .. */
 156 /*     .. Local Scalars .. */
 157 /*     .. */
 158 /*     .. External Functions .. */
 159 /*     .. */
 160 /*     .. Intrinsic Functions .. */
 161 /*     .. */
 162 /*     .. Executable Statements .. */
 163
 164     /* Parameter adjustments */
 165     --work;
 166     --isuppz;
 167     --z__;
 168     --lld;
 169     --ld;
 170     --l;
 171     --d__;
 172
 173     /* Function Body */
 174     eps = dlamch_("Precision");
 175     if (*r__ == 0) {
 176         r1 = *b1;
 177         r2 = *bn;
 178     } else {
 179         r1 = *r__;
 180         r2 = *r__;
 181     }
 182 /*     Storage for LPLUS */
 183     indlpl = 0;
 184 /*     Storage for UMINUS */
 185     indumn = *n;
 186     inds = (*n << 1) + 1;
 187     indp = *n * 3 + 1;
 188     if (*b1 == 1) {
 189         work[inds] = 0.;
 190     } else {
 191         work[inds + *b1 - 1] = lld[*b1 - 1];
 192     }
 193
 194 /*     Compute the stationary transform (using the differential form) */
 195 /*     until the index R2. */
 196
 197     sawnan1 = FALSE_;
 198     neg1 = 0;
 199     s = work[inds + *b1 - 1] - *lambda;
 200     i__1 = r1 - 1;
 201     for (i__ = *b1; i__ <= i__1; ++i__) {
 202         dplus = d__[i__] + s;
 203         work[indlpl + i__] = ld[i__] / dplus;
 204         if (dplus < 0.) {
 205             ++neg1;
 206         }
 207         work[inds + i__] = s * work[indlpl + i__] * l[i__];
 208         s = work[inds + i__] - *lambda;
 209 /* L50: */
 210     }
 211     sawnan1 = disnan_(&s);
 212     if (sawnan1) {
 213         goto L60;
 214     }
 215     i__1 = r2 - 1;
 216     for (i__ = r1; i__ <= i__1; ++i__) {
 217         dplus = d__[i__] + s;
 218         work[indlpl + i__] = ld[i__] / dplus;
 219         work[inds + i__] = s * work[indlpl + i__] * l[i__];
 220         s = work[inds + i__] - *lambda;
 221 /* L51: */
 222     }
 223     sawnan1 = disnan_(&s);
 224
 225 L60:
 226     if (sawnan1) {
 227 /*        Runs a slower version of the above loop if a NaN is detected */
 228         neg1 = 0;
 229         s = work[inds + *b1 - 1] - *lambda;
 230         i__1 = r1 - 1;
 231         for (i__ = *b1; i__ <= i__1; ++i__) {
 232             dplus = d__[i__] + s;
 233             if (abs(dplus) < *pivmin) {
 234                 dplus = -(*pivmin);
 235             }
 236             work[indlpl + i__] = ld[i__] / dplus;
 237             if (dplus < 0.) {
 238                 ++neg1;
 239             }
 240             work[inds + i__] = s * work[indlpl + i__] * l[i__];
 241             if (work[indlpl + i__] == 0.) {
 242                 work[inds + i__] = lld[i__];
 243             }
 244             s = work[inds + i__] - *lambda;
 245 /* L70: */
 246         }
 247         i__1 = r2 - 1;
 248         for (i__ = r1; i__ <= i__1; ++i__) {
 249             dplus = d__[i__] + s;
 250             if (abs(dplus) < *pivmin) {
 251                 dplus = -(*pivmin);
 252             }
 253             work[indlpl + i__] = ld[i__] / dplus;
 254             work[inds + i__] = s * work[indlpl + i__] * l[i__];
 255             if (work[indlpl + i__] == 0.) {
 256                 work[inds + i__] = lld[i__];
 257             }
 258             s = work[inds + i__] - *lambda;
 259 /* L71: */
 260         }
 261     }
 262
 263 /*     Compute the progressive transform (using the differential form) */
 264 /*     until the index R1 */
 265
 266     sawnan2 = FALSE_;
 267     neg2 = 0;
 268     work[indp + *bn - 1] = d__[*bn] - *lambda;
 269     i__1 = r1;
 270     for (i__ = *bn - 1; i__ >= i__1; --i__) {
 271         dminus = lld[i__] + work[indp + i__];
 272         tmp = d__[i__] / dminus;
 273         if (dminus < 0.) {
 274             ++neg2;
 275         }
 276         work[indumn + i__] = l[i__] * tmp;
 277         work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
 278 /* L80: */
 279     }
 280     tmp = work[indp + r1 - 1];
 281     sawnan2 = disnan_(&tmp);
 282     if (sawnan2) {
 283 /*        Runs a slower version of the above loop if a NaN is detected */
 284         neg2 = 0;
 285         i__1 = r1;
 286         for (i__ = *bn - 1; i__ >= i__1; --i__) {
 287             dminus = lld[i__] + work[indp + i__];
 288             if (abs(dminus) < *pivmin) {
 289                 dminus = -(*pivmin);
 290             }
 291             tmp = d__[i__] / dminus;
 292             if (dminus < 0.) {
 293                 ++neg2;
 294             }
 295             work[indumn + i__] = l[i__] * tmp;
 296             work[indp + i__ - 1] = work[indp + i__] * tmp - *lambda;
 297             if (tmp == 0.) {
 298                 work[indp + i__ - 1] = d__[i__] - *lambda;
 299             }
 300 /* L100: */
 301         }
 302     }
 303
 304 /*     Find the index (from R1 to R2) of the largest (in magnitude) */
 305 /*     diagonal element of the inverse */
 306
 307     *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
 308     if (*mingma < 0.) {
 309         ++neg1;
 310     }
 311     if (*wantnc) {
 312         *negcnt = neg1 + neg2;
 313     } else {
 314         *negcnt = -1;
 315     }
 316     if (abs(*mingma) == 0.) {
 317         *mingma = eps * work[inds + r1 - 1];
 318     }
 319     *r__ = r1;
 320     i__1 = r2 - 1;
 321     for (i__ = r1; i__ <= i__1; ++i__) {
 322         tmp = work[inds + i__] + work[indp + i__];
 323         if (tmp == 0.) {
 324             tmp = eps * work[inds + i__];
 325         }
 326         if (abs(tmp) <= abs(*mingma)) {
 327             *mingma = tmp;
 328             *r__ = i__ + 1;
 329         }
 330 /* L110: */
 331     }
 332
 333 /*     Compute the FP vector: solve N^T v = e_r */
 334
 335     isuppz[1] = *b1;
 336     isuppz[2] = *bn;
 337     z__[*r__] = 1.;
 338     *ztz = 1.;
 339
 340 /*     Compute the FP vector upwards from R */
 341
 342     if (! sawnan1 && ! sawnan2) {
 343         i__1 = *b1;
 344         for (i__ = *r__ - 1; i__ >= i__1; --i__) {
 345             z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
 346             if (((d__1 = z__[i__], abs(d__1)) + (d__2 = z__[i__ + 1], abs(
 347                     d__2))) * (d__3 = ld[i__], abs(d__3)) < *gaptol) {
 348                 z__[i__] = 0.;
 349                 isuppz[1] = i__ + 1;
 350                 goto L220;
 351             }
 352             *ztz += z__[i__] * z__[i__];
 353 /* L210: */
 354         }
 355 L220:
 356         ;
 357     } else {
 358 /*        Run slower loop if NaN occurred. */
 359         i__1 = *b1;
 360         for (i__ = *r__ - 1; i__ >= i__1; --i__) {
 361             if (z__[i__ + 1] == 0.) {
 362                 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
 363             } else {
 364                 z__[i__] = -(work[indlpl + i__] * z__[i__ + 1]);
 365             }
 366             if (((d__1 = z__[i__], abs(d__1)) + (d__2 = z__[i__ + 1], abs(
 367                     d__2))) * (d__3 = ld[i__], abs(d__3)) < *gaptol) {
 368                 z__[i__] = 0.;
 369                 isuppz[1] = i__ + 1;
 370                 goto L240;
 371             }
 372             *ztz += z__[i__] * z__[i__];
 373 /* L230: */
 374         }
 375 L240:
 376         ;
 377     }
 378 /*     Compute the FP vector downwards from R in blocks of size BLKSIZ */
 379     if (! sawnan1 && ! sawnan2) {
 380         i__1 = *bn - 1;
 381         for (i__ = *r__; i__ <= i__1; ++i__) {
 382             z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
 383             if (((d__1 = z__[i__], abs(d__1)) + (d__2 = z__[i__ + 1], abs(
 384                     d__2))) * (d__3 = ld[i__], abs(d__3)) < *gaptol) {
 385                 z__[i__ + 1] = 0.;
 386                 isuppz[2] = i__;
 387                 goto L260;
 388             }
 389             *ztz += z__[i__ + 1] * z__[i__ + 1];
 390 /* L250: */
 391         }
 392 L260:
 393         ;
 394     } else {
 395 /*        Run slower loop if NaN occurred. */
 396         i__1 = *bn - 1;
 397         for (i__ = *r__; i__ <= i__1; ++i__) {
 398             if (z__[i__] == 0.) {
 399                 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
 400             } else {
 401                 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
 402             }
 403             if (((d__1 = z__[i__], abs(d__1)) + (d__2 = z__[i__ + 1], abs(
 404                     d__2))) * (d__3 = ld[i__], abs(d__3)) < *gaptol) {
 405                 z__[i__ + 1] = 0.;
 406                 isuppz[2] = i__;
 407                 goto L280;
 408             }
 409             *ztz += z__[i__ + 1] * z__[i__ + 1];
 410 /* L270: */
 411         }
 412 L280:
 413         ;
 414     }
 415
 416 /*     Compute quantities for convergence test */
 417
 418     tmp = 1. / *ztz;
 419     *nrminv = sqrt(tmp);
 420     *resid = abs(*mingma) * *nrminv;
 421     *rqcorr = *mingma * tmp;
 422
 423
 424     return 0;
 425
 426 /*     End of DLAR1V */
 427
 428 } /* dlar1v_ */