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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static real c_b18 = .003f;
   7
   8 /* Subroutine */ int sstemr_(char *jobz, char *range, integer *n, real *d__,
   9         real *e, real *vl, real *vu, integer *il, integer *iu, integer *m,
  10         real *w, real *z__, integer *ldz, integer *nzc, integer *isuppz,
  11         logical *tryrac, real *work, integer *lwork, integer *iwork, integer *
  12         liwork, integer *info)
  13 {
  14     /* System generated locals */
  15     integer z_dim1, z_offset, i__1, i__2;
  16     real r__1, r__2;
  17
  18     /* Builtin functions */
  19     double sqrt(doublereal);
  20
  21     /* Local variables */
  22     integer i__, j;
  23     real r1, r2;
  24     integer jj;
  25     real cs;
  26     integer in;
  27     real sn, wl, wu;
  28     integer iil, iiu;
  29     real eps, tmp;
  30     integer indd, iend, jblk, wend;
  31     real rmin, rmax;
  32     integer itmp;
  33     real tnrm;
  34     integer inde2;
  35     extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)
  36             ;
  37     integer itmp2;
  38     real rtol1, rtol2, scale;
  39     integer indgp;
  40     extern logical lsame_(char *, char *);
  41     integer iinfo;
  42     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
  43     integer iindw, ilast, lwmin;
  44     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  45             integer *), sswap_(integer *, real *, integer *, real *, integer *
  46 );
  47     logical wantz;
  48     extern /* Subroutine */ int slaev2_(real *, real *, real *, real *, real *
  49 , real *, real *);
  50     logical alleig;
  51     integer ibegin;
  52     logical indeig;
  53     integer iindbl;
  54     logical valeig;
  55     extern doublereal slamch_(char *);
  56     integer wbegin;
  57     real safmin;
  58     extern /* Subroutine */ int xerbla_(char *, integer *);
  59     real bignum;
  60     integer inderr, iindwk, indgrs, offset;
  61     extern /* Subroutine */ int slarrc_(char *, integer *, real *, real *,
  62             real *, real *, real *, integer *, integer *, integer *, integer *
  63 ), slarre_(char *, integer *, real *, real *, integer *,
  64             integer *, real *, real *, real *, real *, real *, real *,
  65             integer *, integer *, integer *, real *, real *, real *, integer *
  66 , integer *, real *, real *, real *, integer *, integer *)
  67             ;
  68     real thresh;
  69     integer iinspl, indwrk, ifirst, liwmin, nzcmin;
  70     real pivmin;
  71     extern doublereal slanst_(char *, integer *, real *, real *);
  72     extern /* Subroutine */ int slarrj_(integer *, real *, real *, integer *,
  73             integer *, real *, integer *, real *, real *, real *, integer *,
  74             real *, real *, integer *), slarrr_(integer *, real *, real *,
  75             integer *);
  76     integer nsplit;
  77     extern /* Subroutine */ int slarrv_(integer *, real *, real *, real *,
  78             real *, real *, integer *, integer *, integer *, integer *, real *
  79 , real *, real *, real *, real *, real *, integer *, integer *,
  80             real *, real *, integer *, integer *, real *, integer *, integer *
  81 );
  82     real smlnum;
  83     extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);
  84     logical lquery, zquery;
  85
  86
  87 /*  -- LAPACK computational routine (version 3.1) -- */
  88 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  89 /*     November 2006 */
  90
  91 /*     .. Scalar Arguments .. */
  92 /*     .. */
  93 /*     .. Array Arguments .. */
  94 /*     .. */
  95
  96 /*  Purpose */
  97 /*  ======= */
  98
  99 /*  SSTEMR computes selected eigenvalues and, optionally, eigenvectors */
 100 /*  of a real symmetric tridiagonal matrix T. Any such unreduced matrix has */
 101 /*  a well defined set of pairwise different real eigenvalues, the corresponding */
 102 /*  real eigenvectors are pairwise orthogonal. */
 103
 104 /*  The spectrum may be computed either completely or partially by specifying */
 105 /*  either an interval (VL,VU] or a range of indices IL:IU for the desired */
 106 /*  eigenvalues. */
 107
 108 /*  Depending on the number of desired eigenvalues, these are computed either */
 109 /*  by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are */
 110 /*  computed by the use of various suitable L D L^T factorizations near clusters */
 111 /*  of close eigenvalues (referred to as RRRs, Relatively Robust */
 112 /*  Representations). An informal sketch of the algorithm follows. */
 113
 114 /*  For each unreduced block (submatrix) of T, */
 115 /*     (a) Compute T - sigma I  = L D L^T, so that L and D */
 116 /*         define all the wanted eigenvalues to high relative accuracy. */
 117 /*         This means that small relative changes in the entries of D and L */
 118 /*         cause only small relative changes in the eigenvalues and */
 119 /*         eigenvectors. The standard (unfactored) representation of the */
 120 /*         tridiagonal matrix T does not have this property in general. */
 121 /*     (b) Compute the eigenvalues to suitable accuracy. */
 122 /*         If the eigenvectors are desired, the algorithm attains full */
 123 /*         accuracy of the computed eigenvalues only right before */
 124 /*         the corresponding vectors have to be computed, see steps c) and d). */
 125 /*     (c) For each cluster of close eigenvalues, select a new */
 126 /*         shift close to the cluster, find a new factorization, and refine */
 127 /*         the shifted eigenvalues to suitable accuracy. */
 128 /*     (d) For each eigenvalue with a large enough relative separation compute */
 129 /*         the corresponding eigenvector by forming a rank revealing twisted */
 130 /*         factorization. Go back to (c) for any clusters that remain. */
 131
 132 /*  For more details, see: */
 133 /*  - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations */
 134 /*    to compute orthogonal eigenvectors of symmetric tridiagonal matrices," */
 135 /*    Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. */
 136 /*  - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and */
 137 /*    Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, */
 138 /*    2004.  Also LAPACK Working Note 154. */
 139 /*  - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric */
 140 /*    tridiagonal eigenvalue/eigenvector problem", */
 141 /*    Computer Science Division Technical Report No. UCB/CSD-97-971, */
 142 /*    UC Berkeley, May 1997. */
 143
 144 /*  Notes: */
 145 /*  1.SSTEMR works only on machines which follow IEEE-754 */
 146 /*  floating-point standard in their handling of infinities and NaNs. */
 147 /*  This permits the use of efficient inner loops avoiding a check for */
 148 /*  zero divisors. */
 149
 150 /*  Arguments */
 151 /*  ========= */
 152
 153 /*  JOBZ    (input) CHARACTER*1 */
 154 /*          = 'N':  Compute eigenvalues only; */
 155 /*          = 'V':  Compute eigenvalues and eigenvectors. */
 156
 157 /*  RANGE   (input) CHARACTER*1 */
 158 /*          = 'A': all eigenvalues will be found. */
 159 /*          = 'V': all eigenvalues in the half-open interval (VL,VU] */
 160 /*                 will be found. */
 161 /*          = 'I': the IL-th through IU-th eigenvalues will be found. */
 162
 163 /*  N       (input) INTEGER */
 164 /*          The order of the matrix.  N >= 0. */
 165
 166 /*  D       (input/output) REAL array, dimension (N) */
 167 /*          On entry, the N diagonal elements of the tridiagonal matrix */
 168 /*          T. On exit, D is overwritten. */
 169
 170 /*  E       (input/output) REAL array, dimension (N) */
 171 /*          On entry, the (N-1) subdiagonal elements of the tridiagonal */
 172 /*          matrix T in elements 1 to N-1 of E. E(N) need not be set on */
 173 /*          input, but is used internally as workspace. */
 174 /*          On exit, E is overwritten. */
 175
 176 /*  VL      (input) REAL */
 177 /*  VU      (input) REAL */
 178 /*          If RANGE='V', the lower and upper bounds of the interval to */
 179 /*          be searched for eigenvalues. VL < VU. */
 180 /*          Not referenced if RANGE = 'A' or 'I'. */
 181
 182 /*  IL      (input) INTEGER */
 183 /*  IU      (input) INTEGER */
 184 /*          If RANGE='I', the indices (in ascending order) of the */
 185 /*          smallest and largest eigenvalues to be returned. */
 186 /*          1 <= IL <= IU <= N, if N > 0. */
 187 /*          Not referenced if RANGE = 'A' or 'V'. */
 188
 189 /*  M       (output) INTEGER */
 190 /*          The total number of eigenvalues found.  0 <= M <= N. */
 191 /*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
 192
 193 /*  W       (output) REAL array, dimension (N) */
 194 /*          The first M elements contain the selected eigenvalues in */
 195 /*          ascending order. */
 196
 197 /*  Z       (output) REAL array, dimension (LDZ, max(1,M) ) */
 198 /*          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z */
 199 /*          contain the orthonormal eigenvectors of the matrix T */
 200 /*          corresponding to the selected eigenvalues, with the i-th */
 201 /*          column of Z holding the eigenvector associated with W(i). */
 202 /*          If JOBZ = 'N', then Z is not referenced. */
 203 /*          Note: the user must ensure that at least max(1,M) columns are */
 204 /*          supplied in the array Z; if RANGE = 'V', the exact value of M */
 205 /*          is not known in advance and can be computed with a workspace */
 206 /*          query by setting NZC = -1, see below. */
 207
 208 /*  LDZ     (input) INTEGER */
 209 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
 210 /*          JOBZ = 'V', then LDZ >= max(1,N). */
 211
 212 /*  NZC     (input) INTEGER */
 213 /*          The number of eigenvectors to be held in the array Z. */
 214 /*          If RANGE = 'A', then NZC >= max(1,N). */
 215 /*          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. */
 216 /*          If RANGE = 'I', then NZC >= IU-IL+1. */
 217 /*          If NZC = -1, then a workspace query is assumed; the */
 218 /*          routine calculates the number of columns of the array Z that */
 219 /*          are needed to hold the eigenvectors. */
 220 /*          This value is returned as the first entry of the Z array, and */
 221 /*          no error message related to NZC is issued by XERBLA. */
 222
 223 /*  ISUPPZ  (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */
 224 /*          The support of the eigenvectors in Z, i.e., the indices */
 225 /*          indicating the nonzero elements in Z. The i-th computed eigenvector */
 226 /*          is nonzero only in elements ISUPPZ( 2*i-1 ) through */
 227 /*          ISUPPZ( 2*i ). This is relevant in the case when the matrix */
 228 /*          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. */
 229
 230 /*  TRYRAC  (input/output) LOGICAL */
 231 /*          If TRYRAC.EQ..TRUE., indicates that the code should check whether */
 232 /*          the tridiagonal matrix defines its eigenvalues to high relative */
 233 /*          accuracy.  If so, the code uses relative-accuracy preserving */
 234 /*          algorithms that might be (a bit) slower depending on the matrix. */
 235 /*          If the matrix does not define its eigenvalues to high relative */
 236 /*          accuracy, the code can uses possibly faster algorithms. */
 237 /*          If TRYRAC.EQ..FALSE., the code is not required to guarantee */
 238 /*          relatively accurate eigenvalues and can use the fastest possible */
 239 /*          techniques. */
 240 /*          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix */
 241 /*          does not define its eigenvalues to high relative accuracy. */
 242
 243 /*  WORK    (workspace/output) REAL array, dimension (LWORK) */
 244 /*          On exit, if INFO = 0, WORK(1) returns the optimal */
 245 /*          (and minimal) LWORK. */
 246
 247 /*  LWORK   (input) INTEGER */
 248 /*          The dimension of the array WORK. LWORK >= max(1,18*N) */
 249 /*          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. */
 250 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 251 /*          only calculates the optimal size of the WORK array, returns */
 252 /*          this value as the first entry of the WORK array, and no error */
 253 /*          message related to LWORK is issued by XERBLA. */
 254
 255 /*  IWORK   (workspace/output) INTEGER array, dimension (LIWORK) */
 256 /*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
 257
 258 /*  LIWORK  (input) INTEGER */
 259 /*          The dimension of the array IWORK.  LIWORK >= max(1,10*N) */
 260 /*          if the eigenvectors are desired, and LIWORK >= max(1,8*N) */
 261 /*          if only the eigenvalues are to be computed. */
 262 /*          If LIWORK = -1, then a workspace query is assumed; the */
 263 /*          routine only calculates the optimal size of the IWORK array, */
 264 /*          returns this value as the first entry of the IWORK array, and */
 265 /*          no error message related to LIWORK is issued by XERBLA. */
 266
 267 /*  INFO    (output) INTEGER */
 268 /*          On exit, INFO */
 269 /*          = 0:  successful exit */
 270 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 271 /*          > 0:  if INFO = 1X, internal error in SLARRE, */
 272 /*                if INFO = 2X, internal error in SLARRV. */
 273 /*                Here, the digit X = ABS( IINFO ) < 10, where IINFO is */
 274 /*                the nonzero error code returned by SLARRE or */
 275 /*                SLARRV, respectively. */
 276
 277
 278 /*  Further Details */
 279 /*  =============== */
 280
 281 /*  Based on contributions by */
 282 /*     Beresford Parlett, University of California, Berkeley, USA */
 283 /*     Jim Demmel, University of California, Berkeley, USA */
 284 /*     Inderjit Dhillon, University of Texas, Austin, USA */
 285 /*     Osni Marques, LBNL/NERSC, USA */
 286 /*     Christof Voemel, University of California, Berkeley, USA */
 287
 288 /*  ===================================================================== */
 289
 290 /*     .. Parameters .. */
 291 /*     .. */
 292 /*     .. Local Scalars .. */
 293 /*     .. */
 294 /*     .. */
 295 /*     .. External Functions .. */
 296 /*     .. */
 297 /*     .. External Subroutines .. */
 298 /*     .. */
 299 /*     .. Intrinsic Functions .. */
 300 /*     .. */
 301 /*     .. Executable Statements .. */
 302
 303 /*     Test the input parameters. */
 304
 305     /* Parameter adjustments */
 306     --d__;
 307     --e;
 308     --w;
 309     z_dim1 = *ldz;
 310     z_offset = 1 + z_dim1;
 311     z__ -= z_offset;
 312     --isuppz;
 313     --work;
 314     --iwork;
 315
 316     /* Function Body */
 317     wantz = lsame_(jobz, "V");
 318     alleig = lsame_(range, "A");
 319     valeig = lsame_(range, "V");
 320     indeig = lsame_(range, "I");
 321
 322     lquery = *lwork == -1 || *liwork == -1;
 323     zquery = *nzc == -1;
 324     *tryrac = *info != 0;
 325 /*     SSTEMR needs WORK of size 6*N, IWORK of size 3*N. */
 326 /*     In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N. */
 327 /*     Furthermore, SLARRV needs WORK of size 12*N, IWORK of size 7*N. */
 328     if (wantz) {
 329         lwmin = *n * 18;
 330         liwmin = *n * 10;
 331     } else {
 332 /*        need less workspace if only the eigenvalues are wanted */
 333         lwmin = *n * 12;
 334         liwmin = *n << 3;
 335     }
 336     wl = 0.f;
 337     wu = 0.f;
 338     iil = 0;
 339     iiu = 0;
 340     if (valeig) {
 341 /*        We do not reference VL, VU in the cases RANGE = 'I','A' */
 342 /*        The interval (WL, WU] contains all the wanted eigenvalues. */
 343 /*        It is either given by the user or computed in SLARRE. */
 344         wl = *vl;
 345         wu = *vu;
 346     } else if (indeig) {
 347 /*        We do not reference IL, IU in the cases RANGE = 'V','A' */
 348         iil = *il;
 349         iiu = *iu;
 350     }
 351
 352     *info = 0;
 353     if (! (wantz || lsame_(jobz, "N"))) {
 354         *info = -1;
 355     } else if (! (alleig || valeig || indeig)) {
 356         *info = -2;
 357     } else if (*n < 0) {
 358         *info = -3;
 359     } else if (valeig && *n > 0 && wu <= wl) {
 360         *info = -7;
 361     } else if (indeig && (iil < 1 || iil > *n)) {
 362         *info = -8;
 363     } else if (indeig && (iiu < iil || iiu > *n)) {
 364         *info = -9;
 365     } else if (*ldz < 1 || wantz && *ldz < *n) {
 366         *info = -13;
 367     } else if (*lwork < lwmin && ! lquery) {
 368         *info = -17;
 369     } else if (*liwork < liwmin && ! lquery) {
 370         *info = -19;
 371     }
 372
 373 /*     Get machine constants. */
 374
 375     safmin = slamch_("Safe minimum");
 376     eps = slamch_("Precision");
 377     smlnum = safmin / eps;
 378     bignum = 1.f / smlnum;
 379     rmin = sqrt(smlnum);
 380 /* Computing MIN */
 381     r__1 = sqrt(bignum), r__2 = 1.f / sqrt(sqrt(safmin));
 382     rmax = dmin(r__1,r__2);
 383
 384     if (*info == 0) {
 385         work[1] = (real) lwmin;
 386         iwork[1] = liwmin;
 387
 388         if (wantz && alleig) {
 389             nzcmin = *n;
 390         } else if (wantz && valeig) {
 391             slarrc_("T", n, vl, vu, &d__[1], &e[1], &safmin, &nzcmin, &itmp, &
 392                     itmp2, info);
 393         } else if (wantz && indeig) {
 394             nzcmin = iiu - iil + 1;
 395         } else {
 396 /*           WANTZ .EQ. FALSE. */
 397             nzcmin = 0;
 398         }
 399         if (zquery && *info == 0) {
 400             z__[z_dim1 + 1] = (real) nzcmin;
 401         } else if (*nzc < nzcmin && ! zquery) {
 402             *info = -14;
 403         }
 404     }
 405     if (*info != 0) {
 406
 407         i__1 = -(*info);
 408         xerbla_("SSTEMR", &i__1);
 409
 410         return 0;
 411     } else if (lquery || zquery) {
 412         return 0;
 413     }
 414
 415 /*     Handle N = 0, 1, and 2 cases immediately */
 416
 417     *m = 0;
 418     if (*n == 0) {
 419         return 0;
 420     }
 421
 422     if (*n == 1) {
 423         if (alleig || indeig) {
 424             *m = 1;
 425             w[1] = d__[1];
 426         } else {
 427             if (wl < d__[1] && wu >= d__[1]) {
 428                 *m = 1;
 429                 w[1] = d__[1];
 430             }
 431         }
 432         if (wantz && ! zquery) {
 433             z__[z_dim1 + 1] = 1.f;
 434             isuppz[1] = 1;
 435             isuppz[2] = 1;
 436         }
 437         return 0;
 438     }
 439
 440     if (*n == 2) {
 441         if (! wantz) {
 442             slae2_(&d__[1], &e[1], &d__[2], &r1, &r2);
 443         } else if (wantz && ! zquery) {
 444             slaev2_(&d__[1], &e[1], &d__[2], &r1, &r2, &cs, &sn);
 445         }
 446         if (alleig || valeig && r2 > wl && r2 <= wu || indeig && iil == 1) {
 447             ++(*m);
 448             w[*m] = r2;
 449             if (wantz && ! zquery) {
 450                 z__[*m * z_dim1 + 1] = -sn;
 451                 z__[*m * z_dim1 + 2] = cs;
 452 /*              Note: At most one of SN and CS can be zero. */
 453                 if (sn != 0.f) {
 454                     if (cs != 0.f) {
 455                         isuppz[(*m << 1) - 1] = 1;
 456                         isuppz[(*m << 1) - 1] = 2;
 457                     } else {
 458                         isuppz[(*m << 1) - 1] = 1;
 459                         isuppz[(*m << 1) - 1] = 1;
 460                     }
 461                 } else {
 462                     isuppz[(*m << 1) - 1] = 2;
 463                     isuppz[*m * 2] = 2;
 464                 }
 465             }
 466         }
 467         if (alleig || valeig && r1 > wl && r1 <= wu || indeig && iiu == 2) {
 468             ++(*m);
 469             w[*m] = r1;
 470             if (wantz && ! zquery) {
 471                 z__[*m * z_dim1 + 1] = cs;
 472                 z__[*m * z_dim1 + 2] = sn;
 473 /*              Note: At most one of SN and CS can be zero. */
 474                 if (sn != 0.f) {
 475                     if (cs != 0.f) {
 476                         isuppz[(*m << 1) - 1] = 1;
 477                         isuppz[(*m << 1) - 1] = 2;
 478                     } else {
 479                         isuppz[(*m << 1) - 1] = 1;
 480                         isuppz[(*m << 1) - 1] = 1;
 481                     }
 482                 } else {
 483                     isuppz[(*m << 1) - 1] = 2;
 484                     isuppz[*m * 2] = 2;
 485                 }
 486             }
 487         }
 488         return 0;
 489     }
 490 /*     Continue with general N */
 491     indgrs = 1;
 492     inderr = (*n << 1) + 1;
 493     indgp = *n * 3 + 1;
 494     indd = (*n << 2) + 1;
 495     inde2 = *n * 5 + 1;
 496     indwrk = *n * 6 + 1;
 497
 498     iinspl = 1;
 499     iindbl = *n + 1;
 500     iindw = (*n << 1) + 1;
 501     iindwk = *n * 3 + 1;
 502
 503 /*     Scale matrix to allowable range, if necessary. */
 504 /*     The allowable range is related to the PIVMIN parameter; see the */
 505 /*     comments in SLARRD.  The preference for scaling small values */
 506 /*     up is heuristic; we expect users' matrices not to be close to the */
 507 /*     RMAX threshold. */
 508
 509     scale = 1.f;
 510     tnrm = slanst_("M", n, &d__[1], &e[1]);
 511     if (tnrm > 0.f && tnrm < rmin) {
 512         scale = rmin / tnrm;
 513     } else if (tnrm > rmax) {
 514         scale = rmax / tnrm;
 515     }
 516     if (scale != 1.f) {
 517         sscal_(n, &scale, &d__[1], &c__1);
 518         i__1 = *n - 1;
 519         sscal_(&i__1, &scale, &e[1], &c__1);
 520         tnrm *= scale;
 521         if (valeig) {
 522 /*           If eigenvalues in interval have to be found, */
 523 /*           scale (WL, WU] accordingly */
 524             wl *= scale;
 525             wu *= scale;
 526         }
 527     }
 528
 529 /*     Compute the desired eigenvalues of the tridiagonal after splitting */
 530 /*     into smaller subblocks if the corresponding off-diagonal elements */
 531 /*     are small */
 532 /*     THRESH is the splitting parameter for SLARRE */
 533 /*     A negative THRESH forces the old splitting criterion based on the */
 534 /*     size of the off-diagonal. A positive THRESH switches to splitting */
 535 /*     which preserves relative accuracy. */
 536
 537     if (*tryrac) {
 538 /*        Test whether the matrix warrants the more expensive relative approach. */
 539         slarrr_(n, &d__[1], &e[1], &iinfo);
 540     } else {
 541 /*        The user does not care about relative accurately eigenvalues */
 542         iinfo = -1;
 543     }
 544 /*     Set the splitting criterion */
 545     if (iinfo == 0) {
 546         thresh = eps;
 547     } else {
 548         thresh = -eps;
 549 /*        relative accuracy is desired but T does not guarantee it */
 550         *tryrac = FALSE_;
 551     }
 552
 553     if (*tryrac) {
 554 /*        Copy original diagonal, needed to guarantee relative accuracy */
 555         scopy_(n, &d__[1], &c__1, &work[indd], &c__1);
 556     }
 557 /*     Store the squares of the offdiagonal values of T */
 558     i__1 = *n - 1;
 559     for (j = 1; j <= i__1; ++j) {
 560 /* Computing 2nd power */
 561         r__1 = e[j];
 562         work[inde2 + j - 1] = r__1 * r__1;
 563 /* L5: */
 564     }
 565 /*     Set the tolerance parameters for bisection */
 566     if (! wantz) {
 567 /*        SLARRE computes the eigenvalues to full precision. */
 568         rtol1 = eps * 4.f;
 569         rtol2 = eps * 4.f;
 570     } else {
 571 /*        SLARRE computes the eigenvalues to less than full precision. */
 572 /*        SLARRV will refine the eigenvalue approximations, and we can */
 573 /*        need less accurate initial bisection in SLARRE. */
 574 /*        Note: these settings do only affect the subset case and SLARRE */
 575 /* Computing MAX */
 576         r__1 = sqrt(eps) * .05f, r__2 = eps * 4.f;
 577         rtol1 = dmax(r__1,r__2);
 578 /* Computing MAX */
 579         r__1 = sqrt(eps) * .005f, r__2 = eps * 4.f;
 580         rtol2 = dmax(r__1,r__2);
 581     }
 582     slarre_(range, n, &wl, &wu, &iil, &iiu, &d__[1], &e[1], &work[inde2], &
 583             rtol1, &rtol2, &thresh, &nsplit, &iwork[iinspl], m, &w[1], &work[
 584             inderr], &work[indgp], &iwork[iindbl], &iwork[iindw], &work[
 585             indgrs], &pivmin, &work[indwrk], &iwork[iindwk], &iinfo);
 586     if (iinfo != 0) {
 587         *info = abs(iinfo) + 10;
 588         return 0;
 589     }
 590 /*     Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired */
 591 /*     part of the spectrum. All desired eigenvalues are contained in */
 592 /*     (WL,WU] */
 593     if (wantz) {
 594
 595 /*        Compute the desired eigenvectors corresponding to the computed */
 596 /*        eigenvalues */
 597
 598         slarrv_(n, &wl, &wu, &d__[1], &e[1], &pivmin, &iwork[iinspl], m, &
 599                 c__1, m, &c_b18, &rtol1, &rtol2, &w[1], &work[inderr], &work[
 600                 indgp], &iwork[iindbl], &iwork[iindw], &work[indgrs], &z__[
 601                 z_offset], ldz, &isuppz[1], &work[indwrk], &iwork[iindwk], &
 602                 iinfo);
 603         if (iinfo != 0) {
 604             *info = abs(iinfo) + 20;
 605             return 0;
 606         }
 607     } else {
 608 /*        SLARRE computes eigenvalues of the (shifted) root representation */
 609 /*        SLARRV returns the eigenvalues of the unshifted matrix. */
 610 /*        However, if the eigenvectors are not desired by the user, we need */
 611 /*        to apply the corresponding shifts from SLARRE to obtain the */
 612 /*        eigenvalues of the original matrix. */
 613         i__1 = *m;
 614         for (j = 1; j <= i__1; ++j) {
 615             itmp = iwork[iindbl + j - 1];
 616             w[j] += e[iwork[iinspl + itmp - 1]];
 617 /* L20: */
 618         }
 619     }
 620
 621     if (*tryrac) {
 622 /*        Refine computed eigenvalues so that they are relatively accurate */
 623 /*        with respect to the original matrix T. */
 624         ibegin = 1;
 625         wbegin = 1;
 626         i__1 = iwork[iindbl + *m - 1];
 627         for (jblk = 1; jblk <= i__1; ++jblk) {
 628             iend = iwork[iinspl + jblk - 1];
 629             in = iend - ibegin + 1;
 630             wend = wbegin - 1;
 631 /*           check if any eigenvalues have to be refined in this block */
 632 L36:
 633             if (wend < *m) {
 634                 if (iwork[iindbl + wend] == jblk) {
 635                     ++wend;
 636                     goto L36;
 637                 }
 638             }
 639             if (wend < wbegin) {
 640                 ibegin = iend + 1;
 641                 goto L39;
 642             }
 643             offset = iwork[iindw + wbegin - 1] - 1;
 644             ifirst = iwork[iindw + wbegin - 1];
 645             ilast = iwork[iindw + wend - 1];
 646             rtol2 = eps * 4.f;
 647             slarrj_(&in, &work[indd + ibegin - 1], &work[inde2 + ibegin - 1],
 648                     &ifirst, &ilast, &rtol2, &offset, &w[wbegin], &work[
 649                     inderr + wbegin - 1], &work[indwrk], &iwork[iindwk], &
 650                     pivmin, &tnrm, &iinfo);
 651             ibegin = iend + 1;
 652             wbegin = wend + 1;
 653 L39:
 654             ;
 655         }
 656     }
 657
 658 /*     If matrix was scaled, then rescale eigenvalues appropriately. */
 659
 660     if (scale != 1.f) {
 661         r__1 = 1.f / scale;
 662         sscal_(m, &r__1, &w[1], &c__1);
 663     }
 664
 665 /*     If eigenvalues are not in increasing order, then sort them, */
 666 /*     possibly along with eigenvectors. */
 667
 668     if (nsplit > 1) {
 669         if (! wantz) {
 670             slasrt_("I", m, &w[1], &iinfo);
 671             if (iinfo != 0) {
 672                 *info = 3;
 673                 return 0;
 674             }
 675         } else {
 676             i__1 = *m - 1;
 677             for (j = 1; j <= i__1; ++j) {
 678                 i__ = 0;
 679                 tmp = w[j];
 680                 i__2 = *m;
 681                 for (jj = j + 1; jj <= i__2; ++jj) {
 682                     if (w[jj] < tmp) {
 683                         i__ = jj;
 684                         tmp = w[jj];
 685                     }
 686 /* L50: */
 687                 }
 688                 if (i__ != 0) {
 689                     w[i__] = w[j];
 690                     w[j] = tmp;
 691                     if (wantz) {
 692                         sswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j *
 693                                 z_dim1 + 1], &c__1);
 694                         itmp = isuppz[(i__ << 1) - 1];
 695                         isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
 696                         isuppz[(j << 1) - 1] = itmp;
 697                         itmp = isuppz[i__ * 2];
 698                         isuppz[i__ * 2] = isuppz[j * 2];
 699                         isuppz[j * 2] = itmp;
 700                     }
 701                 }
 702 /* L60: */
 703             }
 704         }
 705     }
 706
 707
 708     work[1] = (real) lwmin;
 709     iwork[1] = liwmin;
 710     return 0;
 711
 712 /*     End of SSTEMR */
 713
 714 } /* sstemr_ */