git.maemo.org Git - opencv/blob - 3rdparty/lapack/dlarrv.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static doublereal c_b5 = 0.;
   6 static integer c__1 = 1;
   7 static integer c__2 = 2;
   8
   9 /* Subroutine */ int dlarrv_(integer *n, doublereal *vl, doublereal *vu,
  10         doublereal *d__, doublereal *l, doublereal *pivmin, integer *isplit,
  11         integer *m, integer *dol, integer *dou, doublereal *minrgp,
  12         doublereal *rtol1, doublereal *rtol2, doublereal *w, doublereal *werr,
  13          doublereal *wgap, integer *iblock, integer *indexw, doublereal *gers,
  14          doublereal *z__, integer *ldz, integer *isuppz, doublereal *work,
  15         integer *iwork, integer *info)
  16 {
  17     /* System generated locals */
  18     integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
  19     doublereal d__1, d__2;
  20     logical L__1;
  21
  22     /* Builtin functions */
  23     double log(doublereal);
  24
  25     /* Local variables */
  26     integer minwsize, i__, j, k, p, q, miniwsize, ii;
  27     doublereal gl;
  28     integer im, in;
  29     doublereal gu, gap, eps, tau, tol, tmp;
  30     integer zto;
  31     doublereal ztz;
  32     integer iend, jblk;
  33     doublereal lgap;
  34     integer done;
  35     doublereal rgap, left;
  36     integer wend, iter;
  37     doublereal bstw;
  38     integer itmp1;
  39     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  40             integer *);
  41     integer indld;
  42     doublereal fudge;
  43     integer idone;
  44     doublereal sigma;
  45     integer iinfo, iindr;
  46     doublereal resid;
  47     logical eskip;
  48     doublereal right;
  49     extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
  50             doublereal *, integer *);
  51     integer nclus, zfrom;
  52     doublereal rqtol;
  53     integer iindc1, iindc2;
  54     extern /* Subroutine */ int dlar1v_(integer *, integer *, integer *,
  55             doublereal *, doublereal *, doublereal *, doublereal *,
  56             doublereal *, doublereal *, doublereal *, doublereal *, logical *,
  57              integer *, doublereal *, doublereal *, integer *, integer *,
  58             doublereal *, doublereal *, doublereal *, doublereal *);
  59     logical stp2ii;
  60     doublereal lambda;
  61     extern doublereal dlamch_(char *);
  62     integer ibegin, indeig;
  63     logical needbs;
  64     integer indlld;
  65     doublereal sgndef, mingma;
  66     extern /* Subroutine */ int dlarrb_(integer *, doublereal *, doublereal *,
  67              integer *, integer *, doublereal *, doublereal *, integer *,
  68             doublereal *, doublereal *, doublereal *, doublereal *, integer *,
  69              doublereal *, doublereal *, integer *, integer *);
  70     integer oldien, oldncl, wbegin;
  71     doublereal spdiam;
  72     integer negcnt;
  73     extern /* Subroutine */ int dlarrf_(integer *, doublereal *, doublereal *,
  74              doublereal *, integer *, integer *, doublereal *, doublereal *,
  75             doublereal *, doublereal *, doublereal *, doublereal *,
  76             doublereal *, doublereal *, doublereal *, doublereal *,
  77             doublereal *, integer *);
  78     integer oldcls;
  79     doublereal savgap;
  80     integer ndepth;
  81     doublereal ssigma;
  82     extern /* Subroutine */ int dlaset_(char *, integer *, integer *,
  83             doublereal *, doublereal *, doublereal *, integer *);
  84     logical usedbs;
  85     integer iindwk, offset;
  86     doublereal gaptol;
  87     integer newcls, oldfst, indwrk, windex, oldlst;
  88     logical usedrq;
  89     integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl;
  90     doublereal bstres;
  91     integer newsiz, zusedu, zusedw;
  92     doublereal nrminv, rqcorr;
  93     logical tryrqc;
  94     integer isupmx;
  95
  96
  97 /*  -- LAPACK auxiliary routine (version 3.1.1) -- */
  98 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  99 /*     November 2006 */
 100
 101 /*     .. Scalar Arguments .. */
 102 /*     .. */
 103 /*     .. Array Arguments .. */
 104 /*     .. */
 105
 106 /*  Purpose */
 107 /*  ======= */
 108
 109 /*  DLARRV computes the eigenvectors of the tridiagonal matrix */
 110 /*  T = L D L^T given L, D and APPROXIMATIONS to the eigenvalues of L D L^T. */
 111 /*  The input eigenvalues should have been computed by DLARRE. */
 112
 113 /*  Arguments */
 114 /*  ========= */
 115
 116 /*  N       (input) INTEGER */
 117 /*          The order of the matrix.  N >= 0. */
 118
 119 /*  VL      (input) DOUBLE PRECISION */
 120 /*  VU      (input) DOUBLE PRECISION */
 121 /*          Lower and upper bounds of the interval that contains the desired */
 122 /*          eigenvalues. VL < VU. Needed to compute gaps on the left or right */
 123 /*          end of the extremal eigenvalues in the desired RANGE. */
 124
 125 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
 126 /*          On entry, the N diagonal elements of the diagonal matrix D. */
 127 /*          On exit, D may be overwritten. */
 128
 129 /*  L       (input/output) DOUBLE PRECISION array, dimension (N) */
 130 /*          On entry, the (N-1) subdiagonal elements of the unit */
 131 /*          bidiagonal matrix L are in elements 1 to N-1 of L */
 132 /*          (if the matrix is not splitted.) At the end of each block */
 133 /*          is stored the corresponding shift as given by DLARRE. */
 134 /*          On exit, L is overwritten. */
 135
 136 /*  PIVMIN  (in) DOUBLE PRECISION */
 137 /*          The minimum pivot allowed in the Sturm sequence. */
 138
 139 /*  ISPLIT  (input) INTEGER array, dimension (N) */
 140 /*          The splitting points, at which T breaks up into blocks. */
 141 /*          The first block consists of rows/columns 1 to */
 142 /*          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 */
 143 /*          through ISPLIT( 2 ), etc. */
 144
 145 /*  M       (input) INTEGER */
 146 /*          The total number of input eigenvalues.  0 <= M <= N. */
 147
 148 /*  DOL     (input) INTEGER */
 149 /*  DOU     (input) INTEGER */
 150 /*          If the user wants to compute only selected eigenvectors from all */
 151 /*          the eigenvalues supplied, he can specify an index range DOL:DOU. */
 152 /*          Or else the setting DOL=1, DOU=M should be applied. */
 153 /*          Note that DOL and DOU refer to the order in which the eigenvalues */
 154 /*          are stored in W. */
 155 /*          If the user wants to compute only selected eigenpairs, then */
 156 /*          the columns DOL-1 to DOU+1 of the eigenvector space Z contain the */
 157 /*          computed eigenvectors. All other columns of Z are set to zero. */
 158
 159 /*  MINRGP  (input) DOUBLE PRECISION */
 160
 161 /*  RTOL1   (input) DOUBLE PRECISION */
 162 /*  RTOL2   (input) DOUBLE PRECISION */
 163 /*           Parameters for bisection. */
 164 /*           An interval [LEFT,RIGHT] has converged if */
 165 /*           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) */
 166
 167 /*  W       (input/output) DOUBLE PRECISION array, dimension (N) */
 168 /*          The first M elements of W contain the APPROXIMATE eigenvalues for */
 169 /*          which eigenvectors are to be computed.  The eigenvalues */
 170 /*          should be grouped by split-off block and ordered from */
 171 /*          smallest to largest within the block ( The output array */
 172 /*          W from DLARRE is expected here ). Furthermore, they are with */
 173 /*          respect to the shift of the corresponding root representation */
 174 /*          for their block. On exit, W holds the eigenvalues of the */
 175 /*          UNshifted matrix. */
 176
 177 /*  WERR    (input/output) DOUBLE PRECISION array, dimension (N) */
 178 /*          The first M elements contain the semiwidth of the uncertainty */
 179 /*          interval of the corresponding eigenvalue in W */
 180
 181 /*  WGAP    (input/output) DOUBLE PRECISION array, dimension (N) */
 182 /*          The separation from the right neighbor eigenvalue in W. */
 183
 184 /*  IBLOCK  (input) INTEGER array, dimension (N) */
 185 /*          The indices of the blocks (submatrices) associated with the */
 186 /*          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue */
 187 /*          W(i) belongs to the first block from the top, =2 if W(i) */
 188 /*          belongs to the second block, etc. */
 189
 190 /*  INDEXW  (input) INTEGER array, dimension (N) */
 191 /*          The indices of the eigenvalues within each block (submatrix); */
 192 /*          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the */
 193 /*          i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. */
 194
 195 /*  GERS    (input) DOUBLE PRECISION array, dimension (2*N) */
 196 /*          The N Gerschgorin intervals (the i-th Gerschgorin interval */
 197 /*          is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should */
 198 /*          be computed from the original UNshifted matrix. */
 199
 200 /*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
 201 /*          If INFO = 0, the first M columns of Z contain the */
 202 /*          orthonormal eigenvectors of the matrix T */
 203 /*          corresponding to the input eigenvalues, with the i-th */
 204 /*          column of Z holding the eigenvector associated with W(i). */
 205 /*          Note: the user must ensure that at least max(1,M) columns are */
 206 /*          supplied in the array Z. */
 207
 208 /*  LDZ     (input) INTEGER */
 209 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
 210 /*          JOBZ = 'V', LDZ >= max(1,N). */
 211
 212 /*  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) ) */
 213 /*          The support of the eigenvectors in Z, i.e., the indices */
 214 /*          indicating the nonzero elements in Z. The I-th eigenvector */
 215 /*          is nonzero only in elements ISUPPZ( 2*I-1 ) through */
 216 /*          ISUPPZ( 2*I ). */
 217
 218 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (12*N) */
 219
 220 /*  IWORK   (workspace) INTEGER array, dimension (7*N) */
 221
 222 /*  INFO    (output) INTEGER */
 223 /*          = 0:  successful exit */
 224
 225 /*          > 0:  A problem occured in DLARRV. */
 226 /*          < 0:  One of the called subroutines signaled an internal problem. */
 227 /*                Needs inspection of the corresponding parameter IINFO */
 228 /*                for further information. */
 229
 230 /*          =-1:  Problem in DLARRB when refining a child's eigenvalues. */
 231 /*          =-2:  Problem in DLARRF when computing the RRR of a child. */
 232 /*                When a child is inside a tight cluster, it can be difficult */
 233 /*                to find an RRR. A partial remedy from the user's point of */
 234 /*                view is to make the parameter MINRGP smaller and recompile. */
 235 /*                However, as the orthogonality of the computed vectors is */
 236 /*                proportional to 1/MINRGP, the user should be aware that */
 237 /*                he might be trading in precision when he decreases MINRGP. */
 238 /*          =-3:  Problem in DLARRB when refining a single eigenvalue */
 239 /*                after the Rayleigh correction was rejected. */
 240 /*          = 5:  The Rayleigh Quotient Iteration failed to converge to */
 241 /*                full accuracy in MAXITR steps. */
 242
 243 /*  Further Details */
 244 /*  =============== */
 245
 246 /*  Based on contributions by */
 247 /*     Beresford Parlett, University of California, Berkeley, USA */
 248 /*     Jim Demmel, University of California, Berkeley, USA */
 249 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 250 /*     Osni Marques, LBNL/NERSC, USA */
 251 /*     Christof Voemel, University of California, Berkeley, USA */
 252
 253 /*  ===================================================================== */
 254
 255 /*     .. Parameters .. */
 256 /*     .. */
 257 /*     .. Local Scalars .. */
 258 /*     .. */
 259 /*     .. External Functions .. */
 260 /*     .. */
 261 /*     .. External Subroutines .. */
 262 /*     .. */
 263 /*     .. Intrinsic Functions .. */
 264 /*     .. */
 265 /*     .. Executable Statements .. */
 266 /*     .. */
 267 /*     The first N entries of WORK are reserved for the eigenvalues */
 268     /* Parameter adjustments */
 269     --d__;
 270     --l;
 271     --isplit;
 272     --w;
 273     --werr;
 274     --wgap;
 275     --iblock;
 276     --indexw;
 277     --gers;
 278     z_dim1 = *ldz;
 279     z_offset = 1 + z_dim1;
 280     z__ -= z_offset;
 281     --isuppz;
 282     --work;
 283     --iwork;
 284
 285     /* Function Body */
 286     indld = *n + 1;
 287     indlld = (*n << 1) + 1;
 288     indwrk = *n * 3 + 1;
 289     minwsize = *n * 12;
 290     i__1 = minwsize;
 291     for (i__ = 1; i__ <= i__1; ++i__) {
 292         work[i__] = 0.;
 293 /* L5: */
 294     }
 295 /*     IWORK(IINDR+1:IINDR+N) hold the twist indices R for the */
 296 /*     factorization used to compute the FP vector */
 297     iindr = 0;
 298 /*     IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current */
 299 /*     layer and the one above. */
 300     iindc1 = *n;
 301     iindc2 = *n << 1;
 302     iindwk = *n * 3 + 1;
 303     miniwsize = *n * 7;
 304     i__1 = miniwsize;
 305     for (i__ = 1; i__ <= i__1; ++i__) {
 306         iwork[i__] = 0;
 307 /* L10: */
 308     }
 309     zusedl = 1;
 310     if (*dol > 1) {
 311 /*        Set lower bound for use of Z */
 312         zusedl = *dol - 1;
 313     }
 314     zusedu = *m;
 315     if (*dou < *m) {
 316 /*        Set lower bound for use of Z */
 317         zusedu = *dou + 1;
 318     }
 319 /*     The width of the part of Z that is used */
 320     zusedw = zusedu - zusedl + 1;
 321     dlaset_("Full", n, &zusedw, &c_b5, &c_b5, &z__[zusedl * z_dim1 + 1], ldz);
 322     eps = dlamch_("Precision");
 323     rqtol = eps * 2.;
 324
 325 /*     Set expert flags for standard code. */
 326     tryrqc = TRUE_;
 327     if (*dol == 1 && *dou == *m) {
 328     } else {
 329 /*        Only selected eigenpairs are computed. Since the other evalues */
 330 /*        are not refined by RQ iteration, bisection has to compute to full */
 331 /*        accuracy. */
 332         *rtol1 = eps * 4.;
 333         *rtol2 = eps * 4.;
 334     }
 335 /*     The entries WBEGIN:WEND in W, WERR, WGAP correspond to the */
 336 /*     desired eigenvalues. The support of the nonzero eigenvector */
 337 /*     entries is contained in the interval IBEGIN:IEND. */
 338 /*     Remark that if k eigenpairs are desired, then the eigenvectors */
 339 /*     are stored in k contiguous columns of Z. */
 340 /*     DONE is the number of eigenvectors already computed */
 341     done = 0;
 342     ibegin = 1;
 343     wbegin = 1;
 344     i__1 = iblock[*m];
 345     for (jblk = 1; jblk <= i__1; ++jblk) {
 346         iend = isplit[jblk];
 347         sigma = l[iend];
 348 /*        Find the eigenvectors of the submatrix indexed IBEGIN */
 349 /*        through IEND. */
 350         wend = wbegin - 1;
 351 L15:
 352         if (wend < *m) {
 353             if (iblock[wend + 1] == jblk) {
 354                 ++wend;
 355                 goto L15;
 356             }
 357         }
 358         if (wend < wbegin) {
 359             ibegin = iend + 1;
 360             goto L170;
 361         } else if (wend < *dol || wbegin > *dou) {
 362             ibegin = iend + 1;
 363             wbegin = wend + 1;
 364             goto L170;
 365         }
 366 /*        Find local spectral diameter of the block */
 367         gl = gers[(ibegin << 1) - 1];
 368         gu = gers[ibegin * 2];
 369         i__2 = iend;
 370         for (i__ = ibegin + 1; i__ <= i__2; ++i__) {
 371 /* Computing MIN */
 372             d__1 = gers[(i__ << 1) - 1];
 373             gl = min(d__1,gl);
 374 /* Computing MAX */
 375             d__1 = gers[i__ * 2];
 376             gu = max(d__1,gu);
 377 /* L20: */
 378         }
 379         spdiam = gu - gl;
 380 /*        OLDIEN is the last index of the previous block */
 381         oldien = ibegin - 1;
 382 /*        Calculate the size of the current block */
 383         in = iend - ibegin + 1;
 384 /*        The number of eigenvalues in the current block */
 385         im = wend - wbegin + 1;
 386 /*        This is for a 1x1 block */
 387         if (ibegin == iend) {
 388             ++done;
 389             z__[ibegin + wbegin * z_dim1] = 1.;
 390             isuppz[(wbegin << 1) - 1] = ibegin;
 391             isuppz[wbegin * 2] = ibegin;
 392             w[wbegin] += sigma;
 393             work[wbegin] = w[wbegin];
 394             ibegin = iend + 1;
 395             ++wbegin;
 396             goto L170;
 397         }
 398 /*        The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) */
 399 /*        Note that these can be approximations, in this case, the corresp. */
 400 /*        entries of WERR give the size of the uncertainty interval. */
 401 /*        The eigenvalue approximations will be refined when necessary as */
 402 /*        high relative accuracy is required for the computation of the */
 403 /*        corresponding eigenvectors. */
 404         dcopy_(&im, &w[wbegin], &c__1, &work[wbegin], &c__1);
 405 /*        We store in W the eigenvalue approximations w.r.t. the original */
 406 /*        matrix T. */
 407         i__2 = im;
 408         for (i__ = 1; i__ <= i__2; ++i__) {
 409             w[wbegin + i__ - 1] += sigma;
 410 /* L30: */
 411         }
 412 /*        NDEPTH is the current depth of the representation tree */
 413         ndepth = 0;
 414 /*        PARITY is either 1 or 0 */
 415         parity = 1;
 416 /*        NCLUS is the number of clusters for the next level of the */
 417 /*        representation tree, we start with NCLUS = 1 for the root */
 418         nclus = 1;
 419         iwork[iindc1 + 1] = 1;
 420         iwork[iindc1 + 2] = im;
 421 /*        IDONE is the number of eigenvectors already computed in the current */
 422 /*        block */
 423         idone = 0;
 424 /*        loop while( IDONE.LT.IM ) */
 425 /*        generate the representation tree for the current block and */
 426 /*        compute the eigenvectors */
 427 L40:
 428         if (idone < im) {
 429 /*           This is a crude protection against infinitely deep trees */
 430             if (ndepth > *m) {
 431                 *info = -2;
 432                 return 0;
 433             }
 434 /*           breadth first processing of the current level of the representation */
 435 /*           tree: OLDNCL = number of clusters on current level */
 436             oldncl = nclus;
 437 /*           reset NCLUS to count the number of child clusters */
 438             nclus = 0;
 439
 440             parity = 1 - parity;
 441             if (parity == 0) {
 442                 oldcls = iindc1;
 443                 newcls = iindc2;
 444             } else {
 445                 oldcls = iindc2;
 446                 newcls = iindc1;
 447             }
 448 /*           Process the clusters on the current level */
 449             i__2 = oldncl;
 450             for (i__ = 1; i__ <= i__2; ++i__) {
 451                 j = oldcls + (i__ << 1);
 452 /*              OLDFST, OLDLST = first, last index of current cluster. */
 453 /*                               cluster indices start with 1 and are relative */
 454 /*                               to WBEGIN when accessing W, WGAP, WERR, Z */
 455                 oldfst = iwork[j - 1];
 456                 oldlst = iwork[j];
 457                 if (ndepth > 0) {
 458 /*                 Retrieve relatively robust representation (RRR) of cluster */
 459 /*                 that has been computed at the previous level */
 460 /*                 The RRR is stored in Z and overwritten once the eigenvectors */
 461 /*                 have been computed or when the cluster is refined */
 462                     if (*dol == 1 && *dou == *m) {
 463 /*                    Get representation from location of the leftmost evalue */
 464 /*                    of the cluster */
 465                         j = wbegin + oldfst - 1;
 466                     } else {
 467                         if (wbegin + oldfst - 1 < *dol) {
 468 /*                       Get representation from the left end of Z array */
 469                             j = *dol - 1;
 470                         } else if (wbegin + oldfst - 1 > *dou) {
 471 /*                       Get representation from the right end of Z array */
 472                             j = *dou;
 473                         } else {
 474                             j = wbegin + oldfst - 1;
 475                         }
 476                     }
 477                     dcopy_(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin]
 478 , &c__1);
 479                     i__3 = in - 1;
 480                     dcopy_(&i__3, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[
 481                             ibegin], &c__1);
 482                     sigma = z__[iend + (j + 1) * z_dim1];
 483 /*                 Set the corresponding entries in Z to zero */
 484                     dlaset_("Full", &in, &c__2, &c_b5, &c_b5, &z__[ibegin + j
 485                             * z_dim1], ldz);
 486                 }
 487 /*              Compute DL and DLL of current RRR */
 488                 i__3 = iend - 1;
 489                 for (j = ibegin; j <= i__3; ++j) {
 490                     tmp = d__[j] * l[j];
 491                     work[indld - 1 + j] = tmp;
 492                     work[indlld - 1 + j] = tmp * l[j];
 493 /* L50: */
 494                 }
 495                 if (ndepth > 0) {
 496 /*                 P and Q are index of the first and last eigenvalue to compute */
 497 /*                 within the current block */
 498                     p = indexw[wbegin - 1 + oldfst];
 499                     q = indexw[wbegin - 1 + oldlst];
 500 /*                 Offset for the arrays WORK, WGAP and WERR, i.e., th P-OFFSET */
 501 /*                 thru' Q-OFFSET elements of these arrays are to be used. */
 502 /*                  OFFSET = P-OLDFST */
 503                     offset = indexw[wbegin] - 1;
 504 /*                 perform limited bisection (if necessary) to get approximate */
 505 /*                 eigenvalues to the precision needed. */
 506                     dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p,
 507                              &q, rtol1, rtol2, &offset, &work[wbegin], &wgap[
 508                             wbegin], &werr[wbegin], &work[indwrk], &iwork[
 509                             iindwk], pivmin, &spdiam, &in, &iinfo);
 510                     if (iinfo != 0) {
 511                         *info = -1;
 512                         return 0;
 513                     }
 514 /*                 We also recompute the extremal gaps. W holds all eigenvalues */
 515 /*                 of the unshifted matrix and must be used for computation */
 516 /*                 of WGAP, the entries of WORK might stem from RRRs with */
 517 /*                 different shifts. The gaps from WBEGIN-1+OLDFST to */
 518 /*                 WBEGIN-1+OLDLST are correctly computed in DLARRB. */
 519 /*                 However, we only allow the gaps to become greater since */
 520 /*                 this is what should happen when we decrease WERR */
 521                     if (oldfst > 1) {
 522 /* Computing MAX */
 523                         d__1 = wgap[wbegin + oldfst - 2], d__2 = w[wbegin +
 524                                 oldfst - 1] - werr[wbegin + oldfst - 1] - w[
 525                                 wbegin + oldfst - 2] - werr[wbegin + oldfst -
 526                                 2];
 527                         wgap[wbegin + oldfst - 2] = max(d__1,d__2);
 528                     }
 529                     if (wbegin + oldlst - 1 < wend) {
 530 /* Computing MAX */
 531                         d__1 = wgap[wbegin + oldlst - 1], d__2 = w[wbegin +
 532                                 oldlst] - werr[wbegin + oldlst] - w[wbegin +
 533                                 oldlst - 1] - werr[wbegin + oldlst - 1];
 534                         wgap[wbegin + oldlst - 1] = max(d__1,d__2);
 535                     }
 536 /*                 Each time the eigenvalues in WORK get refined, we store */
 537 /*                 the newly found approximation with all shifts applied in W */
 538                     i__3 = oldlst;
 539                     for (j = oldfst; j <= i__3; ++j) {
 540                         w[wbegin + j - 1] = work[wbegin + j - 1] + sigma;
 541 /* L53: */
 542                     }
 543                 }
 544 /*              Process the current node. */
 545                 newfst = oldfst;
 546                 i__3 = oldlst;
 547                 for (j = oldfst; j <= i__3; ++j) {
 548                     if (j == oldlst) {
 549 /*                    we are at the right end of the cluster, this is also the */
 550 /*                    boundary of the child cluster */
 551                         newlst = j;
 552                     } else if (wgap[wbegin + j - 1] >= *minrgp * (d__1 = work[
 553                             wbegin + j - 1], abs(d__1))) {
 554 /*                    the right relative gap is big enough, the child cluster */
 555 /*                    (NEWFST,..,NEWLST) is well separated from the following */
 556                         newlst = j;
 557                     } else {
 558 /*                    inside a child cluster, the relative gap is not */
 559 /*                    big enough. */
 560                         goto L140;
 561                     }
 562 /*                 Compute size of child cluster found */
 563                     newsiz = newlst - newfst + 1;
 564 /*                 NEWFTT is the place in Z where the new RRR or the computed */
 565 /*                 eigenvector is to be stored */
 566                     if (*dol == 1 && *dou == *m) {
 567 /*                    Store representation at location of the leftmost evalue */
 568 /*                    of the cluster */
 569                         newftt = wbegin + newfst - 1;
 570                     } else {
 571                         if (wbegin + newfst - 1 < *dol) {
 572 /*                       Store representation at the left end of Z array */
 573                             newftt = *dol - 1;
 574                         } else if (wbegin + newfst - 1 > *dou) {
 575 /*                       Store representation at the right end of Z array */
 576                             newftt = *dou;
 577                         } else {
 578                             newftt = wbegin + newfst - 1;
 579                         }
 580                     }
 581                     if (newsiz > 1) {
 582
 583 /*                    Current child is not a singleton but a cluster. */
 584 /*                    Compute and store new representation of child. */
 585
 586
 587 /*                    Compute left and right cluster gap. */
 588
 589 /*                    LGAP and RGAP are not computed from WORK because */
 590 /*                    the eigenvalue approximations may stem from RRRs */
 591 /*                    different shifts. However, W hold all eigenvalues */
 592 /*                    of the unshifted matrix. Still, the entries in WGAP */
 593 /*                    have to be computed from WORK since the entries */
 594 /*                    in W might be of the same order so that gaps are not */
 595 /*                    exhibited correctly for very close eigenvalues. */
 596                         if (newfst == 1) {
 597 /* Computing MAX */
 598                             d__1 = 0., d__2 = w[wbegin] - werr[wbegin] - *vl;
 599                             lgap = max(d__1,d__2);
 600                         } else {
 601                             lgap = wgap[wbegin + newfst - 2];
 602                         }
 603                         rgap = wgap[wbegin + newlst - 1];
 604
 605 /*                    Compute left- and rightmost eigenvalue of child */
 606 /*                    to high precision in order to shift as close */
 607 /*                    as possible and obtain as large relative gaps */
 608 /*                    as possible */
 609
 610                         for (k = 1; k <= 2; ++k) {
 611                             if (k == 1) {
 612                                 p = indexw[wbegin - 1 + newfst];
 613                             } else {
 614                                 p = indexw[wbegin - 1 + newlst];
 615                             }
 616                             offset = indexw[wbegin] - 1;
 617                             dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin
 618                                     - 1], &p, &p, &rqtol, &rqtol, &offset, &
 619                                     work[wbegin], &wgap[wbegin], &werr[wbegin]
 620 , &work[indwrk], &iwork[iindwk], pivmin, &
 621                                     spdiam, &in, &iinfo);
 622 /* L55: */
 623                         }
 624
 625                         if (wbegin + newlst - 1 < *dol || wbegin + newfst - 1
 626                                 > *dou) {
 627 /*                       if the cluster contains no desired eigenvalues */
 628 /*                       skip the computation of that branch of the rep. tree */
 629
 630 /*                       We could skip before the refinement of the extremal */
 631 /*                       eigenvalues of the child, but then the representation */
 632 /*                       tree could be different from the one when nothing is */
 633 /*                       skipped. For this reason we skip at this place. */
 634                             idone = idone + newlst - newfst + 1;
 635                             goto L139;
 636                         }
 637
 638 /*                    Compute RRR of child cluster. */
 639 /*                    Note that the new RRR is stored in Z */
 640
 641 /*                    DLARRF needs LWORK = 2*N */
 642                         dlarrf_(&in, &d__[ibegin], &l[ibegin], &work[indld +
 643                                 ibegin - 1], &newfst, &newlst, &work[wbegin],
 644                                 &wgap[wbegin], &werr[wbegin], &spdiam, &lgap,
 645                                 &rgap, pivmin, &tau, &z__[ibegin + newftt *
 646                                 z_dim1], &z__[ibegin + (newftt + 1) * z_dim1],
 647                                  &work[indwrk], &iinfo);
 648                         if (iinfo == 0) {
 649 /*                       a new RRR for the cluster was found by DLARRF */
 650 /*                       update shift and store it */
 651                             ssigma = sigma + tau;
 652                             z__[iend + (newftt + 1) * z_dim1] = ssigma;
 653 /*                       WORK() are the midpoints and WERR() the semi-width */
 654 /*                       Note that the entries in W are unchanged. */
 655                             i__4 = newlst;
 656                             for (k = newfst; k <= i__4; ++k) {
 657                                 fudge = eps * 3. * (d__1 = work[wbegin + k -
 658                                         1], abs(d__1));
 659                                 work[wbegin + k - 1] -= tau;
 660                                 fudge += eps * 4. * (d__1 = work[wbegin + k -
 661                                         1], abs(d__1));
 662 /*                          Fudge errors */
 663                                 werr[wbegin + k - 1] += fudge;
 664 /*                          Gaps are not fudged. Provided that WERR is small */
 665 /*                          when eigenvalues are close, a zero gap indicates */
 666 /*                          that a new representation is needed for resolving */
 667 /*                          the cluster. A fudge could lead to a wrong decision */
 668 /*                          of judging eigenvalues 'separated' which in */
 669 /*                          reality are not. This could have a negative impact */
 670 /*                          on the orthogonality of the computed eigenvectors. */
 671 /* L116: */
 672                             }
 673                             ++nclus;
 674                             k = newcls + (nclus << 1);
 675                             iwork[k - 1] = newfst;
 676                             iwork[k] = newlst;
 677                         } else {
 678                             *info = -2;
 679                             return 0;
 680                         }
 681                     } else {
 682
 683 /*                    Compute eigenvector of singleton */
 684
 685                         iter = 0;
 686
 687                         tol = log((doublereal) in) * 4. * eps;
 688
 689                         k = newfst;
 690                         windex = wbegin + k - 1;
 691 /* Computing MAX */
 692                         i__4 = windex - 1;
 693                         windmn = max(i__4,1);
 694 /* Computing MIN */
 695                         i__4 = windex + 1;
 696                         windpl = min(i__4,*m);
 697                         lambda = work[windex];
 698                         ++done;
 699 /*                    Check if eigenvector computation is to be skipped */
 700                         if (windex < *dol || windex > *dou) {
 701                             eskip = TRUE_;
 702                             goto L125;
 703                         } else {
 704                             eskip = FALSE_;
 705                         }
 706                         left = work[windex] - werr[windex];
 707                         right = work[windex] + werr[windex];
 708                         indeig = indexw[windex];
 709 /*                    Note that since we compute the eigenpairs for a child, */
 710 /*                    all eigenvalue approximations are w.r.t the same shift. */
 711 /*                    In this case, the entries in WORK should be used for */
 712 /*                    computing the gaps since they exhibit even very small */
 713 /*                    differences in the eigenvalues, as opposed to the */
 714 /*                    entries in W which might "look" the same. */
 715                         if (k == 1) {
 716 /*                       In the case RANGE='I' and with not much initial */
 717 /*                       accuracy in LAMBDA and VL, the formula */
 718 /*                       LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) */
 719 /*                       can lead to an overestimation of the left gap and */
 720 /*                       thus to inadequately early RQI 'convergence'. */
 721 /*                       Prevent this by forcing a small left gap. */
 722 /* Computing MAX */
 723                             d__1 = abs(left), d__2 = abs(right);
 724                             lgap = eps * max(d__1,d__2);
 725                         } else {
 726                             lgap = wgap[windmn];
 727                         }
 728                         if (k == im) {
 729 /*                       In the case RANGE='I' and with not much initial */
 730 /*                       accuracy in LAMBDA and VU, the formula */
 731 /*                       can lead to an overestimation of the right gap and */
 732 /*                       thus to inadequately early RQI 'convergence'. */
 733 /*                       Prevent this by forcing a small right gap. */
 734 /* Computing MAX */
 735                             d__1 = abs(left), d__2 = abs(right);
 736                             rgap = eps * max(d__1,d__2);
 737                         } else {
 738                             rgap = wgap[windex];
 739                         }
 740                         gap = min(lgap,rgap);
 741                         if (k == 1 || k == im) {
 742 /*                       The eigenvector support can become wrong */
 743 /*                       because significant entries could be cut off due to a */
 744 /*                       large GAPTOL parameter in LAR1V. Prevent this. */
 745                             gaptol = 0.;
 746                         } else {
 747                             gaptol = gap * eps;
 748                         }
 749                         isupmn = in;
 750                         isupmx = 1;
 751 /*                    Update WGAP so that it holds the minimum gap */
 752 /*                    to the left or the right. This is crucial in the */
 753 /*                    case where bisection is used to ensure that the */
 754 /*                    eigenvalue is refined up to the required precision. */
 755 /*                    The correct value is restored afterwards. */
 756                         savgap = wgap[windex];
 757                         wgap[windex] = gap;
 758 /*                    We want to use the Rayleigh Quotient Correction */
 759 /*                    as often as possible since it converges quadratically */
 760 /*                    when we are close enough to the desired eigenvalue. */
 761 /*                    However, the Rayleigh Quotient can have the wrong sign */
 762 /*                    and lead us away from the desired eigenvalue. In this */
 763 /*                    case, the best we can do is to use bisection. */
 764                         usedbs = FALSE_;
 765                         usedrq = FALSE_;
 766 /*                    Bisection is initially turned off unless it is forced */
 767                         needbs = ! tryrqc;
 768 L120:
 769 /*                    Check if bisection should be used to refine eigenvalue */
 770                         if (needbs) {
 771 /*                       Take the bisection as new iterate */
 772                             usedbs = TRUE_;
 773                             itmp1 = iwork[iindr + windex];
 774                             offset = indexw[wbegin] - 1;
 775                             d__1 = eps * 2.;
 776                             dlarrb_(&in, &d__[ibegin], &work[indlld + ibegin
 777                                     - 1], &indeig, &indeig, &c_b5, &d__1, &
 778                                     offset, &work[wbegin], &wgap[wbegin], &
 779                                     werr[wbegin], &work[indwrk], &iwork[
 780                                     iindwk], pivmin, &spdiam, &itmp1, &iinfo);
 781                             if (iinfo != 0) {
 782                                 *info = -3;
 783                                 return 0;
 784                             }
 785                             lambda = work[windex];
 786 /*                       Reset twist index from inaccurate LAMBDA to */
 787 /*                       force computation of true MINGMA */
 788                             iwork[iindr + windex] = 0;
 789                         }
 790 /*                    Given LAMBDA, compute the eigenvector. */
 791                         L__1 = ! usedbs;
 792                         dlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &l[
 793                                 ibegin], &work[indld + ibegin - 1], &work[
 794                                 indlld + ibegin - 1], pivmin, &gaptol, &z__[
 795                                 ibegin + windex * z_dim1], &L__1, &negcnt, &
 796                                 ztz, &mingma, &iwork[iindr + windex], &isuppz[
 797                                 (windex << 1) - 1], &nrminv, &resid, &rqcorr,
 798                                 &work[indwrk]);
 799                         if (iter == 0) {
 800                             bstres = resid;
 801                             bstw = lambda;
 802                         } else if (resid < bstres) {
 803                             bstres = resid;
 804                             bstw = lambda;
 805                         }
 806 /* Computing MIN */
 807                         i__4 = isupmn, i__5 = isuppz[(windex << 1) - 1];
 808                         isupmn = min(i__4,i__5);
 809 /* Computing MAX */
 810                         i__4 = isupmx, i__5 = isuppz[windex * 2];
 811                         isupmx = max(i__4,i__5);
 812                         ++iter;
 813 /*                    sin alpha <= |resid|/gap */
 814 /*                    Note that both the residual and the gap are */
 815 /*                    proportional to the matrix, so ||T|| doesn't play */
 816 /*                    a role in the quotient */
 817
 818 /*                    Convergence test for Rayleigh-Quotient iteration */
 819 /*                    (omitted when Bisection has been used) */
 820
 821                         if (resid > tol * gap && abs(rqcorr) > rqtol * abs(
 822                                 lambda) && ! usedbs) {
 823 /*                       We need to check that the RQCORR update doesn't */
 824 /*                       move the eigenvalue away from the desired one and */
 825 /*                       towards a neighbor. -> protection with bisection */
 826                             if (indeig <= negcnt) {
 827 /*                          The wanted eigenvalue lies to the left */
 828                                 sgndef = -1.;
 829                             } else {
 830 /*                          The wanted eigenvalue lies to the right */
 831                                 sgndef = 1.;
 832                             }
 833 /*                       We only use the RQCORR if it improves the */
 834 /*                       the iterate reasonably. */
 835                             if (rqcorr * sgndef >= 0. && lambda + rqcorr <=
 836                                     right && lambda + rqcorr >= left) {
 837                                 usedrq = TRUE_;
 838 /*                          Store new midpoint of bisection interval in WORK */
 839                                 if (sgndef == 1.) {
 840 /*                             The current LAMBDA is on the left of the true */
 841 /*                             eigenvalue */
 842                                     left = lambda;
 843 /*                             We prefer to assume that the error estimate */
 844 /*                             is correct. We could make the interval not */
 845 /*                             as a bracket but to be modified if the RQCORR */
 846 /*                             chooses to. In this case, the RIGHT side should */
 847 /*                             be modified as follows: */
 848 /*                              RIGHT = MAX(RIGHT, LAMBDA + RQCORR) */
 849                                 } else {
 850 /*                             The current LAMBDA is on the right of the true */
 851 /*                             eigenvalue */
 852                                     right = lambda;
 853 /*                             See comment about assuming the error estimate is */
 854 /*                             correct above. */
 855 /*                              LEFT = MIN(LEFT, LAMBDA + RQCORR) */
 856                                 }
 857                                 work[windex] = (right + left) * .5;
 858 /*                          Take RQCORR since it has the correct sign and */
 859 /*                          improves the iterate reasonably */
 860                                 lambda += rqcorr;
 861 /*                          Update width of error interval */
 862                                 werr[windex] = (right - left) * .5;
 863                             } else {
 864                                 needbs = TRUE_;
 865                             }
 866                             if (right - left < rqtol * abs(lambda)) {
 867 /*                             The eigenvalue is computed to bisection accuracy */
 868 /*                             compute eigenvector and stop */
 869                                 usedbs = TRUE_;
 870                                 goto L120;
 871                             } else if (iter < 10) {
 872                                 goto L120;
 873                             } else if (iter == 10) {
 874                                 needbs = TRUE_;
 875                                 goto L120;
 876                             } else {
 877                                 *info = 5;
 878                                 return 0;
 879                             }
 880                         } else {
 881                             stp2ii = FALSE_;
 882                             if (usedrq && usedbs && bstres <= resid) {
 883                                 lambda = bstw;
 884                                 stp2ii = TRUE_;
 885                             }
 886                             if (stp2ii) {
 887 /*                          improve error angle by second step */
 888                                 L__1 = ! usedbs;
 889                                 dlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin]
 890 , &l[ibegin], &work[indld + ibegin -
 891                                         1], &work[indlld + ibegin - 1],
 892                                         pivmin, &gaptol, &z__[ibegin + windex
 893                                         * z_dim1], &L__1, &negcnt, &ztz, &
 894                                         mingma, &iwork[iindr + windex], &
 895                                         isuppz[(windex << 1) - 1], &nrminv, &
 896                                         resid, &rqcorr, &work[indwrk]);
 897                             }
 898                             work[windex] = lambda;
 899                         }
 900
 901 /*                    Compute FP-vector support w.r.t. whole matrix */
 902
 903                         isuppz[(windex << 1) - 1] += oldien;
 904                         isuppz[windex * 2] += oldien;
 905                         zfrom = isuppz[(windex << 1) - 1];
 906                         zto = isuppz[windex * 2];
 907                         isupmn += oldien;
 908                         isupmx += oldien;
 909 /*                    Ensure vector is ok if support in the RQI has changed */
 910                         if (isupmn < zfrom) {
 911                             i__4 = zfrom - 1;
 912                             for (ii = isupmn; ii <= i__4; ++ii) {
 913                                 z__[ii + windex * z_dim1] = 0.;
 914 /* L122: */
 915                             }
 916                         }
 917                         if (isupmx > zto) {
 918                             i__4 = isupmx;
 919                             for (ii = zto + 1; ii <= i__4; ++ii) {
 920                                 z__[ii + windex * z_dim1] = 0.;
 921 /* L123: */
 922                             }
 923                         }
 924                         i__4 = zto - zfrom + 1;
 925                         dscal_(&i__4, &nrminv, &z__[zfrom + windex * z_dim1],
 926                                 &c__1);
 927 L125:
 928 /*                    Update W */
 929                         w[windex] = lambda + sigma;
 930 /*                    Recompute the gaps on the left and right */
 931 /*                    But only allow them to become larger and not */
 932 /*                    smaller (which can only happen through "bad" */
 933 /*                    cancellation and doesn't reflect the theory */
 934 /*                    where the initial gaps are underestimated due */
 935 /*                    to WERR being too crude.) */
 936                         if (! eskip) {
 937                             if (k > 1) {
 938 /* Computing MAX */
 939                                 d__1 = wgap[windmn], d__2 = w[windex] - werr[
 940                                         windex] - w[windmn] - werr[windmn];
 941                                 wgap[windmn] = max(d__1,d__2);
 942                             }
 943                             if (windex < wend) {
 944 /* Computing MAX */
 945                                 d__1 = savgap, d__2 = w[windpl] - werr[windpl]
 946                                          - w[windex] - werr[windex];
 947                                 wgap[windex] = max(d__1,d__2);
 948                             }
 949                         }
 950                         ++idone;
 951                     }
 952 /*                 here ends the code for the current child */
 953
 954 L139:
 955 /*                 Proceed to any remaining child nodes */
 956                     newfst = j + 1;
 957 L140:
 958                     ;
 959                 }
 960 /* L150: */
 961             }
 962             ++ndepth;
 963             goto L40;
 964         }
 965         ibegin = iend + 1;
 966         wbegin = wend + 1;
 967 L170:
 968         ;
 969     }
 970
 971     return 0;
 972
 973 /*     End of DLARRV */
 974
 975 } /* dlarrv_ */