vcs.maemo.org Git - opencv/blob - 3rdparty/lapack/slarrd.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7 static integer c__3 = 3;
   8 static integer c__2 = 2;
   9 static integer c__0 = 0;
  10
  11 /* Subroutine */ int slarrd_(char *range, char *order, integer *n, real *vl,
  12         real *vu, integer *il, integer *iu, real *gers, real *reltol, real *
  13         d__, real *e, real *e2, real *pivmin, integer *nsplit, integer *
  14         isplit, integer *m, real *w, real *werr, real *wl, real *wu, integer *
  15         iblock, integer *indexw, real *work, integer *iwork, integer *info)
  16 {
  17     /* System generated locals */
  18     integer i__1, i__2, i__3;
  19     real r__1, r__2;
  20
  21     /* Builtin functions */
  22     double log(doublereal);
  23
  24     /* Local variables */
  25     integer i__, j, ib, ie, je, nb;
  26     real gl;
  27     integer im, in;
  28     real gu;
  29     integer iw, jee;
  30     real eps;
  31     integer nwl;
  32     real wlu, wul;
  33     integer nwu;
  34     real tmp1, tmp2;
  35     integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc;
  36     extern logical lsame_(char *, char *);
  37     integer iinfo;
  38     real atoli;
  39     integer iwoff, itmax;
  40     real wkill, rtoli, uflow, tnorm;
  41     integer ibegin, irange, idiscl;
  42     extern doublereal slamch_(char *);
  43     integer idumma[1];
  44     real spdiam;
  45     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  46             integer *, integer *);
  47     integer idiscu;
  48     extern /* Subroutine */ int slaebz_(integer *, integer *, integer *,
  49             integer *, integer *, integer *, real *, real *, real *, real *,
  50             real *, real *, integer *, real *, real *, integer *, integer *,
  51             real *, integer *, integer *);
  52     logical ncnvrg, toofew;
  53
  54
  55 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  56 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  57 /*     November 2006 */
  58
  59 /*     .. Scalar Arguments .. */
  60 /*     .. */
  61 /*     .. Array Arguments .. */
  62 /*     .. */
  63
  64 /*  Purpose */
  65 /*  ======= */
  66
  67 /*  SLARRD computes the eigenvalues of a symmetric tridiagonal */
  68 /*  matrix T to suitable accuracy. This is an auxiliary code to be */
  69 /*  called from SSTEMR. */
  70 /*  The user may ask for all eigenvalues, all eigenvalues */
  71 /*  in the half-open interval (VL, VU], or the IL-th through IU-th */
  72 /*  eigenvalues. */
  73
  74 /*  To avoid overflow, the matrix must be scaled so that its */
  75 /*  largest element is no greater than overflow**(1/2) * */
  76 /*  underflow**(1/4) in absolute value, and for greatest */
  77 /*  accuracy, it should not be much smaller than that. */
  78
  79 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
  80 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
  81 /*  University, July 21, 1966. */
  82
  83 /*  Arguments */
  84 /*  ========= */
  85
  86 /*  RANGE   (input) CHARACTER */
  87 /*          = 'A': ("All")   all eigenvalues will be found. */
  88 /*          = 'V': ("Value") all eigenvalues in the half-open interval */
  89 /*                           (VL, VU] will be found. */
  90 /*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
  91 /*                           entire matrix) will be found. */
  92
  93 /*  ORDER   (input) CHARACTER */
  94 /*          = 'B': ("By Block") the eigenvalues will be grouped by */
  95 /*                              split-off block (see IBLOCK, ISPLIT) and */
  96 /*                              ordered from smallest to largest within */
  97 /*                              the block. */
  98 /*          = 'E': ("Entire matrix") */
  99 /*                              the eigenvalues for the entire matrix */
 100 /*                              will be ordered from smallest to */
 101 /*                              largest. */
 102
 103 /*  N       (input) INTEGER */
 104 /*          The order of the tridiagonal matrix T.  N >= 0. */
 105
 106 /*  VL      (input) REAL */
 107 /*  VU      (input) REAL */
 108 /*          If RANGE='V', the lower and upper bounds of the interval to */
 109 /*          be searched for eigenvalues.  Eigenvalues less than or equal */
 110 /*          to VL, or greater than VU, will not be returned.  VL < VU. */
 111 /*          Not referenced if RANGE = 'A' or 'I'. */
 112
 113 /*  IL      (input) INTEGER */
 114 /*  IU      (input) INTEGER */
 115 /*          If RANGE='I', the indices (in ascending order) of the */
 116 /*          smallest and largest eigenvalues to be returned. */
 117 /*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
 118 /*          Not referenced if RANGE = 'A' or 'V'. */
 119
 120 /*  GERS    (input) REAL             array, dimension (2*N) */
 121 /*          The N Gerschgorin intervals (the i-th Gerschgorin interval */
 122 /*          is (GERS(2*i-1), GERS(2*i)). */
 123
 124 /*  RELTOL  (input) REAL */
 125 /*          The minimum relative width of an interval.  When an interval */
 126 /*          is narrower than RELTOL times the larger (in */
 127 /*          magnitude) endpoint, then it is considered to be */
 128 /*          sufficiently small, i.e., converged.  Note: this should */
 129 /*          always be at least radix*machine epsilon. */
 130
 131 /*  D       (input) REAL             array, dimension (N) */
 132 /*          The n diagonal elements of the tridiagonal matrix T. */
 133
 134 /*  E       (input) REAL             array, dimension (N-1) */
 135 /*          The (n-1) off-diagonal elements of the tridiagonal matrix T. */
 136
 137 /*  E2      (input) REAL             array, dimension (N-1) */
 138 /*          The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
 139
 140 /*  PIVMIN  (input) REAL */
 141 /*          The minimum pivot allowed in the Sturm sequence for T. */
 142
 143 /*  NSPLIT  (input) INTEGER */
 144 /*          The number of diagonal blocks in the matrix T. */
 145 /*          1 <= NSPLIT <= N. */
 146
 147 /*  ISPLIT  (input) INTEGER array, dimension (N) */
 148 /*          The splitting points, at which T breaks up into submatrices. */
 149 /*          The first submatrix consists of rows/columns 1 to ISPLIT(1), */
 150 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
 151 /*          etc., and the NSPLIT-th consists of rows/columns */
 152 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
 153 /*          (Only the first NSPLIT elements will actually be used, but */
 154 /*          since the user cannot know a priori what value NSPLIT will */
 155 /*          have, N words must be reserved for ISPLIT.) */
 156
 157 /*  M       (output) INTEGER */
 158 /*          The actual number of eigenvalues found. 0 <= M <= N. */
 159 /*          (See also the description of INFO=2,3.) */
 160
 161 /*  W       (output) REAL             array, dimension (N) */
 162 /*          On exit, the first M elements of W will contain the */
 163 /*          eigenvalue approximations. SLARRD computes an interval */
 164 /*          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */
 165 /*          approximation is given as the interval midpoint */
 166 /*          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */
 167 /*          WERR(j) = abs( a_j - b_j)/2 */
 168
 169 /*  WERR    (output) REAL             array, dimension (N) */
 170 /*          The error bound on the corresponding eigenvalue approximation */
 171 /*          in W. */
 172
 173 /*  WL      (output) REAL */
 174 /*  WU      (output) REAL */
 175 /*          The interval (WL, WU] contains all the wanted eigenvalues. */
 176 /*          If RANGE='V', then WL=VL and WU=VU. */
 177 /*          If RANGE='A', then WL and WU are the global Gerschgorin bounds */
 178 /*                        on the spectrum. */
 179 /*          If RANGE='I', then WL and WU are computed by SLAEBZ from the */
 180 /*                        index range specified. */
 181
 182 /*  IBLOCK  (output) INTEGER array, dimension (N) */
 183 /*          At each row/column j where E(j) is zero or small, the */
 184 /*          matrix T is considered to split into a block diagonal */
 185 /*          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which */
 186 /*          block (from 1 to the number of blocks) the eigenvalue W(i) */
 187 /*          belongs.  (SLARRD may use the remaining N-M elements as */
 188 /*          workspace.) */
 189
 190 /*  INDEXW  (output) INTEGER array, dimension (N) */
 191 /*          The indices of the eigenvalues within each block (submatrix); */
 192 /*          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */
 193 /*          i-th eigenvalue W(i) is the j-th eigenvalue in block k. */
 194
 195 /*  WORK    (workspace) REAL             array, dimension (4*N) */
 196
 197 /*  IWORK   (workspace) INTEGER array, dimension (3*N) */
 198
 199 /*  INFO    (output) INTEGER */
 200 /*          = 0:  successful exit */
 201 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 202 /*          > 0:  some or all of the eigenvalues failed to converge or */
 203 /*                were not computed: */
 204 /*                =1 or 3: Bisection failed to converge for some */
 205 /*                        eigenvalues; these eigenvalues are flagged by a */
 206 /*                        negative block number.  The effect is that the */
 207 /*                        eigenvalues may not be as accurate as the */
 208 /*                        absolute and relative tolerances.  This is */
 209 /*                        generally caused by unexpectedly inaccurate */
 210 /*                        arithmetic. */
 211 /*                =2 or 3: RANGE='I' only: Not all of the eigenvalues */
 212 /*                        IL:IU were found. */
 213 /*                        Effect: M < IU+1-IL */
 214 /*                        Cause:  non-monotonic arithmetic, causing the */
 215 /*                                Sturm sequence to be non-monotonic. */
 216 /*                        Cure:   recalculate, using RANGE='A', and pick */
 217 /*                                out eigenvalues IL:IU.  In some cases, */
 218 /*                                increasing the PARAMETER "FUDGE" may */
 219 /*                                make things work. */
 220 /*                = 4:    RANGE='I', and the Gershgorin interval */
 221 /*                        initially used was too small.  No eigenvalues */
 222 /*                        were computed. */
 223 /*                        Probable cause: your machine has sloppy */
 224 /*                                        floating-point arithmetic. */
 225 /*                        Cure: Increase the PARAMETER "FUDGE", */
 226 /*                              recompile, and try again. */
 227
 228 /*  Internal Parameters */
 229 /*  =================== */
 230
 231 /*  FUDGE   REAL            , default = 2 */
 232 /*          A "fudge factor" to widen the Gershgorin intervals.  Ideally, */
 233 /*          a value of 1 should work, but on machines with sloppy */
 234 /*          arithmetic, this needs to be larger.  The default for */
 235 /*          publicly released versions should be large enough to handle */
 236 /*          the worst machine around.  Note that this has no effect */
 237 /*          on accuracy of the solution. */
 238
 239 /*  Based on contributions by */
 240 /*     W. Kahan, University of California, Berkeley, USA */
 241 /*     Beresford Parlett, University of California, Berkeley, USA */
 242 /*     Jim Demmel, University of California, Berkeley, USA */
 243 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 244 /*     Osni Marques, LBNL/NERSC, USA */
 245 /*     Christof Voemel, University of California, Berkeley, USA */
 246
 247 /*  ===================================================================== */
 248
 249 /*     .. Parameters .. */
 250 /*     .. */
 251 /*     .. Local Scalars .. */
 252 /*     .. */
 253 /*     .. Local Arrays .. */
 254 /*     .. */
 255 /*     .. External Functions .. */
 256 /*     .. */
 257 /*     .. External Subroutines .. */
 258 /*     .. */
 259 /*     .. Intrinsic Functions .. */
 260 /*     .. */
 261 /*     .. Executable Statements .. */
 262
 263     /* Parameter adjustments */
 264     --iwork;
 265     --work;
 266     --indexw;
 267     --iblock;
 268     --werr;
 269     --w;
 270     --isplit;
 271     --e2;
 272     --e;
 273     --d__;
 274     --gers;
 275
 276     /* Function Body */
 277     *info = 0;
 278
 279 /*     Decode RANGE */
 280
 281     if (lsame_(range, "A")) {
 282         irange = 1;
 283     } else if (lsame_(range, "V")) {
 284         irange = 2;
 285     } else if (lsame_(range, "I")) {
 286         irange = 3;
 287     } else {
 288         irange = 0;
 289     }
 290
 291 /*     Check for Errors */
 292
 293     if (irange <= 0) {
 294         *info = -1;
 295     } else if (! (lsame_(order, "B") || lsame_(order,
 296             "E"))) {
 297         *info = -2;
 298     } else if (*n < 0) {
 299         *info = -3;
 300     } else if (irange == 2) {
 301         if (*vl >= *vu) {
 302             *info = -5;
 303         }
 304     } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
 305         *info = -6;
 306     } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
 307         *info = -7;
 308     }
 309
 310     if (*info != 0) {
 311         return 0;
 312     }
 313 /*     Initialize error flags */
 314     *info = 0;
 315     ncnvrg = FALSE_;
 316     toofew = FALSE_;
 317 /*     Quick return if possible */
 318     *m = 0;
 319     if (*n == 0) {
 320         return 0;
 321     }
 322 /*     Simplification: */
 323     if (irange == 3 && *il == 1 && *iu == *n) {
 324         irange = 1;
 325     }
 326 /*     Get machine constants */
 327     eps = slamch_("P");
 328     uflow = slamch_("U");
 329 /*     Special Case when N=1 */
 330 /*     Treat case of 1x1 matrix for quick return */
 331     if (*n == 1) {
 332         if (irange == 1 || irange == 2 && d__[1] > *vl && d__[1] <= *vu ||
 333                 irange == 3 && *il == 1 && *iu == 1) {
 334             *m = 1;
 335             w[1] = d__[1];
 336 /*           The computation error of the eigenvalue is zero */
 337             werr[1] = 0.f;
 338             iblock[1] = 1;
 339             indexw[1] = 1;
 340         }
 341         return 0;
 342     }
 343 /*     NB is the minimum vector length for vector bisection, or 0 */
 344 /*     if only scalar is to be done. */
 345     nb = ilaenv_(&c__1, "SSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
 346     if (nb <= 1) {
 347         nb = 0;
 348     }
 349 /*     Find global spectral radius */
 350     gl = d__[1];
 351     gu = d__[1];
 352     i__1 = *n;
 353     for (i__ = 1; i__ <= i__1; ++i__) {
 354 /* Computing MIN */
 355         r__1 = gl, r__2 = gers[(i__ << 1) - 1];
 356         gl = dmin(r__1,r__2);
 357 /* Computing MAX */
 358         r__1 = gu, r__2 = gers[i__ * 2];
 359         gu = dmax(r__1,r__2);
 360 /* L5: */
 361     }
 362 /*     Compute global Gerschgorin bounds and spectral diameter */
 363 /* Computing MAX */
 364     r__1 = dabs(gl), r__2 = dabs(gu);
 365     tnorm = dmax(r__1,r__2);
 366     gl = gl - tnorm * 2.f * eps * *n - *pivmin * 4.f;
 367     gu = gu + tnorm * 2.f * eps * *n + *pivmin * 4.f;
 368     spdiam = gu - gl;
 369 /*     Input arguments for SLAEBZ: */
 370 /*     The relative tolerance.  An interval (a,b] lies within */
 371 /*     "relative tolerance" if  b-a < RELTOL*max(|a|,|b|), */
 372     rtoli = *reltol;
 373 /*     Set the absolute tolerance for interval convergence to zero to force */
 374 /*     interval convergence based on relative size of the interval. */
 375 /*     This is dangerous because intervals might not converge when RELTOL is */
 376 /*     small. But at least a very small number should be selected so that for */
 377 /*     strongly graded matrices, the code can get relatively accurate */
 378 /*     eigenvalues. */
 379     atoli = uflow * 4.f + *pivmin * 4.f;
 380     if (irange == 3) {
 381 /*        RANGE='I': Compute an interval containing eigenvalues */
 382 /*        IL through IU. The initial interval [GL,GU] from the global */
 383 /*        Gerschgorin bounds GL and GU is refined by SLAEBZ. */
 384         itmax = (integer) ((log(tnorm + *pivmin) - log(*pivmin)) / log(2.f))
 385                 + 2;
 386         work[*n + 1] = gl;
 387         work[*n + 2] = gl;
 388         work[*n + 3] = gu;
 389         work[*n + 4] = gu;
 390         work[*n + 5] = gl;
 391         work[*n + 6] = gu;
 392         iwork[1] = -1;
 393         iwork[2] = -1;
 394         iwork[3] = *n + 1;
 395         iwork[4] = *n + 1;
 396         iwork[5] = *il - 1;
 397         iwork[6] = *iu;
 398
 399         slaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, &
 400                 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5]
 401 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
 402         if (iinfo != 0) {
 403             *info = iinfo;
 404             return 0;
 405         }
 406 /*        On exit, output intervals may not be ordered by ascending negcount */
 407         if (iwork[6] == *iu) {
 408             *wl = work[*n + 1];
 409             wlu = work[*n + 3];
 410             nwl = iwork[1];
 411             *wu = work[*n + 4];
 412             wul = work[*n + 2];
 413             nwu = iwork[4];
 414         } else {
 415             *wl = work[*n + 2];
 416             wlu = work[*n + 4];
 417             nwl = iwork[2];
 418             *wu = work[*n + 3];
 419             wul = work[*n + 1];
 420             nwu = iwork[3];
 421         }
 422 /*        On exit, the interval [WL, WLU] contains a value with negcount NWL, */
 423 /*        and [WUL, WU] contains a value with negcount NWU. */
 424         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
 425             *info = 4;
 426             return 0;
 427         }
 428     } else if (irange == 2) {
 429         *wl = *vl;
 430         *wu = *vu;
 431     } else if (irange == 1) {
 432         *wl = gl;
 433         *wu = gu;
 434     }
 435 /*     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */
 436 /*     NWL accumulates the number of eigenvalues .le. WL, */
 437 /*     NWU accumulates the number of eigenvalues .le. WU */
 438     *m = 0;
 439     iend = 0;
 440     *info = 0;
 441     nwl = 0;
 442     nwu = 0;
 443
 444     i__1 = *nsplit;
 445     for (jblk = 1; jblk <= i__1; ++jblk) {
 446         ioff = iend;
 447         ibegin = ioff + 1;
 448         iend = isplit[jblk];
 449         in = iend - ioff;
 450
 451         if (in == 1) {
 452 /*           1x1 block */
 453             if (*wl >= d__[ibegin] - *pivmin) {
 454                 ++nwl;
 455             }
 456             if (*wu >= d__[ibegin] - *pivmin) {
 457                 ++nwu;
 458             }
 459             if (irange == 1 || *wl < d__[ibegin] - *pivmin && *wu >= d__[
 460                     ibegin] - *pivmin) {
 461                 ++(*m);
 462                 w[*m] = d__[ibegin];
 463                 werr[*m] = 0.f;
 464 /*              The gap for a single block doesn't matter for the later */
 465 /*              algorithm and is assigned an arbitrary large value */
 466                 iblock[*m] = jblk;
 467                 indexw[*m] = 1;
 468             }
 469 /*        Disabled 2x2 case because of a failure on the following matrix */
 470 /*        RANGE = 'I', IL = IU = 4 */
 471 /*          Original Tridiagonal, d = [ */
 472 /*           -0.150102010615740E+00 */
 473 /*           -0.849897989384260E+00 */
 474 /*           -0.128208148052635E-15 */
 475 /*            0.128257718286320E-15 */
 476 /*          ]; */
 477 /*          e = [ */
 478 /*           -0.357171383266986E+00 */
 479 /*           -0.180411241501588E-15 */
 480 /*           -0.175152352710251E-15 */
 481 /*          ]; */
 482
 483 /*         ELSE IF( IN.EQ.2 ) THEN */
 484 /* *           2x2 block */
 485 /*            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */
 486 /*            TMP1 = HALF*(D(IBEGIN)+D(IEND)) */
 487 /*            L1 = TMP1 - DISC */
 488 /*            IF( WL.GE. L1-PIVMIN ) */
 489 /*     $         NWL = NWL + 1 */
 490 /*            IF( WU.GE. L1-PIVMIN ) */
 491 /*     $         NWU = NWU + 1 */
 492 /*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */
 493 /*     $          L1-PIVMIN ) ) THEN */
 494 /*               M = M + 1 */
 495 /*               W( M ) = L1 */
 496 /* *              The uncertainty of eigenvalues of a 2x2 matrix is very small */
 497 /*               WERR( M ) = EPS * ABS( W( M ) ) * TWO */
 498 /*               IBLOCK( M ) = JBLK */
 499 /*               INDEXW( M ) = 1 */
 500 /*            ENDIF */
 501 /*            L2 = TMP1 + DISC */
 502 /*            IF( WL.GE. L2-PIVMIN ) */
 503 /*     $         NWL = NWL + 1 */
 504 /*            IF( WU.GE. L2-PIVMIN ) */
 505 /*     $         NWU = NWU + 1 */
 506 /*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */
 507 /*     $          L2-PIVMIN ) ) THEN */
 508 /*               M = M + 1 */
 509 /*               W( M ) = L2 */
 510 /* *              The uncertainty of eigenvalues of a 2x2 matrix is very small */
 511 /*               WERR( M ) = EPS * ABS( W( M ) ) * TWO */
 512 /*               IBLOCK( M ) = JBLK */
 513 /*               INDEXW( M ) = 2 */
 514 /*            ENDIF */
 515         } else {
 516 /*           General Case - block of size IN >= 2 */
 517 /*           Compute local Gerschgorin interval and use it as the initial */
 518 /*           interval for SLAEBZ */
 519             gu = d__[ibegin];
 520             gl = d__[ibegin];
 521             tmp1 = 0.f;
 522             i__2 = iend;
 523             for (j = ibegin; j <= i__2; ++j) {
 524 /* Computing MIN */
 525                 r__1 = gl, r__2 = gers[(j << 1) - 1];
 526                 gl = dmin(r__1,r__2);
 527 /* Computing MAX */
 528                 r__1 = gu, r__2 = gers[j * 2];
 529                 gu = dmax(r__1,r__2);
 530 /* L40: */
 531             }
 532             spdiam = gu - gl;
 533             gl = gl - spdiam * 2.f * eps * in - *pivmin * 2.f;
 534             gu = gu + spdiam * 2.f * eps * in + *pivmin * 2.f;
 535
 536             if (irange > 1) {
 537                 if (gu < *wl) {
 538 /*                 the local block contains none of the wanted eigenvalues */
 539                     nwl += in;
 540                     nwu += in;
 541                     goto L70;
 542                 }
 543 /*              refine search interval if possible, only range (WL,WU] matters */
 544                 gl = dmax(gl,*wl);
 545                 gu = dmin(gu,*wu);
 546                 if (gl >= gu) {
 547                     goto L70;
 548                 }
 549             }
 550 /*           Find negcount of initial interval boundaries GL and GU */
 551             work[*n + 1] = gl;
 552             work[*n + in + 1] = gu;
 553             slaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli,
 554                     pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
 555                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
 556                     w[*m + 1], &iblock[*m + 1], &iinfo);
 557             if (iinfo != 0) {
 558                 *info = iinfo;
 559                 return 0;
 560             }
 561
 562             nwl += iwork[1];
 563             nwu += iwork[in + 1];
 564             iwoff = *m - iwork[1];
 565 /*           Compute Eigenvalues */
 566             itmax = (integer) ((log(gu - gl + *pivmin) - log(*pivmin)) / log(
 567                     2.f)) + 2;
 568             slaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli,
 569                     pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
 570                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
 571                      &w[*m + 1], &iblock[*m + 1], &iinfo);
 572             if (iinfo != 0) {
 573                 *info = iinfo;
 574                 return 0;
 575             }
 576
 577 /*           Copy eigenvalues into W and IBLOCK */
 578 /*           Use -JBLK for block number for unconverged eigenvalues. */
 579 /*           Loop over the number of output intervals from SLAEBZ */
 580             i__2 = iout;
 581             for (j = 1; j <= i__2; ++j) {
 582 /*              eigenvalue approximation is middle point of interval */
 583                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5f;
 584 /*              semi length of error interval */
 585                 tmp2 = (r__1 = work[j + *n] - work[j + in + *n], dabs(r__1)) *
 586                          .5f;
 587                 if (j > iout - iinfo) {
 588 /*                 Flag non-convergence. */
 589                     ncnvrg = TRUE_;
 590                     ib = -jblk;
 591                 } else {
 592                     ib = jblk;
 593                 }
 594                 i__3 = iwork[j + in] + iwoff;
 595                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
 596                     w[je] = tmp1;
 597                     werr[je] = tmp2;
 598                     indexw[je] = je - iwoff;
 599                     iblock[je] = ib;
 600 /* L50: */
 601                 }
 602 /* L60: */
 603             }
 604
 605             *m += im;
 606         }
 607 L70:
 608         ;
 609     }
 610 /*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
 611 /*     If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
 612     if (irange == 3) {
 613         idiscl = *il - 1 - nwl;
 614         idiscu = nwu - *iu;
 615
 616         if (idiscl > 0) {
 617             im = 0;
 618             i__1 = *m;
 619             for (je = 1; je <= i__1; ++je) {
 620 /*              Remove some of the smallest eigenvalues from the left so that */
 621 /*              at the end IDISCL =0. Move all eigenvalues up to the left. */
 622                 if (w[je] <= wlu && idiscl > 0) {
 623                     --idiscl;
 624                 } else {
 625                     ++im;
 626                     w[im] = w[je];
 627                     werr[im] = werr[je];
 628                     indexw[im] = indexw[je];
 629                     iblock[im] = iblock[je];
 630                 }
 631 /* L80: */
 632             }
 633             *m = im;
 634         }
 635         if (idiscu > 0) {
 636 /*           Remove some of the largest eigenvalues from the right so that */
 637 /*           at the end IDISCU =0. Move all eigenvalues up to the left. */
 638             im = *m + 1;
 639             for (je = *m; je >= 1; --je) {
 640                 if (w[je] >= wul && idiscu > 0) {
 641                     --idiscu;
 642                 } else {
 643                     --im;
 644                     w[im] = w[je];
 645                     werr[im] = werr[je];
 646                     indexw[im] = indexw[je];
 647                     iblock[im] = iblock[je];
 648                 }
 649 /* L81: */
 650             }
 651             jee = 0;
 652             i__1 = *m;
 653             for (je = im; je <= i__1; ++je) {
 654                 ++jee;
 655                 w[jee] = w[je];
 656                 werr[jee] = werr[je];
 657                 indexw[jee] = indexw[je];
 658                 iblock[jee] = iblock[je];
 659 /* L82: */
 660             }
 661             *m = *m - im + 1;
 662         }
 663         if (idiscl > 0 || idiscu > 0) {
 664 /*           Code to deal with effects of bad arithmetic. (If N(w) is */
 665 /*           monotone non-decreasing, this should never happen.) */
 666 /*           Some low eigenvalues to be discarded are not in (WL,WLU], */
 667 /*           or high eigenvalues to be discarded are not in (WUL,WU] */
 668 /*           so just kill off the smallest IDISCL/largest IDISCU */
 669 /*           eigenvalues, by marking the corresponding IBLOCK = 0 */
 670             if (idiscl > 0) {
 671                 wkill = *wu;
 672                 i__1 = idiscl;
 673                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 674                     iw = 0;
 675                     i__2 = *m;
 676                     for (je = 1; je <= i__2; ++je) {
 677                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
 678                             iw = je;
 679                             wkill = w[je];
 680                         }
 681 /* L90: */
 682                     }
 683                     iblock[iw] = 0;
 684 /* L100: */
 685                 }
 686             }
 687             if (idiscu > 0) {
 688                 wkill = *wl;
 689                 i__1 = idiscu;
 690                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 691                     iw = 0;
 692                     i__2 = *m;
 693                     for (je = 1; je <= i__2; ++je) {
 694                         if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) {
 695                             iw = je;
 696                             wkill = w[je];
 697                         }
 698 /* L110: */
 699                     }
 700                     iblock[iw] = 0;
 701 /* L120: */
 702                 }
 703             }
 704 /*           Now erase all eigenvalues with IBLOCK set to zero */
 705             im = 0;
 706             i__1 = *m;
 707             for (je = 1; je <= i__1; ++je) {
 708                 if (iblock[je] != 0) {
 709                     ++im;
 710                     w[im] = w[je];
 711                     werr[im] = werr[je];
 712                     indexw[im] = indexw[je];
 713                     iblock[im] = iblock[je];
 714                 }
 715 /* L130: */
 716             }
 717             *m = im;
 718         }
 719         if (idiscl < 0 || idiscu < 0) {
 720             toofew = TRUE_;
 721         }
 722     }
 723
 724     if (irange == 1 && *m != *n || irange == 3 && *m != *iu - *il + 1) {
 725         toofew = TRUE_;
 726     }
 727 /*     If ORDER='B', do nothing the eigenvalues are already sorted by */
 728 /*        block. */
 729 /*     If ORDER='E', sort the eigenvalues from smallest to largest */
 730     if (lsame_(order, "E") && *nsplit > 1) {
 731         i__1 = *m - 1;
 732         for (je = 1; je <= i__1; ++je) {
 733             ie = 0;
 734             tmp1 = w[je];
 735             i__2 = *m;
 736             for (j = je + 1; j <= i__2; ++j) {
 737                 if (w[j] < tmp1) {
 738                     ie = j;
 739                     tmp1 = w[j];
 740                 }
 741 /* L140: */
 742             }
 743             if (ie != 0) {
 744                 tmp2 = werr[ie];
 745                 itmp1 = iblock[ie];
 746                 itmp2 = indexw[ie];
 747                 w[ie] = w[je];
 748                 werr[ie] = werr[je];
 749                 iblock[ie] = iblock[je];
 750                 indexw[ie] = indexw[je];
 751                 w[je] = tmp1;
 752                 werr[je] = tmp2;
 753                 iblock[je] = itmp1;
 754                 indexw[je] = itmp2;
 755             }
 756 /* L150: */
 757         }
 758     }
 759
 760     *info = 0;
 761     if (ncnvrg) {
 762         ++(*info);
 763     }
 764     if (toofew) {
 765         *info += 2;
 766     }
 767     return 0;
 768
 769 /*     End of SLARRD */
 770
 771 } /* slarrd_ */