template_lapack_larrk.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_LARRK_HEADER
00036 #define TEMPLATE_LAPACK_LARRK_HEADER
00037 
00038 template<class Treal>
00039 int template_lapack_larrk(integer *n, integer *iw, Treal *gl, 
00040         Treal *gu, Treal *d__, Treal *e2, Treal *pivmin, 
00041         Treal *reltol, Treal *w, Treal *werr, integer *info)
00042 {
00043     /* System generated locals */
00044     integer i__1;
00045     Treal d__1, d__2;
00046 
00047 
00048     /* Local variables */
00049     integer i__, it;
00050     Treal mid, eps, tmp1, tmp2, left, atoli, right;
00051     integer itmax;
00052     Treal rtoli, tnorm;
00053     integer negcnt;
00054 
00055 
00056 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00057 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00058 /*     November 2006 */
00059 
00060 /*     .. Scalar Arguments .. */
00061 /*     .. */
00062 /*     .. Array Arguments .. */
00063 /*     .. */
00064 
00065 /*  Purpose */
00066 /*  ======= */
00067 
00068 /*  DLARRK computes one eigenvalue of a symmetric tridiagonal */
00069 /*  matrix T to suitable accuracy. This is an auxiliary code to be */
00070 /*  called from DSTEMR. */
00071 
00072 /*  To avoid overflow, the matrix must be scaled so that its */
00073 /*  largest element is no greater than overflow**(1/2) * */
00074 /*  underflow**(1/4) in absolute value, and for greatest */
00075 /*  accuracy, it should not be much smaller than that. */
00076 
00077 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
00078 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
00079 /*  University, July 21, 1966. */
00080 
00081 /*  Arguments */
00082 /*  ========= */
00083 
00084 /*  N       (input) INTEGER */
00085 /*          The order of the tridiagonal matrix T.  N >= 0. */
00086 
00087 /*  IW      (input) INTEGER */
00088 /*          The index of the eigenvalues to be returned. */
00089 
00090 /*  GL      (input) DOUBLE PRECISION */
00091 /*  GU      (input) DOUBLE PRECISION */
00092 /*          An upper and a lower bound on the eigenvalue. */
00093 
00094 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
00095 /*          The n diagonal elements of the tridiagonal matrix T. */
00096 
00097 /*  E2      (input) DOUBLE PRECISION array, dimension (N-1) */
00098 /*          The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */
00099 
00100 /*  PIVMIN  (input) DOUBLE PRECISION */
00101 /*          The minimum pivot allowed in the Sturm sequence for T. */
00102 
00103 /*  RELTOL  (input) DOUBLE PRECISION */
00104 /*          The minimum relative width of an interval.  When an interval */
00105 /*          is narrower than RELTOL times the larger (in */
00106 /*          magnitude) endpoint, then it is considered to be */
00107 /*          sufficiently small, i.e., converged.  Note: this should */
00108 /*          always be at least radix*machine epsilon. */
00109 
00110 /*  W       (output) DOUBLE PRECISION */
00111 
00112 /*  WERR    (output) DOUBLE PRECISION */
00113 /*          The error bound on the corresponding eigenvalue approximation */
00114 /*          in W. */
00115 
00116 /*  INFO    (output) INTEGER */
00117 /*          = 0:       Eigenvalue converged */
00118 /*          = -1:      Eigenvalue did NOT converge */
00119 
00120 /*  Internal Parameters */
00121 /*  =================== */
00122 
00123 /*  FUDGE   DOUBLE PRECISION, default = 2 */
00124 /*          A "fudge factor" to widen the Gershgorin intervals. */
00125 
00126 /*  ===================================================================== */
00127 
00128 /*     .. Parameters .. */
00129 /*     .. */
00130 /*     .. Local Scalars .. */
00131 /*     .. */
00132 /*     .. External Functions .. */
00133 /*     .. */
00134 /*     .. Intrinsic Functions .. */
00135 /*     .. */
00136 /*     .. Executable Statements .. */
00137 
00138 /*     Get machine constants */
00139     /* Parameter adjustments */
00140     --e2;
00141     --d__;
00142 
00143     /* Function Body */
00144     eps = template_lapack_lamch("P", (Treal)0);
00145 /* Computing MAX */
00146     d__1 = absMACRO(*gl), d__2 = absMACRO(*gu);
00147     tnorm = maxMACRO(d__1,d__2);
00148     rtoli = *reltol;
00149     atoli = *pivmin * 4.;
00150     itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 2;
00151     *info = -1;
00152     left = *gl - tnorm * 2. * eps * *n - *pivmin * 4.;
00153     right = *gu + tnorm * 2. * eps * *n + *pivmin * 4.;
00154     it = 0;
00155 L10:
00156 
00157 /*     Check if interval converged or maximum number of iterations reached */
00158 
00159     tmp1 = (d__1 = right - left, absMACRO(d__1));
00160 /* Computing MAX */
00161     d__1 = absMACRO(right), d__2 = absMACRO(left);
00162     tmp2 = maxMACRO(d__1,d__2);
00163 /* Computing MAX */
00164     d__1 = maxMACRO(atoli,*pivmin), d__2 = rtoli * tmp2;
00165     if (tmp1 < maxMACRO(d__1,d__2)) {
00166         *info = 0;
00167         goto L30;
00168     }
00169     if (it > itmax) {
00170         goto L30;
00171     }
00172 
00173 /*     Count number of negative pivots for mid-point */
00174 
00175     ++it;
00176     mid = (left + right) * .5;
00177     negcnt = 0;
00178     tmp1 = d__[1] - mid;
00179     if (absMACRO(tmp1) < *pivmin) {
00180         tmp1 = -(*pivmin);
00181     }
00182     if (tmp1 <= 0.) {
00183         ++negcnt;
00184     }
00185 
00186     i__1 = *n;
00187     for (i__ = 2; i__ <= i__1; ++i__) {
00188         tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid;
00189         if (absMACRO(tmp1) < *pivmin) {
00190             tmp1 = -(*pivmin);
00191         }
00192         if (tmp1 <= 0.) {
00193             ++negcnt;
00194         }
00195 /* L20: */
00196     }
00197     if (negcnt >= *iw) {
00198         right = mid;
00199     } else {
00200         left = mid;
00201     }
00202     goto L10;
00203 L30:
00204 
00205 /*     Converged or maximum number of iterations reached */
00206 
00207     *w = (left + right) * .5;
00208     *werr = (d__1 = right - left, absMACRO(d__1)) * .5;
00209     return 0;
00210 
00211 /*     End of DLARRK */
00212 
00213 } /* dlarrk_ */
00214 
00215 #endif

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