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

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