git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlarre.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
   8 /* Subroutine */ int dlarre_(char *range, integer *n, doublereal *vl,
   9         doublereal *vu, integer *il, integer *iu, doublereal *d__, doublereal
  10         *e, doublereal *e2, doublereal *rtol1, doublereal *rtol2, doublereal *
  11         spltol, integer *nsplit, integer *isplit, integer *m, doublereal *w,
  12         doublereal *werr, doublereal *wgap, integer *iblock, integer *indexw,
  13         doublereal *gers, doublereal *pivmin, doublereal *work, integer *
  14         iwork, integer *info)
  15 {
  16     /* System generated locals */
  17     integer i__1, i__2;
  18     doublereal d__1, d__2, d__3;
  19
  20     /* Builtin functions */
  21     double sqrt(doublereal), log(doublereal);
  22
  23     /* Local variables */
  24     integer i__, j;
  25     doublereal s1, s2;
  26     integer mb;
  27     doublereal gl;
  28     integer in, mm;
  29     doublereal gu;
  30     integer cnt;
  31     doublereal eps, tau, tmp, rtl;
  32     integer cnt1, cnt2;
  33     doublereal tmp1, eabs;
  34     integer iend, jblk;
  35     doublereal eold;
  36     integer indl;
  37     doublereal dmax__, emax;
  38     integer wend, idum, indu;
  39     doublereal rtol;
  40     integer iseed[4];
  41     doublereal avgap, sigma;
  42     extern logical lsame_(char *, char *);
  43     integer iinfo;
  44     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  45             doublereal *, integer *);
  46     logical norep;
  47     extern /* Subroutine */ int dlasq2_(integer *, doublereal *, integer *);
  48     extern doublereal dlamch_(char *);
  49     integer ibegin;
  50     logical forceb;
  51     integer irange;
  52     doublereal sgndef;
  53     extern /* Subroutine */ int dlarra_(integer *, doublereal *, doublereal *,
  54              doublereal *, doublereal *, doublereal *, integer *, integer *,
  55             integer *), dlarrb_(integer *, doublereal *, doublereal *,
  56             integer *, integer *, doublereal *, doublereal *, integer *,
  57             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  58              doublereal *, doublereal *, integer *, integer *), dlarrc_(char *
  59 , integer *, doublereal *, doublereal *, doublereal *, doublereal
  60             *, doublereal *, integer *, integer *, integer *, integer *);
  61     integer wbegin;
  62     extern /* Subroutine */ int dlarrd_(char *, char *, integer *, doublereal
  63             *, doublereal *, integer *, integer *, doublereal *, doublereal *,
  64              doublereal *, doublereal *, doublereal *, doublereal *, integer *
  65 , integer *, integer *, doublereal *, doublereal *, doublereal *,
  66             doublereal *, integer *, integer *, doublereal *, integer *,
  67             integer *);
  68     doublereal safmin, spdiam;
  69     extern /* Subroutine */ int dlarrk_(integer *, integer *, doublereal *,
  70             doublereal *, doublereal *, doublereal *, doublereal *,
  71             doublereal *, doublereal *, doublereal *, integer *);
  72     logical usedqd;
  73     doublereal clwdth, isleft;
  74     extern /* Subroutine */ int dlarnv_(integer *, integer *, integer *,
  75             doublereal *);
  76     doublereal isrght, bsrtol, dpivot;
  77
  78
  79 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  80 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  81 /*     November 2006 */
  82
  83 /*     .. Scalar Arguments .. */
  84 /*     .. */
  85 /*     .. Array Arguments .. */
  86 /*     .. */
  87
  88 /*  Purpose */
  89 /*  ======= */
  90
  91 /*  To find the desired eigenvalues of a given real symmetric */
  92 /*  tridiagonal matrix T, DLARRE sets any "small" off-diagonal */
  93 /*  elements to zero, and for each unreduced block T_i, it finds */
  94 /*  (a) a suitable shift at one end of the block's spectrum, */
  95 /*  (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and */
  96 /*  (c) eigenvalues of each L_i D_i L_i^T. */
  97 /*  The representations and eigenvalues found are then used by */
  98 /*  DSTEMR to compute the eigenvectors of T. */
  99 /*  The accuracy varies depending on whether bisection is used to */
 100 /*  find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to */
 101 /*  conpute all and then discard any unwanted one. */
 102 /*  As an added benefit, DLARRE also outputs the n */
 103 /*  Gerschgorin intervals for the matrices L_i D_i L_i^T. */
 104
 105 /*  Arguments */
 106 /*  ========= */
 107
 108 /*  RANGE   (input) CHARACTER */
 109 /*          = 'A': ("All")   all eigenvalues will be found. */
 110 /*          = 'V': ("Value") all eigenvalues in the half-open interval */
 111 /*                           (VL, VU] will be found. */
 112 /*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
 113 /*                           entire matrix) will be found. */
 114
 115 /*  N       (input) INTEGER */
 116 /*          The order of the matrix. N > 0. */
 117
 118 /*  VL      (input/output) DOUBLE PRECISION */
 119 /*  VU      (input/output) DOUBLE PRECISION */
 120 /*          If RANGE='V', the lower and upper bounds for the eigenvalues. */
 121 /*          Eigenvalues less than or equal to VL, or greater than VU, */
 122 /*          will not be returned.  VL < VU. */
 123 /*          If RANGE='I' or ='A', DLARRE computes bounds on the desired */
 124 /*          part of the spectrum. */
 125
 126 /*  IL      (input) INTEGER */
 127 /*  IU      (input) INTEGER */
 128 /*          If RANGE='I', the indices (in ascending order) of the */
 129 /*          smallest and largest eigenvalues to be returned. */
 130 /*          1 <= IL <= IU <= N. */
 131
 132 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
 133 /*          On entry, the N diagonal elements of the tridiagonal */
 134 /*          matrix T. */
 135 /*          On exit, the N diagonal elements of the diagonal */
 136 /*          matrices D_i. */
 137
 138 /*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
 139 /*          On entry, the first (N-1) entries contain the subdiagonal */
 140 /*          elements of the tridiagonal matrix T; E(N) need not be set. */
 141 /*          On exit, E contains the subdiagonal elements of the unit */
 142 /*          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), */
 143 /*          1 <= I <= NSPLIT, contain the base points sigma_i on output. */
 144
 145 /*  E2      (input/output) DOUBLE PRECISION array, dimension (N) */
 146 /*          On entry, the first (N-1) entries contain the SQUARES of the */
 147 /*          subdiagonal elements of the tridiagonal matrix T; */
 148 /*          E2(N) need not be set. */
 149 /*          On exit, the entries E2( ISPLIT( I ) ), */
 150 /*          1 <= I <= NSPLIT, have been set to zero */
 151
 152 /*  RTOL1   (input) DOUBLE PRECISION */
 153 /*  RTOL2   (input) DOUBLE PRECISION */
 154 /*           Parameters for bisection. */
 155 /*           An interval [LEFT,RIGHT] has converged if */
 156 /*           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
 157
 158 /*  SPLTOL (input) DOUBLE PRECISION */
 159 /*          The threshold for splitting. */
 160
 161 /*  NSPLIT  (output) INTEGER */
 162 /*          The number of blocks T splits into. 1 <= NSPLIT <= N. */
 163
 164 /*  ISPLIT  (output) INTEGER array, dimension (N) */
 165 /*          The splitting points, at which T breaks up into blocks. */
 166 /*          The first block consists of rows/columns 1 to ISPLIT(1), */
 167 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
 168 /*          etc., and the NSPLIT-th consists of rows/columns */
 169 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
 170
 171 /*  M       (output) INTEGER */
 172 /*          The total number of eigenvalues (of all L_i D_i L_i^T) */
 173 /*          found. */
 174
 175 /*  W       (output) DOUBLE PRECISION array, dimension (N) */
 176 /*          The first M elements contain the eigenvalues. The */
 177 /*          eigenvalues of each of the blocks, L_i D_i L_i^T, are */
 178 /*          sorted in ascending order ( DLARRE may use the */
 179 /*          remaining N-M elements as workspace). */
 180
 181 /*  WERR    (output) DOUBLE PRECISION array, dimension (N) */
 182 /*          The error bound on the corresponding eigenvalue in W. */
 183
 184 /*  WGAP    (output) DOUBLE PRECISION array, dimension (N) */
 185 /*          The separation from the right neighbor eigenvalue in W. */
 186 /*          The gap is only with respect to the eigenvalues of the same block */
 187 /*          as each block has its own representation tree. */
 188 /*          Exception: at the right end of a block we store the left gap */
 189
 190 /*  IBLOCK  (output) INTEGER array, dimension (N) */
 191 /*          The indices of the blocks (submatrices) associated with the */
 192 /*          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue */
 193 /*          W(i) belongs to the first block from the top, =2 if W(i) */
 194 /*          belongs to the second block, etc. */
 195
 196 /*  INDEXW  (output) INTEGER array, dimension (N) */
 197 /*          The indices of the eigenvalues within each block (submatrix); */
 198 /*          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the */
 199 /*          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2 */
 200
 201 /*  GERS    (output) DOUBLE PRECISION array, dimension (2*N) */
 202 /*          The N Gerschgorin intervals (the i-th Gerschgorin interval */
 203 /*          is (GERS(2*i-1), GERS(2*i)). */
 204
 205 /*  PIVMIN  (output) DOUBLE PRECISION */
 206 /*          The minimum pivot in the Sturm sequence for T. */
 207
 208 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (6*N) */
 209 /*          Workspace. */
 210
 211 /*  IWORK   (workspace) INTEGER array, dimension (5*N) */
 212 /*          Workspace. */
 213
 214 /*  INFO    (output) INTEGER */
 215 /*          = 0:  successful exit */
 216 /*          > 0:  A problem occured in DLARRE. */
 217 /*          < 0:  One of the called subroutines signaled an internal problem. */
 218 /*                Needs inspection of the corresponding parameter IINFO */
 219 /*                for further information. */
 220
 221 /*          =-1:  Problem in DLARRD. */
 222 /*          = 2:  No base representation could be found in MAXTRY iterations. */
 223 /*                Increasing MAXTRY and recompilation might be a remedy. */
 224 /*          =-3:  Problem in DLARRB when computing the refined root */
 225 /*                representation for DLASQ2. */
 226 /*          =-4:  Problem in DLARRB when preforming bisection on the */
 227 /*                desired part of the spectrum. */
 228 /*          =-5:  Problem in DLASQ2. */
 229 /*          =-6:  Problem in DLASQ2. */
 230
 231 /*  Further Details */
 232 /*  The base representations are required to suffer very little */
 233 /*  element growth and consequently define all their eigenvalues to */
 234 /*  high relative accuracy. */
 235 /*  =============== */
 236
 237 /*  Based on contributions by */
 238 /*     Beresford Parlett, University of California, Berkeley, USA */
 239 /*     Jim Demmel, University of California, Berkeley, USA */
 240 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 241 /*     Osni Marques, LBNL/NERSC, USA */
 242 /*     Christof Voemel, University of California, Berkeley, USA */
 243
 244 /*  ===================================================================== */
 245
 246 /*     .. Parameters .. */
 247 /*     .. */
 248 /*     .. Local Scalars .. */
 249 /*     .. */
 250 /*     .. Local Arrays .. */
 251 /*     .. */
 252 /*     .. External Functions .. */
 253 /*     .. */
 254 /*     .. External Subroutines .. */
 255 /*     .. */
 256 /*     .. Intrinsic Functions .. */
 257 /*     .. */
 258 /*     .. Executable Statements .. */
 259
 260     /* Parameter adjustments */
 261     --iwork;
 262     --work;
 263     --gers;
 264     --indexw;
 265     --iblock;
 266     --wgap;
 267     --werr;
 268     --w;
 269     --isplit;
 270     --e2;
 271     --e;
 272     --d__;
 273
 274     /* Function Body */
 275     *info = 0;
 276
 277 /*     Decode RANGE */
 278
 279     if (lsame_(range, "A")) {
 280         irange = 1;
 281     } else if (lsame_(range, "V")) {
 282         irange = 3;
 283     } else if (lsame_(range, "I")) {
 284         irange = 2;
 285     }
 286     *m = 0;
 287 /*     Get machine constants */
 288     safmin = dlamch_("S");
 289     eps = dlamch_("P");
 290 /*     Set parameters */
 291     rtl = sqrt(eps);
 292     bsrtol = sqrt(eps);
 293 /*     Treat case of 1x1 matrix for quick return */
 294     if (*n == 1) {
 295         if (irange == 1 || irange == 3 && d__[1] > *vl && d__[1] <= *vu ||
 296                 irange == 2 && *il == 1 && *iu == 1) {
 297             *m = 1;
 298             w[1] = d__[1];
 299 /*           The computation error of the eigenvalue is zero */
 300             werr[1] = 0.;
 301             wgap[1] = 0.;
 302             iblock[1] = 1;
 303             indexw[1] = 1;
 304             gers[1] = d__[1];
 305             gers[2] = d__[1];
 306         }
 307 /*        store the shift for the initial RRR, which is zero in this case */
 308         e[1] = 0.;
 309         return 0;
 310     }
 311 /*     General case: tridiagonal matrix of order > 1 */
 312
 313 /*     Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter. */
 314 /*     Compute maximum off-diagonal entry and pivmin. */
 315     gl = d__[1];
 316     gu = d__[1];
 317     eold = 0.;
 318     emax = 0.;
 319     e[*n] = 0.;
 320     i__1 = *n;
 321     for (i__ = 1; i__ <= i__1; ++i__) {
 322         werr[i__] = 0.;
 323         wgap[i__] = 0.;
 324         eabs = (d__1 = e[i__], abs(d__1));
 325         if (eabs >= emax) {
 326             emax = eabs;
 327         }
 328         tmp1 = eabs + eold;
 329         gers[(i__ << 1) - 1] = d__[i__] - tmp1;
 330 /* Computing MIN */
 331         d__1 = gl, d__2 = gers[(i__ << 1) - 1];
 332         gl = min(d__1,d__2);
 333         gers[i__ * 2] = d__[i__] + tmp1;
 334 /* Computing MAX */
 335         d__1 = gu, d__2 = gers[i__ * 2];
 336         gu = max(d__1,d__2);
 337         eold = eabs;
 338 /* L5: */
 339     }
 340 /*     The minimum pivot allowed in the Sturm sequence for T */
 341 /* Computing MAX */
 342 /* Computing 2nd power */
 343     d__3 = emax;
 344     d__1 = 1., d__2 = d__3 * d__3;
 345     *pivmin = safmin * max(d__1,d__2);
 346 /*     Compute spectral diameter. The Gerschgorin bounds give an */
 347 /*     estimate that is wrong by at most a factor of SQRT(2) */
 348     spdiam = gu - gl;
 349 /*     Compute splitting points */
 350     dlarra_(n, &d__[1], &e[1], &e2[1], spltol, &spdiam, nsplit, &isplit[1], &
 351             iinfo);
 352 /*     Can force use of bisection instead of faster DQDS. */
 353 /*     Option left in the code for future multisection work. */
 354     forceb = FALSE_;
 355     if (irange == 1 && ! forceb) {
 356 /*        Set interval [VL,VU] that contains all eigenvalues */
 357         *vl = gl;
 358         *vu = gu;
 359     } else {
 360 /*        We call DLARRD to find crude approximations to the eigenvalues */
 361 /*        in the desired range. In case IRANGE = INDRNG, we also obtain the */
 362 /*        interval (VL,VU] that contains all the wanted eigenvalues. */
 363 /*        An interval [LEFT,RIGHT] has converged if */
 364 /*        RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT)) */
 365 /*        DLARRD needs a WORK of size 4*N, IWORK of size 3*N */
 366         dlarrd_(range, "B", n, vl, vu, il, iu, &gers[1], &bsrtol, &d__[1], &e[
 367                 1], &e2[1], pivmin, nsplit, &isplit[1], &mm, &w[1], &werr[1],
 368                 vl, vu, &iblock[1], &indexw[1], &work[1], &iwork[1], &iinfo);
 369         if (iinfo != 0) {
 370             *info = -1;
 371             return 0;
 372         }
 373 /*        Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0 */
 374         i__1 = *n;
 375         for (i__ = mm + 1; i__ <= i__1; ++i__) {
 376             w[i__] = 0.;
 377             werr[i__] = 0.;
 378             iblock[i__] = 0;
 379             indexw[i__] = 0;
 380 /* L14: */
 381         }
 382     }
 383 /* ** */
 384 /*     Loop over unreduced blocks */
 385     ibegin = 1;
 386     wbegin = 1;
 387     i__1 = *nsplit;
 388     for (jblk = 1; jblk <= i__1; ++jblk) {
 389         iend = isplit[jblk];
 390         in = iend - ibegin + 1;
 391 /*        1 X 1 block */
 392         if (in == 1) {
 393             if (irange == 1 || irange == 3 && d__[ibegin] > *vl && d__[ibegin]
 394                      <= *vu || irange == 2 && iblock[wbegin] == jblk) {
 395                 ++(*m);
 396                 w[*m] = d__[ibegin];
 397                 werr[*m] = 0.;
 398 /*              The gap for a single block doesn't matter for the later */
 399 /*              algorithm and is assigned an arbitrary large value */
 400                 wgap[*m] = 0.;
 401                 iblock[*m] = jblk;
 402                 indexw[*m] = 1;
 403                 ++wbegin;
 404             }
 405 /*           E( IEND ) holds the shift for the initial RRR */
 406             e[iend] = 0.;
 407             ibegin = iend + 1;
 408             goto L170;
 409         }
 410
 411 /*        Blocks of size larger than 1x1 */
 412
 413 /*        E( IEND ) will hold the shift for the initial RRR, for now set it =0 */
 414         e[iend] = 0.;
 415
 416 /*        Find local outer bounds GL,GU for the block */
 417         gl = d__[ibegin];
 418         gu = d__[ibegin];
 419         i__2 = iend;
 420         for (i__ = ibegin; i__ <= i__2; ++i__) {
 421 /* Computing MIN */
 422             d__1 = gers[(i__ << 1) - 1];
 423             gl = min(d__1,gl);
 424 /* Computing MAX */
 425             d__1 = gers[i__ * 2];
 426             gu = max(d__1,gu);
 427 /* L15: */
 428         }
 429         spdiam = gu - gl;
 430         if (! (irange == 1 && ! forceb)) {
 431 /*           Count the number of eigenvalues in the current block. */
 432             mb = 0;
 433             i__2 = mm;
 434             for (i__ = wbegin; i__ <= i__2; ++i__) {
 435                 if (iblock[i__] == jblk) {
 436                     ++mb;
 437                 } else {
 438                     goto L21;
 439                 }
 440 /* L20: */
 441             }
 442 L21:
 443             if (mb == 0) {
 444 /*              No eigenvalue in the current block lies in the desired range */
 445 /*              E( IEND ) holds the shift for the initial RRR */
 446                 e[iend] = 0.;
 447                 ibegin = iend + 1;
 448                 goto L170;
 449             } else {
 450 /*              Decide whether dqds or bisection is more efficient */
 451                 usedqd = (doublereal) mb > in * .5 && ! forceb;
 452                 wend = wbegin + mb - 1;
 453 /*              Calculate gaps for the current block */
 454 /*              In later stages, when representations for individual */
 455 /*              eigenvalues are different, we use SIGMA = E( IEND ). */
 456                 sigma = 0.;
 457                 i__2 = wend - 1;
 458                 for (i__ = wbegin; i__ <= i__2; ++i__) {
 459 /* Computing MAX */
 460                     d__1 = 0., d__2 = w[i__ + 1] - werr[i__ + 1] - (w[i__] +
 461                             werr[i__]);
 462                     wgap[i__] = max(d__1,d__2);
 463 /* L30: */
 464                 }
 465 /* Computing MAX */
 466                 d__1 = 0., d__2 = *vu - sigma - (w[wend] + werr[wend]);
 467                 wgap[wend] = max(d__1,d__2);
 468 /*              Find local index of the first and last desired evalue. */
 469                 indl = indexw[wbegin];
 470                 indu = indexw[wend];
 471             }
 472         }
 473         if (irange == 1 && ! forceb || usedqd) {
 474 /*           Case of DQDS */
 475 /*           Find approximations to the extremal eigenvalues of the block */
 476             dlarrk_(&in, &c__1, &gl, &gu, &d__[ibegin], &e2[ibegin], pivmin, &
 477                     rtl, &tmp, &tmp1, &iinfo);
 478             if (iinfo != 0) {
 479                 *info = -1;
 480                 return 0;
 481             }
 482 /* Computing MAX */
 483             d__2 = gl, d__3 = tmp - tmp1 - eps * 100. * (d__1 = tmp - tmp1,
 484                     abs(d__1));
 485             isleft = max(d__2,d__3);
 486             dlarrk_(&in, &in, &gl, &gu, &d__[ibegin], &e2[ibegin], pivmin, &
 487                     rtl, &tmp, &tmp1, &iinfo);
 488             if (iinfo != 0) {
 489                 *info = -1;
 490                 return 0;
 491             }
 492 /* Computing MIN */
 493             d__2 = gu, d__3 = tmp + tmp1 + eps * 100. * (d__1 = tmp + tmp1,
 494                     abs(d__1));
 495             isrght = min(d__2,d__3);
 496 /*           Improve the estimate of the spectral diameter */
 497             spdiam = isrght - isleft;
 498         } else {
 499 /*           Case of bisection */
 500 /*           Find approximations to the wanted extremal eigenvalues */
 501 /* Computing MAX */
 502             d__2 = gl, d__3 = w[wbegin] - werr[wbegin] - eps * 100. * (d__1 =
 503                     w[wbegin] - werr[wbegin], abs(d__1));
 504             isleft = max(d__2,d__3);
 505 /* Computing MIN */
 506             d__2 = gu, d__3 = w[wend] + werr[wend] + eps * 100. * (d__1 = w[
 507                     wend] + werr[wend], abs(d__1));
 508             isrght = min(d__2,d__3);
 509         }
 510 /*        Decide whether the base representation for the current block */
 511 /*        L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I */
 512 /*        should be on the left or the right end of the current block. */
 513 /*        The strategy is to shift to the end which is "more populated" */
 514 /*        Furthermore, decide whether to use DQDS for the computation of */
 515 /*        the eigenvalue approximations at the end of DLARRE or bisection. */
 516 /*        dqds is chosen if all eigenvalues are desired or the number of */
 517 /*        eigenvalues to be computed is large compared to the blocksize. */
 518         if (irange == 1 && ! forceb) {
 519 /*           If all the eigenvalues have to be computed, we use dqd */
 520             usedqd = TRUE_;
 521 /*           INDL is the local index of the first eigenvalue to compute */
 522             indl = 1;
 523             indu = in;
 524 /*           MB =  number of eigenvalues to compute */
 525             mb = in;
 526             wend = wbegin + mb - 1;
 527 /*           Define 1/4 and 3/4 points of the spectrum */
 528             s1 = isleft + spdiam * .25;
 529             s2 = isrght - spdiam * .25;
 530         } else {
 531 /*           DLARRD has computed IBLOCK and INDEXW for each eigenvalue */
 532 /*           approximation. */
 533 /*           choose sigma */
 534             if (usedqd) {
 535                 s1 = isleft + spdiam * .25;
 536                 s2 = isrght - spdiam * .25;
 537             } else {
 538                 tmp = min(isrght,*vu) - max(isleft,*vl);
 539                 s1 = max(isleft,*vl) + tmp * .25;
 540                 s2 = min(isrght,*vu) - tmp * .25;
 541             }
 542         }
 543 /*        Compute the negcount at the 1/4 and 3/4 points */
 544         if (mb > 1) {
 545             dlarrc_("T", &in, &s1, &s2, &d__[ibegin], &e[ibegin], pivmin, &
 546                     cnt, &cnt1, &cnt2, &iinfo);
 547         }
 548         if (mb == 1) {
 549             sigma = gl;
 550             sgndef = 1.;
 551         } else if (cnt1 - indl >= indu - cnt2) {
 552             if (irange == 1 && ! forceb) {
 553                 sigma = max(isleft,gl);
 554             } else if (usedqd) {
 555 /*              use Gerschgorin bound as shift to get pos def matrix */
 556 /*              for dqds */
 557                 sigma = isleft;
 558             } else {
 559 /*              use approximation of the first desired eigenvalue of the */
 560 /*              block as shift */
 561                 sigma = max(isleft,*vl);
 562             }
 563             sgndef = 1.;
 564         } else {
 565             if (irange == 1 && ! forceb) {
 566                 sigma = min(isrght,gu);
 567             } else if (usedqd) {
 568 /*              use Gerschgorin bound as shift to get neg def matrix */
 569 /*              for dqds */
 570                 sigma = isrght;
 571             } else {
 572 /*              use approximation of the first desired eigenvalue of the */
 573 /*              block as shift */
 574                 sigma = min(isrght,*vu);
 575             }
 576             sgndef = -1.;
 577         }
 578 /*        An initial SIGMA has been chosen that will be used for computing */
 579 /*        T - SIGMA I = L D L^T */
 580 /*        Define the increment TAU of the shift in case the initial shift */
 581 /*        needs to be refined to obtain a factorization with not too much */
 582 /*        element growth. */
 583         if (usedqd) {
 584 /*           The initial SIGMA was to the outer end of the spectrum */
 585 /*           the matrix is definite and we need not retreat. */
 586             tau = spdiam * eps * *n + *pivmin * 2.;
 587         } else {
 588             if (mb > 1) {
 589                 clwdth = w[wend] + werr[wend] - w[wbegin] - werr[wbegin];
 590                 avgap = (d__1 = clwdth / (doublereal) (wend - wbegin), abs(
 591                         d__1));
 592                 if (sgndef == 1.) {
 593 /* Computing MAX */
 594                     d__1 = wgap[wbegin];
 595                     tau = max(d__1,avgap) * .5;
 596 /* Computing MAX */
 597                     d__1 = tau, d__2 = werr[wbegin];
 598                     tau = max(d__1,d__2);
 599                 } else {
 600 /* Computing MAX */
 601                     d__1 = wgap[wend - 1];
 602                     tau = max(d__1,avgap) * .5;
 603 /* Computing MAX */
 604                     d__1 = tau, d__2 = werr[wend];
 605                     tau = max(d__1,d__2);
 606                 }
 607             } else {
 608                 tau = werr[wbegin];
 609             }
 610         }
 611
 612         for (idum = 1; idum <= 6; ++idum) {
 613 /*           Compute L D L^T factorization of tridiagonal matrix T - sigma I. */
 614 /*           Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of */
 615 /*           pivots in WORK(2*IN+1:3*IN) */
 616             dpivot = d__[ibegin] - sigma;
 617             work[1] = dpivot;
 618             dmax__ = abs(work[1]);
 619             j = ibegin;
 620             i__2 = in - 1;
 621             for (i__ = 1; i__ <= i__2; ++i__) {
 622                 work[(in << 1) + i__] = 1. / work[i__];
 623                 tmp = e[j] * work[(in << 1) + i__];
 624                 work[in + i__] = tmp;
 625                 dpivot = d__[j + 1] - sigma - tmp * e[j];
 626                 work[i__ + 1] = dpivot;
 627 /* Computing MAX */
 628                 d__1 = dmax__, d__2 = abs(dpivot);
 629                 dmax__ = max(d__1,d__2);
 630                 ++j;
 631 /* L70: */
 632             }
 633 /*           check for element growth */
 634             if (dmax__ > spdiam * 64.) {
 635                 norep = TRUE_;
 636             } else {
 637                 norep = FALSE_;
 638             }
 639             if (usedqd && ! norep) {
 640 /*              Ensure the definiteness of the representation */
 641 /*              All entries of D (of L D L^T) must have the same sign */
 642                 i__2 = in;
 643                 for (i__ = 1; i__ <= i__2; ++i__) {
 644                     tmp = sgndef * work[i__];
 645                     if (tmp < 0.) {
 646                         norep = TRUE_;
 647                     }
 648 /* L71: */
 649                 }
 650             }
 651             if (norep) {
 652 /*              Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin */
 653 /*              shift which makes the matrix definite. So we should end up */
 654 /*              here really only in the case of IRANGE = VALRNG or INDRNG. */
 655                 if (idum == 5) {
 656                     if (sgndef == 1.) {
 657 /*                    The fudged Gerschgorin shift should succeed */
 658                         sigma = gl - spdiam * 2. * eps * *n - *pivmin * 4.;
 659                     } else {
 660                         sigma = gu + spdiam * 2. * eps * *n + *pivmin * 4.;
 661                     }
 662                 } else {
 663                     sigma -= sgndef * tau;
 664                     tau *= 2.;
 665                 }
 666             } else {
 667 /*              an initial RRR is found */
 668                 goto L83;
 669             }
 670 /* L80: */
 671         }
 672 /*        if the program reaches this point, no base representation could be */
 673 /*        found in MAXTRY iterations. */
 674         *info = 2;
 675         return 0;
 676 L83:
 677 /*        At this point, we have found an initial base representation */
 678 /*        T - SIGMA I = L D L^T with not too much element growth. */
 679 /*        Store the shift. */
 680         e[iend] = sigma;
 681 /*        Store D and L. */
 682         dcopy_(&in, &work[1], &c__1, &d__[ibegin], &c__1);
 683         i__2 = in - 1;
 684         dcopy_(&i__2, &work[in + 1], &c__1, &e[ibegin], &c__1);
 685         if (mb > 1) {
 686
 687 /*           Perturb each entry of the base representation by a small */
 688 /*           (but random) relative amount to overcome difficulties with */
 689 /*           glued matrices. */
 690
 691             for (i__ = 1; i__ <= 4; ++i__) {
 692                 iseed[i__ - 1] = 1;
 693 /* L122: */
 694             }
 695             i__2 = (in << 1) - 1;
 696             dlarnv_(&c__2, iseed, &i__2, &work[1]);
 697             i__2 = in - 1;
 698             for (i__ = 1; i__ <= i__2; ++i__) {
 699                 d__[ibegin + i__ - 1] *= eps * 8. * work[i__] + 1.;
 700                 e[ibegin + i__ - 1] *= eps * 8. * work[in + i__] + 1.;
 701 /* L125: */
 702             }
 703             d__[iend] *= eps * 4. * work[in] + 1.;
 704
 705         }
 706
 707 /*        Don't update the Gerschgorin intervals because keeping track */
 708 /*        of the updates would be too much work in DLARRV. */
 709 /*        We update W instead and use it to locate the proper Gerschgorin */
 710 /*        intervals. */
 711 /*        Compute the required eigenvalues of L D L' by bisection or dqds */
 712         if (! usedqd) {
 713 /*           If DLARRD has been used, shift the eigenvalue approximations */
 714 /*           according to their representation. This is necessary for */
 715 /*           a uniform DLARRV since dqds computes eigenvalues of the */
 716 /*           shifted representation. In DLARRV, W will always hold the */
 717 /*           UNshifted eigenvalue approximation. */
 718             i__2 = wend;
 719             for (j = wbegin; j <= i__2; ++j) {
 720                 w[j] -= sigma;
 721                 werr[j] += (d__1 = w[j], abs(d__1)) * eps;
 722 /* L134: */
 723             }
 724 /*           call DLARRB to reduce eigenvalue error of the approximations */
 725 /*           from DLARRD */
 726             i__2 = iend - 1;
 727             for (i__ = ibegin; i__ <= i__2; ++i__) {
 728 /* Computing 2nd power */
 729                 d__1 = e[i__];
 730                 work[i__] = d__[i__] * (d__1 * d__1);
 731 /* L135: */
 732             }
 733 /*           use bisection to find EV from INDL to INDU */
 734             i__2 = indl - 1;
 735             dlarrb_(&in, &d__[ibegin], &work[ibegin], &indl, &indu, rtol1,
 736                     rtol2, &i__2, &w[wbegin], &wgap[wbegin], &werr[wbegin], &
 737                     work[(*n << 1) + 1], &iwork[1], pivmin, &spdiam, &in, &
 738                     iinfo);
 739             if (iinfo != 0) {
 740                 *info = -4;
 741                 return 0;
 742             }
 743 /*           DLARRB computes all gaps correctly except for the last one */
 744 /*           Record distance to VU/GU */
 745 /* Computing MAX */
 746             d__1 = 0., d__2 = *vu - sigma - (w[wend] + werr[wend]);
 747             wgap[wend] = max(d__1,d__2);
 748             i__2 = indu;
 749             for (i__ = indl; i__ <= i__2; ++i__) {
 750                 ++(*m);
 751                 iblock[*m] = jblk;
 752                 indexw[*m] = i__;
 753 /* L138: */
 754             }
 755         } else {
 756 /*           Call dqds to get all eigs (and then possibly delete unwanted */
 757 /*           eigenvalues). */
 758 /*           Note that dqds finds the eigenvalues of the L D L^T representation */
 759 /*           of T to high relative accuracy. High relative accuracy */
 760 /*           might be lost when the shift of the RRR is subtracted to obtain */
 761 /*           the eigenvalues of T. However, T is not guaranteed to define its */
 762 /*           eigenvalues to high relative accuracy anyway. */
 763 /*           Set RTOL to the order of the tolerance used in DLASQ2 */
 764 /*           This is an ESTIMATED error, the worst case bound is 4*N*EPS */
 765 /*           which is usually too large and requires unnecessary work to be */
 766 /*           done by bisection when computing the eigenvectors */
 767             rtol = log((doublereal) in) * 4. * eps;
 768             j = ibegin;
 769             i__2 = in - 1;
 770             for (i__ = 1; i__ <= i__2; ++i__) {
 771                 work[(i__ << 1) - 1] = (d__1 = d__[j], abs(d__1));
 772                 work[i__ * 2] = e[j] * e[j] * work[(i__ << 1) - 1];
 773                 ++j;
 774 /* L140: */
 775             }
 776             work[(in << 1) - 1] = (d__1 = d__[iend], abs(d__1));
 777             work[in * 2] = 0.;
 778             dlasq2_(&in, &work[1], &iinfo);
 779             if (iinfo != 0) {
 780 /*              If IINFO = -5 then an index is part of a tight cluster */
 781 /*              and should be changed. The index is in IWORK(1) and the */
 782 /*              gap is in WORK(N+1) */
 783                 *info = -5;
 784                 return 0;
 785             } else {
 786 /*              Test that all eigenvalues are positive as expected */
 787                 i__2 = in;
 788                 for (i__ = 1; i__ <= i__2; ++i__) {
 789                     if (work[i__] < 0.) {
 790                         *info = -6;
 791                         return 0;
 792                     }
 793 /* L149: */
 794                 }
 795             }
 796             if (sgndef > 0.) {
 797                 i__2 = indu;
 798                 for (i__ = indl; i__ <= i__2; ++i__) {
 799                     ++(*m);
 800                     w[*m] = work[in - i__ + 1];
 801                     iblock[*m] = jblk;
 802                     indexw[*m] = i__;
 803 /* L150: */
 804                 }
 805             } else {
 806                 i__2 = indu;
 807                 for (i__ = indl; i__ <= i__2; ++i__) {
 808                     ++(*m);
 809                     w[*m] = -work[i__];
 810                     iblock[*m] = jblk;
 811                     indexw[*m] = i__;
 812 /* L160: */
 813                 }
 814             }
 815             i__2 = *m;
 816             for (i__ = *m - mb + 1; i__ <= i__2; ++i__) {
 817 /*              the value of RTOL below should be the tolerance in DLASQ2 */
 818                 werr[i__] = rtol * (d__1 = w[i__], abs(d__1));
 819 /* L165: */
 820             }
 821             i__2 = *m - 1;
 822             for (i__ = *m - mb + 1; i__ <= i__2; ++i__) {
 823 /*              compute the right gap between the intervals */
 824 /* Computing MAX */
 825                 d__1 = 0., d__2 = w[i__ + 1] - werr[i__ + 1] - (w[i__] + werr[
 826                         i__]);
 827                 wgap[i__] = max(d__1,d__2);
 828 /* L166: */
 829             }
 830 /* Computing MAX */
 831             d__1 = 0., d__2 = *vu - sigma - (w[*m] + werr[*m]);
 832             wgap[*m] = max(d__1,d__2);
 833         }
 834 /*        proceed with next block */
 835         ibegin = iend + 1;
 836         wbegin = wend + 1;
 837 L170:
 838         ;
 839     }
 840
 841     return 0;
 842
 843 /*     end of DLARRE */
 844
 845 } /* dlarre_ */