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_LARRC_HEADER 00038 #define TEMPLATE_LAPACK_LARRC_HEADER 00039 00040 template<class Treal> 00041 int template_lapack_larrc(const char *jobt, const integer *n, const Treal *vl, 00042 const Treal *vu, Treal *d__, Treal *e, Treal *pivmin, 00043 integer *eigcnt, integer *lcnt, integer *rcnt, integer *info) 00044 { 00045 /* System generated locals */ 00046 integer i__1; 00047 Treal d__1; 00048 00049 /* Local variables */ 00050 integer i__; 00051 Treal sl, su, tmp, tmp2; 00052 logical matt; 00053 Treal lpivot, rpivot; 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 /* Find the number of eigenvalues of the symmetric tridiagonal matrix T */ 00069 /* that are in the interval (VL,VU] if JOBT = 'T', and of L D L^T */ 00070 /* if JOBT = 'L'. */ 00071 00072 /* Arguments */ 00073 /* ========= */ 00074 00075 /* JOBT (input) CHARACTER*1 */ 00076 /* = 'T': Compute Sturm count for matrix T. */ 00077 /* = 'L': Compute Sturm count for matrix L D L^T. */ 00078 00079 /* N (input) INTEGER */ 00080 /* The order of the matrix. N > 0. */ 00081 00082 /* VL (input) DOUBLE PRECISION */ 00083 /* VU (input) DOUBLE PRECISION */ 00084 /* The lower and upper bounds for the eigenvalues. */ 00085 00086 /* D (input) DOUBLE PRECISION array, dimension (N) */ 00087 /* JOBT = 'T': The N diagonal elements of the tridiagonal matrix T. */ 00088 /* JOBT = 'L': The N diagonal elements of the diagonal matrix D. */ 00089 00090 /* E (input) DOUBLE PRECISION array, dimension (N) */ 00091 /* JOBT = 'T': The N-1 offdiagonal elements of the matrix T. */ 00092 /* JOBT = 'L': The N-1 offdiagonal elements of the matrix L. */ 00093 00094 /* PIVMIN (input) DOUBLE PRECISION */ 00095 /* The minimum pivot in the Sturm sequence for T. */ 00096 00097 /* EIGCNT (output) INTEGER */ 00098 /* The number of eigenvalues of the symmetric tridiagonal matrix T */ 00099 /* that are in the interval (VL,VU] */ 00100 00101 /* LCNT (output) INTEGER */ 00102 /* RCNT (output) INTEGER */ 00103 /* The left and right negcounts of the interval. */ 00104 00105 /* INFO (output) INTEGER */ 00106 00107 /* Further Details */ 00108 /* =============== */ 00109 00110 /* Based on contributions by */ 00111 /* Beresford Parlett, University of California, Berkeley, USA */ 00112 /* Jim Demmel, University of California, Berkeley, USA */ 00113 /* Inderjit Dhillon, University of Texas, Austin, USA */ 00114 /* Osni Marques, LBNL/NERSC, USA */ 00115 /* Christof Voemel, University of California, Berkeley, USA */ 00116 00117 /* ===================================================================== */ 00118 00119 /* .. Parameters .. */ 00120 /* .. */ 00121 /* .. Local Scalars .. */ 00122 /* .. */ 00123 /* .. External Functions .. */ 00124 /* .. */ 00125 /* .. Executable Statements .. */ 00126 00127 /* Parameter adjustments */ 00128 --e; 00129 --d__; 00130 00131 /* Function Body */ 00132 *info = 0; 00133 *lcnt = 0; 00134 *rcnt = 0; 00135 *eigcnt = 0; 00136 matt = template_blas_lsame(jobt, "T"); 00137 if (matt) { 00138 /* Sturm sequence count on T */ 00139 lpivot = d__[1] - *vl; 00140 rpivot = d__[1] - *vu; 00141 if (lpivot <= 0.) { 00142 ++(*lcnt); 00143 } 00144 if (rpivot <= 0.) { 00145 ++(*rcnt); 00146 } 00147 i__1 = *n - 1; 00148 for (i__ = 1; i__ <= i__1; ++i__) { 00149 /* Computing 2nd power */ 00150 d__1 = e[i__]; 00151 tmp = d__1 * d__1; 00152 lpivot = d__[i__ + 1] - *vl - tmp / lpivot; 00153 rpivot = d__[i__ + 1] - *vu - tmp / rpivot; 00154 if (lpivot <= 0.) { 00155 ++(*lcnt); 00156 } 00157 if (rpivot <= 0.) { 00158 ++(*rcnt); 00159 } 00160 /* L10: */ 00161 } 00162 } else { 00163 /* Sturm sequence count on L D L^T */ 00164 sl = -(*vl); 00165 su = -(*vu); 00166 i__1 = *n - 1; 00167 for (i__ = 1; i__ <= i__1; ++i__) { 00168 lpivot = d__[i__] + sl; 00169 rpivot = d__[i__] + su; 00170 if (lpivot <= 0.) { 00171 ++(*lcnt); 00172 } 00173 if (rpivot <= 0.) { 00174 ++(*rcnt); 00175 } 00176 tmp = e[i__] * d__[i__] * e[i__]; 00177 00178 tmp2 = tmp / lpivot; 00179 if (tmp2 == 0.) { 00180 sl = tmp - *vl; 00181 } else { 00182 sl = sl * tmp2 - *vl; 00183 } 00184 00185 tmp2 = tmp / rpivot; 00186 if (tmp2 == 0.) { 00187 su = tmp - *vu; 00188 } else { 00189 su = su * tmp2 - *vu; 00190 } 00191 /* L20: */ 00192 } 00193 lpivot = d__[*n] + sl; 00194 rpivot = d__[*n] + su; 00195 if (lpivot <= 0.) { 00196 ++(*lcnt); 00197 } 00198 if (rpivot <= 0.) { 00199 ++(*rcnt); 00200 } 00201 } 00202 *eigcnt = *rcnt - *lcnt; 00203 return 0; 00204 00205 /* end of DLARRC */ 00206 00207 } /* dlarrc_ */ 00208 00209 #endif