ergo
template_lapack_stebz.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_STEBZ_HEADER
00036 #define TEMPLATE_LAPACK_STEBZ_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal 
00041         *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, 
00042         const Treal *d__, const Treal *e, integer *m, integer *nsplit, 
00043         Treal *w, integer *iblock, integer *isplit, Treal *work, 
00044         integer *iwork, integer *info)
00045 {
00046 /*  -- LAPACK routine (version 3.0) --   
00047        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00048        Courant Institute, Argonne National Lab, and Rice University   
00049        June 30, 1999   
00050 
00051 
00052     Purpose   
00053     =======   
00054 
00055     DSTEBZ computes the eigenvalues of a symmetric tridiagonal   
00056     matrix T.  The user may ask for all eigenvalues, all eigenvalues   
00057     in the half-open interval (VL, VU], or the IL-th through IU-th   
00058     eigenvalues.   
00059 
00060     To avoid overflow, the matrix must be scaled so that its   
00061     largest element is no greater than overflow**(1/2) *   
00062     underflow**(1/4) in absolute value, and for greatest   
00063     accuracy, it should not be much smaller than that.   
00064 
00065     See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal   
00066     Matrix", Report CS41, Computer Science Dept., Stanford   
00067     University, July 21, 1966.   
00068 
00069     Arguments   
00070     =========   
00071 
00072     RANGE   (input) CHARACTER   
00073             = 'A': ("All")   all eigenvalues will be found.   
00074             = 'V': ("Value") all eigenvalues in the half-open interval   
00075                              (VL, VU] will be found.   
00076             = 'I': ("Index") the IL-th through IU-th eigenvalues (of the   
00077                              entire matrix) will be found.   
00078 
00079     ORDER   (input) CHARACTER   
00080             = 'B': ("By Block") the eigenvalues will be grouped by   
00081                                 split-off block (see IBLOCK, ISPLIT) and   
00082                                 ordered from smallest to largest within   
00083                                 the block.   
00084             = 'E': ("Entire matrix")   
00085                                 the eigenvalues for the entire matrix   
00086                                 will be ordered from smallest to   
00087                                 largest.   
00088 
00089     N       (input) INTEGER   
00090             The order of the tridiagonal matrix T.  N >= 0.   
00091 
00092     VL      (input) DOUBLE PRECISION   
00093     VU      (input) DOUBLE PRECISION   
00094             If RANGE='V', the lower and upper bounds of the interval to   
00095             be searched for eigenvalues.  Eigenvalues less than or equal   
00096             to VL, or greater than VU, will not be returned.  VL < VU.   
00097             Not referenced if RANGE = 'A' or 'I'.   
00098 
00099     IL      (input) INTEGER   
00100     IU      (input) INTEGER   
00101             If RANGE='I', the indices (in ascending order) of the   
00102             smallest and largest eigenvalues to be returned.   
00103             1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.   
00104             Not referenced if RANGE = 'A' or 'V'.   
00105 
00106     ABSTOL  (input) DOUBLE PRECISION   
00107             The absolute tolerance for the eigenvalues.  An eigenvalue   
00108             (or cluster) is considered to be located if it has been   
00109             determined to lie in an interval whose width is ABSTOL or   
00110             less.  If ABSTOL is less than or equal to zero, then ULP*|T|   
00111             will be used, where |T| means the 1-norm of T.   
00112 
00113             Eigenvalues will be computed most accurately when ABSTOL is   
00114             set to twice the underflow threshold 2*DLAMCH('S'), not zero.   
00115 
00116     D       (input) DOUBLE PRECISION array, dimension (N)   
00117             The n diagonal elements of the tridiagonal matrix T.   
00118 
00119     E       (input) DOUBLE PRECISION array, dimension (N-1)   
00120             The (n-1) off-diagonal elements of the tridiagonal matrix T.   
00121 
00122     M       (output) INTEGER   
00123             The actual number of eigenvalues found. 0 <= M <= N.   
00124             (See also the description of INFO=2,3.)   
00125 
00126     NSPLIT  (output) INTEGER   
00127             The number of diagonal blocks in the matrix T.   
00128             1 <= NSPLIT <= N.   
00129 
00130     W       (output) DOUBLE PRECISION array, dimension (N)   
00131             On exit, the first M elements of W will contain the   
00132             eigenvalues.  (DSTEBZ may use the remaining N-M elements as   
00133             workspace.)   
00134 
00135     IBLOCK  (output) INTEGER array, dimension (N)   
00136             At each row/column j where E(j) is zero or small, the   
00137             matrix T is considered to split into a block diagonal   
00138             matrix.  On exit, if INFO = 0, IBLOCK(i) specifies to which   
00139             block (from 1 to the number of blocks) the eigenvalue W(i)   
00140             belongs.  (DSTEBZ may use the remaining N-M elements as   
00141             workspace.)   
00142 
00143     ISPLIT  (output) INTEGER array, dimension (N)   
00144             The splitting points, at which T breaks up into submatrices.   
00145             The first submatrix consists of rows/columns 1 to ISPLIT(1),   
00146             the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),   
00147             etc., and the NSPLIT-th consists of rows/columns   
00148             ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.   
00149             (Only the first NSPLIT elements will actually be used, but   
00150             since the user cannot know a priori what value NSPLIT will   
00151             have, N words must be reserved for ISPLIT.)   
00152 
00153     WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)   
00154 
00155     IWORK   (workspace) INTEGER array, dimension (3*N)   
00156 
00157     INFO    (output) INTEGER   
00158             = 0:  successful exit   
00159             < 0:  if INFO = -i, the i-th argument had an illegal value   
00160             > 0:  some or all of the eigenvalues failed to converge or   
00161                   were not computed:   
00162                   =1 or 3: Bisection failed to converge for some   
00163                           eigenvalues; these eigenvalues are flagged by a   
00164                           negative block number.  The effect is that the   
00165                           eigenvalues may not be as accurate as the   
00166                           absolute and relative tolerances.  This is   
00167                           generally caused by unexpectedly inaccurate   
00168                           arithmetic.   
00169                   =2 or 3: RANGE='I' only: Not all of the eigenvalues   
00170                           IL:IU were found.   
00171                           Effect: M < IU+1-IL   
00172                           Cause:  non-monotonic arithmetic, causing the   
00173                                   Sturm sequence to be non-monotonic.   
00174                           Cure:   recalculate, using RANGE='A', and pick   
00175                                   out eigenvalues IL:IU.  In some cases,   
00176                                   increasing the PARAMETER "FUDGE" may   
00177                                   make things work.   
00178                   = 4:    RANGE='I', and the Gershgorin interval   
00179                           initially used was too small.  No eigenvalues   
00180                           were computed.   
00181                           Probable cause: your machine has sloppy   
00182                                           floating-point arithmetic.   
00183                           Cure: Increase the PARAMETER "FUDGE",   
00184                                 recompile, and try again.   
00185 
00186     Internal Parameters   
00187     ===================   
00188 
00189     RELFAC  DOUBLE PRECISION, default = 2.0e0   
00190             The relative tolerance.  An interval (a,b] lies within   
00191             "relative tolerance" if  b-a < RELFAC*ulp*max(|a|,|b|),   
00192             where "ulp" is the machine precision (distance from 1 to   
00193             the next larger floating point number.)   
00194 
00195     FUDGE   DOUBLE PRECISION, default = 2   
00196             A "fudge factor" to widen the Gershgorin intervals.  Ideally,   
00197             a value of 1 should work, but on machines with sloppy   
00198             arithmetic, this needs to be larger.  The default for   
00199             publicly released versions should be large enough to handle   
00200             the worst machine around.  Note that this has no effect   
00201             on accuracy of the solution.   
00202 
00203     =====================================================================   
00204 
00205 
00206        Parameter adjustments */
00207     /* Table of constant values */
00208      integer c__1 = 1;
00209      integer c_n1 = -1;
00210      integer c__3 = 3;
00211      integer c__2 = 2;
00212      integer c__0 = 0;
00213     
00214     /* System generated locals */
00215     integer i__1, i__2, i__3;
00216     Treal d__1, d__2, d__3, d__4, d__5;
00217     /* Local variables */
00218      integer iend, ioff, iout, itmp1, j, jdisc;
00219      integer iinfo;
00220      Treal atoli;
00221      integer iwoff;
00222      Treal bnorm;
00223      integer itmax;
00224      Treal wkill, rtoli, tnorm;
00225      integer ib, jb, ie, je, nb;
00226      Treal gl;
00227      integer im, in;
00228      integer ibegin;
00229      Treal gu;
00230      integer iw;
00231      Treal wl;
00232      integer irange, idiscl;
00233      Treal safemn, wu;
00234      integer idumma[1];
00235      integer idiscu, iorder;
00236      logical ncnvrg;
00237      Treal pivmin;
00238      logical toofew;
00239      integer nwl;
00240      Treal ulp, wlu, wul;
00241      integer nwu;
00242      Treal tmp1, tmp2;
00243 
00244 
00245     --iwork;
00246     --work;
00247     --isplit;
00248     --iblock;
00249     --w;
00250     --e;
00251     --d__;
00252 
00253     /* Initialization added by Elias to get rid of compiler warnings. */
00254     wlu = wul = 0;
00255     /* Function Body */
00256     *info = 0;
00257 
00258 /*     Decode RANGE */
00259 
00260     if (template_blas_lsame(range, "A")) {
00261         irange = 1;
00262     } else if (template_blas_lsame(range, "V")) {
00263         irange = 2;
00264     } else if (template_blas_lsame(range, "I")) {
00265         irange = 3;
00266     } else {
00267         irange = 0;
00268     }
00269 
00270 /*     Decode ORDER */
00271 
00272     if (template_blas_lsame(order, "B")) {
00273         iorder = 2;
00274     } else if (template_blas_lsame(order, "E")) {
00275         iorder = 1;
00276     } else {
00277         iorder = 0;
00278     }
00279 
00280 /*     Check for Errors */
00281 
00282     if (irange <= 0) {
00283         *info = -1;
00284     } else if (iorder <= 0) {
00285         *info = -2;
00286     } else if (*n < 0) {
00287         *info = -3;
00288     } else if (irange == 2) {
00289         if (*vl >= *vu) {
00290             *info = -5;
00291         }
00292     } else if (irange == 3 && (*il < 1 || *il > maxMACRO(1,*n))) {
00293         *info = -6;
00294     } else if (irange == 3 && (*iu < minMACRO(*n,*il) || *iu > *n)) {
00295         *info = -7;
00296     }
00297 
00298     if (*info != 0) {
00299         i__1 = -(*info);
00300         template_blas_erbla("STEBZ ", &i__1);
00301         return 0;
00302     }
00303 
00304 /*     Initialize error flags */
00305 
00306     *info = 0;
00307     ncnvrg = FALSE_;
00308     toofew = FALSE_;
00309 
00310 /*     Quick return if possible */
00311 
00312     *m = 0;
00313     if (*n == 0) {
00314         return 0;
00315     }
00316 
00317 /*     Simplifications: */
00318 
00319     if (irange == 3 && *il == 1 && *iu == *n) {
00320         irange = 1;
00321     }
00322 
00323 /*     Get machine constants   
00324        NB is the minimum vector length for vector bisection, or 0   
00325        if only scalar is to be done. */
00326 
00327     safemn = template_lapack_lamch("S", (Treal)0);
00328     ulp = template_lapack_lamch("P", (Treal)0);
00329     rtoli = ulp * 2.;
00330     nb = template_lapack_ilaenv(&c__1, "DSTEBZ", " ", n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
00331             ftnlen)1);
00332     if (nb <= 1) {
00333         nb = 0;
00334     }
00335 
00336 /*     Special Case when N=1 */
00337 
00338     if (*n == 1) {
00339         *nsplit = 1;
00340         isplit[1] = 1;
00341         if (irange == 2 && (*vl >= d__[1] || *vu < d__[1])) {
00342             *m = 0;
00343         } else {
00344             w[1] = d__[1];
00345             iblock[1] = 1;
00346             *m = 1;
00347         }
00348         return 0;
00349     }
00350 
00351 /*     Compute Splitting Points */
00352 
00353     *nsplit = 1;
00354     work[*n] = 0.;
00355     pivmin = 1.;
00356 
00357 /* DIR$ NOVECTOR */
00358     i__1 = *n;
00359     for (j = 2; j <= i__1; ++j) {
00360 /* Computing 2nd power */
00361         d__1 = e[j - 1];
00362         tmp1 = d__1 * d__1;
00363 /* Computing 2nd power */
00364         d__2 = ulp;
00365         if ((d__1 = d__[j] * d__[j - 1], absMACRO(d__1)) * (d__2 * d__2) + safemn 
00366                 > tmp1) {
00367             isplit[*nsplit] = j - 1;
00368             ++(*nsplit);
00369             work[j - 1] = 0.;
00370         } else {
00371             work[j - 1] = tmp1;
00372             pivmin = maxMACRO(pivmin,tmp1);
00373         }
00374 /* L10: */
00375     }
00376     isplit[*nsplit] = *n;
00377     pivmin *= safemn;
00378 
00379 /*     Compute Interval and ATOLI */
00380 
00381     if (irange == 3) {
00382 
00383 /*        RANGE='I': Compute the interval containing eigenvalues   
00384                      IL through IU.   
00385 
00386           Compute Gershgorin interval for entire (split) matrix   
00387           and use it as the initial interval */
00388 
00389         gu = d__[1];
00390         gl = d__[1];
00391         tmp1 = 0.;
00392 
00393         i__1 = *n - 1;
00394         for (j = 1; j <= i__1; ++j) {
00395             tmp2 = template_blas_sqrt(work[j]);
00396 /* Computing MAX */
00397             d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
00398             gu = maxMACRO(d__1,d__2);
00399 /* Computing MIN */
00400             d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
00401             gl = minMACRO(d__1,d__2);
00402             tmp1 = tmp2;
00403 /* L20: */
00404         }
00405 
00406 /* Computing MAX */
00407         d__1 = gu, d__2 = d__[*n] + tmp1;
00408         gu = maxMACRO(d__1,d__2);
00409 /* Computing MIN */
00410         d__1 = gl, d__2 = d__[*n] - tmp1;
00411         gl = minMACRO(d__1,d__2);
00412 /* Computing MAX */
00413         d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00414         tnorm = maxMACRO(d__1,d__2);
00415         gl = gl - tnorm * 2. * ulp * *n - pivmin * 4.;
00416         gu = gu + tnorm * 2. * ulp * *n + pivmin * 2.;
00417 
00418 /*        Compute Iteration parameters */
00419 
00420         itmax = (integer) ((template_blas_log(tnorm + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)) + 2;
00421         if (*abstol <= 0.) {
00422             atoli = ulp * tnorm;
00423         } else {
00424             atoli = *abstol;
00425         }
00426 
00427         work[*n + 1] = gl;
00428         work[*n + 2] = gl;
00429         work[*n + 3] = gu;
00430         work[*n + 4] = gu;
00431         work[*n + 5] = gl;
00432         work[*n + 6] = gu;
00433         iwork[1] = -1;
00434         iwork[2] = -1;
00435         iwork[3] = *n + 1;
00436         iwork[4] = *n + 1;
00437         iwork[5] = *il - 1;
00438         iwork[6] = *iu;
00439 
00440         template_lapack_laebz(&c__3, &itmax, n, &c__2, &c__2, &nb, &atoli, &rtoli, &pivmin, 
00441                 &d__[1], &e[1], &work[1], &iwork[5], &work[*n + 1], &work[*n 
00442                 + 5], &iout, &iwork[1], &w[1], &iblock[1], &iinfo);
00443 
00444         if (iwork[6] == *iu) {
00445             wl = work[*n + 1];
00446             wlu = work[*n + 3];
00447             nwl = iwork[1];
00448             wu = work[*n + 4];
00449             wul = work[*n + 2];
00450             nwu = iwork[4];
00451         } else {
00452             wl = work[*n + 2];
00453             wlu = work[*n + 4];
00454             nwl = iwork[2];
00455             wu = work[*n + 3];
00456             wul = work[*n + 1];
00457             nwu = iwork[3];
00458         }
00459 
00460         if (nwl < 0 || nwl >= *n || nwu < 1 || nwu > *n) {
00461             *info = 4;
00462             return 0;
00463         }
00464     } else {
00465 
00466 /*        RANGE='A' or 'V' -- Set ATOLI   
00467 
00468    Computing MAX */
00469         d__3 = absMACRO(d__[1]) + absMACRO(e[1]), d__4 = (d__1 = d__[*n], absMACRO(d__1)) + (
00470                 d__2 = e[*n - 1], absMACRO(d__2));
00471         tnorm = maxMACRO(d__3,d__4);
00472 
00473         i__1 = *n - 1;
00474         for (j = 2; j <= i__1; ++j) {
00475 /* Computing MAX */
00476             d__4 = tnorm, d__5 = (d__1 = d__[j], absMACRO(d__1)) + (d__2 = e[j - 1]
00477                     , absMACRO(d__2)) + (d__3 = e[j], absMACRO(d__3));
00478             tnorm = maxMACRO(d__4,d__5);
00479 /* L30: */
00480         }
00481 
00482         if (*abstol <= 0.) {
00483             atoli = ulp * tnorm;
00484         } else {
00485             atoli = *abstol;
00486         }
00487 
00488         if (irange == 2) {
00489             wl = *vl;
00490             wu = *vu;
00491         } else {
00492             wl = 0.;
00493             wu = 0.;
00494         }
00495     }
00496 
00497 /*     Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.   
00498        NWL accumulates the number of eigenvalues .le. WL,   
00499        NWU accumulates the number of eigenvalues .le. WU */
00500 
00501     *m = 0;
00502     iend = 0;
00503     *info = 0;
00504     nwl = 0;
00505     nwu = 0;
00506 
00507     i__1 = *nsplit;
00508     for (jb = 1; jb <= i__1; ++jb) {
00509         ioff = iend;
00510         ibegin = ioff + 1;
00511         iend = isplit[jb];
00512         in = iend - ioff;
00513 
00514         if (in == 1) {
00515 
00516 /*           Special Case -- IN=1 */
00517 
00518             if (irange == 1 || wl >= d__[ibegin] - pivmin) {
00519                 ++nwl;
00520             }
00521             if (irange == 1 || wu >= d__[ibegin] - pivmin) {
00522                 ++nwu;
00523             }
00524             if (irange == 1 || ( wl < d__[ibegin] - pivmin && wu >= d__[ibegin] 
00525                                  - pivmin ) ) {
00526                 ++(*m);
00527                 w[*m] = d__[ibegin];
00528                 iblock[*m] = jb;
00529             }
00530         } else {
00531 
00532 /*           General Case -- IN > 1   
00533 
00534              Compute Gershgorin Interval   
00535              and use it as the initial interval */
00536 
00537             gu = d__[ibegin];
00538             gl = d__[ibegin];
00539             tmp1 = 0.;
00540 
00541             i__2 = iend - 1;
00542             for (j = ibegin; j <= i__2; ++j) {
00543                 tmp2 = (d__1 = e[j], absMACRO(d__1));
00544 /* Computing MAX */
00545                 d__1 = gu, d__2 = d__[j] + tmp1 + tmp2;
00546                 gu = maxMACRO(d__1,d__2);
00547 /* Computing MIN */
00548                 d__1 = gl, d__2 = d__[j] - tmp1 - tmp2;
00549                 gl = minMACRO(d__1,d__2);
00550                 tmp1 = tmp2;
00551 /* L40: */
00552             }
00553 
00554 /* Computing MAX */
00555             d__1 = gu, d__2 = d__[iend] + tmp1;
00556             gu = maxMACRO(d__1,d__2);
00557 /* Computing MIN */
00558             d__1 = gl, d__2 = d__[iend] - tmp1;
00559             gl = minMACRO(d__1,d__2);
00560 /* Computing MAX */
00561             d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00562             bnorm = maxMACRO(d__1,d__2);
00563             gl = gl - bnorm * 2. * ulp * in - pivmin * 2.;
00564             gu = gu + bnorm * 2. * ulp * in + pivmin * 2.;
00565 
00566 /*           Compute ATOLI for the current submatrix */
00567 
00568             if (*abstol <= 0.) {
00569 /* Computing MAX */
00570                 d__1 = absMACRO(gl), d__2 = absMACRO(gu);
00571                 atoli = ulp * maxMACRO(d__1,d__2);
00572             } else {
00573                 atoli = *abstol;
00574             }
00575 
00576             if (irange > 1) {
00577                 if (gu < wl) {
00578                     nwl += in;
00579                     nwu += in;
00580                     goto L70;
00581                 }
00582                 gl = maxMACRO(gl,wl);
00583                 gu = minMACRO(gu,wu);
00584                 if (gl >= gu) {
00585                     goto L70;
00586                 }
00587             }
00588 
00589 /*           Set Up Initial Interval */
00590 
00591             work[*n + 1] = gl;
00592             work[*n + in + 1] = gu;
00593             template_lapack_laebz(&c__1, &c__0, &in, &in, &c__1, &nb, &atoli, &rtoli, &
00594                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
00595                     work[*n + 1], &work[*n + (in << 1) + 1], &im, &iwork[1], &
00596                     w[*m + 1], &iblock[*m + 1], &iinfo);
00597 
00598             nwl += iwork[1];
00599             nwu += iwork[in + 1];
00600             iwoff = *m - iwork[1];
00601 
00602 /*           Compute Eigenvalues */
00603 
00604             itmax = (integer) ((template_blas_log(gu - gl + pivmin) - template_blas_log(pivmin)) / template_blas_log(2.)
00605                     ) + 2;
00606             template_lapack_laebz(&c__2, &itmax, &in, &in, &c__1, &nb, &atoli, &rtoli, &
00607                     pivmin, &d__[ibegin], &e[ibegin], &work[ibegin], idumma, &
00608                     work[*n + 1], &work[*n + (in << 1) + 1], &iout, &iwork[1],
00609                      &w[*m + 1], &iblock[*m + 1], &iinfo);
00610 
00611 /*           Copy Eigenvalues Into W and IBLOCK   
00612              Use -JB for block number for unconverged eigenvalues. */
00613 
00614             i__2 = iout;
00615             for (j = 1; j <= i__2; ++j) {
00616                 tmp1 = (work[j + *n] + work[j + in + *n]) * .5;
00617 
00618 /*              Flag non-convergence. */
00619 
00620                 if (j > iout - iinfo) {
00621                     ncnvrg = TRUE_;
00622                     ib = -jb;
00623                 } else {
00624                     ib = jb;
00625                 }
00626                 i__3 = iwork[j + in] + iwoff;
00627                 for (je = iwork[j] + 1 + iwoff; je <= i__3; ++je) {
00628                     w[je] = tmp1;
00629                     iblock[je] = ib;
00630 /* L50: */
00631                 }
00632 /* L60: */
00633             }
00634 
00635             *m += im;
00636         }
00637 L70:
00638         ;
00639     }
00640 
00641 /*     If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU   
00642        If NWL+1 < IL or NWU > IU, discard extra eigenvalues. */
00643 
00644     if (irange == 3) {
00645         im = 0;
00646         idiscl = *il - 1 - nwl;
00647         idiscu = nwu - *iu;
00648 
00649         if (idiscl > 0 || idiscu > 0) {
00650             i__1 = *m;
00651             for (je = 1; je <= i__1; ++je) {
00652                 if (w[je] <= wlu && idiscl > 0) {
00653                     --idiscl;
00654                 } else if (w[je] >= wul && idiscu > 0) {
00655                     --idiscu;
00656                 } else {
00657                     ++im;
00658                     w[im] = w[je];
00659                     iblock[im] = iblock[je];
00660                 }
00661 /* L80: */
00662             }
00663             *m = im;
00664         }
00665         if (idiscl > 0 || idiscu > 0) {
00666 
00667 /*           Code to deal with effects of bad arithmetic:   
00668              Some low eigenvalues to be discarded are not in (WL,WLU],   
00669              or high eigenvalues to be discarded are not in (WUL,WU]   
00670              so just kill off the smallest IDISCL/largest IDISCU   
00671              eigenvalues, by simply finding the smallest/largest   
00672              eigenvalue(s).   
00673 
00674              (If N(w) is monotone non-decreasing, this should never   
00675                  happen.) */
00676 
00677             if (idiscl > 0) {
00678                 wkill = wu;
00679                 i__1 = idiscl;
00680                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00681                     iw = 0;
00682                     i__2 = *m;
00683                     for (je = 1; je <= i__2; ++je) {
00684                         if (iblock[je] != 0 && (w[je] < wkill || iw == 0)) {
00685                             iw = je;
00686                             wkill = w[je];
00687                         }
00688 /* L90: */
00689                     }
00690                     iblock[iw] = 0;
00691 /* L100: */
00692                 }
00693             }
00694             if (idiscu > 0) {
00695 
00696                 wkill = wl;
00697                 i__1 = idiscu;
00698                 for (jdisc = 1; jdisc <= i__1; ++jdisc) {
00699                     iw = 0;
00700                     i__2 = *m;
00701                     for (je = 1; je <= i__2; ++je) {
00702                         if (iblock[je] != 0 && (w[je] > wkill || iw == 0)) {
00703                             iw = je;
00704                             wkill = w[je];
00705                         }
00706 /* L110: */
00707                     }
00708                     iblock[iw] = 0;
00709 /* L120: */
00710                 }
00711             }
00712             im = 0;
00713             i__1 = *m;
00714             for (je = 1; je <= i__1; ++je) {
00715                 if (iblock[je] != 0) {
00716                     ++im;
00717                     w[im] = w[je];
00718                     iblock[im] = iblock[je];
00719                 }
00720 /* L130: */
00721             }
00722             *m = im;
00723         }
00724         if (idiscl < 0 || idiscu < 0) {
00725             toofew = TRUE_;
00726         }
00727     }
00728 
00729 /*     If ORDER='B', do nothing -- the eigenvalues are already sorted   
00730           by block.   
00731        If ORDER='E', sort the eigenvalues from smallest to largest */
00732 
00733     if (iorder == 1 && *nsplit > 1) {
00734         i__1 = *m - 1;
00735         for (je = 1; je <= i__1; ++je) {
00736             ie = 0;
00737             tmp1 = w[je];
00738             i__2 = *m;
00739             for (j = je + 1; j <= i__2; ++j) {
00740                 if (w[j] < tmp1) {
00741                     ie = j;
00742                     tmp1 = w[j];
00743                 }
00744 /* L140: */
00745             }
00746 
00747             if (ie != 0) {
00748                 itmp1 = iblock[ie];
00749                 w[ie] = w[je];
00750                 iblock[ie] = iblock[je];
00751                 w[je] = tmp1;
00752                 iblock[je] = itmp1;
00753             }
00754 /* L150: */
00755         }
00756     }
00757 
00758     *info = 0;
00759     if (ncnvrg) {
00760         ++(*info);
00761     }
00762     if (toofew) {
00763         *info += 2;
00764     }
00765     return 0;
00766 
00767 /*     End of DSTEBZ */
00768 
00769 } /* dstebz_ */
00770 
00771 #endif