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

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