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