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