00001 /* Ergo, version 3.7, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2018 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, 00004 * and Anastasia Kruchinina. 00005 * 00006 * This program is free software: you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation, either version 3 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00018 * 00019 * Primary academic reference: 00020 * Ergo: An open-source program for linear-scaling electronic structure 00021 * calculations, 00022 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia 00023 * Kruchinina, 00024 * SoftwareX 7, 107 (2018), 00025 * <http://dx.doi.org/10.1016/j.softx.2018.03.005> 00026 * 00027 * For further information about Ergo, see <http://www.ergoscf.org>. 00028 */ 00029 00030 /* This file belongs to the template_lapack part of the Ergo source 00031 * code. The source files in the template_lapack directory are modified 00032 * versions of files originally distributed as CLAPACK, see the 00033 * Copyright/license notice in the file template_lapack/COPYING. 00034 */ 00035 00036 00037 #ifndef TEMPLATE_LAPACK_LARRK_HEADER 00038 #define TEMPLATE_LAPACK_LARRK_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_larrk(integer *n, integer *iw, Treal *gl, 00042 Treal *gu, Treal *d__, Treal *e2, Treal *pivmin, 00043 Treal *reltol, Treal *w, Treal *werr, integer *info) 00044 { 00045 /* System generated locals */ 00046 integer i__1; 00047 Treal d__1, d__2; 00048 00049 00050 /* Local variables */ 00051 integer i__, it; 00052 Treal mid, eps, tmp1, tmp2, left, atoli, right; 00053 integer itmax; 00054 Treal rtoli, tnorm; 00055 integer negcnt; 00056 00057 00058 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00059 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00060 /* November 2006 */ 00061 00062 /* .. Scalar Arguments .. */ 00063 /* .. */ 00064 /* .. Array Arguments .. */ 00065 /* .. */ 00066 00067 /* Purpose */ 00068 /* ======= */ 00069 00070 /* DLARRK computes one eigenvalue of a symmetric tridiagonal */ 00071 /* matrix T to suitable accuracy. This is an auxiliary code to be */ 00072 /* called from DSTEMR. */ 00073 00074 /* To avoid overflow, the matrix must be scaled so that its */ 00075 /* largest element is no greater than overflow**(1/2) * */ 00076 /* underflow**(1/4) in absolute value, and for greatest */ 00077 /* accuracy, it should not be much smaller than that. */ 00078 00079 /* See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */ 00080 /* Matrix", Report CS41, Computer Science Dept., Stanford */ 00081 /* University, July 21, 1966. */ 00082 00083 /* Arguments */ 00084 /* ========= */ 00085 00086 /* N (input) INTEGER */ 00087 /* The order of the tridiagonal matrix T. N >= 0. */ 00088 00089 /* IW (input) INTEGER */ 00090 /* The index of the eigenvalues to be returned. */ 00091 00092 /* GL (input) DOUBLE PRECISION */ 00093 /* GU (input) DOUBLE PRECISION */ 00094 /* An upper and a lower bound on the eigenvalue. */ 00095 00096 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00097 /* The n diagonal elements of the tridiagonal matrix T. */ 00098 00099 /* E2 (input) DOUBLE PRECISION array, dimension (N-1) */ 00100 /* The (n-1) squared off-diagonal elements of the tridiagonal matrix T. */ 00101 00102 /* PIVMIN (input) DOUBLE PRECISION */ 00103 /* The minimum pivot allowed in the Sturm sequence for T. */ 00104 00105 /* RELTOL (input) DOUBLE PRECISION */ 00106 /* The minimum relative width of an interval. When an interval */ 00107 /* is narrower than RELTOL times the larger (in */ 00108 /* magnitude) endpoint, then it is considered to be */ 00109 /* sufficiently small, i.e., converged. Note: this should */ 00110 /* always be at least radix*machine epsilon. */ 00111 00112 /* W (output) DOUBLE PRECISION */ 00113 00114 /* WERR (output) DOUBLE PRECISION */ 00115 /* The error bound on the corresponding eigenvalue approximation */ 00116 /* in W. */ 00117 00118 /* INFO (output) INTEGER */ 00119 /* = 0: Eigenvalue converged */ 00120 /* = -1: Eigenvalue did NOT converge */ 00121 00122 /* Internal Parameters */ 00123 /* =================== */ 00124 00125 /* FUDGE DOUBLE PRECISION, default = 2 */ 00126 /* A "fudge factor" to widen the Gershgorin intervals. */ 00127 00128 /* ===================================================================== */ 00129 00130 /* .. Parameters .. */ 00131 /* .. */ 00132 /* .. Local Scalars .. */ 00133 /* .. */ 00134 /* .. External Functions .. */ 00135 /* .. */ 00136 /* .. Intrinsic Functions .. */ 00137 /* .. */ 00138 /* .. Executable Statements .. */ 00139 00140 /* Get machine constants */ 00141 /* Parameter adjustments */ 00142 --e2; 00143 --d__; 00144 00145 /* Function Body */ 00146 eps = template_lapack_lamch("P", (Treal)0); 00147 /* Computing MAX */ 00148 d__1 = absMACRO(*gl), d__2 = absMACRO(*gu); 00149 tnorm = maxMACRO(d__1,d__2); 00150 rtoli = *reltol; 00151 atoli = *pivmin * 4.; 00152 itmax = (integer) ((template_blas_log(tnorm + *pivmin) - template_blas_log(*pivmin)) / template_blas_log(2.)) + 2; 00153 *info = -1; 00154 left = *gl - tnorm * 2. * eps * *n - *pivmin * 4.; 00155 right = *gu + tnorm * 2. * eps * *n + *pivmin * 4.; 00156 it = 0; 00157 L10: 00158 00159 /* Check if interval converged or maximum number of iterations reached */ 00160 00161 tmp1 = (d__1 = right - left, absMACRO(d__1)); 00162 /* Computing MAX */ 00163 d__1 = absMACRO(right), d__2 = absMACRO(left); 00164 tmp2 = maxMACRO(d__1,d__2); 00165 /* Computing MAX */ 00166 d__1 = maxMACRO(atoli,*pivmin), d__2 = rtoli * tmp2; 00167 if (tmp1 < maxMACRO(d__1,d__2)) { 00168 *info = 0; 00169 goto L30; 00170 } 00171 if (it > itmax) { 00172 goto L30; 00173 } 00174 00175 /* Count number of negative pivots for mid-point */ 00176 00177 ++it; 00178 mid = (left + right) * .5; 00179 negcnt = 0; 00180 tmp1 = d__[1] - mid; 00181 if (absMACRO(tmp1) < *pivmin) { 00182 tmp1 = -(*pivmin); 00183 } 00184 if (tmp1 <= 0.) { 00185 ++negcnt; 00186 } 00187 00188 i__1 = *n; 00189 for (i__ = 2; i__ <= i__1; ++i__) { 00190 tmp1 = d__[i__] - e2[i__ - 1] / tmp1 - mid; 00191 if (absMACRO(tmp1) < *pivmin) { 00192 tmp1 = -(*pivmin); 00193 } 00194 if (tmp1 <= 0.) { 00195 ++negcnt; 00196 } 00197 /* L20: */ 00198 } 00199 if (negcnt >= *iw) { 00200 right = mid; 00201 } else { 00202 left = mid; 00203 } 00204 goto L10; 00205 L30: 00206 00207 /* Converged or maximum number of iterations reached */ 00208 00209 *w = (left + right) * .5; 00210 *werr = (d__1 = right - left, absMACRO(d__1)) * .5; 00211 return 0; 00212 00213 /* End of DLARRK */ 00214 00215 } /* dlarrk_ */ 00216 00217 #endif