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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 /* Subroutine */ int slarrf_(integer *n, real *d__, real *l, real *ld,
   8         integer *clstrt, integer *clend, real *w, real *wgap, real *werr,
   9         real *spdiam, real *clgapl, real *clgapr, real *pivmin, real *sigma,
  10         real *dplus, real *lplus, real *work, integer *info)
  11 {
  12     /* System generated locals */
  13     integer i__1;
  14     real r__1, r__2, r__3;
  15
  16     /* Builtin functions */
  17     double sqrt(doublereal);
  18
  19     /* Local variables */
  20     integer i__;
  21     real s, bestshift, smlgrowth, eps, tmp, max1, max2, rrr1, rrr2, znm2,
  22             growthbound, fail, fact, oldp;
  23     integer indx;
  24     real prod;
  25     integer ktry;
  26     real fail2, avgap, ldmax, rdmax;
  27     integer shift;
  28     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  29             integer *);
  30     logical dorrr1;
  31     real ldelta;
  32     extern doublereal slamch_(char *);
  33     logical nofail;
  34     real mingap, lsigma, rdelta;
  35     logical forcer;
  36     real rsigma, clwdth;
  37     extern logical sisnan_(real *);
  38     logical sawnan1, sawnan2, tryrrr1;
  39
  40
  41 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  42 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  43 /*     November 2006 */
  44 /* * */
  45 /*     .. Scalar Arguments .. */
  46 /*     .. */
  47 /*     .. Array Arguments .. */
  48 /*     .. */
  49
  50 /*  Purpose */
  51 /*  ======= */
  52
  53 /*  Given the initial representation L D L^T and its cluster of close */
  54 /*  eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ... */
  55 /*  W( CLEND ), SLARRF finds a new relatively robust representation */
  56 /*  L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the */
  57 /*  eigenvalues of L(+) D(+) L(+)^T is relatively isolated. */
  58
  59 /*  Arguments */
  60 /*  ========= */
  61
  62 /*  N       (input) INTEGER */
  63 /*          The order of the matrix (subblock, if the matrix splitted). */
  64
  65 /*  D       (input) REAL             array, dimension (N) */
  66 /*          The N diagonal elements of the diagonal matrix D. */
  67
  68 /*  L       (input) REAL             array, dimension (N-1) */
  69 /*          The (N-1) subdiagonal elements of the unit bidiagonal */
  70 /*          matrix L. */
  71
  72 /*  LD      (input) REAL             array, dimension (N-1) */
  73 /*          The (N-1) elements L(i)*D(i). */
  74
  75 /*  CLSTRT  (input) INTEGER */
  76 /*          The index of the first eigenvalue in the cluster. */
  77
  78 /*  CLEND   (input) INTEGER */
  79 /*          The index of the last eigenvalue in the cluster. */
  80
  81 /*  W       (input) REAL             array, dimension >=  (CLEND-CLSTRT+1) */
  82 /*          The eigenvalue APPROXIMATIONS of L D L^T in ascending order. */
  83 /*          W( CLSTRT ) through W( CLEND ) form the cluster of relatively */
  84 /*          close eigenalues. */
  85
  86 /*  WGAP    (input/output) REAL             array, dimension >=  (CLEND-CLSTRT+1) */
  87 /*          The separation from the right neighbor eigenvalue in W. */
  88
  89 /*  WERR    (input) REAL             array, dimension >=  (CLEND-CLSTRT+1) */
  90 /*          WERR contain the semiwidth of the uncertainty */
  91 /*          interval of the corresponding eigenvalue APPROXIMATION in W */
  92
  93 /*  SPDIAM (input) estimate of the spectral diameter obtained from the */
  94 /*          Gerschgorin intervals */
  95
  96 /*  CLGAPL, CLGAPR (input) absolute gap on each end of the cluster. */
  97 /*          Set by the calling routine to protect against shifts too close */
  98 /*          to eigenvalues outside the cluster. */
  99
 100 /*  PIVMIN  (input) DOUBLE PRECISION */
 101 /*          The minimum pivot allowed in the Sturm sequence. */
 102
 103 /*  SIGMA   (output) REAL */
 104 /*          The shift used to form L(+) D(+) L(+)^T. */
 105
 106 /*  DPLUS   (output) REAL             array, dimension (N) */
 107 /*          The N diagonal elements of the diagonal matrix D(+). */
 108
 109 /*  LPLUS   (output) REAL             array, dimension (N-1) */
 110 /*          The first (N-1) elements of LPLUS contain the subdiagonal */
 111 /*          elements of the unit bidiagonal matrix L(+). */
 112
 113 /*  WORK    (workspace) REAL             array, dimension (2*N) */
 114 /*          Workspace. */
 115
 116 /*  Further Details */
 117 /*  =============== */
 118
 119 /*  Based on contributions by */
 120 /*     Beresford Parlett, University of California, Berkeley, USA */
 121 /*     Jim Demmel, University of California, Berkeley, USA */
 122 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 123 /*     Osni Marques, LBNL/NERSC, USA */
 124 /*     Christof Voemel, University of California, Berkeley, USA */
 125
 126 /*  ===================================================================== */
 127
 128 /*     .. Parameters .. */
 129 /*     .. */
 130 /*     .. Local Scalars .. */
 131 /*     .. */
 132 /*     .. External Functions .. */
 133 /*     .. */
 134 /*     .. External Subroutines .. */
 135 /*     .. */
 136 /*     .. Intrinsic Functions .. */
 137 /*     .. */
 138 /*     .. Executable Statements .. */
 139
 140     /* Parameter adjustments */
 141     --work;
 142     --lplus;
 143     --dplus;
 144     --werr;
 145     --wgap;
 146     --w;
 147     --ld;
 148     --l;
 149     --d__;
 150
 151     /* Function Body */
 152     *info = 0;
 153     fact = 2.f;
 154     eps = slamch_("Precision");
 155     shift = 0;
 156     forcer = FALSE_;
 157 /*     Note that we cannot guarantee that for any of the shifts tried, */
 158 /*     the factorization has a small or even moderate element growth. */
 159 /*     There could be Ritz values at both ends of the cluster and despite */
 160 /*     backing off, there are examples where all factorizations tried */
 161 /*     (in IEEE mode, allowing zero pivots & infinities) have INFINITE */
 162 /*     element growth. */
 163 /*     For this reason, we should use PIVMIN in this subroutine so that at */
 164 /*     least the L D L^T factorization exists. It can be checked afterwards */
 165 /*     whether the element growth caused bad residuals/orthogonality. */
 166 /*     Decide whether the code should accept the best among all */
 167 /*     representations despite large element growth or signal INFO=1 */
 168     nofail = TRUE_;
 169
 170 /*     Compute the average gap length of the cluster */
 171     clwdth = (r__1 = w[*clend] - w[*clstrt], dabs(r__1)) + werr[*clend] +
 172             werr[*clstrt];
 173     avgap = clwdth / (real) (*clend - *clstrt);
 174     mingap = dmin(*clgapl,*clgapr);
 175 /*     Initial values for shifts to both ends of cluster */
 176 /* Computing MIN */
 177     r__1 = w[*clstrt], r__2 = w[*clend];
 178     lsigma = dmin(r__1,r__2) - werr[*clstrt];
 179 /* Computing MAX */
 180     r__1 = w[*clstrt], r__2 = w[*clend];
 181     rsigma = dmax(r__1,r__2) + werr[*clend];
 182 /*     Use a small fudge to make sure that we really shift to the outside */
 183     lsigma -= dabs(lsigma) * 2.f * eps;
 184     rsigma += dabs(rsigma) * 2.f * eps;
 185 /*     Compute upper bounds for how much to back off the initial shifts */
 186     ldmax = mingap * .25f + *pivmin * 2.f;
 187     rdmax = mingap * .25f + *pivmin * 2.f;
 188 /* Computing MAX */
 189     r__1 = avgap, r__2 = wgap[*clstrt];
 190     ldelta = dmax(r__1,r__2) / fact;
 191 /* Computing MAX */
 192     r__1 = avgap, r__2 = wgap[*clend - 1];
 193     rdelta = dmax(r__1,r__2) / fact;
 194
 195 /*     Initialize the record of the best representation found */
 196
 197     s = slamch_("S");
 198     smlgrowth = 1.f / s;
 199     fail = (real) (*n - 1) * mingap / (*spdiam * eps);
 200     fail2 = (real) (*n - 1) * mingap / (*spdiam * sqrt(eps));
 201     bestshift = lsigma;
 202
 203 /*     while (KTRY <= KTRYMAX) */
 204     ktry = 0;
 205     growthbound = *spdiam * 8.f;
 206 L5:
 207     sawnan1 = FALSE_;
 208     sawnan2 = FALSE_;
 209 /*     Ensure that we do not back off too much of the initial shifts */
 210     ldelta = dmin(ldmax,ldelta);
 211     rdelta = dmin(rdmax,rdelta);
 212 /*     Compute the element growth when shifting to both ends of the cluster */
 213 /*     accept the shift if there is no element growth at one of the two ends */
 214 /*     Left end */
 215     s = -lsigma;
 216     dplus[1] = d__[1] + s;
 217     if (dabs(dplus[1]) < *pivmin) {
 218         dplus[1] = -(*pivmin);
 219 /*        Need to set SAWNAN1 because refined RRR test should not be used */
 220 /*        in this case */
 221         sawnan1 = TRUE_;
 222     }
 223     max1 = dabs(dplus[1]);
 224     i__1 = *n - 1;
 225     for (i__ = 1; i__ <= i__1; ++i__) {
 226         lplus[i__] = ld[i__] / dplus[i__];
 227         s = s * lplus[i__] * l[i__] - lsigma;
 228         dplus[i__ + 1] = d__[i__ + 1] + s;
 229         if ((r__1 = dplus[i__ + 1], dabs(r__1)) < *pivmin) {
 230             dplus[i__ + 1] = -(*pivmin);
 231 /*           Need to set SAWNAN1 because refined RRR test should not be used */
 232 /*           in this case */
 233             sawnan1 = TRUE_;
 234         }
 235 /* Computing MAX */
 236         r__2 = max1, r__3 = (r__1 = dplus[i__ + 1], dabs(r__1));
 237         max1 = dmax(r__2,r__3);
 238 /* L6: */
 239     }
 240     sawnan1 = sawnan1 || sisnan_(&max1);
 241     if (forcer || max1 <= growthbound && ! sawnan1) {
 242         *sigma = lsigma;
 243         shift = 1;
 244         goto L100;
 245     }
 246 /*     Right end */
 247     s = -rsigma;
 248     work[1] = d__[1] + s;
 249     if (dabs(work[1]) < *pivmin) {
 250         work[1] = -(*pivmin);
 251 /*        Need to set SAWNAN2 because refined RRR test should not be used */
 252 /*        in this case */
 253         sawnan2 = TRUE_;
 254     }
 255     max2 = dabs(work[1]);
 256     i__1 = *n - 1;
 257     for (i__ = 1; i__ <= i__1; ++i__) {
 258         work[*n + i__] = ld[i__] / work[i__];
 259         s = s * work[*n + i__] * l[i__] - rsigma;
 260         work[i__ + 1] = d__[i__ + 1] + s;
 261         if ((r__1 = work[i__ + 1], dabs(r__1)) < *pivmin) {
 262             work[i__ + 1] = -(*pivmin);
 263 /*           Need to set SAWNAN2 because refined RRR test should not be used */
 264 /*           in this case */
 265             sawnan2 = TRUE_;
 266         }
 267 /* Computing MAX */
 268         r__2 = max2, r__3 = (r__1 = work[i__ + 1], dabs(r__1));
 269         max2 = dmax(r__2,r__3);
 270 /* L7: */
 271     }
 272     sawnan2 = sawnan2 || sisnan_(&max2);
 273     if (forcer || max2 <= growthbound && ! sawnan2) {
 274         *sigma = rsigma;
 275         shift = 2;
 276         goto L100;
 277     }
 278 /*     If we are at this point, both shifts led to too much element growth */
 279 /*     Record the better of the two shifts (provided it didn't lead to NaN) */
 280     if (sawnan1 && sawnan2) {
 281 /*        both MAX1 and MAX2 are NaN */
 282         goto L50;
 283     } else {
 284         if (! sawnan1) {
 285             indx = 1;
 286             if (max1 <= smlgrowth) {
 287                 smlgrowth = max1;
 288                 bestshift = lsigma;
 289             }
 290         }
 291         if (! sawnan2) {
 292             if (sawnan1 || max2 <= max1) {
 293                 indx = 2;
 294             }
 295             if (max2 <= smlgrowth) {
 296                 smlgrowth = max2;
 297                 bestshift = rsigma;
 298             }
 299         }
 300     }
 301 /*     If we are here, both the left and the right shift led to */
 302 /*     element growth. If the element growth is moderate, then */
 303 /*     we may still accept the representation, if it passes a */
 304 /*     refined test for RRR. This test supposes that no NaN occurred. */
 305 /*     Moreover, we use the refined RRR test only for isolated clusters. */
 306     if (clwdth < mingap / 128.f && dmin(max1,max2) < fail2 && ! sawnan1 && !
 307             sawnan2) {
 308         dorrr1 = TRUE_;
 309     } else {
 310         dorrr1 = FALSE_;
 311     }
 312     tryrrr1 = TRUE_;
 313     if (tryrrr1 && dorrr1) {
 314         if (indx == 1) {
 315             tmp = (r__1 = dplus[*n], dabs(r__1));
 316             znm2 = 1.f;
 317             prod = 1.f;
 318             oldp = 1.f;
 319             for (i__ = *n - 1; i__ >= 1; --i__) {
 320                 if (prod <= eps) {
 321                     prod = dplus[i__ + 1] * work[*n + i__ + 1] / (dplus[i__] *
 322                              work[*n + i__]) * oldp;
 323                 } else {
 324                     prod *= (r__1 = work[*n + i__], dabs(r__1));
 325                 }
 326                 oldp = prod;
 327 /* Computing 2nd power */
 328                 r__1 = prod;
 329                 znm2 += r__1 * r__1;
 330 /* Computing MAX */
 331                 r__2 = tmp, r__3 = (r__1 = dplus[i__] * prod, dabs(r__1));
 332                 tmp = dmax(r__2,r__3);
 333 /* L15: */
 334             }
 335             rrr1 = tmp / (*spdiam * sqrt(znm2));
 336             if (rrr1 <= 8.f) {
 337                 *sigma = lsigma;
 338                 shift = 1;
 339                 goto L100;
 340             }
 341         } else if (indx == 2) {
 342             tmp = (r__1 = work[*n], dabs(r__1));
 343             znm2 = 1.f;
 344             prod = 1.f;
 345             oldp = 1.f;
 346             for (i__ = *n - 1; i__ >= 1; --i__) {
 347                 if (prod <= eps) {
 348                     prod = work[i__ + 1] * lplus[i__ + 1] / (work[i__] *
 349                             lplus[i__]) * oldp;
 350                 } else {
 351                     prod *= (r__1 = lplus[i__], dabs(r__1));
 352                 }
 353                 oldp = prod;
 354 /* Computing 2nd power */
 355                 r__1 = prod;
 356                 znm2 += r__1 * r__1;
 357 /* Computing MAX */
 358                 r__2 = tmp, r__3 = (r__1 = work[i__] * prod, dabs(r__1));
 359                 tmp = dmax(r__2,r__3);
 360 /* L16: */
 361             }
 362             rrr2 = tmp / (*spdiam * sqrt(znm2));
 363             if (rrr2 <= 8.f) {
 364                 *sigma = rsigma;
 365                 shift = 2;
 366                 goto L100;
 367             }
 368         }
 369     }
 370 L50:
 371     if (ktry < 1) {
 372 /*        If we are here, both shifts failed also the RRR test. */
 373 /*        Back off to the outside */
 374 /* Computing MAX */
 375         r__1 = lsigma - ldelta, r__2 = lsigma - ldmax;
 376         lsigma = dmax(r__1,r__2);
 377 /* Computing MIN */
 378         r__1 = rsigma + rdelta, r__2 = rsigma + rdmax;
 379         rsigma = dmin(r__1,r__2);
 380         ldelta *= 2.f;
 381         rdelta *= 2.f;
 382         ++ktry;
 383         goto L5;
 384     } else {
 385 /*        None of the representations investigated satisfied our */
 386 /*        criteria. Take the best one we found. */
 387         if (smlgrowth < fail || nofail) {
 388             lsigma = bestshift;
 389             rsigma = bestshift;
 390             forcer = TRUE_;
 391             goto L5;
 392         } else {
 393             *info = 1;
 394             return 0;
 395         }
 396     }
 397 L100:
 398     if (shift == 1) {
 399     } else if (shift == 2) {
 400 /*        store new L and D back into DPLUS, LPLUS */
 401         scopy_(n, &work[1], &c__1, &dplus[1], &c__1);
 402         i__1 = *n - 1;
 403         scopy_(&i__1, &work[*n + 1], &c__1, &lplus[1], &c__1);
 404     }
 405     return 0;
 406
 407 /*     End of SLARRF */
 408
 409 } /* slarrf_ */