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

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