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