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_LARRR_HEADER 00038 #define TEMPLATE_LAPACK_LARRR_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, 00042 integer *info) 00043 { 00044 /* System generated locals */ 00045 integer i__1; 00046 Treal d__1; 00047 00048 00049 /* Local variables */ 00050 integer i__; 00051 Treal eps, tmp, tmp2, rmin; 00052 Treal offdig, safmin; 00053 logical yesrel; 00054 Treal smlnum, offdig2; 00055 00056 00057 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00058 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00059 /* November 2006 */ 00060 00061 /* .. Scalar Arguments .. */ 00062 /* .. */ 00063 /* .. Array Arguments .. */ 00064 /* .. */ 00065 00066 00067 /* Purpose */ 00068 /* ======= */ 00069 00070 /* Perform tests to decide whether the symmetric tridiagonal matrix T */ 00071 /* warrants expensive computations which guarantee high relative accuracy */ 00072 /* in the eigenvalues. */ 00073 00074 /* Arguments */ 00075 /* ========= */ 00076 00077 /* N (input) INTEGER */ 00078 /* The order of the matrix. N > 0. */ 00079 00080 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00081 /* The N diagonal elements of the tridiagonal matrix T. */ 00082 00083 /* E (input/output) DOUBLE PRECISION array, dimension (N) */ 00084 /* On entry, the first (N-1) entries contain the subdiagonal */ 00085 /* elements of the tridiagonal matrix T; E(N) is set to ZERO. */ 00086 00087 /* INFO (output) INTEGER */ 00088 /* INFO = 0(default) : the matrix warrants computations preserving */ 00089 /* relative accuracy. */ 00090 /* INFO = 1 : the matrix warrants computations guaranteeing */ 00091 /* only absolute accuracy. */ 00092 00093 /* Further Details */ 00094 /* =============== */ 00095 00096 /* Based on contributions by */ 00097 /* Beresford Parlett, University of California, Berkeley, USA */ 00098 /* Jim Demmel, University of California, Berkeley, USA */ 00099 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00100 /* Osni Marques, LBNL/NERSC, USA */ 00101 /* Christof Voemel, University of California, Berkeley, USA */ 00102 00103 /* ===================================================================== */ 00104 00105 /* .. Parameters .. */ 00106 /* .. */ 00107 /* .. Local Scalars .. */ 00108 /* .. */ 00109 /* .. External Functions .. */ 00110 /* .. */ 00111 /* .. Intrinsic Functions .. */ 00112 /* .. */ 00113 /* .. Executable Statements .. */ 00114 00115 /* As a default, do NOT go for relative-accuracy preserving computations. */ 00116 /* Parameter adjustments */ 00117 --e; 00118 --d__; 00119 00120 /* Function Body */ 00121 *info = 1; 00122 safmin = template_lapack_lamch("Safe minimum", (Treal)0); 00123 eps = template_lapack_lamch("Precision", (Treal)0); 00124 smlnum = safmin / eps; 00125 rmin = template_blas_sqrt(smlnum); 00126 /* Tests for relative accuracy */ 00127 00128 /* Test for scaled diagonal dominance */ 00129 /* Scale the diagonal entries to one and check whether the sum of the */ 00130 /* off-diagonals is less than one */ 00131 00132 /* The sdd relative error bounds have a 1/(1- 2*x) factor in them, */ 00133 /* x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */ 00134 /* accuracy is promised. In the notation of the code fragment below, */ 00135 /* 1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */ 00136 /* We don't think it is worth going into "sdd mode" unless the relative */ 00137 /* condition number is reasonable, not 1/macheps. */ 00138 /* The threshold should be compatible with other thresholds used in the */ 00139 /* code. We set OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */ 00140 /* to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */ 00141 /* instead of the current OFFDIG + OFFDIG2 < 1 */ 00142 00143 yesrel = TRUE_; 00144 offdig = 0.; 00145 tmp = template_blas_sqrt((absMACRO(d__[1]))); 00146 if (tmp < rmin) { 00147 yesrel = FALSE_; 00148 } 00149 if (! yesrel) { 00150 goto L11; 00151 } 00152 i__1 = *n; 00153 for (i__ = 2; i__ <= i__1; ++i__) { 00154 tmp2 = template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1))); 00155 if (tmp2 < rmin) { 00156 yesrel = FALSE_; 00157 } 00158 if (! yesrel) { 00159 goto L11; 00160 } 00161 offdig2 = (d__1 = e[i__ - 1], absMACRO(d__1)) / (tmp * tmp2); 00162 if (offdig + offdig2 >= .999) { 00163 yesrel = FALSE_; 00164 } 00165 if (! yesrel) { 00166 goto L11; 00167 } 00168 tmp = tmp2; 00169 offdig = offdig2; 00170 /* L10: */ 00171 } 00172 L11: 00173 if (yesrel) { 00174 *info = 0; 00175 return 0; 00176 } else { 00177 } 00178 00179 00180 /* *** MORE TO BE IMPLEMENTED *** */ 00181 00182 00183 /* Test if the lower bidiagonal matrix L from T = L D L^T */ 00184 /* (zero shift facto) is well conditioned */ 00185 00186 00187 /* Test if the upper bidiagonal matrix U from T = U D U^T */ 00188 /* (zero shift facto) is well conditioned. */ 00189 /* In this case, the matrix needs to be flipped and, at the end */ 00190 /* of the eigenvector computation, the flip needs to be applied */ 00191 /* to the computed eigenvectors (and the support) */ 00192 00193 00194 return 0; 00195 00196 /* END OF DLARRR */ 00197 00198 } /* dlarrr_ */ 00199 00200 #endif