ergo
template_lapack_sytd2.h
Go to the documentation of this file.
1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2019 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 
40 #include "template_lapack_common.h"
41 
42 template<class Treal>
43 int 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
template_blas_dot
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
template_lapack_larfg
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43
minMACRO
#define minMACRO(a, b)
Definition: template_blas_common.h:46
template_lapack_sytd2
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
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_common.h
template_blas_erbla
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
template_blas_symv
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
template_blas_lsame
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
template_blas_axpy
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
integer
int integer
Definition: template_blas_common.h:40
template_blas_syr2
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
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
a_ref
#define a_ref(a_1, a_2)