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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__10 = 10;
   6 static integer c__1 = 1;
   7 static integer c__2 = 2;
   8 static integer c__3 = 3;
   9 static integer c__4 = 4;
  10 static integer c_n1 = -1;
  11
  12 /* Subroutine */ int ssyevr_(char *jobz, char *range, char *uplo, integer *n,
  13         real *a, integer *lda, real *vl, real *vu, integer *il, integer *iu,
  14         real *abstol, integer *m, real *w, real *z__, integer *ldz, integer *
  15         isuppz, real *work, integer *lwork, integer *iwork, integer *liwork,
  16         integer *info)
  17 {
  18     /* System generated locals */
  19     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
  20     real r__1, r__2;
  21
  22     /* Builtin functions */
  23     double sqrt(doublereal);
  24
  25     /* Local variables */
  26     integer i__, j, nb, jj;
  27     real eps, vll, vuu, tmp1;
  28     integer indd, inde;
  29     real anrm;
  30     integer imax;
  31     real rmin, rmax;
  32     logical test;
  33     integer inddd, indee;
  34     real sigma;
  35     extern logical lsame_(char *, char *);
  36     integer iinfo;
  37     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *);
  38     char order[1];
  39     integer indwk, lwmin;
  40     logical lower;
  41     extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
  42             integer *), sswap_(integer *, real *, integer *, real *, integer *
  43 );
  44     logical wantz, alleig, indeig;
  45     integer iscale, ieeeok, indibl, indifl;
  46     logical valeig;
  47     extern doublereal slamch_(char *);
  48     real safmin;
  49     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  50             integer *, integer *);
  51     extern /* Subroutine */ int xerbla_(char *, integer *);
  52     real abstll, bignum;
  53     integer indtau, indisp, indiwo, indwkn, liwmin;
  54     logical tryrac;
  55     extern /* Subroutine */ int sstein_(integer *, real *, real *, integer *,
  56             real *, integer *, integer *, real *, integer *, real *, integer *
  57 , integer *, integer *), ssterf_(integer *, real *, real *,
  58             integer *);
  59     integer llwrkn, llwork, nsplit;
  60     real smlnum;
  61     extern doublereal slansy_(char *, char *, integer *, real *, integer *,
  62             real *);
  63     extern /* Subroutine */ int sstebz_(char *, char *, integer *, real *,
  64             real *, integer *, integer *, real *, real *, real *, integer *,
  65             integer *, real *, integer *, integer *, real *, integer *,
  66             integer *), sstemr_(char *, char *, integer *,
  67             real *, real *, real *, real *, integer *, integer *, integer *,
  68             real *, real *, integer *, integer *, integer *, logical *, real *
  69 , integer *, integer *, integer *, integer *);
  70     integer lwkopt;
  71     logical lquery;
  72     extern /* Subroutine */ int sormtr_(char *, char *, char *, integer *,
  73             integer *, real *, integer *, real *, real *, integer *, real *,
  74             integer *, integer *), ssytrd_(char *,
  75             integer *, real *, integer *, real *, real *, real *, real *,
  76             integer *, integer *);
  77
  78
  79 /*  -- LAPACK driver routine (version 3.1) -- */
  80 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  81 /*     November 2006 */
  82
  83 /*     .. Scalar Arguments .. */
  84 /*     .. */
  85 /*     .. Array Arguments .. */
  86 /*     .. */
  87
  88 /*  Purpose */
  89 /*  ======= */
  90
  91 /*  SSYEVR computes selected eigenvalues and, optionally, eigenvectors */
  92 /*  of a real symmetric matrix A.  Eigenvalues and eigenvectors can be */
  93 /*  selected by specifying either a range of values or a range of */
  94 /*  indices for the desired eigenvalues. */
  95
  96 /*  SSYEVR first reduces the matrix A to tridiagonal form T with a call */
  97 /*  to SSYTRD.  Then, whenever possible, SSYEVR calls SSTEMR to compute */
  98 /*  the eigenspectrum using Relatively Robust Representations.  SSTEMR */
  99 /*  computes eigenvalues by the dqds algorithm, while orthogonal */
 100 /*  eigenvectors are computed from various "good" L D L^T representations */
 101 /*  (also known as Relatively Robust Representations). Gram-Schmidt */
 102 /*  orthogonalization is avoided as far as possible. More specifically, */
 103 /*  the various steps of the algorithm are as follows. */
 104
 105 /*  For each unreduced block (submatrix) of T, */
 106 /*     (a) Compute T - sigma I  = L D L^T, so that L and D */
 107 /*         define all the wanted eigenvalues to high relative accuracy. */
 108 /*         This means that small relative changes in the entries of D and L */
 109 /*         cause only small relative changes in the eigenvalues and */
 110 /*         eigenvectors. The standard (unfactored) representation of the */
 111 /*         tridiagonal matrix T does not have this property in general. */
 112 /*     (b) Compute the eigenvalues to suitable accuracy. */
 113 /*         If the eigenvectors are desired, the algorithm attains full */
 114 /*         accuracy of the computed eigenvalues only right before */
 115 /*         the corresponding vectors have to be computed, see steps c) and d). */
 116 /*     (c) For each cluster of close eigenvalues, select a new */
 117 /*         shift close to the cluster, find a new factorization, and refine */
 118 /*         the shifted eigenvalues to suitable accuracy. */
 119 /*     (d) For each eigenvalue with a large enough relative separation compute */
 120 /*         the corresponding eigenvector by forming a rank revealing twisted */
 121 /*         factorization. Go back to (c) for any clusters that remain. */
 122
 123 /*  The desired accuracy of the output can be specified by the input */
 124 /*  parameter ABSTOL. */
 125
 126 /*  For more details, see SSTEMR's documentation and: */
 127 /*  - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations */
 128 /*    to compute orthogonal eigenvectors of symmetric tridiagonal matrices," */
 129 /*    Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. */
 130 /*  - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and */
 131 /*    Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, */
 132 /*    2004.  Also LAPACK Working Note 154. */
 133 /*  - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric */
 134 /*    tridiagonal eigenvalue/eigenvector problem", */
 135 /*    Computer Science Division Technical Report No. UCB/CSD-97-971, */
 136 /*    UC Berkeley, May 1997. */
 137
 138
 139 /*  Note 1 : SSYEVR calls SSTEMR when the full spectrum is requested */
 140 /*  on machines which conform to the ieee-754 floating point standard. */
 141 /*  SSYEVR calls SSTEBZ and SSTEIN on non-ieee machines and */
 142 /*  when partial spectrum requests are made. */
 143
 144 /*  Normal execution of SSTEMR may create NaNs and infinities and */
 145 /*  hence may abort due to a floating point exception in environments */
 146 /*  which do not handle NaNs and infinities in the ieee standard default */
 147 /*  manner. */
 148
 149 /*  Arguments */
 150 /*  ========= */
 151
 152 /*  JOBZ    (input) CHARACTER*1 */
 153 /*          = 'N':  Compute eigenvalues only; */
 154 /*          = 'V':  Compute eigenvalues and eigenvectors. */
 155
 156 /*  RANGE   (input) CHARACTER*1 */
 157 /*          = 'A': all eigenvalues will be found. */
 158 /*          = 'V': all eigenvalues in the half-open interval (VL,VU] */
 159 /*                 will be found. */
 160 /*          = 'I': the IL-th through IU-th eigenvalues will be found. */
 161 /* ********* For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and */
 162 /* ********* SSTEIN are called */
 163
 164 /*  UPLO    (input) CHARACTER*1 */
 165 /*          = 'U':  Upper triangle of A is stored; */
 166 /*          = 'L':  Lower triangle of A is stored. */
 167
 168 /*  N       (input) INTEGER */
 169 /*          The order of the matrix A.  N >= 0. */
 170
 171 /*  A       (input/output) REAL array, dimension (LDA, N) */
 172 /*          On entry, the symmetric matrix A.  If UPLO = 'U', the */
 173 /*          leading N-by-N upper triangular part of A contains the */
 174 /*          upper triangular part of the matrix A.  If UPLO = 'L', */
 175 /*          the leading N-by-N lower triangular part of A contains */
 176 /*          the lower triangular part of the matrix A. */
 177 /*          On exit, the lower triangle (if UPLO='L') or the upper */
 178 /*          triangle (if UPLO='U') of A, including the diagonal, is */
 179 /*          destroyed. */
 180
 181 /*  LDA     (input) INTEGER */
 182 /*          The leading dimension of the array A.  LDA >= max(1,N). */
 183
 184 /*  VL      (input) REAL */
 185 /*  VU      (input) REAL */
 186 /*          If RANGE='V', the lower and upper bounds of the interval to */
 187 /*          be searched for eigenvalues. VL < VU. */
 188 /*          Not referenced if RANGE = 'A' or 'I'. */
 189
 190 /*  IL      (input) INTEGER */
 191 /*  IU      (input) INTEGER */
 192 /*          If RANGE='I', the indices (in ascending order) of the */
 193 /*          smallest and largest eigenvalues to be returned. */
 194 /*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
 195 /*          Not referenced if RANGE = 'A' or 'V'. */
 196
 197 /*  ABSTOL  (input) REAL */
 198 /*          The absolute error tolerance for the eigenvalues. */
 199 /*          An approximate eigenvalue is accepted as converged */
 200 /*          when it is determined to lie in an interval [a,b] */
 201 /*          of width less than or equal to */
 202
 203 /*                  ABSTOL + EPS *   max( |a|,|b| ) , */
 204
 205 /*          where EPS is the machine precision.  If ABSTOL is less than */
 206 /*          or equal to zero, then  EPS*|T|  will be used in its place, */
 207 /*          where |T| is the 1-norm of the tridiagonal matrix obtained */
 208 /*          by reducing A to tridiagonal form. */
 209
 210 /*          See "Computing Small Singular Values of Bidiagonal Matrices */
 211 /*          with Guaranteed High Relative Accuracy," by Demmel and */
 212 /*          Kahan, LAPACK Working Note #3. */
 213
 214 /*          If high relative accuracy is important, set ABSTOL to */
 215 /*          SLAMCH( 'Safe minimum' ).  Doing so will guarantee that */
 216 /*          eigenvalues are computed to high relative accuracy when */
 217 /*          possible in future releases.  The current code does not */
 218 /*          make any guarantees about high relative accuracy, but */
 219 /*          future releases will. See J. Barlow and J. Demmel, */
 220 /*          "Computing Accurate Eigensystems of Scaled Diagonally */
 221 /*          Dominant Matrices", LAPACK Working Note #7, for a discussion */
 222 /*          of which matrices define their eigenvalues to high relative */
 223 /*          accuracy. */
 224
 225 /*  M       (output) INTEGER */
 226 /*          The total number of eigenvalues found.  0 <= M <= N. */
 227 /*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
 228
 229 /*  W       (output) REAL array, dimension (N) */
 230 /*          The first M elements contain the selected eigenvalues in */
 231 /*          ascending order. */
 232
 233 /*  Z       (output) REAL array, dimension (LDZ, max(1,M)) */
 234 /*          If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
 235 /*          contain the orthonormal eigenvectors of the matrix A */
 236 /*          corresponding to the selected eigenvalues, with the i-th */
 237 /*          column of Z holding the eigenvector associated with W(i). */
 238 /*          If JOBZ = 'N', then Z is not referenced. */
 239 /*          Note: the user must ensure that at least max(1,M) columns are */
 240 /*          supplied in the array Z; if RANGE = 'V', the exact value of M */
 241 /*          is not known in advance and an upper bound must be used. */
 242 /*          Supplying N columns is always safe. */
 243
 244 /*  LDZ     (input) INTEGER */
 245 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
 246 /*          JOBZ = 'V', LDZ >= max(1,N). */
 247
 248 /*  ISUPPZ  (output) INTEGER array, dimension ( 2*max(1,M) ) */
 249 /*          The support of the eigenvectors in Z, i.e., the indices */
 250 /*          indicating the nonzero elements in Z. The i-th eigenvector */
 251 /*          is nonzero only in elements ISUPPZ( 2*i-1 ) through */
 252 /*          ISUPPZ( 2*i ). */
 253 /* ********* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 */
 254
 255 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
 256 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 257
 258 /*  LWORK   (input) INTEGER */
 259 /*          The dimension of the array WORK.  LWORK >= max(1,26*N). */
 260 /*          For optimal efficiency, LWORK >= (NB+6)*N, */
 261 /*          where NB is the max of the blocksize for SSYTRD and SORMTR */
 262 /*          returned by ILAENV. */
 263
 264 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 265 /*          only calculates the optimal sizes of the WORK and IWORK */
 266 /*          arrays, returns these values as the first entries of the WORK */
 267 /*          and IWORK arrays, and no error message related to LWORK or */
 268 /*          LIWORK is issued by XERBLA. */
 269
 270 /*  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */
 271 /*          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK. */
 272
 273 /*  LIWORK  (input) INTEGER */
 274 /*          The dimension of the array IWORK.  LIWORK >= max(1,10*N). */
 275
 276 /*          If LIWORK = -1, then a workspace query is assumed; the */
 277 /*          routine only calculates the optimal sizes of the WORK and */
 278 /*          IWORK arrays, returns these values as the first entries of */
 279 /*          the WORK and IWORK arrays, and no error message related to */
 280 /*          LWORK or LIWORK is issued by XERBLA. */
 281
 282 /*  INFO    (output) INTEGER */
 283 /*          = 0:  successful exit */
 284 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 285 /*          > 0:  Internal error */
 286
 287 /*  Further Details */
 288 /*  =============== */
 289
 290 /*  Based on contributions by */
 291 /*     Inderjit Dhillon, IBM Almaden, USA */
 292 /*     Osni Marques, LBNL/NERSC, USA */
 293 /*     Ken Stanley, Computer Science Division, University of */
 294 /*       California at Berkeley, USA */
 295 /*     Jason Riedy, Computer Science Division, University of */
 296 /*       California at Berkeley, USA */
 297
 298 /* ===================================================================== */
 299
 300 /*     .. Parameters .. */
 301 /*     .. */
 302 /*     .. Local Scalars .. */
 303 /*     .. */
 304 /*     .. External Functions .. */
 305 /*     .. */
 306 /*     .. External Subroutines .. */
 307 /*     .. */
 308 /*     .. Intrinsic Functions .. */
 309 /*     .. */
 310 /*     .. Executable Statements .. */
 311
 312 /*     Test the input parameters. */
 313
 314     /* Parameter adjustments */
 315     a_dim1 = *lda;
 316     a_offset = 1 + a_dim1;
 317     a -= a_offset;
 318     --w;
 319     z_dim1 = *ldz;
 320     z_offset = 1 + z_dim1;
 321     z__ -= z_offset;
 322     --isuppz;
 323     --work;
 324     --iwork;
 325
 326     /* Function Body */
 327     ieeeok = ilaenv_(&c__10, "SSYEVR", "N", &c__1, &c__2, &c__3, &c__4);
 328
 329     lower = lsame_(uplo, "L");
 330     wantz = lsame_(jobz, "V");
 331     alleig = lsame_(range, "A");
 332     valeig = lsame_(range, "V");
 333     indeig = lsame_(range, "I");
 334
 335     lquery = *lwork == -1 || *liwork == -1;
 336
 337 /* Computing MAX */
 338     i__1 = 1, i__2 = *n * 26;
 339     lwmin = max(i__1,i__2);
 340 /* Computing MAX */
 341     i__1 = 1, i__2 = *n * 10;
 342     liwmin = max(i__1,i__2);
 343
 344     *info = 0;
 345     if (! (wantz || lsame_(jobz, "N"))) {
 346         *info = -1;
 347     } else if (! (alleig || valeig || indeig)) {
 348         *info = -2;
 349     } else if (! (lower || lsame_(uplo, "U"))) {
 350         *info = -3;
 351     } else if (*n < 0) {
 352         *info = -4;
 353     } else if (*lda < max(1,*n)) {
 354         *info = -6;
 355     } else {
 356         if (valeig) {
 357             if (*n > 0 && *vu <= *vl) {
 358                 *info = -8;
 359             }
 360         } else if (indeig) {
 361             if (*il < 1 || *il > max(1,*n)) {
 362                 *info = -9;
 363             } else if (*iu < min(*n,*il) || *iu > *n) {
 364                 *info = -10;
 365             }
 366         }
 367     }
 368     if (*info == 0) {
 369         if (*ldz < 1 || wantz && *ldz < *n) {
 370             *info = -15;
 371         }
 372     }
 373
 374     if (*info == 0) {
 375         nb = ilaenv_(&c__1, "SSYTRD", uplo, n, &c_n1, &c_n1, &c_n1);
 376 /* Computing MAX */
 377         i__1 = nb, i__2 = ilaenv_(&c__1, "SORMTR", uplo, n, &c_n1, &c_n1, &
 378                 c_n1);
 379         nb = max(i__1,i__2);
 380 /* Computing MAX */
 381         i__1 = (nb + 1) * *n;
 382         lwkopt = max(i__1,lwmin);
 383         work[1] = (real) lwkopt;
 384         iwork[1] = liwmin;
 385
 386         if (*lwork < lwmin && ! lquery) {
 387             *info = -18;
 388         } else if (*liwork < liwmin && ! lquery) {
 389             *info = -20;
 390         }
 391     }
 392
 393     if (*info != 0) {
 394         i__1 = -(*info);
 395         xerbla_("SSYEVR", &i__1);
 396         return 0;
 397     } else if (lquery) {
 398         return 0;
 399     }
 400
 401 /*     Quick return if possible */
 402
 403     *m = 0;
 404     if (*n == 0) {
 405         work[1] = 1.f;
 406         return 0;
 407     }
 408
 409     if (*n == 1) {
 410         work[1] = 26.f;
 411         if (alleig || indeig) {
 412             *m = 1;
 413             w[1] = a[a_dim1 + 1];
 414         } else {
 415             if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
 416                 *m = 1;
 417                 w[1] = a[a_dim1 + 1];
 418             }
 419         }
 420         if (wantz) {
 421             z__[z_dim1 + 1] = 1.f;
 422         }
 423         return 0;
 424     }
 425
 426 /*     Get machine constants. */
 427
 428     safmin = slamch_("Safe minimum");
 429     eps = slamch_("Precision");
 430     smlnum = safmin / eps;
 431     bignum = 1.f / smlnum;
 432     rmin = sqrt(smlnum);
 433 /* Computing MIN */
 434     r__1 = sqrt(bignum), r__2 = 1.f / sqrt(sqrt(safmin));
 435     rmax = dmin(r__1,r__2);
 436
 437 /*     Scale matrix to allowable range, if necessary. */
 438
 439     iscale = 0;
 440     abstll = *abstol;
 441     if (valeig) {
 442         vll = *vl;
 443         vuu = *vu;
 444     }
 445     anrm = slansy_("M", uplo, n, &a[a_offset], lda, &work[1]);
 446     if (anrm > 0.f && anrm < rmin) {
 447         iscale = 1;
 448         sigma = rmin / anrm;
 449     } else if (anrm > rmax) {
 450         iscale = 1;
 451         sigma = rmax / anrm;
 452     }
 453     if (iscale == 1) {
 454         if (lower) {
 455             i__1 = *n;
 456             for (j = 1; j <= i__1; ++j) {
 457                 i__2 = *n - j + 1;
 458                 sscal_(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
 459 /* L10: */
 460             }
 461         } else {
 462             i__1 = *n;
 463             for (j = 1; j <= i__1; ++j) {
 464                 sscal_(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
 465 /* L20: */
 466             }
 467         }
 468         if (*abstol > 0.f) {
 469             abstll = *abstol * sigma;
 470         }
 471         if (valeig) {
 472             vll = *vl * sigma;
 473             vuu = *vu * sigma;
 474         }
 475     }
 476 /*     Initialize indices into workspaces.  Note: The IWORK indices are */
 477 /*     used only if SSTERF or SSTEMR fail. */
 478 /*     WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the */
 479 /*     elementary reflectors used in SSYTRD. */
 480     indtau = 1;
 481 /*     WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries. */
 482     indd = indtau + *n;
 483 /*     WORK(INDE:INDE+N-1) stores the off-diagonal entries of the */
 484 /*     tridiagonal matrix from SSYTRD. */
 485     inde = indd + *n;
 486 /*     WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over */
 487 /*     -written by SSTEMR (the SSTERF path copies the diagonal to W). */
 488     inddd = inde + *n;
 489 /*     WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over */
 490 /*     -written while computing the eigenvalues in SSTERF and SSTEMR. */
 491     indee = inddd + *n;
 492 /*     INDWK is the starting offset of the left-over workspace, and */
 493 /*     LLWORK is the remaining workspace size. */
 494     indwk = indee + *n;
 495     llwork = *lwork - indwk + 1;
 496 /*     IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and */
 497 /*     stores the block indices of each of the M<=N eigenvalues. */
 498     indibl = 1;
 499 /*     IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and */
 500 /*     stores the starting and finishing indices of each block. */
 501     indisp = indibl + *n;
 502 /*     IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors */
 503 /*     that corresponding to eigenvectors that fail to converge in */
 504 /*     SSTEIN.  This information is discarded; if any fail, the driver */
 505 /*     returns INFO > 0. */
 506     indifl = indisp + *n;
 507 /*     INDIWO is the offset of the remaining integer workspace. */
 508     indiwo = indisp + *n;
 509
 510 /*     Call SSYTRD to reduce symmetric matrix to tridiagonal form. */
 511
 512     ssytrd_(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
 513             indtau], &work[indwk], &llwork, &iinfo);
 514
 515 /*     If all eigenvalues are desired */
 516 /*     then call SSTERF or SSTEMR and SORMTR. */
 517
 518     test = FALSE_;
 519     if (indeig) {
 520         if (*il == 1 && *iu == *n) {
 521             test = TRUE_;
 522         }
 523     }
 524     if ((alleig || test) && ieeeok == 1) {
 525         if (! wantz) {
 526             scopy_(n, &work[indd], &c__1, &w[1], &c__1);
 527             i__1 = *n - 1;
 528             scopy_(&i__1, &work[inde], &c__1, &work[indee], &c__1);
 529             ssterf_(n, &w[1], &work[indee], info);
 530         } else {
 531             i__1 = *n - 1;
 532             scopy_(&i__1, &work[inde], &c__1, &work[indee], &c__1);
 533             scopy_(n, &work[indd], &c__1, &work[inddd], &c__1);
 534
 535             if (*abstol <= *n * 2.f * eps) {
 536                 tryrac = TRUE_;
 537             } else {
 538                 tryrac = FALSE_;
 539             }
 540             sstemr_(jobz, "A", n, &work[inddd], &work[indee], vl, vu, il, iu,
 541                     m, &w[1], &z__[z_offset], ldz, n, &isuppz[1], &tryrac, &
 542                     work[indwk], lwork, &iwork[1], liwork, info);
 543
 544
 545
 546 /*        Apply orthogonal matrix used in reduction to tridiagonal */
 547 /*        form to eigenvectors returned by SSTEIN. */
 548
 549             if (wantz && *info == 0) {
 550                 indwkn = inde;
 551                 llwrkn = *lwork - indwkn + 1;
 552                 sormtr_("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau]
 553 , &z__[z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
 554             }
 555         }
 556
 557
 558         if (*info == 0) {
 559 /*           Everything worked.  Skip SSTEBZ/SSTEIN.  IWORK(:) are */
 560 /*           undefined. */
 561             *m = *n;
 562             goto L30;
 563         }
 564         *info = 0;
 565     }
 566
 567 /*     Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN. */
 568 /*     Also call SSTEBZ and SSTEIN if SSTEMR fails. */
 569
 570     if (wantz) {
 571         *(unsigned char *)order = 'B';
 572     } else {
 573         *(unsigned char *)order = 'E';
 574     }
 575     sstebz_(range, order, n, &vll, &vuu, il, iu, &abstll, &work[indd], &work[
 576             inde], m, &nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[
 577             indwk], &iwork[indiwo], info);
 578
 579     if (wantz) {
 580         sstein_(n, &work[indd], &work[inde], m, &w[1], &iwork[indibl], &iwork[
 581                 indisp], &z__[z_offset], ldz, &work[indwk], &iwork[indiwo], &
 582                 iwork[indifl], info);
 583
 584 /*        Apply orthogonal matrix used in reduction to tridiagonal */
 585 /*        form to eigenvectors returned by SSTEIN. */
 586
 587         indwkn = inde;
 588         llwrkn = *lwork - indwkn + 1;
 589         sormtr_("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau], &z__[
 590                 z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
 591     }
 592
 593 /*     If matrix was scaled, then rescale eigenvalues appropriately. */
 594
 595 /*  Jump here if SSTEMR/SSTEIN succeeded. */
 596 L30:
 597     if (iscale == 1) {
 598         if (*info == 0) {
 599             imax = *m;
 600         } else {
 601             imax = *info - 1;
 602         }
 603         r__1 = 1.f / sigma;
 604         sscal_(&imax, &r__1, &w[1], &c__1);
 605     }
 606
 607 /*     If eigenvalues are not in order, then sort them, along with */
 608 /*     eigenvectors.  Note: We do not sort the IFAIL portion of IWORK. */
 609 /*     It may not be initialized (if SSTEMR/SSTEIN succeeded), and we do */
 610 /*     not return this detailed information to the user. */
 611
 612     if (wantz) {
 613         i__1 = *m - 1;
 614         for (j = 1; j <= i__1; ++j) {
 615             i__ = 0;
 616             tmp1 = w[j];
 617             i__2 = *m;
 618             for (jj = j + 1; jj <= i__2; ++jj) {
 619                 if (w[jj] < tmp1) {
 620                     i__ = jj;
 621                     tmp1 = w[jj];
 622                 }
 623 /* L40: */
 624             }
 625
 626             if (i__ != 0) {
 627                 w[i__] = w[j];
 628                 w[j] = tmp1;
 629                 sswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
 630                          &c__1);
 631             }
 632 /* L50: */
 633         }
 634     }
 635
 636 /*     Set WORK(1) to optimal workspace size. */
 637
 638     work[1] = (real) lwkopt;
 639     iwork[1] = liwmin;
 640
 641     return 0;
 642
 643 /*     End of SSYEVR */
 644
 645 } /* ssyevr_ */