git.maemo.org Git - opencv/blob - 3rdparty/lapack/dstebz.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 dstebz_(char *range, char *order, integer *n, doublereal
  12         *vl, doublereal *vu, integer *il, integer *iu, doublereal *abstol,
  13         doublereal *d__, doublereal *e, integer *m, integer *nsplit,
  14         doublereal *w, integer *iblock, integer *isplit, doublereal *work,
  15         integer *iwork, integer *info)
  16 {
  17     /* System generated locals */
  18     integer i__1, i__2, i__3;
  19     doublereal d__1, d__2, d__3, d__4, d__5;
  20
  21     /* Builtin functions */
  22     double sqrt(doublereal), log(doublereal);
  23
  24     /* Local variables */
  25     integer j, ib, jb, ie, je, nb;
  26     doublereal gl;
  27     integer im, in;
  28     doublereal gu;
  29     integer iw;
  30     doublereal wl, wu;
  31     integer nwl;
  32     doublereal ulp, wlu, wul;
  33     integer nwu;
  34     doublereal tmp1, tmp2;
  35     integer iend, ioff, iout, itmp1, jdisc;
  36     extern logical lsame_(char *, char *);
  37     integer iinfo;
  38     doublereal atoli;
  39     integer iwoff;
  40     doublereal bnorm;
  41     integer itmax;
  42     doublereal wkill, rtoli, tnorm;
  43     extern doublereal dlamch_(char *);
  44     integer ibegin;
  45     extern /* Subroutine */ int dlaebz_(integer *, integer *, integer *,
  46             integer *, integer *, integer *, doublereal *, doublereal *,
  47             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  48              doublereal *, doublereal *, integer *, integer *, doublereal *,
  49             integer *, integer *);
  50     integer irange, idiscl;
  51     doublereal safemn;
  52     integer idumma[1];
  53     extern /* Subroutine */ int xerbla_(char *, integer *);
  54     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  55             integer *, integer *);
  56     integer idiscu, iorder;
  57     logical ncnvrg;
  58     doublereal pivmin;
  59     logical toofew;
  60
  61
  62 /*  -- LAPACK routine (version 3.1) -- */
  63 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  64 /*     November 2006 */
  65 /*     8-18-00:  Increase FUDGE factor for T3E (eca) */
  66
  67 /*     .. Scalar Arguments .. */
  68 /*     .. */
  69 /*     .. Array Arguments .. */
  70 /*     .. */
  71
  72 /*  Purpose */
  73 /*  ======= */
  74
  75 /*  DSTEBZ computes the eigenvalues of a symmetric tridiagonal */
  76 /*  matrix T.  The user may ask for all eigenvalues, all eigenvalues */
  77 /*  in the half-open interval (VL, VU], or the IL-th through IU-th */
  78 /*  eigenvalues. */
  79
  80 /*  To avoid overflow, the matrix must be scaled so that its */
  81 /*  largest element is no greater than overflow**(1/2) * */
  82 /*  underflow**(1/4) in absolute value, and for greatest */
  83 /*  accuracy, it should not be much smaller than that. */
  84
  85 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
  86 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
  87 /*  University, July 21, 1966. */
  88
  89 /*  Arguments */
  90 /*  ========= */
  91
  92 /*  RANGE   (input) CHARACTER*1 */
  93 /*          = 'A': ("All")   all eigenvalues will be found. */
  94 /*          = 'V': ("Value") all eigenvalues in the half-open interval */
  95 /*                           (VL, VU] will be found. */
  96 /*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
  97 /*                           entire matrix) will be found. */
  98
  99 /*  ORDER   (input) CHARACTER*1 */
 100 /*          = 'B': ("By Block") the eigenvalues will be grouped by */
 101 /*                              split-off block (see IBLOCK, ISPLIT) and */
 102 /*                              ordered from smallest to largest within */
 103 /*                              the block. */
 104 /*          = 'E': ("Entire matrix") */
 105 /*                              the eigenvalues for the entire matrix */
 106 /*                              will be ordered from smallest to */
 107 /*                              largest. */
 108
 109 /*  N       (input) INTEGER */
 110 /*          The order of the tridiagonal matrix T.  N >= 0. */
 111
 112 /*  VL      (input) DOUBLE PRECISION */
 113 /*  VU      (input) DOUBLE PRECISION */
 114 /*          If RANGE='V', the lower and upper bounds of the interval to */
 115 /*          be searched for eigenvalues.  Eigenvalues less than or equal */
 116 /*          to VL, or greater than VU, will not be returned.  VL < VU. */
 117 /*          Not referenced if RANGE = 'A' or 'I'. */
 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, if N > 0; IL = 1 and IU = 0 if N = 0. */
 124 /*          Not referenced if RANGE = 'A' or 'V'. */
 125
 126 /*  ABSTOL  (input) DOUBLE PRECISION */
 127 /*          The absolute tolerance for the eigenvalues.  An eigenvalue */
 128 /*          (or cluster) is considered to be located if it has been */
 129 /*          determined to lie in an interval whose width is ABSTOL or */
 130 /*          less.  If ABSTOL is less than or equal to zero, then ULP*|T| */
 131 /*          will be used, where |T| means the 1-norm of T. */
 132
 133 /*          Eigenvalues will be computed most accurately when ABSTOL is */
 134 /*          set to twice the underflow threshold 2*DLAMCH('S'), not zero. */
 135
 136 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
 137 /*          The n diagonal elements of the tridiagonal matrix T. */
 138
 139 /*  E       (input) DOUBLE PRECISION array, dimension (N-1) */
 140 /*          The (n-1) off-diagonal elements of the tridiagonal matrix T. */
 141
 142 /*  M       (output) INTEGER */
 143 /*          The actual number of eigenvalues found. 0 <= M <= N. */
 144 /*          (See also the description of INFO=2,3.) */
 145
 146 /*  NSPLIT  (output) INTEGER */
 147 /*          The number of diagonal blocks in the matrix T. */
 148 /*          1 <= NSPLIT <= N. */
 149
 150 /*  W       (output) DOUBLE PRECISION array, dimension (N) */
 151 /*          On exit, the first M elements of W will contain the */
 152 /*          eigenvalues.  (DSTEBZ may use the remaining N-M elements as */
 153 /*          workspace.) */
 154
 155 /*  IBLOCK  (output) INTEGER array, dimension (N) */
 156 /*          At each row/column j where E(j) is zero or small, the */
 157 /*          matrix T is considered to split into a block diagonal */
 158 /*          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which */
 159 /*          block (from 1 to the number of blocks) the eigenvalue W(i) */
 160 /*          belongs.  (DSTEBZ may use the remaining N-M elements as */
 161 /*          workspace.) */
 162
 163 /*  ISPLIT  (output) INTEGER array, dimension (N) */
 164 /*          The splitting points, at which T breaks up into submatrices. */
 165 /*          The first submatrix consists of rows/columns 1 to ISPLIT(1), */
 166 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
 167 /*          etc., and the NSPLIT-th consists of rows/columns */
 168 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
 169 /*          (Only the first NSPLIT elements will actually be used, but */
 170 /*          since the user cannot know a priori what value NSPLIT will */
 171 /*          have, N words must be reserved for ISPLIT.) */
 172
 173 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N) */
 174
 175 /*  IWORK   (workspace) INTEGER array, dimension (3*N) */
 176
 177 /*  INFO    (output) INTEGER */
 178 /*          = 0:  successful exit */
 179 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 180 /*          > 0:  some or all of the eigenvalues failed to converge or */
 181 /*                were not computed: */
 182 /*                =1 or 3: Bisection failed to converge for some */
 183 /*                        eigenvalues; these eigenvalues are flagged by a */
 184 /*                        negative block number.  The effect is that the */
 185 /*                        eigenvalues may not be as accurate as the */
 186 /*                        absolute and relative tolerances.  This is */
 187 /*                        generally caused by unexpectedly inaccurate */
 188 /*                        arithmetic. */
 189 /*                =2 or 3: RANGE='I' only: Not all of the eigenvalues */
 190 /*                        IL:IU were found. */
 191 /*                        Effect: M < IU+1-IL */
 192 /*                        Cause:  non-monotonic arithmetic, causing the */
 193 /*                                Sturm sequence to be non-monotonic. */
 194 /*                        Cure:   recalculate, using RANGE='A', and pick */
 195 /*                                out eigenvalues IL:IU.  In some cases, */
 196 /*                                increasing the PARAMETER "FUDGE" may */
 197 /*                                make things work. */
 198 /*                = 4:    RANGE='I', and the Gershgorin interval */
 199 /*                        initially used was too small.  No eigenvalues */
 200 /*                        were computed. */
 201 /*                        Probable cause: your machine has sloppy */
 202 /*                                        floating-point arithmetic. */
 203 /*                        Cure: Increase the PARAMETER "FUDGE", */
 204 /*                              recompile, and try again. */
 205
 206 /*  Internal Parameters */
 207 /*  =================== */
 208
 209 /*  RELFAC  DOUBLE PRECISION, default = 2.0e0 */
 210 /*          The relative tolerance.  An interval (a,b] lies within */
 211 /*          "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|), */
 212 /*          where "ulp" is the machine precision (distance from 1 to */
 213 /*          the next larger floating point number.) */
 214
 215 /*  FUDGE   DOUBLE PRECISION, default = 2 */
 216 /*          A "fudge factor" to widen the Gershgorin intervals.  Ideally, */
 217 /*          a value of 1 should work, but on machines with sloppy */
 218 /*          arithmetic, this needs to be larger.  The default for */
 219 /*          publicly released versions should be large enough to handle */
 220 /*          the worst machine around.  Note that this has no effect */
 221 /*          on accuracy of the solution. */
 222
 223 /*  ===================================================================== */
 224
 225 /*     .. Parameters .. */
 226 /*     .. */
 227 /*     .. Local Scalars .. */
 228 /*     .. */
 229 /*     .. Local Arrays .. */
 230 /*     .. */
 231 /*     .. External Functions .. */
 232 /*     .. */
 233 /*     .. External Subroutines .. */
 234 /*     .. */
 235 /*     .. Intrinsic Functions .. */
 236 /*     .. */
 237 /*     .. Executable Statements .. */
 238
 239     /* Parameter adjustments */
 240     --iwork;
 241     --work;
 242     --isplit;
 243     --iblock;
 244     --w;
 245     --e;
 246     --d__;
 247
 248     /* Function Body */
 249     *info = 0;
 250
 251 /*     Decode RANGE */
 252
 253     if (lsame_(range, "A")) {
 254         irange = 1;
 255     } else if (lsame_(range, "V")) {
 256         irange = 2;
 257     } else if (lsame_(range, "I")) {
 258         irange = 3;
 259     } else {
 260         irange = 0;
 261     }
 262
 263 /*     Decode ORDER */
 264
 265     if (lsame_(order, "B")) {
 266         iorder = 2;
 267     } else if (lsame_(order, "E")) {
 268         iorder = 1;
 269     } else {
 270         iorder = 0;
 271     }
 272
 273 /*     Check for Errors */
 274
 275     if (irange <= 0) {
 276         *info = -1;
 277     } else if (iorder <= 0) {
 278         *info = -2;
 279     } else if (*n < 0) {
 280         *info = -3;
 281     } else if (irange == 2) {
 282         if (*vl >= *vu) {
 283             *info = -5;
 284         }
 285     } else if (irange == 3 && (*il < 1 || *il > max(1,*n))) {
 286         *info = -6;
 287     } else if (irange == 3 && (*iu < min(*n,*il) || *iu > *n)) {
 288         *info = -7;
 289     }
 290
 291     if (*info != 0) {
 292         i__1 = -(*info);
 293         xerbla_("DSTEBZ", &i__1);
 294         return 0;
 295     }
 296
 297 /*     Initialize error flags */
 298
 299     *info = 0;
 300     ncnvrg = FALSE_;
 301     toofew = FALSE_;
 302
 303 /*     Quick return if possible */
 304
 305     *m = 0;
 306     if (*n == 0) {
 307         return 0;
 308     }
 309
 310 /*     Simplifications: */
 311
 312     if (irange == 3 && *il == 1 && *iu == *n) {
 313         irange = 1;
 314     }
 315
 316 /*     Get machine constants */
 317 /*     NB is the minimum vector length for vector bisection, or 0 */
 318 /*     if only scalar is to be done. */
 319
 320     safemn = dlamch_("S");
 321     ulp = dlamch_("P");
 322     rtoli = ulp * 2.;
 323     nb = ilaenv_(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1);
 324     if (nb <= 1) {
 325         nb = 0;
 326     }
 327
 328 /*     Special Case when N=1 */
 329
 330     if (*n == 1) {
 331         *nsplit = 1;
 332         isplit[1] = 1;
 333         if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
 334             *m = 0;
 335         } else {
 336             w[1] = d__[1];
 337             iblock[1] = 1;
 338             *m = 1;
 339         }
 340         return 0;
 341     }
 342
 343 /*     Compute Splitting Points */
 344
 345     *nsplit = 1;
 346     work[*n] = 0.;
 347     pivmin = 1.;
 348
 349 /* DIR$ NOVECTOR */
 350     i__1 = *n;
 351     for (j = 2; j <= i__1; ++j) {
 352 /* Computing 2nd power */
 353         d__1 = e[j - 1];
 354         tmp1 = d__1 * d__1;
 355 /* Computing 2nd power */
 356         d__2 = ulp;
 357         if ((d__1 = d__[j] * d__[j - 1], abs(d__1)) * (d__2 * d__2) + safemn
 358                 > tmp1) {
 359             isplit[*nsplit] = j - 1;
 360             ++(*nsplit);
 361             work[j - 1] = 0.;
 362         } else {
 363             work[j - 1] = tmp1;
 364             pivmin = max(pivmin,tmp1);
 365         }
 366 /* L10: */
 367     }
 368     isplit[*nsplit] = *n;
 369     pivmin *= safemn;
 370
 371 /*     Compute Interval and ATOLI */
 372
 373     if (irange == 3) {
 374
 375 /*        RANGE='I': Compute the interval containing eigenvalues */
 376 /*                   IL through IU. */
 377
 378 /*        Compute Gershgorin interval for entire (split) matrix */
 379 /*        and use it as the initial interval */
 380
 381         gu = d__[1];
 382         gl = d__[1];
 383         tmp1 = 0.;
 384
 385         i__1 = *n - 1;
 386         for (j = 1; j <= i__1; ++j) {
 387             tmp2 = sqrt(work[j]);
 388 /* Computing MAX */
 389             d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 390             gu = max(d__1,d__2);
 391 /* Computing MIN */
 392             d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 393             gl = min(d__1,d__2);
 394             tmp1 = tmp2;
 395 /* L20: */
 396         }
 397
 398 /* Computing MAX */
 399         d__1 = gu, d__2 = d__[*n] + tmp1;
 400         gu = max(d__1,d__2);
 401 /* Computing MIN */
 402         d__1 = gl, d__2 = d__[*n] - tmp1;
 403         gl = min(d__1,d__2);
 404 /* Computing MAX */
 405         d__1 = abs(gl), d__2 = abs(gu);
 406         tnorm = max(d__1,d__2);
 407         gl = gl - tnorm * 2.1 * ulp * *n - pivmin * 4.2000000000000002;
 408         gu = gu + tnorm * 2.1 * ulp * *n + pivmin * 2.1;
 409
 410 /*        Compute Iteration parameters */
 411
 412         itmax = (integer) ((log(tnorm + pivmin) - log(pivmin)) / log(2.)) + 2;
 413         if (*abstol <= 0.) {
 414             atoli = ulp * tnorm;
 415         } else {
 416             atoli = *abstol;
 417         }
 418
 419         work[*n + 1] = gl;
 420         work[*n + 2] = gl;
 421         work[*n + 3] = gu;
 422         work[*n + 4] = gu;
 423         work[*n + 5] = gl;
 424         work[*n + 6] = gu;
 425         iwork[1] = -1;
 426         iwork[2] = -1;
 427         iwork[3] = *n + 1;
 428         iwork[4] = *n + 1;
 429         iwork[5] = *il - 1;
 430         iwork[6] = *iu;
 431
 432         dlaebz_(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin,
 433                 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n
 434                 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
 435
 436         if (iwork[6] == *iu) {
 437             wl = work[*n + 1];
 438             wlu = work[*n + 3];
 439             nwl = iwork[1];
 440             wu = work[*n + 4];
 441             wul = work[*n + 2];
 442             nwu = iwork[4];
 443         } else {
 444             wl = work[*n + 2];
 445             wlu = work[*n + 4];
 446             nwl = iwork[2];
 447             wu = work[*n + 3];
 448             wul = work[*n + 1];
 449             nwu = iwork[3];
 450         }
 451
 452         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
 453             *info = 4;
 454             return 0;
 455         }
 456     } else {
 457
 458 /*        RANGE='A' or 'V' -- Set ATOLI */
 459
 460 /* Computing MAX */
 461         d__3 = abs(d__[1]) + abs(e[1]), d__4 = (d__1 = d__[*n], abs(d__1)) + (
 462                 d__2 = e[*n - 1], abs(d__2));
 463         tnorm = max(d__3,d__4);
 464
 465         i__1 = *n - 1;
 466         for (j = 2; j <= i__1; ++j) {
 467 /* Computing MAX */
 468             d__4 = tnorm, d__5 = (d__1 = d__[j], abs(d__1)) + (d__2 = e[j - 1]
 469                     , abs(d__2)) + (d__3 = e[j], abs(d__3));
 470             tnorm = max(d__4,d__5);
 471 /* L30: */
 472         }
 473
 474         if (*abstol <= 0.) {
 475             atoli = ulp * tnorm;
 476         } else {
 477             atoli = *abstol;
 478         }
 479
 480         if (irange == 2) {
 481             wl = *vl;
 482             wu = *vu;
 483         } else {
 484             wl = 0.;
 485             wu = 0.;
 486         }
 487     }
 488
 489 /*     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU. */
 490 /*     NWL accumulates the number of eigenvalues .le. WL, */
 491 /*     NWU accumulates the number of eigenvalues .le. WU */
 492
 493     *m = 0;
 494     iend = 0;
 495     *info = 0;
 496     nwl = 0;
 497     nwu = 0;
 498
 499     i__1 = *nsplit;
 500     for (jb = 1; jb <= i__1; ++jb) {
 501         ioff = iend;
 502         ibegin = ioff + 1;
 503         iend = isplit[jb];
 504         in = iend - ioff;
 505
 506         if (in == 1) {
 507
 508 /*           Special Case -- IN=1 */
 509
 510             if (irange == 1 || wl >= d__[ibegin] - pivmin) {
 511                 ++nwl;
 512             }
 513             if (irange == 1 || wu >= d__[ibegin] - pivmin) {
 514                 ++nwu;
 515             }
 516             if (irange == 1 || wl < d__[ibegin] - pivmin && wu >= d__[ibegin]
 517                     - pivmin) {
 518                 ++(*m);
 519                 w[*m] = d__[ibegin];
 520                 iblock[*m] = jb;
 521             }
 522         } else {
 523
 524 /*           General Case -- IN > 1 */
 525
 526 /*           Compute Gershgorin Interval */
 527 /*           and use it as the initial interval */
 528
 529             gu = d__[ibegin];
 530             gl = d__[ibegin];
 531             tmp1 = 0.;
 532
 533             i__2 = iend - 1;
 534             for (j = ibegin; j <= i__2; ++j) {
 535                 tmp2 = (d__1 = e[j], abs(d__1));
 536 /* Computing MAX */
 537                 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
 538                 gu = max(d__1,d__2);
 539 /* Computing MIN */
 540                 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
 541                 gl = min(d__1,d__2);
 542                 tmp1 = tmp2;
 543 /* L40: */
 544             }
 545
 546 /* Computing MAX */
 547             d__1 = gu, d__2 = d__[iend] + tmp1;
 548             gu = max(d__1,d__2);
 549 /* Computing MIN */
 550             d__1 = gl, d__2 = d__[iend] - tmp1;
 551             gl = min(d__1,d__2);
 552 /* Computing MAX */
 553             d__1 = abs(gl), d__2 = abs(gu);
 554             bnorm = max(d__1,d__2);
 555             gl = gl - bnorm * 2.1 * ulp * in - pivmin * 2.1;
 556             gu = gu + bnorm * 2.1 * ulp * in + pivmin * 2.1;
 557
 558 /*           Compute ATOLI for the current submatrix */
 559
 560             if (*abstol <= 0.) {
 561 /* Computing MAX */
 562                 d__1 = abs(gl), d__2 = abs(gu);
 563                 atoli = ulp * max(d__1,d__2);
 564             } else {
 565                 atoli = *abstol;
 566             }
 567
 568             if (irange > 1) {
 569                 if (gu < wl) {
 570                     nwl += in;
 571                     nwu += in;
 572                     goto L70;
 573                 }
 574                 gl = max(gl,wl);
 575                 gu = min(gu,wu);
 576                 if (gl >= gu) {
 577                     goto L70;
 578                 }
 579             }
 580
 581 /*           Set Up Initial Interval */
 582
 583             work[*n + 1] = gl;
 584             work[*n + in + 1] = gu;
 585             dlaebz_(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 586                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 587                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
 588                     w[*m + 1], &iblock[*m + 1], &iinfo);
 589
 590             nwl += iwork[1];
 591             nwu += iwork[in + 1];
 592             iwoff = *m - iwork[1];
 593
 594 /*           Compute Eigenvalues */
 595
 596             itmax = (integer) ((log(gu - gl + pivmin) - log(pivmin)) / log(2.)
 597                     ) + 2;
 598             dlaebz_(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
 599                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
 600                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
 601                      &w[*m + 1], &iblock[*m + 1], &iinfo);
 602
 603 /*           Copy Eigenvalues Into W and IBLOCK */
 604 /*           Use -JB for block number for unconverged eigenvalues. */
 605
 606             i__2 = iout;
 607             for (j = 1; j <= i__2; ++j) {
 608                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
 609
 610 /*              Flag non-convergence. */
 611
 612                 if (j > iout - iinfo) {
 613                     ncnvrg = TRUE_;
 614                     ib = -jb;
 615                 } else {
 616                     ib = jb;
 617                 }
 618                 i__3 = iwork[j + in] + iwoff;
 619                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
 620                     w[je] = tmp1;
 621                     iblock[je] = ib;
 622 /* L50: */
 623                 }
 624 /* L60: */
 625             }
 626
 627             *m += im;
 628         }
 629 L70:
 630         ;
 631     }
 632
 633 /*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
 634 /*     If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
 635
 636     if (irange == 3) {
 637         im = 0;
 638         idiscl = *il - 1 - nwl;
 639         idiscu = nwu - *iu;
 640
 641         if (idiscl > 0 || idiscu > 0) {
 642             i__1 = *m;
 643             for (je = 1; je <= i__1; ++je) {
 644                 if (w[je] <= wlu && idiscl > 0) {
 645                     --idiscl;
 646                 } else if (w[je] >= wul && idiscu > 0) {
 647                     --idiscu;
 648                 } else {
 649                     ++im;
 650                     w[im] = w[je];
 651                     iblock[im] = iblock[je];
 652                 }
 653 /* L80: */
 654             }
 655             *m = im;
 656         }
 657         if (idiscl > 0 || idiscu > 0) {
 658
 659 /*           Code to deal with effects of bad arithmetic: */
 660 /*           Some low eigenvalues to be discarded are not in (WL,WLU], */
 661 /*           or high eigenvalues to be discarded are not in (WUL,WU] */
 662 /*           so just kill off the smallest IDISCL/largest IDISCU */
 663 /*           eigenvalues, by simply finding the smallest/largest */
 664 /*           eigenvalue(s). */
 665
 666 /*           (If N(w) is monotone non-decreasing, this should never */
 667 /*               happen.) */
 668
 669             if (idiscl > 0) {
 670                 wkill = wu;
 671                 i__1 = idiscl;
 672                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
 673                     iw = 0;
 674                     i__2 = *m;
 675                     for (je = 1; je <= i__2; ++je) {
 676                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
 677                             iw = je;
 678                             wkill = w[je];
 679                         }
 680 /* L90: */
 681                     }
 682                     iblock[iw] = 0;
 683 /* L100: */
 684                 }
 685             }
 686             if (idiscu > 0) {
 687
 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             im = 0;
 705             i__1 = *m;
 706             for (je = 1; je <= i__1; ++je) {
 707                 if (iblock[je] != 0) {
 708                     ++im;
 709                     w[im] = w[je];
 710                     iblock[im] = iblock[je];
 711                 }
 712 /* L130: */
 713             }
 714             *m = im;
 715         }
 716         if (idiscl < 0 || idiscu < 0) {
 717             toofew = TRUE_;
 718         }
 719     }
 720
 721 /*     If ORDER='B', do nothing -- the eigenvalues are already sorted */
 722 /*        by block. */
 723 /*     If ORDER='E', sort the eigenvalues from smallest to largest */
 724
 725     if (iorder == 1 && *nsplit > 1) {
 726         i__1 = *m - 1;
 727         for (je = 1; je <= i__1; ++je) {
 728             ie = 0;
 729             tmp1 = w[je];
 730             i__2 = *m;
 731             for (j = je + 1; j <= i__2; ++j) {
 732                 if (w[j] < tmp1) {
 733                     ie = j;
 734                     tmp1 = w[j];
 735                 }
 736 /* L140: */
 737             }
 738
 739             if (ie != 0) {
 740                 itmp1 = iblock[ie];
 741                 w[ie] = w[je];
 742                 iblock[ie] = iblock[je];
 743                 w[je] = tmp1;
 744                 iblock[je] = itmp1;
 745             }
 746 /* L150: */
 747         }
 748     }
 749
 750     *info = 0;
 751     if (ncnvrg) {
 752         ++(*info);
 753     }
 754     if (toofew) {
 755         *info += 2;
 756     }
 757     return 0;
 758
 759 /*     End of DSTEBZ */
 760
 761 } /* dstebz_ */