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