template_lapack_laebz.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_LAEBZ_HEADER
00036 #define TEMPLATE_LAPACK_LAEBZ_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_laebz(const integer *ijob, const integer *nitmax, const integer *n, 
00041         const integer *mmax, const integer *minp, const integer *nbmin, const Treal *abstol, 
00042         const Treal *reltol, const Treal *pivmin, const Treal *d__, const Treal *
00043         e, const Treal *e2, integer *nval, Treal *ab, Treal *c__, 
00044         integer *mout, integer *nab, Treal *work, integer *iwork, 
00045         integer *info)
00046 {
00047 /*  -- LAPACK auxiliary routine (version 3.0) --   
00048        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00049        Courant Institute, Argonne National Lab, and Rice University   
00050        June 30, 1999   
00051 
00052 
00053     Purpose   
00054     =======   
00055 
00056     DLAEBZ contains the iteration loops which compute and use the   
00057     function N(w), which is the count of eigenvalues of a symmetric   
00058     tridiagonal matrix T less than or equal to its argument  w.  It   
00059     performs a choice of two types of loops:   
00060 
00061     IJOB=1, followed by   
00062     IJOB=2: It takes as input a list of intervals and returns a list of   
00063             sufficiently small intervals whose union contains the same   
00064             eigenvalues as the union of the original intervals.   
00065             The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.   
00066             The output interval (AB(j,1),AB(j,2)] will contain   
00067             eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.   
00068 
00069     IJOB=3: It performs a binary search in each input interval   
00070             (AB(j,1),AB(j,2)] for a point  w(j)  such that   
00071             N(w(j))=NVAL(j), and uses  C(j)  as the starting point of   
00072             the search.  If such a w(j) is found, then on output   
00073             AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output   
00074             (AB(j,1),AB(j,2)] will be a small interval containing the   
00075             point where N(w) jumps through NVAL(j), unless that point   
00076             lies outside the initial interval.   
00077 
00078     Note that the intervals are in all cases half-open intervals,   
00079     i.e., of the form  (a,b] , which includes  b  but not  a .   
00080 
00081     To avoid underflow, the matrix should be scaled so that its largest   
00082     element is no greater than  overflow**(1/2) * underflow**(1/4)   
00083     in absolute value.  To assure the most accurate computation   
00084     of small eigenvalues, the matrix should be scaled to be   
00085     not much smaller than that, either.   
00086 
00087     See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal   
00088     Matrix", Report CS41, Computer Science Dept., Stanford   
00089     University, July 21, 1966   
00090 
00091     Note: the arguments are, in general, *not* checked for unreasonable   
00092     values.   
00093 
00094     Arguments   
00095     =========   
00096 
00097     IJOB    (input) INTEGER   
00098             Specifies what is to be done:   
00099             = 1:  Compute NAB for the initial intervals.   
00100             = 2:  Perform bisection iteration to find eigenvalues of T.   
00101             = 3:  Perform bisection iteration to invert N(w), i.e.,   
00102                   to find a point which has a specified number of   
00103                   eigenvalues of T to its left.   
00104             Other values will cause DLAEBZ to return with INFO=-1.   
00105 
00106     NITMAX  (input) INTEGER   
00107             The maximum number of "levels" of bisection to be   
00108             performed, i.e., an interval of width W will not be made   
00109             smaller than 2^(-NITMAX) * W.  If not all intervals   
00110             have converged after NITMAX iterations, then INFO is set   
00111             to the number of non-converged intervals.   
00112 
00113     N       (input) INTEGER   
00114             The dimension n of the tridiagonal matrix T.  It must be at   
00115             least 1.   
00116 
00117     MMAX    (input) INTEGER   
00118             The maximum number of intervals.  If more than MMAX intervals   
00119             are generated, then DLAEBZ will quit with INFO=MMAX+1.   
00120 
00121     MINP    (input) INTEGER   
00122             The initial number of intervals.  It may not be greater than   
00123             MMAX.   
00124 
00125     NBMIN   (input) INTEGER   
00126             The smallest number of intervals that should be processed   
00127             using a vector loop.  If zero, then only the scalar loop   
00128             will be used.   
00129 
00130     ABSTOL  (input) DOUBLE PRECISION   
00131             The minimum (absolute) width of an interval.  When an   
00132             interval is narrower than ABSTOL, or than RELTOL times the   
00133             larger (in magnitude) endpoint, then it is considered to be   
00134             sufficiently small, i.e., converged.  This must be at least   
00135             zero.   
00136 
00137     RELTOL  (input) DOUBLE PRECISION   
00138             The minimum relative width of an interval.  When an interval   
00139             is narrower than ABSTOL, or than RELTOL times the larger (in   
00140             magnitude) endpoint, then it is considered to be   
00141             sufficiently small, i.e., converged.  Note: this should   
00142             always be at least radix*machine epsilon.   
00143 
00144     PIVMIN  (input) DOUBLE PRECISION   
00145             The minimum absolute value of a "pivot" in the Sturm   
00146             sequence loop.  This *must* be at least  max |e(j)**2| *   
00147             safe_min  and at least safe_min, where safe_min is at least   
00148             the smallest number that can divide one without overflow.   
00149 
00150     D       (input) DOUBLE PRECISION array, dimension (N)   
00151             The diagonal elements of the tridiagonal matrix T.   
00152 
00153     E       (input) DOUBLE PRECISION array, dimension (N)   
00154             The offdiagonal elements of the tridiagonal matrix T in   
00155             positions 1 through N-1.  E(N) is arbitrary.   
00156 
00157     E2      (input) DOUBLE PRECISION array, dimension (N)   
00158             The squares of the offdiagonal elements of the tridiagonal   
00159             matrix T.  E2(N) is ignored.   
00160 
00161     NVAL    (input/output) INTEGER array, dimension (MINP)   
00162             If IJOB=1 or 2, not referenced.   
00163             If IJOB=3, the desired values of N(w).  The elements of NVAL   
00164             will be reordered to correspond with the intervals in AB.   
00165             Thus, NVAL(j) on output will not, in general be the same as   
00166             NVAL(j) on input, but it will correspond with the interval   
00167             (AB(j,1),AB(j,2)] on output.   
00168 
00169     AB      (input/output) DOUBLE PRECISION array, dimension (MMAX,2)   
00170             The endpoints of the intervals.  AB(j,1) is  a(j), the left   
00171             endpoint of the j-th interval, and AB(j,2) is b(j), the   
00172             right endpoint of the j-th interval.  The input intervals   
00173             will, in general, be modified, split, and reordered by the   
00174             calculation.   
00175 
00176     C       (input/output) DOUBLE PRECISION array, dimension (MMAX)   
00177             If IJOB=1, ignored.   
00178             If IJOB=2, workspace.   
00179             If IJOB=3, then on input C(j) should be initialized to the   
00180             first search point in the binary search.   
00181 
00182     MOUT    (output) INTEGER   
00183             If IJOB=1, the number of eigenvalues in the intervals.   
00184             If IJOB=2 or 3, the number of intervals output.   
00185             If IJOB=3, MOUT will equal MINP.   
00186 
00187     NAB     (input/output) INTEGER array, dimension (MMAX,2)   
00188             If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).   
00189             If IJOB=2, then on input, NAB(i,j) should be set.  It must   
00190                satisfy the condition:   
00191                N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)),   
00192                which means that in interval i only eigenvalues   
00193                NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,   
00194                NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with   
00195                IJOB=1.   
00196                On output, NAB(i,j) will contain   
00197                max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of   
00198                the input interval that the output interval   
00199                (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the   
00200                the input values of NAB(k,1) and NAB(k,2).   
00201             If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),   
00202                unless N(w) > NVAL(i) for all search points  w , in which   
00203                case NAB(i,1) will not be modified, i.e., the output   
00204                value will be the same as the input value (modulo   
00205                reorderings -- see NVAL and AB), or unless N(w) < NVAL(i)   
00206                for all search points  w , in which case NAB(i,2) will   
00207                not be modified.  Normally, NAB should be set to some   
00208                distinctive value(s) before DLAEBZ is called.   
00209 
00210     WORK    (workspace) DOUBLE PRECISION array, dimension (MMAX)   
00211             Workspace.   
00212 
00213     IWORK   (workspace) INTEGER array, dimension (MMAX)   
00214             Workspace.   
00215 
00216     INFO    (output) INTEGER   
00217             = 0:       All intervals converged.   
00218             = 1--MMAX: The last INFO intervals did not converge.   
00219             = MMAX+1:  More than MMAX intervals were generated.   
00220 
00221     Further Details   
00222     ===============   
00223 
00224         This routine is intended to be called only by other LAPACK   
00225     routines, thus the interface is less user-friendly.  It is intended   
00226     for two purposes:   
00227 
00228     (a) finding eigenvalues.  In this case, DLAEBZ should have one or   
00229         more initial intervals set up in AB, and DLAEBZ should be called   
00230         with IJOB=1.  This sets up NAB, and also counts the eigenvalues.   
00231         Intervals with no eigenvalues would usually be thrown out at   
00232         this point.  Also, if not all the eigenvalues in an interval i   
00233         are desired, NAB(i,1) can be increased or NAB(i,2) decreased.   
00234         For example, set NAB(i,1)=NAB(i,2)-1 to get the largest   
00235         eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX   
00236         no smaller than the value of MOUT returned by the call with   
00237         IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1   
00238         through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the   
00239         tolerance specified by ABSTOL and RELTOL.   
00240 
00241     (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).   
00242         In this case, start with a Gershgorin interval  (a,b).  Set up   
00243         AB to contain 2 search intervals, both initially (a,b).  One   
00244         NVAL element should contain  f-1  and the other should contain  l   
00245         , while C should contain a and b, resp.  NAB(i,1) should be -1   
00246         and NAB(i,2) should be N+1, to flag an error if the desired   
00247         interval does not lie in (a,b).  DLAEBZ is then called with   
00248         IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals --   
00249         j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while   
00250         if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r   
00251         >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and   
00252         N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and   
00253         w(l-r)=...=w(l+k) are handled similarly.   
00254 
00255     =====================================================================   
00256 
00257 
00258        Check for Errors   
00259 
00260        Parameter adjustments */
00261     /* System generated locals */
00262     integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4, 
00263             i__5, i__6;
00264     Treal d__1, d__2, d__3, d__4;
00265     /* Local variables */
00266      integer itmp1, itmp2, j, kfnew, klnew, kf, ji, kl, jp, jit;
00267      Treal tmp1, tmp2;
00268 #define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1]
00269 #define nab_ref(a_1,a_2) nab[(a_2)*nab_dim1 + a_1]
00270 
00271     nab_dim1 = *mmax;
00272     nab_offset = 1 + nab_dim1 * 1;
00273     nab -= nab_offset;
00274     ab_dim1 = *mmax;
00275     ab_offset = 1 + ab_dim1 * 1;
00276     ab -= ab_offset;
00277     --d__;
00278     --e;
00279     --e2;
00280     --nval;
00281     --c__;
00282     --work;
00283     --iwork;
00284 
00285     /* Function Body */
00286     *info = 0;
00287     if (*ijob < 1 || *ijob > 3) {
00288         *info = -1;
00289         return 0;
00290     }
00291 
00292 /*     Initialize NAB */
00293 
00294     if (*ijob == 1) {
00295 
00296 /*        Compute the number of eigenvalues in the initial intervals. */
00297 
00298         *mout = 0;
00299 /* DIR$ NOVECTOR */
00300         i__1 = *minp;
00301         for (ji = 1; ji <= i__1; ++ji) {
00302             for (jp = 1; jp <= 2; ++jp) {
00303                 tmp1 = d__[1] - ab_ref(ji, jp);
00304                 if (absMACRO(tmp1) < *pivmin) {
00305                     tmp1 = -(*pivmin);
00306                 }
00307                 nab_ref(ji, jp) = 0;
00308                 if (tmp1 <= 0.) {
00309                     nab_ref(ji, jp) = 1;
00310                 }
00311 
00312                 i__2 = *n;
00313                 for (j = 2; j <= i__2; ++j) {
00314                     tmp1 = d__[j] - e2[j - 1] / tmp1 - ab_ref(ji, jp);
00315                     if (absMACRO(tmp1) < *pivmin) {
00316                         tmp1 = -(*pivmin);
00317                     }
00318                     if (tmp1 <= 0.) {
00319                         nab_ref(ji, jp) = nab_ref(ji, jp) + 1;
00320                     }
00321 /* L10: */
00322                 }
00323 /* L20: */
00324             }
00325             *mout = *mout + nab_ref(ji, 2) - nab_ref(ji, 1);
00326 /* L30: */
00327         }
00328         return 0;
00329     }
00330 
00331 /*     Initialize for loop   
00332 
00333        KF and KL have the following meaning:   
00334           Intervals 1,...,KF-1 have converged.   
00335           Intervals KF,...,KL  still need to be refined. */
00336 
00337     kf = 1;
00338     kl = *minp;
00339 
00340 /*     If IJOB=2, initialize C.   
00341        If IJOB=3, use the user-supplied starting point. */
00342 
00343     if (*ijob == 2) {
00344         i__1 = *minp;
00345         for (ji = 1; ji <= i__1; ++ji) {
00346             c__[ji] = (ab_ref(ji, 1) + ab_ref(ji, 2)) * .5;
00347 /* L40: */
00348         }
00349     }
00350 
00351 /*     Iteration loop */
00352 
00353     i__1 = *nitmax;
00354     for (jit = 1; jit <= i__1; ++jit) {
00355 
00356 /*        Loop over intervals */
00357 
00358         if (kl - kf + 1 >= *nbmin && *nbmin > 0) {
00359 
00360 /*           Begin of Parallel Version of the loop */
00361 
00362             i__2 = kl;
00363             for (ji = kf; ji <= i__2; ++ji) {
00364 
00365 /*              Compute N(c), the number of eigenvalues less than c */
00366 
00367                 work[ji] = d__[1] - c__[ji];
00368                 iwork[ji] = 0;
00369                 if (work[ji] <= *pivmin) {
00370                     iwork[ji] = 1;
00371 /* Computing MIN */
00372                     d__1 = work[ji], d__2 = -(*pivmin);
00373                     work[ji] = minMACRO(d__1,d__2);
00374                 }
00375 
00376                 i__3 = *n;
00377                 for (j = 2; j <= i__3; ++j) {
00378                     work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji];
00379                     if (work[ji] <= *pivmin) {
00380                         ++iwork[ji];
00381 /* Computing MIN */
00382                         d__1 = work[ji], d__2 = -(*pivmin);
00383                         work[ji] = minMACRO(d__1,d__2);
00384                     }
00385 /* L50: */
00386                 }
00387 /* L60: */
00388             }
00389 
00390             if (*ijob <= 2) {
00391 
00392 /*              IJOB=2: Choose all intervals containing eigenvalues. */
00393 
00394                 klnew = kl;
00395                 i__2 = kl;
00396                 for (ji = kf; ji <= i__2; ++ji) {
00397 
00398 /*                 Insure that N(w) is monotone   
00399 
00400    Computing MIN   
00401    Computing MAX */
00402                     i__5 = nab_ref(ji, 1), i__6 = iwork[ji];
00403                     i__3 = nab_ref(ji, 2), i__4 = maxMACRO(i__5,i__6);
00404                     iwork[ji] = minMACRO(i__3,i__4);
00405 
00406 /*                 Update the Queue -- add intervals if both halves   
00407                    contain eigenvalues. */
00408 
00409                     if (iwork[ji] == nab_ref(ji, 2)) {
00410 
00411 /*                    No eigenvalue in the upper interval:   
00412                       just use the lower interval. */
00413 
00414                         ab_ref(ji, 2) = c__[ji];
00415 
00416                     } else if (iwork[ji] == nab_ref(ji, 1)) {
00417 
00418 /*                    No eigenvalue in the lower interval:   
00419                       just use the upper interval. */
00420 
00421                         ab_ref(ji, 1) = c__[ji];
00422                     } else {
00423                         ++klnew;
00424                         if (klnew <= *mmax) {
00425 
00426 /*                       Eigenvalue in both intervals -- add upper to   
00427                          queue. */
00428 
00429                             ab_ref(klnew, 2) = ab_ref(ji, 2);
00430                             nab_ref(klnew, 2) = nab_ref(ji, 2);
00431                             ab_ref(klnew, 1) = c__[ji];
00432                             nab_ref(klnew, 1) = iwork[ji];
00433                             ab_ref(ji, 2) = c__[ji];
00434                             nab_ref(ji, 2) = iwork[ji];
00435                         } else {
00436                             *info = *mmax + 1;
00437                         }
00438                     }
00439 /* L70: */
00440                 }
00441                 if (*info != 0) {
00442                     return 0;
00443                 }
00444                 kl = klnew;
00445             } else {
00446 
00447 /*              IJOB=3: Binary search.  Keep only the interval containing   
00448                         w   s.t. N(w) = NVAL */
00449 
00450                 i__2 = kl;
00451                 for (ji = kf; ji <= i__2; ++ji) {
00452                     if (iwork[ji] <= nval[ji]) {
00453                         ab_ref(ji, 1) = c__[ji];
00454                         nab_ref(ji, 1) = iwork[ji];
00455                     }
00456                     if (iwork[ji] >= nval[ji]) {
00457                         ab_ref(ji, 2) = c__[ji];
00458                         nab_ref(ji, 2) = iwork[ji];
00459                     }
00460 /* L80: */
00461                 }
00462             }
00463 
00464         } else {
00465 
00466 /*           End of Parallel Version of the loop   
00467 
00468              Begin of Serial Version of the loop */
00469 
00470             klnew = kl;
00471             i__2 = kl;
00472             for (ji = kf; ji <= i__2; ++ji) {
00473 
00474 /*              Compute N(w), the number of eigenvalues less than w */
00475 
00476                 tmp1 = c__[ji];
00477                 tmp2 = d__[1] - tmp1;
00478                 itmp1 = 0;
00479                 if (tmp2 <= *pivmin) {
00480                     itmp1 = 1;
00481 /* Computing MIN */
00482                     d__1 = tmp2, d__2 = -(*pivmin);
00483                     tmp2 = minMACRO(d__1,d__2);
00484                 }
00485 
00486 /*              A series of compiler directives to defeat vectorization   
00487                 for the next loop   
00488 
00489    $PL$ CMCHAR=' '   
00490    DIR$          NEXTSCALAR   
00491    $DIR          SCALAR   
00492    DIR$          NEXT SCALAR   
00493    VD$L          NOVECTOR   
00494    DEC$          NOVECTOR   
00495    VD$           NOVECTOR   
00496    VDIR          NOVECTOR   
00497    VOCL          LOOP,SCALAR   
00498    IBM           PREFER SCALAR   
00499    $PL$ CMCHAR='*' */
00500 
00501                 i__3 = *n;
00502                 for (j = 2; j <= i__3; ++j) {
00503                     tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
00504                     if (tmp2 <= *pivmin) {
00505                         ++itmp1;
00506 /* Computing MIN */
00507                         d__1 = tmp2, d__2 = -(*pivmin);
00508                         tmp2 = minMACRO(d__1,d__2);
00509                     }
00510 /* L90: */
00511                 }
00512 
00513                 if (*ijob <= 2) {
00514 
00515 /*                 IJOB=2: Choose all intervals containing eigenvalues.   
00516 
00517                    Insure that N(w) is monotone   
00518 
00519    Computing MIN   
00520    Computing MAX */
00521                     i__5 = nab_ref(ji, 1);
00522                     i__3 = nab_ref(ji, 2), i__4 = maxMACRO(i__5,itmp1);
00523                     itmp1 = minMACRO(i__3,i__4);
00524 
00525 /*                 Update the Queue -- add intervals if both halves   
00526                    contain eigenvalues. */
00527 
00528                     if (itmp1 == nab_ref(ji, 2)) {
00529 
00530 /*                    No eigenvalue in the upper interval:   
00531                       just use the lower interval. */
00532 
00533                         ab_ref(ji, 2) = tmp1;
00534 
00535                     } else if (itmp1 == nab_ref(ji, 1)) {
00536 
00537 /*                    No eigenvalue in the lower interval:   
00538                       just use the upper interval. */
00539 
00540                         ab_ref(ji, 1) = tmp1;
00541                     } else if (klnew < *mmax) {
00542 
00543 /*                    Eigenvalue in both intervals -- add upper to queue. */
00544 
00545                         ++klnew;
00546                         ab_ref(klnew, 2) = ab_ref(ji, 2);
00547                         nab_ref(klnew, 2) = nab_ref(ji, 2);
00548                         ab_ref(klnew, 1) = tmp1;
00549                         nab_ref(klnew, 1) = itmp1;
00550                         ab_ref(ji, 2) = tmp1;
00551                         nab_ref(ji, 2) = itmp1;
00552                     } else {
00553                         *info = *mmax + 1;
00554                         return 0;
00555                     }
00556                 } else {
00557 
00558 /*                 IJOB=3: Binary search.  Keep only the interval   
00559                            containing  w  s.t. N(w) = NVAL */
00560 
00561                     if (itmp1 <= nval[ji]) {
00562                         ab_ref(ji, 1) = tmp1;
00563                         nab_ref(ji, 1) = itmp1;
00564                     }
00565                     if (itmp1 >= nval[ji]) {
00566                         ab_ref(ji, 2) = tmp1;
00567                         nab_ref(ji, 2) = itmp1;
00568                     }
00569                 }
00570 /* L100: */
00571             }
00572             kl = klnew;
00573 
00574 /*           End of Serial Version of the loop */
00575 
00576         }
00577 
00578 /*        Check for convergence */
00579 
00580         kfnew = kf;
00581         i__2 = kl;
00582         for (ji = kf; ji <= i__2; ++ji) {
00583             tmp1 = (d__1 = ab_ref(ji, 2) - ab_ref(ji, 1), absMACRO(d__1));
00584 /* Computing MAX */
00585             d__3 = (d__1 = ab_ref(ji, 2), absMACRO(d__1)), d__4 = (d__2 = ab_ref(
00586                     ji, 1), absMACRO(d__2));
00587             tmp2 = maxMACRO(d__3,d__4);
00588 /* Computing MAX */
00589             d__1 = maxMACRO(*abstol,*pivmin), d__2 = *reltol * tmp2;
00590             if (tmp1 < maxMACRO(d__1,d__2) || nab_ref(ji, 1) >= nab_ref(ji, 2)) {
00591 
00592 /*              Converged -- Swap with position KFNEW,   
00593                              then increment KFNEW */
00594 
00595                 if (ji > kfnew) {
00596                     tmp1 = ab_ref(ji, 1);
00597                     tmp2 = ab_ref(ji, 2);
00598                     itmp1 = nab_ref(ji, 1);
00599                     itmp2 = nab_ref(ji, 2);
00600                     ab_ref(ji, 1) = ab_ref(kfnew, 1);
00601                     ab_ref(ji, 2) = ab_ref(kfnew, 2);
00602                     nab_ref(ji, 1) = nab_ref(kfnew, 1);
00603                     nab_ref(ji, 2) = nab_ref(kfnew, 2);
00604                     ab_ref(kfnew, 1) = tmp1;
00605                     ab_ref(kfnew, 2) = tmp2;
00606                     nab_ref(kfnew, 1) = itmp1;
00607                     nab_ref(kfnew, 2) = itmp2;
00608                     if (*ijob == 3) {
00609                         itmp1 = nval[ji];
00610                         nval[ji] = nval[kfnew];
00611                         nval[kfnew] = itmp1;
00612                     }
00613                 }
00614                 ++kfnew;
00615             }
00616 /* L110: */
00617         }
00618         kf = kfnew;
00619 
00620 /*        Choose Midpoints */
00621 
00622         i__2 = kl;
00623         for (ji = kf; ji <= i__2; ++ji) {
00624             c__[ji] = (ab_ref(ji, 1) + ab_ref(ji, 2)) * .5;
00625 /* L120: */
00626         }
00627 
00628 /*        If no more intervals to refine, quit. */
00629 
00630         if (kf > kl) {
00631             goto L140;
00632         }
00633 /* L130: */
00634     }
00635 
00636 /*     Converged */
00637 
00638 L140:
00639 /* Computing MAX */
00640     i__1 = kl + 1 - kf;
00641     *info = maxMACRO(i__1,0);
00642     *mout = kl;
00643 
00644     return 0;
00645 
00646 /*     End of DLAEBZ */
00647 
00648 } /* dlaebz_ */
00649 
00650 #undef nab_ref
00651 #undef ab_ref
00652 
00653 
00654 #endif

Generated on Mon Sep 17 14:32:56 2012 for ergo by  doxygen 1.4.7