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_SYTD2_HEADER 00036 #define TEMPLATE_LAPACK_SYTD2_HEADER 00037 00038 #include "template_lapack_common.h" 00039 00040 template<class Treal> 00041 int template_lapack_sytd2(const char *uplo, const integer *n, Treal *a, const integer * 00042 lda, Treal *d__, Treal *e, Treal *tau, integer *info) 00043 { 00044 /* -- LAPACK routine (version 3.0) -- 00045 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00046 Courant Institute, Argonne National Lab, and Rice University 00047 October 31, 1992 00048 00049 00050 Purpose 00051 ======= 00052 00053 DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal 00054 form T by an orthogonal similarity transformation: Q' * A * Q = T. 00055 00056 Arguments 00057 ========= 00058 00059 UPLO (input) CHARACTER*1 00060 Specifies whether the upper or lower triangular part of the 00061 symmetric matrix A is stored: 00062 = 'U': Upper triangular 00063 = 'L': Lower triangular 00064 00065 N (input) INTEGER 00066 The order of the matrix A. N >= 0. 00067 00068 A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 00069 On entry, the symmetric matrix A. If UPLO = 'U', the leading 00070 n-by-n upper triangular part of A contains the upper 00071 triangular part of the matrix A, and the strictly lower 00072 triangular part of A is not referenced. If UPLO = 'L', the 00073 leading n-by-n lower triangular part of A contains the lower 00074 triangular part of the matrix A, and the strictly upper 00075 triangular part of A is not referenced. 00076 On exit, if UPLO = 'U', the diagonal and first superdiagonal 00077 of A are overwritten by the corresponding elements of the 00078 tridiagonal matrix T, and the elements above the first 00079 superdiagonal, with the array TAU, represent the orthogonal 00080 matrix Q as a product of elementary reflectors; if UPLO 00081 = 'L', the diagonal and first subdiagonal of A are over- 00082 written by the corresponding elements of the tridiagonal 00083 matrix T, and the elements below the first subdiagonal, with 00084 the array TAU, represent the orthogonal matrix Q as a product 00085 of elementary reflectors. See Further Details. 00086 00087 LDA (input) INTEGER 00088 The leading dimension of the array A. LDA >= max(1,N). 00089 00090 D (output) DOUBLE PRECISION array, dimension (N) 00091 The diagonal elements of the tridiagonal matrix T: 00092 D(i) = A(i,i). 00093 00094 E (output) DOUBLE PRECISION array, dimension (N-1) 00095 The off-diagonal elements of the tridiagonal matrix T: 00096 E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 00097 00098 TAU (output) DOUBLE PRECISION array, dimension (N-1) 00099 The scalar factors of the elementary reflectors (see Further 00100 Details). 00101 00102 INFO (output) INTEGER 00103 = 0: successful exit 00104 < 0: if INFO = -i, the i-th argument had an illegal value. 00105 00106 Further Details 00107 =============== 00108 00109 If UPLO = 'U', the matrix Q is represented as a product of elementary 00110 reflectors 00111 00112 Q = H(n-1) . . . H(2) H(1). 00113 00114 Each H(i) has the form 00115 00116 H(i) = I - tau * v * v' 00117 00118 where tau is a real scalar, and v is a real vector with 00119 v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in 00120 A(1:i-1,i+1), and tau in TAU(i). 00121 00122 If UPLO = 'L', the matrix Q is represented as a product of elementary 00123 reflectors 00124 00125 Q = H(1) H(2) . . . H(n-1). 00126 00127 Each H(i) has the form 00128 00129 H(i) = I - tau * v * v' 00130 00131 where tau is a real scalar, and v is a real vector with 00132 v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i), 00133 and tau in TAU(i). 00134 00135 The contents of A on exit are illustrated by the following examples 00136 with n = 5: 00137 00138 if UPLO = 'U': if UPLO = 'L': 00139 00140 ( d e v2 v3 v4 ) ( d ) 00141 ( d e v3 v4 ) ( e d ) 00142 ( d e v4 ) ( v1 e d ) 00143 ( d e ) ( v1 v2 e d ) 00144 ( d ) ( v1 v2 v3 e d ) 00145 00146 where d and e denote diagonal and off-diagonal elements of T, and vi 00147 denotes an element of the vector defining H(i). 00148 00149 ===================================================================== 00150 00151 00152 Test the input parameters 00153 00154 Parameter adjustments */ 00155 /* Table of constant values */ 00156 integer c__1 = 1; 00157 Treal c_b8 = 0.; 00158 Treal c_b14 = -1.; 00159 00160 /* System generated locals */ 00161 integer a_dim1, a_offset, i__1, i__2, i__3; 00162 /* Local variables */ 00163 Treal taui; 00164 integer i__; 00165 Treal alpha; 00166 logical upper; 00167 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00168 00169 00170 a_dim1 = *lda; 00171 a_offset = 1 + a_dim1 * 1; 00172 a -= a_offset; 00173 --d__; 00174 --e; 00175 --tau; 00176 00177 /* Function Body */ 00178 *info = 0; 00179 upper = template_blas_lsame(uplo, "U"); 00180 if (! upper && ! template_blas_lsame(uplo, "L")) { 00181 *info = -1; 00182 } else if (*n < 0) { 00183 *info = -2; 00184 } else if (*lda < maxMACRO(1,*n)) { 00185 *info = -4; 00186 } 00187 if (*info != 0) { 00188 i__1 = -(*info); 00189 template_blas_erbla("SYTD2 ", &i__1); 00190 return 0; 00191 } 00192 00193 /* Quick return if possible */ 00194 00195 if (*n <= 0) { 00196 return 0; 00197 } 00198 00199 if (upper) { 00200 00201 /* Reduce the upper triangle of A */ 00202 00203 for (i__ = *n - 1; i__ >= 1; --i__) { 00204 00205 /* Generate elementary reflector H(i) = I - tau * v * v' 00206 to annihilate A(1:i-1,i+1) */ 00207 00208 template_lapack_larfg(&i__, &a_ref(i__, i__ + 1), &a_ref(1, i__ + 1), &c__1, & 00209 taui); 00210 e[i__] = a_ref(i__, i__ + 1); 00211 00212 if (taui != 0.) { 00213 00214 /* Apply H(i) from both sides to A(1:i,1:i) */ 00215 00216 a_ref(i__, i__ + 1) = 1.; 00217 00218 /* Compute x := tau * A * v storing x in TAU(1:i) */ 00219 00220 template_blas_symv(uplo, &i__, &taui, &a[a_offset], lda, &a_ref(1, i__ + 00221 1), &c__1, &c_b8, &tau[1], &c__1); 00222 00223 /* Compute w := x - 1/2 * tau * (x'*v) * v */ 00224 00225 alpha = taui * -.5 * template_blas_dot(&i__, &tau[1], &c__1, &a_ref(1, 00226 i__ + 1), &c__1); 00227 template_blas_axpy(&i__, &alpha, &a_ref(1, i__ + 1), &c__1, &tau[1], & 00228 c__1); 00229 00230 /* Apply the transformation as a rank-2 update: 00231 A := A - v * w' - w * v' */ 00232 00233 template_blas_syr2(uplo, &i__, &c_b14, &a_ref(1, i__ + 1), &c__1, &tau[1], 00234 &c__1, &a[a_offset], lda); 00235 00236 a_ref(i__, i__ + 1) = e[i__]; 00237 } 00238 d__[i__ + 1] = a_ref(i__ + 1, i__ + 1); 00239 tau[i__] = taui; 00240 /* L10: */ 00241 } 00242 d__[1] = a_ref(1, 1); 00243 } else { 00244 00245 /* Reduce the lower triangle of A */ 00246 00247 i__1 = *n - 1; 00248 for (i__ = 1; i__ <= i__1; ++i__) { 00249 00250 /* Generate elementary reflector H(i) = I - tau * v * v' 00251 to annihilate A(i+2:n,i) 00252 00253 Computing MIN */ 00254 i__2 = i__ + 2; 00255 i__3 = *n - i__; 00256 template_lapack_larfg(&i__3, &a_ref(i__ + 1, i__), &a_ref(minMACRO(i__2,*n), i__), & 00257 c__1, &taui); 00258 e[i__] = a_ref(i__ + 1, i__); 00259 00260 if (taui != 0.) { 00261 00262 /* Apply H(i) from both sides to A(i+1:n,i+1:n) */ 00263 00264 a_ref(i__ + 1, i__) = 1.; 00265 00266 /* Compute x := tau * A * v storing y in TAU(i:n-1) */ 00267 00268 i__2 = *n - i__; 00269 template_blas_symv(uplo, &i__2, &taui, &a_ref(i__ + 1, i__ + 1), lda, & 00270 a_ref(i__ + 1, i__), &c__1, &c_b8, &tau[i__], &c__1); 00271 00272 /* Compute w := x - 1/2 * tau * (x'*v) * v */ 00273 00274 i__2 = *n - i__; 00275 alpha = taui * -.5 * template_blas_dot(&i__2, &tau[i__], &c__1, &a_ref( 00276 i__ + 1, i__), &c__1); 00277 i__2 = *n - i__; 00278 template_blas_axpy(&i__2, &alpha, &a_ref(i__ + 1, i__), &c__1, &tau[i__], 00279 &c__1); 00280 00281 /* Apply the transformation as a rank-2 update: 00282 A := A - v * w' - w * v' */ 00283 00284 i__2 = *n - i__; 00285 template_blas_syr2(uplo, &i__2, &c_b14, &a_ref(i__ + 1, i__), &c__1, &tau[ 00286 i__], &c__1, &a_ref(i__ + 1, i__ + 1), lda) 00287 ; 00288 00289 a_ref(i__ + 1, i__) = e[i__]; 00290 } 00291 d__[i__] = a_ref(i__, i__); 00292 tau[i__] = taui; 00293 /* L20: */ 00294 } 00295 d__[*n] = a_ref(*n, *n); 00296 } 00297 00298 return 0; 00299 00300 /* End of DSYTD2 */ 00301 00302 } /* dsytd2_ */ 00303 00304 #undef a_ref 00305 00306 00307 #endif