ergo
template_lapack_sytd2.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_SYTD2_HEADER
38#define TEMPLATE_LAPACK_SYTD2_HEADER
39
41
42template<class Treal>
43int template_lapack_sytd2(const char *uplo, const integer *n, Treal *a, const integer *
44 lda, Treal *d__, Treal *e, Treal *tau, integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 October 31, 1992
50
51
52 Purpose
53 =======
54
55 DSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal
56 form T by an orthogonal similarity transformation: Q' * A * Q = T.
57
58 Arguments
59 =========
60
61 UPLO (input) CHARACTER*1
62 Specifies whether the upper or lower triangular part of the
63 symmetric matrix A is stored:
64 = 'U': Upper triangular
65 = 'L': Lower triangular
66
67 N (input) INTEGER
68 The order of the matrix A. N >= 0.
69
70 A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
71 On entry, the symmetric matrix A. If UPLO = 'U', the leading
72 n-by-n upper triangular part of A contains the upper
73 triangular part of the matrix A, and the strictly lower
74 triangular part of A is not referenced. If UPLO = 'L', the
75 leading n-by-n lower triangular part of A contains the lower
76 triangular part of the matrix A, and the strictly upper
77 triangular part of A is not referenced.
78 On exit, if UPLO = 'U', the diagonal and first superdiagonal
79 of A are overwritten by the corresponding elements of the
80 tridiagonal matrix T, and the elements above the first
81 superdiagonal, with the array TAU, represent the orthogonal
82 matrix Q as a product of elementary reflectors; if UPLO
83 = 'L', the diagonal and first subdiagonal of A are over-
84 written by the corresponding elements of the tridiagonal
85 matrix T, and the elements below the first subdiagonal, with
86 the array TAU, represent the orthogonal matrix Q as a product
87 of elementary reflectors. See Further Details.
88
89 LDA (input) INTEGER
90 The leading dimension of the array A. LDA >= max(1,N).
91
92 D (output) DOUBLE PRECISION array, dimension (N)
93 The diagonal elements of the tridiagonal matrix T:
94 D(i) = A(i,i).
95
96 E (output) DOUBLE PRECISION array, dimension (N-1)
97 The off-diagonal elements of the tridiagonal matrix T:
98 E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
99
100 TAU (output) DOUBLE PRECISION array, dimension (N-1)
101 The scalar factors of the elementary reflectors (see Further
102 Details).
103
104 INFO (output) INTEGER
105 = 0: successful exit
106 < 0: if INFO = -i, the i-th argument had an illegal value.
107
108 Further Details
109 ===============
110
111 If UPLO = 'U', the matrix Q is represented as a product of elementary
112 reflectors
113
114 Q = H(n-1) . . . H(2) H(1).
115
116 Each H(i) has the form
117
118 H(i) = I - tau * v * v'
119
120 where tau is a real scalar, and v is a real vector with
121 v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in
122 A(1:i-1,i+1), and tau in TAU(i).
123
124 If UPLO = 'L', the matrix Q is represented as a product of elementary
125 reflectors
126
127 Q = H(1) H(2) . . . H(n-1).
128
129 Each H(i) has the form
130
131 H(i) = I - tau * v * v'
132
133 where tau is a real scalar, and v is a real vector with
134 v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in A(i+2:n,i),
135 and tau in TAU(i).
136
137 The contents of A on exit are illustrated by the following examples
138 with n = 5:
139
140 if UPLO = 'U': if UPLO = 'L':
141
142 ( d e v2 v3 v4 ) ( d )
143 ( d e v3 v4 ) ( e d )
144 ( d e v4 ) ( v1 e d )
145 ( d e ) ( v1 v2 e d )
146 ( d ) ( v1 v2 v3 e d )
147
148 where d and e denote diagonal and off-diagonal elements of T, and vi
149 denotes an element of the vector defining H(i).
150
151 =====================================================================
152
153
154 Test the input parameters
155
156 Parameter adjustments */
157 /* Table of constant values */
158 integer c__1 = 1;
159 Treal c_b8 = 0.;
160 Treal c_b14 = -1.;
161
162 /* System generated locals */
163 integer a_dim1, a_offset, i__1, i__2, i__3;
164 /* Local variables */
165 Treal taui;
166 integer i__;
167 Treal alpha;
168 logical upper;
169#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
170
171
172 a_dim1 = *lda;
173 a_offset = 1 + a_dim1 * 1;
174 a -= a_offset;
175 --d__;
176 --e;
177 --tau;
178
179 /* Function Body */
180 *info = 0;
181 upper = template_blas_lsame(uplo, "U");
182 if (! upper && ! template_blas_lsame(uplo, "L")) {
183 *info = -1;
184 } else if (*n < 0) {
185 *info = -2;
186 } else if (*lda < maxMACRO(1,*n)) {
187 *info = -4;
188 }
189 if (*info != 0) {
190 i__1 = -(*info);
191 template_blas_erbla("SYTD2 ", &i__1);
192 return 0;
193 }
194
195/* Quick return if possible */
196
197 if (*n <= 0) {
198 return 0;
199 }
200
201 if (upper) {
202
203/* Reduce the upper triangle of A */
204
205 for (i__ = *n - 1; i__ >= 1; --i__) {
206
207/* Generate elementary reflector H(i) = I - tau * v * v'
208 to annihilate A(1:i-1,i+1) */
209
210 template_lapack_larfg(&i__, &a_ref(i__, i__ + 1), &a_ref(1, i__ + 1), &c__1, &
211 taui);
212 e[i__] = a_ref(i__, i__ + 1);
213
214 if (taui != 0.) {
215
216/* Apply H(i) from both sides to A(1:i,1:i) */
217
218 a_ref(i__, i__ + 1) = 1.;
219
220/* Compute x := tau * A * v storing x in TAU(1:i) */
221
222 template_blas_symv(uplo, &i__, &taui, &a[a_offset], lda, &a_ref(1, i__ +
223 1), &c__1, &c_b8, &tau[1], &c__1);
224
225/* Compute w := x - 1/2 * tau * (x'*v) * v */
226
227 alpha = taui * -.5 * template_blas_dot(&i__, &tau[1], &c__1, &a_ref(1,
228 i__ + 1), &c__1);
229 template_blas_axpy(&i__, &alpha, &a_ref(1, i__ + 1), &c__1, &tau[1], &
230 c__1);
231
232/* Apply the transformation as a rank-2 update:
233 A := A - v * w' - w * v' */
234
235 template_blas_syr2(uplo, &i__, &c_b14, &a_ref(1, i__ + 1), &c__1, &tau[1],
236 &c__1, &a[a_offset], lda);
237
238 a_ref(i__, i__ + 1) = e[i__];
239 }
240 d__[i__ + 1] = a_ref(i__ + 1, i__ + 1);
241 tau[i__] = taui;
242/* L10: */
243 }
244 d__[1] = a_ref(1, 1);
245 } else {
246
247/* Reduce the lower triangle of A */
248
249 i__1 = *n - 1;
250 for (i__ = 1; i__ <= i__1; ++i__) {
251
252/* Generate elementary reflector H(i) = I - tau * v * v'
253 to annihilate A(i+2:n,i)
254
255 Computing MIN */
256 i__2 = i__ + 2;
257 i__3 = *n - i__;
258 template_lapack_larfg(&i__3, &a_ref(i__ + 1, i__), &a_ref(minMACRO(i__2,*n), i__), &
259 c__1, &taui);
260 e[i__] = a_ref(i__ + 1, i__);
261
262 if (taui != 0.) {
263
264/* Apply H(i) from both sides to A(i+1:n,i+1:n) */
265
266 a_ref(i__ + 1, i__) = 1.;
267
268/* Compute x := tau * A * v storing y in TAU(i:n-1) */
269
270 i__2 = *n - i__;
271 template_blas_symv(uplo, &i__2, &taui, &a_ref(i__ + 1, i__ + 1), lda, &
272 a_ref(i__ + 1, i__), &c__1, &c_b8, &tau[i__], &c__1);
273
274/* Compute w := x - 1/2 * tau * (x'*v) * v */
275
276 i__2 = *n - i__;
277 alpha = taui * -.5 * template_blas_dot(&i__2, &tau[i__], &c__1, &a_ref(
278 i__ + 1, i__), &c__1);
279 i__2 = *n - i__;
280 template_blas_axpy(&i__2, &alpha, &a_ref(i__ + 1, i__), &c__1, &tau[i__],
281 &c__1);
282
283/* Apply the transformation as a rank-2 update:
284 A := A - v * w' - w * v' */
285
286 i__2 = *n - i__;
287 template_blas_syr2(uplo, &i__2, &c_b14, &a_ref(i__ + 1, i__), &c__1, &tau[
288 i__], &c__1, &a_ref(i__ + 1, i__ + 1), lda)
289 ;
290
291 a_ref(i__ + 1, i__) = e[i__];
292 }
293 d__[i__] = a_ref(i__, i__);
294 tau[i__] = taui;
295/* L20: */
296 }
297 d__[*n] = a_ref(*n, *n);
298 }
299
300 return 0;
301
302/* End of DSYTD2 */
303
304} /* dsytd2_ */
305
306#undef a_ref
307
308
309#endif
int template_blas_axpy(const integer *n, const Treal *da, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_axpy.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
Treal template_blas_dot(const integer *n, const Treal *dx, const integer *incx, const Treal *dy, const integer *incy)
Definition: template_blas_dot.h:43
int template_blas_symv(const char *uplo, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, const Treal *x, const integer *incx, const Treal *beta, Treal *y, const integer *incy)
Definition: template_blas_symv.h:42
int template_blas_syr2(const char *uplo, const integer *n, const Treal *alpha, const Treal *x, const integer *incx, const Treal *y, const integer *incy, Treal *a, const integer *lda)
Definition: template_blas_syr2.h:42
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43
#define a_ref(a_1, a_2)
int template_lapack_sytd2(const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *d__, Treal *e, Treal *tau, integer *info)
Definition: template_lapack_sytd2.h:43