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