template_lapack_larrd.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_LARRD_HEADER
00036 #define TEMPLATE_LAPACK_LARRD_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_larrd(const char *range, const char *order, const integer *n, Treal 
00040         *vl, Treal *vu, integer *il, integer *iu, Treal *gers, 
00041         Treal *reltol, Treal *d__, Treal *e, Treal *e2, 
00042         Treal *pivmin, integer *nsplit, integer *isplit, integer *m, 
00043         Treal *w, Treal *werr, Treal *wl, Treal *wu, 
00044         integer *iblock, integer *indexw, Treal *work, integer *iwork, 
00045         integer *info)
00046 {
00047     /* System generated locals */
00048     integer i__1, i__2, i__3;
00049     Treal d__1, d__2;
00050 
00051 
00052     /* Local variables */
00053     integer i__, j, ib, ie, je, nb;
00054     Treal gl;
00055     integer im, in;
00056     Treal gu;
00057     integer iw, jee;
00058     Treal eps;
00059     integer nwl;
00060     Treal wlu = 0;
00061     Treal wul = 0;
00062     integer nwu;
00063     Treal tmp1, tmp2;
00064     integer iend, jblk, ioff, iout, itmp1, itmp2, jdisc;
00065     integer iinfo;
00066     Treal atoli;
00067     integer iwoff, itmax;
00068     Treal wkill, rtoli, uflow, tnorm;
00069     integer ibegin;
00070     integer irange, idiscl, idumma[1];
00071     integer idiscu;
00072     logical ncnvrg, toofew;
00073 
00074 
00075 /*  -- LAPACK auxiliary routine (version 3.2.1)                        -- */
00076 /*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
00077 /*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
00078 /*  -- April 2009                                                      -- */
00079 
00080 /*     .. Scalar Arguments .. */
00081 /*     .. */
00082 /*     .. Array Arguments .. */
00083 /*     .. */
00084 
00085 /*  Purpose */
00086 /*  ======= */
00087 
00088 /*  DLARRD computes the eigenvalues of a symmetric tridiagonal */
00089 /*  matrix T to suitable accuracy. This is an auxiliary code to be */
00090 /*  called from DSTEMR. */
00091 /*  The user may ask for all eigenvalues, all eigenvalues */
00092 /*  in the half-open interval (VL, VU], or the IL-th through IU-th */
00093 /*  eigenvalues. */
00094 
00095 /*  To avoid overflow, the matrix must be scaled so that its */
00096 /*  largest element is no greater than overflow**(1/2) * */
00097 /*  underflow**(1/4) in absolute value, and for greatest */
00098 /*  accuracy, it should not be much smaller than that. */
00099 
00100 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
00101 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
00102 /*  University, July 21, 1966. */
00103 
00104 /*  Arguments */
00105 /*  ========= */
00106 
00107 /*  RANGE   (input) CHARACTER */
00108 /*          = 'A': ("All")   all eigenvalues will be found. */
00109 /*          = 'V': ("Value") all eigenvalues in the half-open interval */
00110 /*                           (VL, VU] will be found. */
00111 /*          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the */
00112 /*                           entire matrix) will be found. */
00113 
00114 /*  ORDER   (input) CHARACTER */
00115 /*          = 'B': ("By Block") the eigenvalues will be grouped by */
00116 /*                              split-off block (see IBLOCK, ISPLIT) and */
00117 /*                              ordered from smallest to largest within */
00118 /*                              the block. */
00119 /*          = 'E': ("Entire matrix") */
00120 /*                              the eigenvalues for the entire matrix */
00121 /*                              will be ordered from smallest to */
00122 /*                              largest. */
00123 
00124 /*  N       (input) INTEGER */
00125 /*          The order of the tridiagonal matrix T.  N >= 0. */
00126 
00127 /*  VL      (input) DOUBLE PRECISION */
00128 /*  VU      (input) DOUBLE PRECISION */
00129 /*          If RANGE='V', the lower and upper bounds of the interval to */
00130 /*          be searched for eigenvalues.  Eigenvalues less than or equal */
00131 /*          to VL, or greater than VU, will not be returned.  VL < VU. */
00132 /*          Not referenced if RANGE = 'A' or 'I'. */
00133 
00134 /*  IL      (input) INTEGER */
00135 /*  IU      (input) INTEGER */
00136 /*          If RANGE='I', the indices (in ascending order) of the */
00137 /*          smallest and largest eigenvalues to be returned. */
00138 /*          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
00139 /*          Not referenced if RANGE = 'A' or 'V'. */
00140 
00141 /*  GERS    (input) DOUBLE PRECISION array, dimension (2*N) */
00142 /*          The N Gerschgorin intervals (the i-th Gerschgorin interval */
00143 /*          is (GERS(2*i-1), GERS(2*i)). */
00144 
00145 /*  RELTOL  (input) DOUBLE PRECISION */
00146 /*          The minimum relative width of an interval.  When an interval */
00147 /*          is narrower than RELTOL times the larger (in */
00148 /*          magnitude) endpoint, then it is considered to be */
00149 /*          sufficiently small, i.e., converged.  Note: this should */
00150 /*          always be at least radix*machine epsilon. */
00151 
00152 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
00153 /*          The n diagonal elements of the tridiagonal matrix T. */
00154 
00155 /*  E       (input) DOUBLE PRECISION array, dimension (N-1) */
00156 /*          The (n-1) off-diagonal elements of the tridiagonal matrix T. */
00157 
00158 /*  E2      (input) DOUBLE PRECISION array, dimension (N-1) */
00159 /*          The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
00160 
00161 /*  PIVMIN  (input) DOUBLE PRECISION */
00162 /*          The minimum pivot allowed in the Sturm sequence for T. */
00163 
00164 /*  NSPLIT  (input) INTEGER */
00165 /*          The number of diagonal blocks in the matrix T. */
00166 /*          1 <= NSPLIT <= N. */
00167 
00168 /*  ISPLIT  (input) INTEGER array, dimension (N) */
00169 /*          The splitting points, at which T breaks up into submatrices. */
00170 /*          The first submatrix consists of rows/columns 1 to ISPLIT(1), */
00171 /*          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), */
00172 /*          etc., and the NSPLIT-th consists of rows/columns */
00173 /*          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. */
00174 /*          (Only the first NSPLIT elements will actually be used, but */
00175 /*          since the user cannot know a priori what value NSPLIT will */
00176 /*          have, N words must be reserved for ISPLIT.) */
00177 
00178 /*  M       (output) INTEGER */
00179 /*          The actual number of eigenvalues found. 0 <= M <= N. */
00180 /*          (See also the description of INFO=2,3.) */
00181 
00182 /*  W       (output) DOUBLE PRECISION array, dimension (N) */
00183 /*          On exit, the first M elements of W will contain the */
00184 /*          eigenvalue approximations. DLARRD computes an interval */
00185 /*          I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue */
00186 /*          approximation is given as the interval midpoint */
00187 /*          W(j)= ( a_j + b_j)/2. The corresponding error is bounded by */
00188 /*          WERR(j) = abs( a_j - b_j)/2 */
00189 
00190 /*  WERR    (output) DOUBLE PRECISION array, dimension (N) */
00191 /*          The error bound on the corresponding eigenvalue approximation */
00192 /*          in W. */
00193 
00194 /*  WL      (output) DOUBLE PRECISION */
00195 /*  WU      (output) DOUBLE PRECISION */
00196 /*          The interval (WL, WU] contains all the wanted eigenvalues. */
00197 /*          If RANGE='V', then WL=VL and WU=VU. */
00198 /*          If RANGE='A', then WL and WU are the global Gerschgorin bounds */
00199 /*                        on the spectrum. */
00200 /*          If RANGE='I', then WL and WU are computed by DLAEBZ from the */
00201 /*                        index range specified. */
00202 
00203 /*  IBLOCK  (output) INTEGER array, dimension (N) */
00204 /*          At each row/column j where E(j) is zero or small, the */
00205 /*          matrix T is considered to split into a block diagonal */
00206 /*          matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which */
00207 /*          block (from 1 to the number of blocks) the eigenvalue W(i) */
00208 /*          belongs.  (DLARRD may use the remaining N-M elements as */
00209 /*          workspace.) */
00210 
00211 /*  INDEXW  (output) INTEGER array, dimension (N) */
00212 /*          The indices of the eigenvalues within each block (submatrix); */
00213 /*          for example, INDEXW(i)= j and IBLOCK(i)=k imply that the */
00214 /*          i-th eigenvalue W(i) is the j-th eigenvalue in block k. */
00215 
00216 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N) */
00217 
00218 /*  IWORK   (workspace) INTEGER array, dimension (3*N) */
00219 
00220 /*  INFO    (output) INTEGER */
00221 /*          = 0:  successful exit */
00222 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00223 /*          > 0:  some or all of the eigenvalues failed to converge or */
00224 /*                were not computed: */
00225 /*                =1 or 3: Bisection failed to converge for some */
00226 /*                        eigenvalues; these eigenvalues are flagged by a */
00227 /*                        negative block number.  The effect is that the */
00228 /*                        eigenvalues may not be as accurate as the */
00229 /*                        absolute and relative tolerances.  This is */
00230 /*                        generally caused by unexpectedly inaccurate */
00231 /*                        arithmetic. */
00232 /*                =2 or 3: RANGE='I' only: Not all of the eigenvalues */
00233 /*                        IL:IU were found. */
00234 /*                        Effect: M < IU+1-IL */
00235 /*                        Cause:  non-monotonic arithmetic, causing the */
00236 /*                                Sturm sequence to be non-monotonic. */
00237 /*                        Cure:   recalculate, using RANGE='A', and pick */
00238 /*                                out eigenvalues IL:IU.  In some cases, */
00239 /*                                increasing the PARAMETER "FUDGE" may */
00240 /*                                make things work. */
00241 /*                = 4:    RANGE='I', and the Gershgorin interval */
00242 /*                        initially used was too small.  No eigenvalues */
00243 /*                        were computed. */
00244 /*                        Probable cause: your machine has sloppy */
00245 /*                                        floating-point arithmetic. */
00246 /*                        Cure: Increase the PARAMETER "FUDGE", */
00247 /*                              recompile, and try again. */
00248 
00249 /*  Internal Parameters */
00250 /*  =================== */
00251 
00252 /*  FUDGE   DOUBLE PRECISION, default = 2 */
00253 /*          A "fudge factor" to widen the Gershgorin intervals.  Ideally, */
00254 /*          a value of 1 should work, but on machines with sloppy */
00255 /*          arithmetic, this needs to be larger.  The default for */
00256 /*          publicly released versions should be large enough to handle */
00257 /*          the worst machine around.  Note that this has no effect */
00258 /*          on accuracy of the solution. */
00259 
00260 /*  Based on contributions by */
00261 /*     W. Kahan, University of California, Berkeley, USA */
00262 /*     Beresford Parlett, University of California, Berkeley, USA */
00263 /*     Jim Demmel, University of California, Berkeley, USA */
00264 /*     Inderjit Dhillon, University of Texas, Austin, USA */
00265 /*     Osni Marques, LBNL/NERSC, USA */
00266 /*     Christof Voemel, University of California, Berkeley, USA */
00267 
00268 /*  ===================================================================== */
00269 
00270 /*     .. Parameters .. */
00271 /*     .. */
00272 /*     .. Local Scalars .. */
00273 /*     .. */
00274 /*     .. Local Arrays .. */
00275 /*     .. */
00276 /*     .. External Functions .. */
00277 /*     .. */
00278 /*     .. External Subroutines .. */
00279 /*     .. */
00280 /*     .. Intrinsic Functions .. */
00281 /*     .. */
00282 /*     .. Executable Statements .. */
00283 
00284 /* Table of constant values */
00285 
00286     integer c__1 = 1;
00287     integer c_n1 = -1;
00288     integer c__3 = 3;
00289     integer c__2 = 2;
00290     integer c__0 = 0;
00291 
00292     /* Parameter adjustments */
00293     --iwork;
00294     --work;
00295     --indexw;
00296     --iblock;
00297     --werr;
00298     --w;
00299     --isplit;
00300     --e2;
00301     --e;
00302     --d__;
00303     --gers;
00304 
00305     /* Function Body */
00306     *info = 0;
00307 
00308 /*     Decode RANGE */
00309 
00310     if (template_blas_lsame(range, "A")) {
00311         irange = 1;
00312     } else if (template_blas_lsame(range, "V")) {
00313         irange = 2;
00314     } else if (template_blas_lsame(range, "I")) {
00315         irange = 3;
00316     } else {
00317         irange = 0;
00318     }
00319 
00320 /*     Check for Errors */
00321 
00322     if (irange <= 0) {
00323         *info = -1;
00324     } else if (! (template_blas_lsame(order, "B") || template_blas_lsame(order, 
00325             "E"))) {
00326         *info = -2;
00327     } else if (*n < 0) {
00328         *info = -3;
00329     } else if (irange == 2) {
00330         if (*vl >= *vu) {
00331             *info = -5;
00332         }
00333     } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) {
00334         *info = -6;
00335     } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) {
00336         *info = -7;
00337     }
00338 
00339     if (*info != 0) {
00340         return 0;
00341     }
00342 /*     Initialize error flags */
00343     *info = 0;
00344     ncnvrg = FALSE_;
00345     toofew = FALSE_;
00346 /*     Quick return if possible */
00347     *m = 0;
00348     if (*n == 0) {
00349         return 0;
00350     }
00351 /*     Simplification: */
00352     if (irange == 3 && *il == 1 && *iu == *n) {
00353         irange = 1;
00354     }
00355 /*     Get machine constants */
00356     eps = template_lapack_lamch("P",(Treal)0);
00357     uflow = template_lapack_lamch("U",(Treal)0);
00358 /*     Special Case when N=1 */
00359 /*     Treat case of 1x1 matrix for quick return */
00360     if (*n == 1) {
00361       if ( irange == 1 || ( irange == 2 && d__[1] > *vl && d__[1] <= *vu ) || 
00362            ( irange == 3 && *il == 1 && *iu == 1 ) ) {
00363             *m = 1;
00364             w[1] = d__[1];
00365 /*           The computation error of the eigenvalue is zero */
00366             werr[1] = 0.;
00367             iblock[1] = 1;
00368             indexw[1] = 1;
00369         }
00370         return 0;
00371     }
00372 /*     NB is the minimum vector length for vector bisection, or 0 */
00373 /*     if only scalar is to be done. */
00374     nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1);
00375     if (nb <= 1) {
00376         nb = 0;
00377     }
00378 /*     Find global spectral radius */
00379     gl = d__[1];
00380     gu = d__[1];
00381     i__1 = *n;
00382     for (i__ = 1; i__ <= i__1; ++i__) {
00383 /* Computing MIN */
00384         d__1 = gl, d__2 = gers[(i__ << 1) - 1];
00385         gl = minMACRO(d__1,d__2);
00386 /* Computing MAX */
00387         d__1 = gu, d__2 = gers[i__ * 2];
00388         gu = maxMACRO(d__1,d__2);
00389 /* L5: */
00390     }
00391 /*     Compute global Gerschgorin bounds and spectral diameter */
00392 /* Computing MAX */
00393     d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00394     tnorm = maxMACRO(d__1,d__2);
00395     gl = gl - tnorm * 2. * eps * *n - *pivmin * 4.;
00396     gu = gu + tnorm * 2. * eps * *n + *pivmin * 4.;
00397 /*     [JAN/28/2009] remove the line below since SPDIAM variable not use */
00398 /*     SPDIAM = GU - GL */
00399 /*     Input arguments for DLAEBZ: */
00400 /*     The relative tolerance.  An interval (a,b] lies within */
00401 /*     "relative tolerance" if  b-a < RELTOL*max(|a|,|b|), */
00402     rtoli = *reltol;
00403 /*     Set the absolute tolerance for interval convergence to zero to force */
00404 /*     interval convergence based on relative size of the interval. */
00405 /*     This is dangerous because intervals might not converge when RELTOL is */
00406 /*     small. But at least a very small number should be selected so that for */
00407 /*     strongly graded matrices, the code can get relatively accurate */
00408 /*     eigenvalues. */
00409     atoli = uflow * 4. + *pivmin * 4.;
00410     if (irange == 3) {
00411 /*        RANGE='I': Compute an interval containing eigenvalues */
00412 /*        IL through IU. The initial interval [GL,GU] from the global */
00413 /*        Gerschgorin bounds GL and GU is refined by DLAEBZ. */
00414         itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 
00415                 2;
00416         work[*n + 1] = gl;
00417         work[*n + 2] = gl;
00418         work[*n + 3] = gu;
00419         work[*n + 4] = gu;
00420         work[*n + 5] = gl;
00421         work[*n + 6] = gu;
00422         iwork[1] = -1;
00423         iwork[2] = -1;
00424         iwork[3] = *n + 1;
00425         iwork[4] = *n + 1;
00426         iwork[5] = *il - 1;
00427         iwork[6] = *iu;
00428 
00429         template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, pivmin, &
00430                 d__[1], &e[1], &e2[1], &iwork[5], &work[*n + 1], &work[*n + 5]
00431 , &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
00432         if (iinfo != 0) {
00433             *info = iinfo;
00434             return 0;
00435         }
00436 /*        On exit, output intervals may not be ordered by ascending negcount */
00437         if (iwork[6] == *iu) {
00438             *wl = work[*n + 1];
00439             wlu = work[*n + 3];
00440             nwl = iwork[1];
00441             *wu = work[*n + 4];
00442             wul = work[*n + 2];
00443             nwu = iwork[4];
00444         } else {
00445             *wl = work[*n + 2];
00446             wlu = work[*n + 4];
00447             nwl = iwork[2];
00448             *wu = work[*n + 3];
00449             wul = work[*n + 1];
00450             nwu = iwork[3];
00451         }
00452 /*        On exit, the interval [WL, WLU] contains a value with negcount NWL, */
00453 /*        and [WUL, WU] contains a value with negcount NWU. */
00454         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
00455             *info = 4;
00456             return 0;
00457         }
00458     } else if (irange == 2) {
00459         *wl = *vl;
00460         *wu = *vu;
00461     } else if (irange == 1) {
00462         *wl = gl;
00463         *wu = gu;
00464     }
00465 /*     Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. */
00466 /*     NWL accumulates the number of eigenvalues .le. WL, */
00467 /*     NWU accumulates the number of eigenvalues .le. WU */
00468     *m = 0;
00469     iend = 0;
00470     *info = 0;
00471     nwl = 0;
00472     nwu = 0;
00473 
00474     i__1 = *nsplit;
00475     for (jblk = 1; jblk <= i__1; ++jblk) {
00476         ioff = iend;
00477         ibegin = ioff + 1;
00478         iend = isplit[jblk];
00479         in = iend - ioff;
00480 
00481         if (in == 1) {
00482 /*           1x1 block */
00483             if (*wl >= d__[ibegin] - *pivmin) {
00484                 ++nwl;
00485             }
00486             if (*wu >= d__[ibegin] - *pivmin) {
00487                 ++nwu;
00488             }
00489             if (irange == 1 || ( *wl < d__[ibegin] - *pivmin && *wu >= d__[
00490                                                                            ibegin] - *pivmin ) ) {
00491                 ++(*m);
00492                 w[*m] = d__[ibegin];
00493                 werr[*m] = 0.;
00494 /*              The gap for a single block doesn't matter for the later */
00495 /*              algorithm and is assigned an arbitrary large value */
00496                 iblock[*m] = jblk;
00497                 indexw[*m] = 1;
00498             }
00499 /*        Disabled 2x2 case because of a failure on the following matrix */
00500 /*        RANGE = 'I', IL = IU = 4 */
00501 /*          Original Tridiagonal, d = [ */
00502 /*           -0.150102010615740E+00 */
00503 /*           -0.849897989384260E+00 */
00504 /*           -0.128208148052635E-15 */
00505 /*            0.128257718286320E-15 */
00506 /*          ]; */
00507 /*          e = [ */
00508 /*           -0.357171383266986E+00 */
00509 /*           -0.180411241501588E-15 */
00510 /*           -0.175152352710251E-15 */
00511 /*          ]; */
00512 
00513 /*         ELSE IF( IN.EQ.2 ) THEN */
00514 /* *           2x2 block */
00515 /*            DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) */
00516 /*            TMP1 = HALF*(D(IBEGIN)+D(IEND)) */
00517 /*            L1 = TMP1 - DISC */
00518 /*            IF( WL.GE. L1-PIVMIN ) */
00519 /*     $         NWL = NWL + 1 */
00520 /*            IF( WU.GE. L1-PIVMIN ) */
00521 /*     $         NWU = NWU + 1 */
00522 /*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. */
00523 /*     $          L1-PIVMIN ) ) THEN */
00524 /*               M = M + 1 */
00525 /*               W( M ) = L1 */
00526 /* *              The uncertainty of eigenvalues of a 2x2 matrix is very small */
00527 /*               WERR( M ) = EPS * ABS( W( M ) ) * TWO */
00528 /*               IBLOCK( M ) = JBLK */
00529 /*               INDEXW( M ) = 1 */
00530 /*            ENDIF */
00531 /*            L2 = TMP1 + DISC */
00532 /*            IF( WL.GE. L2-PIVMIN ) */
00533 /*     $         NWL = NWL + 1 */
00534 /*            IF( WU.GE. L2-PIVMIN ) */
00535 /*     $         NWU = NWU + 1 */
00536 /*            IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. */
00537 /*     $          L2-PIVMIN ) ) THEN */
00538 /*               M = M + 1 */
00539 /*               W( M ) = L2 */
00540 /* *              The uncertainty of eigenvalues of a 2x2 matrix is very small */
00541 /*               WERR( M ) = EPS * ABS( W( M ) ) * TWO */
00542 /*               IBLOCK( M ) = JBLK */
00543 /*               INDEXW( M ) = 2 */
00544 /*            ENDIF */
00545         } else {
00546 /*           General Case - block of size IN >= 2 */
00547 /*           Compute local Gerschgorin interval and use it as the initial */
00548 /*           interval for DLAEBZ */
00549             gu = d__[ibegin];
00550             gl = d__[ibegin];
00551             tmp1 = 0.;
00552             i__2 = iend;
00553             for (j = ibegin; j <= i__2; ++j) {
00554 /* Computing MIN */
00555                 d__1 = gl, d__2 = gers[(j << 1) - 1];
00556                 gl = minMACRO(d__1,d__2);
00557 /* Computing MAX */
00558                 d__1 = gu, d__2 = gers[j * 2];
00559                 gu = maxMACRO(d__1,d__2);
00560 /* L40: */
00561             }
00562 /*           [JAN/28/2009] */
00563 /*           change SPDIAM by TNORM in lines 2 and 3 thereafter */
00564 /*           line 1: remove computation of SPDIAM (not useful anymore) */
00565 /*           SPDIAM = GU - GL */
00566 /*           GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN */
00567 /*           GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN */
00568             gl = gl - tnorm * 2. * eps * in - *pivmin * 2.;
00569             gu = gu + tnorm * 2. * eps * in + *pivmin * 2.;
00570 
00571             if (irange > 1) {
00572                 if (gu < *wl) {
00573 /*                 the local block contains none of the wanted eigenvalues */
00574                     nwl += in;
00575                     nwu += in;
00576                     goto L70;
00577                 }
00578 /*              refine search interval if possible, only range (WL,WU] matters */
00579                 gl = maxMACRO(gl,*wl);
00580                 gu = minMACRO(gu,*wu);
00581                 if (gl >= gu) {
00582                     goto L70;
00583                 }
00584             }
00585 /*           Find negcount of initial interval boundaries GL and GU */
00586             work[*n + 1] = gl;
00587             work[*n + in + 1] = gu;
00588             template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, 
00589                     pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
00590                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
00591                     w[*m + 1], &iblock[*m + 1], &iinfo);
00592             if (iinfo != 0) {
00593                 *info = iinfo;
00594                 return 0;
00595             }
00596 
00597             nwl += iwork[1];
00598             nwu += iwork[in + 1];
00599             iwoff = *m - iwork[1];
00600 /*           Compute Eigenvalues */
00601             itmax = (integer) ((template_blas_log(gu - gl + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(
00602                     2.)) + 2;
00603             template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, 
00604                     pivmin, &d__[ibegin], &e[ibegin], &e2[ibegin], idumma, &
00605                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1], 
00606                      &w[*m + 1], &iblock[*m + 1], &iinfo);
00607             if (iinfo != 0) {
00608                 *info = iinfo;
00609                 return 0;
00610             }
00611 
00612 /*           Copy eigenvalues into W and IBLOCK */
00613 /*           Use -JBLK for block number for unconverged eigenvalues. */
00614 /*           Loop over the number of output intervals from DLAEBZ */
00615             i__2 = iout;
00616             for (j = 1; j <= i__2; ++j) {
00617 /*              eigenvalue approximation is middle point of interval */
00618                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
00619 /*              semi length of error interval */
00620                 tmp2 = (d__1 = work[j + *n] - work[j + in + *n], absMACRO(d__1)) * 
00621                         .5;
00622                 if (j > iout - iinfo) {
00623 /*                 Flag non-convergence. */
00624                     ncnvrg = TRUE_;
00625                     ib = -jblk;
00626                 } else {
00627                     ib = jblk;
00628                 }
00629                 i__3 = iwork[j + in] + iwoff;
00630                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
00631                     w[je] = tmp1;
00632                     werr[je] = tmp2;
00633                     indexw[je] = je - iwoff;
00634                     iblock[je] = ib;
00635 /* L50: */
00636                 }
00637 /* L60: */
00638             }
00639 
00640             *m += im;
00641         }
00642 L70:
00643         ;
00644     }
00645 /*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU */
00646 /*     If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
00647     if (irange == 3) {
00648         idiscl = *il - 1 - nwl;
00649         idiscu = nwu - *iu;
00650 
00651         if (idiscl > 0) {
00652             im = 0;
00653             i__1 = *m;
00654             for (je = 1; je <= i__1; ++je) {
00655 /*              Remove some of the smallest eigenvalues from the left so that */
00656 /*              at the end IDISCL =0. Move all eigenvalues up to the left. */
00657                 if (w[je] <= wlu && idiscl > 0) {
00658                     --idiscl;
00659                 } else {
00660                     ++im;
00661                     w[im] = w[je];
00662                     werr[im] = werr[je];
00663                     indexw[im] = indexw[je];
00664                     iblock[im] = iblock[je];
00665                 }
00666 /* L80: */
00667             }
00668             *m = im;
00669         }
00670         if (idiscu > 0) {
00671 /*           Remove some of the largest eigenvalues from the right so that */
00672 /*           at the end IDISCU =0. Move all eigenvalues up to the left. */
00673             im = *m + 1;
00674             for (je = *m; je >= 1; --je) {
00675                 if (w[je] >= wul && idiscu > 0) {
00676                     --idiscu;
00677                 } else {
00678                     --im;
00679                     w[im] = w[je];
00680                     werr[im] = werr[je];
00681                     indexw[im] = indexw[je];
00682                     iblock[im] = iblock[je];
00683                 }
00684 /* L81: */
00685             }
00686             jee = 0;
00687             i__1 = *m;
00688             for (je = im; je <= i__1; ++je) {
00689                 ++jee;
00690                 w[jee] = w[je];
00691                 werr[jee] = werr[je];
00692                 indexw[jee] = indexw[je];
00693                 iblock[jee] = iblock[je];
00694 /* L82: */
00695             }
00696             *m = *m - im + 1;
00697         }
00698         if (idiscl > 0 || idiscu > 0) {
00699 /*           Code to deal with effects of bad arithmetic. (If N(w) is */
00700 /*           monotone non-decreasing, this should never happen.) */
00701 /*           Some low eigenvalues to be discarded are not in (WL,WLU], */
00702 /*           or high eigenvalues to be discarded are not in (WUL,WU] */
00703 /*           so just kill off the smallest IDISCL/largest IDISCU */
00704 /*           eigenvalues, by marking the corresponding IBLOCK = 0 */
00705             if (idiscl > 0) {
00706                 wkill = *wu;
00707                 i__1 = idiscl;
00708                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00709                     iw = 0;
00710                     i__2 = *m;
00711                     for (je = 1; je <= i__2; ++je) {
00712                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
00713                             iw = je;
00714                             wkill = w[je];
00715                         }
00716 /* L90: */
00717                     }
00718                     iblock[iw] = 0;
00719 /* L100: */
00720                 }
00721             }
00722             if (idiscu > 0) {
00723                 wkill = *wl;
00724                 i__1 = idiscu;
00725                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00726                     iw = 0;
00727                     i__2 = *m;
00728                     for (je = 1; je <= i__2; ++je) {
00729                         if (iblock[je] != 0 && (w[je] >= wkill || iw == 0)) {
00730                             iw = je;
00731                             wkill = w[je];
00732                         }
00733 /* L110: */
00734                     }
00735                     iblock[iw] = 0;
00736 /* L120: */
00737                 }
00738             }
00739 /*           Now erase all eigenvalues with IBLOCK set to zero */
00740             im = 0;
00741             i__1 = *m;
00742             for (je = 1; je <= i__1; ++je) {
00743                 if (iblock[je] != 0) {
00744                     ++im;
00745                     w[im] = w[je];
00746                     werr[im] = werr[je];
00747                     indexw[im] = indexw[je];
00748                     iblock[im] = iblock[je];
00749                 }
00750 /* L130: */
00751             }
00752             *m = im;
00753         }
00754         if (idiscl < 0 || idiscu < 0) {
00755             toofew = TRUE_;
00756         }
00757     }
00758 
00759     if ( ( irange == 1 && *m != *n ) || ( irange == 3 && *m != *iu - *il + 1 ) ) {
00760         toofew = TRUE_;
00761     }
00762 /*     If ORDER='B', do nothing the eigenvalues are already sorted by */
00763 /*        block. */
00764 /*     If ORDER='E', sort the eigenvalues from smallest to largest */
00765     if (template_blas_lsame(order, "E") && *nsplit > 1) {
00766         i__1 = *m - 1;
00767         for (je = 1; je <= i__1; ++je) {
00768             ie = 0;
00769             tmp1 = w[je];
00770             i__2 = *m;
00771             for (j = je + 1; j <= i__2; ++j) {
00772                 if (w[j] < tmp1) {
00773                     ie = j;
00774                     tmp1 = w[j];
00775                 }
00776 /* L140: */
00777             }
00778             if (ie != 0) {
00779                 tmp2 = werr[ie];
00780                 itmp1 = iblock[ie];
00781                 itmp2 = indexw[ie];
00782                 w[ie] = w[je];
00783                 werr[ie] = werr[je];
00784                 iblock[ie] = iblock[je];
00785                 indexw[ie] = indexw[je];
00786                 w[je] = tmp1;
00787                 werr[je] = tmp2;
00788                 iblock[je] = itmp1;
00789                 indexw[je] = itmp2;
00790             }
00791 /* L150: */
00792         }
00793     }
00794 
00795     *info = 0;
00796     if (ncnvrg) {
00797         ++(*info);
00798     }
00799     if (toofew) {
00800         *info += 2;
00801     }
00802     return 0;
00803 
00804 /*     End of DLARRD */
00805 
00806 } /* dlarrd_ */
00807 
00808 #endif

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1