template_lapack_stemr.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027  
00028  /* This file belongs to the template_lapack part of the Ergo source 
00029   * code. The source files in the template_lapack directory are modified
00030   * versions of files originally distributed as CLAPACK, see the
00031   * Copyright/license notice in the file template_lapack/COPYING.
00032   */
00033  
00034 
00035 #ifndef TEMPLATE_LAPACK_STEMR_HEADER
00036 #define TEMPLATE_LAPACK_STEMR_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_stemr(const char *jobz, const char *range, const integer *n, Treal *
00040         d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il, 
00041         const integer *iu, integer *m, Treal *w, Treal *z__, const integer *ldz, 
00042          const integer *nzc, integer *isuppz, logical *tryrac, Treal *work, 
00043         integer *lwork, integer *iwork, integer *liwork, integer *info)
00044 {
00045     /* System generated locals */
00046     integer z_dim1, z_offset, i__1, i__2;
00047     Treal d__1, d__2;
00048 
00049     /* Builtin functions */
00050 
00051     /* Local variables */
00052     integer i__, j;
00053     Treal r1, r2;
00054     integer jj;
00055     Treal cs;
00056     integer in;
00057     Treal sn, wl, wu;
00058     integer iil, iiu;
00059     Treal eps, tmp;
00060     integer indd, iend, jblk, wend;
00061     Treal rmin, rmax;
00062     integer itmp;
00063     Treal tnrm;
00064     integer inde2, itmp2;
00065     Treal rtol1, rtol2;
00066     Treal scale;
00067     integer indgp;
00068     integer iinfo, iindw, ilast;
00069     integer lwmin;
00070     logical wantz;
00071     logical alleig;
00072     integer ibegin;
00073     logical indeig;
00074     integer iindbl;
00075     logical valeig;
00076     integer wbegin;
00077     Treal safmin;
00078     Treal bignum;
00079     integer inderr, iindwk, indgrs, offset;
00080     Treal thresh;
00081     integer iinspl, ifirst, indwrk, liwmin, nzcmin;
00082     Treal pivmin;
00083     integer nsplit;
00084     Treal smlnum;
00085     logical lquery, zquery;
00086 
00087 
00088 /*  -- LAPACK computational routine (version 3.2) -- */
00089 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00090 /*     November 2006 */
00091 
00092 /*     .. Scalar Arguments .. */
00093 /*     .. */
00094 /*     .. Array Arguments .. */
00095 /*     .. */
00096 
00097 /*  Purpose */
00098 /*  ======= */
00099 
00100 /*  DSTEMR computes selected eigenvalues and, optionally, eigenvectors */
00101 /*  of a real symmetric tridiagonal matrix T. Any such unreduced matrix has */
00102 /*  a well defined set of pairwise different real eigenvalues, the corresponding */
00103 /*  real eigenvectors are pairwise orthogonal. */
00104 
00105 /*  The spectrum may be computed either completely or partially by specifying */
00106 /*  either an interval (VL,VU] or a range of indices IL:IU for the desired */
00107 /*  eigenvalues. */
00108 
00109 /*  Depending on the number of desired eigenvalues, these are computed either */
00110 /*  by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are */
00111 /*  computed by the use of various suitable L D L^T factorizations near clusters */
00112 /*  of close eigenvalues (referred to as RRRs, Relatively Robust */
00113 /*  Representations). An informal sketch of the algorithm follows. */
00114 
00115 /*  For each unreduced block (submatrix) of T, */
00116 /*     (a) Compute T - sigma I  = L D L^T, so that L and D */
00117 /*         define all the wanted eigenvalues to high relative accuracy. */
00118 /*         This means that small relative changes in the entries of D and L */
00119 /*         cause only small relative changes in the eigenvalues and */
00120 /*         eigenvectors. The standard (unfactored) representation of the */
00121 /*         tridiagonal matrix T does not have this property in general. */
00122 /*     (b) Compute the eigenvalues to suitable accuracy. */
00123 /*         If the eigenvectors are desired, the algorithm attains full */
00124 /*         accuracy of the computed eigenvalues only right before */
00125 /*         the corresponding vectors have to be computed, see steps c) and d). */
00126 /*     (c) For each cluster of close eigenvalues, select a new */
00127 /*         shift close to the cluster, find a new factorization, and refine */
00128 /*         the shifted eigenvalues to suitable accuracy. */
00129 /*     (d) For each eigenvalue with a large enough relative separation compute */
00130 /*         the corresponding eigenvector by forming a rank revealing twisted */
00131 /*         factorization. Go back to (c) for any clusters that remain. */
00132 
00133 /*  For more details, see: */
00134 /*  - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations */
00135 /*    to compute orthogonal eigenvectors of symmetric tridiagonal matrices," */
00136 /*    Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. */
00137 /*  - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and */
00138 /*    Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, */
00139 /*    2004.  Also LAPACK Working Note 154. */
00140 /*  - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric */
00141 /*    tridiagonal eigenvalue/eigenvector problem", */
00142 /*    Computer Science Division Technical Report No. UCB/CSD-97-971, */
00143 /*    UC Berkeley, May 1997. */
00144 
00145 /*  Notes: */
00146 /*  1.DSTEMR works only on machines which follow IEEE-754 */
00147 /*  floating-point standard in their handling of infinities and NaNs. */
00148 /*  This permits the use of efficient inner loops avoiding a check for */
00149 /*  zero divisors. */
00150 
00151 /*  Arguments */
00152 /*  ========= */
00153 
00154 /*  JOBZ    (input) CHARACTER*1 */
00155 /*          = 'N':  Compute eigenvalues only; */
00156 /*          = 'V':  Compute eigenvalues and eigenvectors. */
00157 
00158 /*  RANGE   (input) CHARACTER*1 */
00159 /*          = 'A': all eigenvalues will be found. */
00160 /*          = 'V': all eigenvalues in the half-open interval (VL,VU] */
00161 /*                 will be found. */
00162 /*          = 'I': the IL-th through IU-th eigenvalues will be found. */
00163 
00164 /*  N       (input) INTEGER */
00165 /*          The order of the matrix.  N >= 0. */
00166 
00167 /*  D       (input/output) DOUBLE PRECISION array, dimension (N) */
00168 /*          On entry, the N diagonal elements of the tridiagonal matrix */
00169 /*          T. On exit, D is overwritten. */
00170 
00171 /*  E       (input/output) DOUBLE PRECISION array, dimension (N) */
00172 /*          On entry, the (N-1) subdiagonal elements of the tridiagonal */
00173 /*          matrix T in elements 1 to N-1 of E. E(N) need not be set on */
00174 /*          input, but is used internally as workspace. */
00175 /*          On exit, E is overwritten. */
00176 
00177 /*  VL      (input) DOUBLE PRECISION */
00178 /*  VU      (input) DOUBLE PRECISION */
00179 /*          If RANGE='V', the lower and upper bounds of the interval to */
00180 /*          be searched for eigenvalues. VL < VU. */
00181 /*          Not referenced if RANGE = 'A' or 'I'. */
00182 
00183 /*  IL      (input) INTEGER */
00184 /*  IU      (input) INTEGER */
00185 /*          If RANGE='I', the indices (in ascending order) of the */
00186 /*          smallest and largest eigenvalues to be returned. */
00187 /*          1 <= IL <= IU <= N, if N > 0. */
00188 /*          Not referenced if RANGE = 'A' or 'V'. */
00189 
00190 /*  M       (output) INTEGER */
00191 /*          The total number of eigenvalues found.  0 <= M <= N. */
00192 /*          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
00193 
00194 /*  W       (output) DOUBLE PRECISION array, dimension (N) */
00195 /*          The first M elements contain the selected eigenvalues in */
00196 /*          ascending order. */
00197 
00198 /*  Z       (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
00199 /*          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z */
00200 /*          contain the orthonormal eigenvectors of the matrix T */
00201 /*          corresponding to the selected eigenvalues, with the i-th */
00202 /*          column of Z holding the eigenvector associated with W(i). */
00203 /*          If JOBZ = 'N', then Z is not referenced. */
00204 /*          Note: the user must ensure that at least max(1,M) columns are */
00205 /*          supplied in the array Z; if RANGE = 'V', the exact value of M */
00206 /*          is not known in advance and can be computed with a workspace */
00207 /*          query by setting NZC = -1, see below. */
00208 
00209 /*  LDZ     (input) INTEGER */
00210 /*          The leading dimension of the array Z.  LDZ >= 1, and if */
00211 /*          JOBZ = 'V', then LDZ >= max(1,N). */
00212 
00213 /*  NZC     (input) INTEGER */
00214 /*          The number of eigenvectors to be held in the array Z. */
00215 /*          If RANGE = 'A', then NZC >= max(1,N). */
00216 /*          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. */
00217 /*          If RANGE = 'I', then NZC >= IU-IL+1. */
00218 /*          If NZC = -1, then a workspace query is assumed; the */
00219 /*          routine calculates the number of columns of the array Z that */
00220 /*          are needed to hold the eigenvectors. */
00221 /*          This value is returned as the first entry of the Z array, and */
00222 /*          no error message related to NZC is issued by XERBLA. */
00223 
00224 /*  ISUPPZ  (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) */
00225 /*          The support of the eigenvectors in Z, i.e., the indices */
00226 /*          indicating the nonzero elements in Z. The i-th computed eigenvector */
00227 /*          is nonzero only in elements ISUPPZ( 2*i-1 ) through */
00228 /*          ISUPPZ( 2*i ). This is relevant in the case when the matrix */
00229 /*          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. */
00230 
00231 /*  TRYRAC  (input/output) LOGICAL */
00232 /*          If TRYRAC.EQ..TRUE., indicates that the code should check whether */
00233 /*          the tridiagonal matrix defines its eigenvalues to high relative */
00234 /*          accuracy.  If so, the code uses relative-accuracy preserving */
00235 /*          algorithms that might be (a bit) slower depending on the matrix. */
00236 /*          If the matrix does not define its eigenvalues to high relative */
00237 /*          accuracy, the code can uses possibly faster algorithms. */
00238 /*          If TRYRAC.EQ..FALSE., the code is not required to guarantee */
00239 /*          relatively accurate eigenvalues and can use the fastest possible */
00240 /*          techniques. */
00241 /*          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix */
00242 /*          does not define its eigenvalues to high relative accuracy. */
00243 
00244 /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
00245 /*          On exit, if INFO = 0, WORK(1) returns the optimal */
00246 /*          (and minimal) LWORK. */
00247 
00248 /*  LWORK   (input) INTEGER */
00249 /*          The dimension of the array WORK. LWORK >= max(1,18*N) */
00250 /*          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. */
00251 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00252 /*          only calculates the optimal size of the WORK array, returns */
00253 /*          this value as the first entry of the WORK array, and no error */
00254 /*          message related to LWORK is issued by XERBLA. */
00255 
00256 /*  IWORK   (workspace/output) INTEGER array, dimension (LIWORK) */
00257 /*          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
00258 
00259 /*  LIWORK  (input) INTEGER */
00260 /*          The dimension of the array IWORK.  LIWORK >= max(1,10*N) */
00261 /*          if the eigenvectors are desired, and LIWORK >= max(1,8*N) */
00262 /*          if only the eigenvalues are to be computed. */
00263 /*          If LIWORK = -1, then a workspace query is assumed; the */
00264 /*          routine only calculates the optimal size of the IWORK array, */
00265 /*          returns this value as the first entry of the IWORK array, and */
00266 /*          no error message related to LIWORK is issued by XERBLA. */
00267 
00268 /*  INFO    (output) INTEGER */
00269 /*          On exit, INFO */
00270 /*          = 0:  successful exit */
00271 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00272 /*          > 0:  if INFO = 1X, internal error in DLARRE, */
00273 /*                if INFO = 2X, internal error in DLARRV. */
00274 /*                Here, the digit X = ABS( IINFO ) < 10, where IINFO is */
00275 /*                the nonzero error code returned by DLARRE or */
00276 /*                DLARRV, respectively. */
00277 
00278 
00279 /*  Further Details */
00280 /*  =============== */
00281 
00282 /*  Based on contributions by */
00283 /*     Beresford Parlett, University of California, Berkeley, USA */
00284 /*     Jim Demmel, University of California, Berkeley, USA */
00285 /*     Inderjit Dhillon, University of Texas, Austin, USA */
00286 /*     Osni Marques, LBNL/NERSC, USA */
00287 /*     Christof Voemel, University of California, Berkeley, USA */
00288 
00289 /*  ===================================================================== */
00290 
00291 /*     .. Parameters .. */
00292 /*     .. */
00293 /*     .. Local Scalars .. */
00294 /*     .. */
00295 /*     .. */
00296 /*     .. External Functions .. */
00297 /*     .. */
00298 /*     .. External Subroutines .. */
00299 /*     .. */
00300 /*     .. Intrinsic Functions .. */
00301 /*     .. */
00302 /*     .. Executable Statements .. */
00303 
00304 /*     Test the input parameters. */
00305 
00306     /* Parameter adjustments */
00307     /* Table of constant values */
00308     integer c__1 = 1;
00309     Treal c_b18 = .001;
00310 
00311     --d__;
00312     --e;
00313     --w;
00314     z_dim1 = *ldz;
00315     z_offset = 1 + z_dim1;
00316     z__ -= z_offset;
00317     --isuppz;
00318     --work;
00319     --iwork;
00320 
00321     /* Function Body */
00322     wantz = template_blas_lsame(jobz, "V");
00323     alleig = template_blas_lsame(range, "A");
00324     valeig = template_blas_lsame(range, "V");
00325     indeig = template_blas_lsame(range, "I");
00326 
00327     lquery = *lwork == -1 || *liwork == -1;
00328     zquery = *nzc == -1;
00329 /*     DSTEMR needs WORK of size 6*N, IWORK of size 3*N. */
00330 /*     In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N. */
00331 /*     Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N. */
00332     if (wantz) {
00333         lwmin = *n * 18;
00334         liwmin = *n * 10;
00335     } else {
00336 /*        need less workspace if only the eigenvalues are wanted */
00337         lwmin = *n * 12;
00338         liwmin = *n << 3;
00339     }
00340     wl = 0.;
00341     wu = 0.;
00342     iil = 0;
00343     iiu = 0;
00344     if (valeig) {
00345 /*        We do not reference VL, VU in the cases RANGE = 'I','A' */
00346 /*        The interval (WL, WU] contains all the wanted eigenvalues. */
00347 /*        It is either given by the user or computed in DLARRE. */
00348         wl = *vl;
00349         wu = *vu;
00350     } else if (indeig) {
00351 /*        We do not reference IL, IU in the cases RANGE = 'V','A' */
00352         iil = *il;
00353         iiu = *iu;
00354     }
00355 
00356     *info = 0;
00357     if (! (wantz || template_blas_lsame(jobz, "N"))) {
00358         *info = -1;
00359     } else if (! (alleig || valeig || indeig)) {
00360         *info = -2;
00361     } else if (*n < 0) {
00362         *info = -3;
00363     } else if (valeig && *n > 0 && wu <= wl) {
00364         *info = -7;
00365     } else if (indeig && (iil < 1 || iil > *n)) {
00366         *info = -8;
00367     } else if (indeig && (iiu < iil || iiu > *n)) {
00368         *info = -9;
00369     } else if (*ldz < 1 || ( wantz && *ldz < *n ) ) {
00370         *info = -13;
00371     } else if (*lwork < lwmin && ! lquery) {
00372         *info = -17;
00373     } else if (*liwork < liwmin && ! lquery) {
00374         *info = -19;
00375     }
00376 
00377 /*     Get machine constants. */
00378 
00379     safmin = template_lapack_lamch("Safe minimum", (Treal)0);
00380     eps = template_lapack_lamch("Precision", (Treal)0);
00381     smlnum = safmin / eps;
00382     bignum = 1. / smlnum;
00383     rmin = template_blas_sqrt(smlnum);
00384 /* Computing MIN */
00385     d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
00386     rmax = minMACRO(d__1,d__2);
00387 
00388     if (*info == 0) {
00389         work[1] = (Treal) lwmin;
00390         iwork[1] = liwmin;
00391 
00392         if (wantz && alleig) {
00393             nzcmin = *n;
00394         } else if (wantz && valeig) {
00395             template_lapack_larrc("T", n, vl, vu, &d__[1], &e[1], &safmin, &nzcmin, &itmp, &
00396                     itmp2, info);
00397         } else if (wantz && indeig) {
00398             nzcmin = iiu - iil + 1;
00399         } else {
00400 /*           WANTZ .EQ. FALSE. */
00401             nzcmin = 0;
00402         }
00403         if (zquery && *info == 0) {
00404             z__[z_dim1 + 1] = (Treal) nzcmin;
00405         } else if (*nzc < nzcmin && ! zquery) {
00406             *info = -14;
00407         }
00408     }
00409     if (*info != 0) {
00410 
00411         i__1 = -(*info);
00412         template_blas_erbla("DSTEMR", &i__1);
00413 
00414         return 0;
00415     } else if (lquery || zquery) {
00416         return 0;
00417     }
00418 
00419 /*     Handle N = 0, 1, and 2 cases immediately */
00420 
00421     *m = 0;
00422     if (*n == 0) {
00423         return 0;
00424     }
00425 
00426     if (*n == 1) {
00427         if (alleig || indeig) {
00428             *m = 1;
00429             w[1] = d__[1];
00430         } else {
00431             if (wl < d__[1] && wu >= d__[1]) {
00432                 *m = 1;
00433                 w[1] = d__[1];
00434             }
00435         }
00436         if (wantz && ! zquery) {
00437             z__[z_dim1 + 1] = 1.;
00438             isuppz[1] = 1;
00439             isuppz[2] = 1;
00440         }
00441         return 0;
00442     }
00443 
00444     if (*n == 2) {
00445         if (! wantz) {
00446             template_lapack_lae2(&d__[1], &e[1], &d__[2], &r1, &r2);
00447         } else if (wantz && ! zquery) {
00448             template_lapack_laev2(&d__[1], &e[1], &d__[2], &r1, &r2, &cs, &sn);
00449         }
00450         if (alleig || ( valeig && r2 > wl && r2 <= wu ) || ( indeig && iil == 1 ) ) {
00451             ++(*m);
00452             w[*m] = r2;
00453             if (wantz && ! zquery) {
00454                 z__[*m * z_dim1 + 1] = -sn;
00455                 z__[*m * z_dim1 + 2] = cs;
00456 /*              Note: At most one of SN and CS can be zero. */
00457                 if (sn != 0.) {
00458                     if (cs != 0.) {
00459                         isuppz[(*m << 1) - 1] = 1;
00460                         isuppz[(*m << 1) - 1] = 2;
00461                     } else {
00462                         isuppz[(*m << 1) - 1] = 1;
00463                         isuppz[(*m << 1) - 1] = 1;
00464                     }
00465                 } else {
00466                     isuppz[(*m << 1) - 1] = 2;
00467                     isuppz[*m * 2] = 2;
00468                 }
00469             }
00470         }
00471         if (alleig || ( valeig && r1 > wl && r1 <= wu ) || ( indeig && iiu == 2 ) ) {
00472             ++(*m);
00473             w[*m] = r1;
00474             if (wantz && ! zquery) {
00475                 z__[*m * z_dim1 + 1] = cs;
00476                 z__[*m * z_dim1 + 2] = sn;
00477 /*              Note: At most one of SN and CS can be zero. */
00478                 if (sn != 0.) {
00479                     if (cs != 0.) {
00480                         isuppz[(*m << 1) - 1] = 1;
00481                         isuppz[(*m << 1) - 1] = 2;
00482                     } else {
00483                         isuppz[(*m << 1) - 1] = 1;
00484                         isuppz[(*m << 1) - 1] = 1;
00485                     }
00486                 } else {
00487                     isuppz[(*m << 1) - 1] = 2;
00488                     isuppz[*m * 2] = 2;
00489                 }
00490             }
00491         }
00492         return 0;
00493     }
00494 /*     Continue with general N */
00495     indgrs = 1;
00496     inderr = (*n << 1) + 1;
00497     indgp = *n * 3 + 1;
00498     indd = (*n << 2) + 1;
00499     inde2 = *n * 5 + 1;
00500     indwrk = *n * 6 + 1;
00501 
00502     iinspl = 1;
00503     iindbl = *n + 1;
00504     iindw = (*n << 1) + 1;
00505     iindwk = *n * 3 + 1;
00506 
00507 /*     Scale matrix to allowable range, if necessary. */
00508 /*     The allowable range is related to the PIVMIN parameter; see the */
00509 /*     comments in DLARRD.  The preference for scaling small values */
00510 /*     up is heuristic; we expect users' matrices not to be close to the */
00511 /*     RMAX threshold. */
00512 
00513     scale = 1.;
00514     tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
00515     if (tnrm > 0. && tnrm < rmin) {
00516         scale = rmin / tnrm;
00517     } else if (tnrm > rmax) {
00518         scale = rmax / tnrm;
00519     }
00520     if (scale != 1.) {
00521         template_blas_scal(n, &scale, &d__[1], &c__1);
00522         i__1 = *n - 1;
00523         template_blas_scal(&i__1, &scale, &e[1], &c__1);
00524         tnrm *= scale;
00525         if (valeig) {
00526 /*           If eigenvalues in interval have to be found, */
00527 /*           scale (WL, WU] accordingly */
00528             wl *= scale;
00529             wu *= scale;
00530         }
00531     }
00532 
00533 /*     Compute the desired eigenvalues of the tridiagonal after splitting */
00534 /*     into smaller subblocks if the corresponding off-diagonal elements */
00535 /*     are small */
00536 /*     THRESH is the splitting parameter for DLARRE */
00537 /*     A negative THRESH forces the old splitting criterion based on the */
00538 /*     size of the off-diagonal. A positive THRESH switches to splitting */
00539 /*     which preserves relative accuracy. */
00540 
00541     if (*tryrac) {
00542 /*        Test whether the matrix warrants the more expensive relative approach. */
00543         template_lapack_larrr(n, &d__[1], &e[1], &iinfo);
00544     } else {
00545 /*        The user does not care about relative accurately eigenvalues */
00546         iinfo = -1;
00547     }
00548 /*     Set the splitting criterion */
00549     if (iinfo == 0) {
00550         thresh = eps;
00551     } else {
00552         thresh = -eps;
00553 /*        relative accuracy is desired but T does not guarantee it */
00554         *tryrac = FALSE_;
00555     }
00556 
00557     if (*tryrac) {
00558 /*        Copy original diagonal, needed to guarantee relative accuracy */
00559         template_blas_copy(n, &d__[1], &c__1, &work[indd], &c__1);
00560     }
00561 /*     Store the squares of the offdiagonal values of T */
00562     i__1 = *n - 1;
00563     for (j = 1; j <= i__1; ++j) {
00564 /* Computing 2nd power */
00565         d__1 = e[j];
00566         work[inde2 + j - 1] = d__1 * d__1;
00567 /* L5: */
00568     }
00569 /*     Set the tolerance parameters for bisection */
00570     if (! wantz) {
00571 /*        DLARRE computes the eigenvalues to full precision. */
00572         rtol1 = eps * 4.;
00573         rtol2 = eps * 4.;
00574     } else {
00575 /*        DLARRE computes the eigenvalues to less than full precision. */
00576 /*        DLARRV will refine the eigenvalue approximations, and we can */
00577 /*        need less accurate initial bisection in DLARRE. */
00578 /*        Note: these settings do only affect the subset case and DLARRE */
00579         rtol1 = template_blas_sqrt(eps);
00580 /* Computing MAX */
00581         d__1 = template_blas_sqrt(eps) * .005, d__2 = eps * 4.;
00582         rtol2 = maxMACRO(d__1,d__2);
00583     }
00584     template_lapack_larre(range, n, &wl, &wu, &iil, &iiu, &d__[1], &e[1], &work[inde2], &
00585             rtol1, &rtol2, &thresh, &nsplit, &iwork[iinspl], m, &w[1], &work[
00586             inderr], &work[indgp], &iwork[iindbl], &iwork[iindw], &work[
00587             indgrs], &pivmin, &work[indwrk], &iwork[iindwk], &iinfo);
00588     if (iinfo != 0) {
00589         *info = absMACRO(iinfo) + 10;
00590         return 0;
00591     }
00592 /*     Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired */
00593 /*     part of the spectrum. All desired eigenvalues are contained in */
00594 /*     (WL,WU] */
00595     if (wantz) {
00596 
00597 /*        Compute the desired eigenvectors corresponding to the computed */
00598 /*        eigenvalues */
00599 
00600         template_lapack_larrv(n, &wl, &wu, &d__[1], &e[1], &pivmin, &iwork[iinspl], m, &
00601                 c__1, m, &c_b18, &rtol1, &rtol2, &w[1], &work[inderr], &work[
00602                 indgp], &iwork[iindbl], &iwork[iindw], &work[indgrs], &z__[
00603                 z_offset], ldz, &isuppz[1], &work[indwrk], &iwork[iindwk], &
00604                 iinfo);
00605         if (iinfo != 0) {
00606             *info = absMACRO(iinfo) + 20;
00607             return 0;
00608         }
00609     } else {
00610 /*        DLARRE computes eigenvalues of the (shifted) root representation */
00611 /*        DLARRV returns the eigenvalues of the unshifted matrix. */
00612 /*        However, if the eigenvectors are not desired by the user, we need */
00613 /*        to apply the corresponding shifts from DLARRE to obtain the */
00614 /*        eigenvalues of the original matrix. */
00615         i__1 = *m;
00616         for (j = 1; j <= i__1; ++j) {
00617             itmp = iwork[iindbl + j - 1];
00618             w[j] += e[iwork[iinspl + itmp - 1]];
00619 /* L20: */
00620         }
00621     }
00622 
00623     if (*tryrac) {
00624 /*        Refine computed eigenvalues so that they are relatively accurate */
00625 /*        with respect to the original matrix T. */
00626         ibegin = 1;
00627         wbegin = 1;
00628         i__1 = iwork[iindbl + *m - 1];
00629         for (jblk = 1; jblk <= i__1; ++jblk) {
00630             iend = iwork[iinspl + jblk - 1];
00631             in = iend - ibegin + 1;
00632             wend = wbegin - 1;
00633 /*           check if any eigenvalues have to be refined in this block */
00634 L36:
00635             if (wend < *m) {
00636                 if (iwork[iindbl + wend] == jblk) {
00637                     ++wend;
00638                     goto L36;
00639                 }
00640             }
00641             if (wend < wbegin) {
00642                 ibegin = iend + 1;
00643                 goto L39;
00644             }
00645             offset = iwork[iindw + wbegin - 1] - 1;
00646             ifirst = iwork[iindw + wbegin - 1];
00647             ilast = iwork[iindw + wend - 1];
00648             rtol2 = eps * 4.;
00649             template_lapack_larrj(&in, &work[indd + ibegin - 1], &work[inde2 + ibegin - 1], 
00650                     &ifirst, &ilast, &rtol2, &offset, &w[wbegin], &work[
00651                     inderr + wbegin - 1], &work[indwrk], &iwork[iindwk], &
00652                     pivmin, &tnrm, &iinfo);
00653             ibegin = iend + 1;
00654             wbegin = wend + 1;
00655 L39:
00656             ;
00657         }
00658     }
00659 
00660 /*     If matrix was scaled, then rescale eigenvalues appropriately. */
00661 
00662     if (scale != 1.) {
00663         d__1 = 1. / scale;
00664         template_blas_scal(m, &d__1, &w[1], &c__1);
00665     }
00666 
00667 /*     If eigenvalues are not in increasing order, then sort them, */
00668 /*     possibly along with eigenvectors. */
00669 
00670     if (nsplit > 1) {
00671         if (! wantz) {
00672             template_lapack_lasrt("I", m, &w[1], &iinfo);
00673             if (iinfo != 0) {
00674                 *info = 3;
00675                 return 0;
00676             }
00677         } else {
00678             i__1 = *m - 1;
00679             for (j = 1; j <= i__1; ++j) {
00680                 i__ = 0;
00681                 tmp = w[j];
00682                 i__2 = *m;
00683                 for (jj = j + 1; jj <= i__2; ++jj) {
00684                     if (w[jj] < tmp) {
00685                         i__ = jj;
00686                         tmp = w[jj];
00687                     }
00688 /* L50: */
00689                 }
00690                 if (i__ != 0) {
00691                     w[i__] = w[j];
00692                     w[j] = tmp;
00693                     if (wantz) {
00694                         template_blas_swap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * 
00695                                 z_dim1 + 1], &c__1);
00696                         itmp = isuppz[(i__ << 1) - 1];
00697                         isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
00698                         isuppz[(j << 1) - 1] = itmp;
00699                         itmp = isuppz[i__ * 2];
00700                         isuppz[i__ * 2] = isuppz[j * 2];
00701                         isuppz[j * 2] = itmp;
00702                     }
00703                 }
00704 /* L60: */
00705             }
00706         }
00707     }
00708 
00709 
00710     work[1] = (Treal) lwmin;
00711     iwork[1] = liwmin;
00712     return 0;
00713 
00714 /*     End of DSTEMR */
00715 
00716 } /* dstemr_ */
00717 
00718 
00719 #endif

Generated on Wed Nov 21 09:32:37 2012 for ergo by  doxygen 1.4.7