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_LANEG_HEADER 00038 #define TEMPLATE_LAPACK_LANEG_HEADER 00039 00040 00041 #include "template_lapack_isnan.h" 00042 00043 00044 template<class Treal> 00045 int template_lapack_laneg(integer *n, Treal *d__, Treal *lld, Treal * 00046 sigma, Treal *pivmin, integer *r__) 00047 { 00048 /* System generated locals */ 00049 integer ret_val, i__1, i__2, i__3, i__4; 00050 00051 /* Local variables */ 00052 integer j; 00053 Treal p, t; 00054 integer bj; 00055 Treal tmp; 00056 integer neg1, neg2; 00057 Treal bsav, gamma, dplus; 00058 integer negcnt; 00059 logical sawnan; 00060 Treal dminus; 00061 00062 00063 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00064 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00065 /* November 2006 */ 00066 00067 /* .. Scalar Arguments .. */ 00068 /* .. */ 00069 /* .. Array Arguments .. */ 00070 /* .. */ 00071 00072 /* Purpose */ 00073 /* ======= */ 00074 00075 /* DLANEG computes the Sturm count, the number of negative pivots */ 00076 /* encountered while factoring tridiagonal T - sigma I = L D L^T. */ 00077 /* This implementation works directly on the factors without forming */ 00078 /* the tridiagonal matrix T. The Sturm count is also the number of */ 00079 /* eigenvalues of T less than sigma. */ 00080 00081 /* This routine is called from DLARRB. */ 00082 00083 /* The current routine does not use the PIVMIN parameter but rather */ 00084 /* requires IEEE-754 propagation of Infinities and NaNs. This */ 00085 /* routine also has no input range restrictions but does require */ 00086 /* default exception handling such that x/0 produces Inf when x is */ 00087 /* non-zero, and Inf/Inf produces NaN. For more information, see: */ 00088 00089 /* Marques, Riedy, and Voemel, "Benefits of IEEE-754 Features in */ 00090 /* Modern Symmetric Tridiagonal Eigensolvers," SIAM Journal on */ 00091 /* Scientific Computing, v28, n5, 2006. DOI 10.1137/050641624 */ 00092 /* (Tech report version in LAWN 172 with the same title.) */ 00093 00094 /* Arguments */ 00095 /* ========= */ 00096 00097 /* N (input) INTEGER */ 00098 /* The order of the matrix. */ 00099 00100 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00101 /* The N diagonal elements of the diagonal matrix D. */ 00102 00103 /* LLD (input) DOUBLE PRECISION array, dimension (N-1) */ 00104 /* The (N-1) elements L(i)*L(i)*D(i). */ 00105 00106 /* SIGMA (input) DOUBLE PRECISION */ 00107 /* Shift amount in T - sigma I = L D L^T. */ 00108 00109 /* PIVMIN (input) DOUBLE PRECISION */ 00110 /* The minimum pivot in the Sturm sequence. May be used */ 00111 /* when zero pivots are encountered on non-IEEE-754 */ 00112 /* architectures. */ 00113 00114 /* R (input) INTEGER */ 00115 /* The twist index for the twisted factorization that is used */ 00116 /* for the negcount. */ 00117 00118 /* Further Details */ 00119 /* =============== */ 00120 00121 /* Based on contributions by */ 00122 /* Osni Marques, LBNL/NERSC, USA */ 00123 /* Christof Voemel, University of California, Berkeley, USA */ 00124 /* Jason Riedy, University of California, Berkeley, USA */ 00125 00126 /* ===================================================================== */ 00127 00128 /* .. Parameters .. */ 00129 /* Some architectures propagate Infinities and NaNs very slowly, so */ 00130 /* the code computes counts in BLKLEN chunks. Then a NaN can */ 00131 /* propagate at most BLKLEN columns before being detected. This is */ 00132 /* not a general tuning parameter; it needs only to be just large */ 00133 /* enough that the overhead is tiny in common cases. */ 00134 /* .. */ 00135 /* .. Local Scalars .. */ 00136 /* .. */ 00137 /* .. Intrinsic Functions .. */ 00138 /* .. */ 00139 /* .. External Functions .. */ 00140 /* .. */ 00141 /* .. Executable Statements .. */ 00142 /* Parameter adjustments */ 00143 --lld; 00144 --d__; 00145 00146 /* Function Body */ 00147 negcnt = 0; 00148 /* I) upper part: L D L^T - SIGMA I = L+ D+ L+^T */ 00149 t = -(*sigma); 00150 i__1 = *r__ - 1; 00151 for (bj = 1; bj <= i__1; bj += 128) { 00152 neg1 = 0; 00153 bsav = t; 00154 /* Computing MIN */ 00155 i__3 = bj + 127, i__4 = *r__ - 1; 00156 i__2 = minMACRO(i__3,i__4); 00157 for (j = bj; j <= i__2; ++j) { 00158 dplus = d__[j] + t; 00159 if (dplus < 0.) { 00160 ++neg1; 00161 } 00162 tmp = t / dplus; 00163 t = tmp * lld[j] - *sigma; 00164 /* L21: */ 00165 } 00166 sawnan = template_lapack_isnan(&t); 00167 /* Run a slower version of the above loop if a NaN is detected. */ 00168 /* A NaN should occur only with a zero pivot after an infinite */ 00169 /* pivot. In that case, substituting 1 for T/DPLUS is the */ 00170 /* correct limit. */ 00171 if (sawnan) { 00172 neg1 = 0; 00173 t = bsav; 00174 /* Computing MIN */ 00175 i__3 = bj + 127, i__4 = *r__ - 1; 00176 i__2 = minMACRO(i__3,i__4); 00177 for (j = bj; j <= i__2; ++j) { 00178 dplus = d__[j] + t; 00179 if (dplus < 0.) { 00180 ++neg1; 00181 } 00182 tmp = t / dplus; 00183 if (template_lapack_isnan(&tmp)) { 00184 tmp = 1.; 00185 } 00186 t = tmp * lld[j] - *sigma; 00187 /* L22: */ 00188 } 00189 } 00190 negcnt += neg1; 00191 /* L210: */ 00192 } 00193 00194 /* II) lower part: L D L^T - SIGMA I = U- D- U-^T */ 00195 p = d__[*n] - *sigma; 00196 i__1 = *r__; 00197 for (bj = *n - 1; bj >= i__1; bj += -128) { 00198 neg2 = 0; 00199 bsav = p; 00200 /* Computing MAX */ 00201 i__3 = bj - 127; 00202 i__2 = maxMACRO(i__3,*r__); 00203 for (j = bj; j >= i__2; --j) { 00204 dminus = lld[j] + p; 00205 if (dminus < 0.) { 00206 ++neg2; 00207 } 00208 tmp = p / dminus; 00209 p = tmp * d__[j] - *sigma; 00210 /* L23: */ 00211 } 00212 sawnan = template_lapack_isnan(&p); 00213 /* As above, run a slower version that substitutes 1 for Inf/Inf. */ 00214 00215 if (sawnan) { 00216 neg2 = 0; 00217 p = bsav; 00218 /* Computing MAX */ 00219 i__3 = bj - 127; 00220 i__2 = maxMACRO(i__3,*r__); 00221 for (j = bj; j >= i__2; --j) { 00222 dminus = lld[j] + p; 00223 if (dminus < 0.) { 00224 ++neg2; 00225 } 00226 tmp = p / dminus; 00227 if (template_lapack_isnan(&tmp)) { 00228 tmp = 1.; 00229 } 00230 p = tmp * d__[j] - *sigma; 00231 /* L24: */ 00232 } 00233 } 00234 negcnt += neg2; 00235 /* L230: */ 00236 } 00237 00238 /* III) Twist index */ 00239 /* T was shifted by SIGMA initially. */ 00240 gamma = t + *sigma + p; 00241 if (gamma < 0.) { 00242 ++negcnt; 00243 } 00244 ret_val = negcnt; 00245 return ret_val; 00246 } /* dlaneg_ */ 00247 00248 #endif