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