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