ergo
template_lapack_sygs2.h
Go to the documentation of this file.
1 /* Ergo, version 3.4, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2014 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_SYGS2_HEADER
36 #define TEMPLATE_LAPACK_SYGS2_HEADER
37 
38 #include "template_lapack_common.h"
39 
40 template<class Treal>
41 int template_lapack_sygs2(const integer *itype, const char *uplo, const integer *n,
42  Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *
43  info)
44 {
45 /* -- LAPACK routine (version 3.0) --
46  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47  Courant Institute, Argonne National Lab, and Rice University
48  February 29, 1992
49 
50 
51  Purpose
52  =======
53 
54  DSYGS2 reduces a real symmetric-definite generalized eigenproblem
55  to standard form.
56 
57  If ITYPE = 1, the problem is A*x = lambda*B*x,
58  and A is overwritten by inv(U')*A*inv(U) or inv(L)*A*inv(L')
59 
60  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
61  B*A*x = lambda*x, and A is overwritten by U*A*U` or L'*A*L.
62 
63  B must have been previously factorized as U'*U or L*L' by DPOTRF.
64 
65  Arguments
66  =========
67 
68  ITYPE (input) INTEGER
69  = 1: compute inv(U')*A*inv(U) or inv(L)*A*inv(L');
70  = 2 or 3: compute U*A*U' or L'*A*L.
71 
72  UPLO (input) CHARACTER
73  Specifies whether the upper or lower triangular part of the
74  symmetric matrix A is stored, and how B has been factorized.
75  = 'U': Upper triangular
76  = 'L': Lower triangular
77 
78  N (input) INTEGER
79  The order of the matrices A and B. N >= 0.
80 
81  A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
82  On entry, the symmetric matrix A. If UPLO = 'U', the leading
83  n by n upper triangular part of A contains the upper
84  triangular part of the matrix A, and the strictly lower
85  triangular part of A is not referenced. If UPLO = 'L', the
86  leading n by n lower triangular part of A contains the lower
87  triangular part of the matrix A, and the strictly upper
88  triangular part of A is not referenced.
89 
90  On exit, if INFO = 0, the transformed matrix, stored in the
91  same format as A.
92 
93  LDA (input) INTEGER
94  The leading dimension of the array A. LDA >= max(1,N).
95 
96  B (input) DOUBLE PRECISION array, dimension (LDB,N)
97  The triangular factor from the Cholesky factorization of B,
98  as returned by DPOTRF.
99 
100  LDB (input) INTEGER
101  The leading dimension of the array B. LDB >= max(1,N).
102 
103  INFO (output) INTEGER
104  = 0: successful exit.
105  < 0: if INFO = -i, the i-th argument had an illegal value.
106 
107  =====================================================================
108 
109 
110  Test the input parameters.
111 
112  Parameter adjustments */
113  /* Table of constant values */
114  Treal c_b6 = -1.;
115  integer c__1 = 1;
116  Treal c_b27 = 1.;
117 
118  /* System generated locals */
119  integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
120  Treal d__1;
121  /* Local variables */
122  integer k;
123  logical upper;
124  Treal ct;
125  Treal akk, bkk;
126 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
127 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
128 
129 
130  a_dim1 = *lda;
131  a_offset = 1 + a_dim1 * 1;
132  a -= a_offset;
133  b_dim1 = *ldb;
134  b_offset = 1 + b_dim1 * 1;
135  b -= b_offset;
136 
137  /* Function Body */
138  *info = 0;
139  upper = template_blas_lsame(uplo, "U");
140  if (*itype < 1 || *itype > 3) {
141  *info = -1;
142  } else if (! upper && ! template_blas_lsame(uplo, "L")) {
143  *info = -2;
144  } else if (*n < 0) {
145  *info = -3;
146  } else if (*lda < maxMACRO(1,*n)) {
147  *info = -5;
148  } else if (*ldb < maxMACRO(1,*n)) {
149  *info = -7;
150  }
151  if (*info != 0) {
152  i__1 = -(*info);
153  template_blas_erbla("SYGS2 ", &i__1);
154  return 0;
155  }
156 
157  if (*itype == 1) {
158  if (upper) {
159 
160 /* Compute inv(U')*A*inv(U) */
161 
162  i__1 = *n;
163  for (k = 1; k <= i__1; ++k) {
164 
165 /* Update the upper triangle of A(k:n,k:n) */
166 
167  akk = a_ref(k, k);
168  bkk = b_ref(k, k);
169 /* Computing 2nd power */
170  d__1 = bkk;
171  akk /= d__1 * d__1;
172  a_ref(k, k) = akk;
173  if (k < *n) {
174  i__2 = *n - k;
175  d__1 = 1. / bkk;
176  template_blas_scal(&i__2, &d__1, &a_ref(k, k + 1), lda);
177  ct = akk * -.5;
178  i__2 = *n - k;
179  template_blas_axpy(&i__2, &ct, &b_ref(k, k + 1), ldb, &a_ref(k, k + 1)
180  , lda);
181  i__2 = *n - k;
182  template_blas_syr2(uplo, &i__2, &c_b6, &a_ref(k, k + 1),
183  lda, &b_ref(k, k + 1), ldb,
184  &a_ref(k + 1, k + 1), lda);
185  i__2 = *n - k;
186  template_blas_axpy(&i__2, &ct, &b_ref(k, k + 1), ldb, &a_ref(k, k + 1)
187  , lda);
188  i__2 = *n - k;
189  template_blas_trsv(uplo, "Transpose", "Non-unit", &i__2, &b_ref(k + 1,
190  k + 1), ldb, &a_ref(k, k + 1), lda);
191  }
192 /* L10: */
193  }
194  } else {
195 
196 /* Compute inv(L)*A*inv(L') */
197 
198  i__1 = *n;
199  for (k = 1; k <= i__1; ++k) {
200 
201 /* Update the lower triangle of A(k:n,k:n) */
202 
203  akk = a_ref(k, k);
204  bkk = b_ref(k, k);
205 /* Computing 2nd power */
206  d__1 = bkk;
207  akk /= d__1 * d__1;
208  a_ref(k, k) = akk;
209  if (k < *n) {
210  i__2 = *n - k;
211  d__1 = 1. / bkk;
212  template_blas_scal(&i__2, &d__1, &a_ref(k + 1, k), &c__1);
213  ct = akk * -.5;
214  i__2 = *n - k;
215  template_blas_axpy(&i__2, &ct, &b_ref(k + 1, k), &c__1, &a_ref(k + 1,
216  k), &c__1);
217  i__2 = *n - k;
218  template_blas_syr2(uplo, &i__2, &c_b6, &a_ref(k + 1, k),
219  &c__1, &b_ref(k + 1, k),
220  &c__1, &a_ref(k + 1, k + 1),
221  lda);
222 
223  i__2 = *n - k;
224  template_blas_axpy(&i__2, &ct, &b_ref(k + 1, k), &c__1, &a_ref(k + 1,
225  k), &c__1);
226  i__2 = *n - k;
227  template_blas_trsv(uplo, "No transpose", "Non-unit", &i__2, &b_ref(k
228  + 1, k + 1), ldb, &a_ref(k + 1, k), &c__1);
229  }
230 /* L20: */
231  }
232  }
233  } else {
234  if (upper) {
235 
236 /* Compute U*A*U' */
237 
238  i__1 = *n;
239  for (k = 1; k <= i__1; ++k) {
240 
241 /* Update the upper triangle of A(1:k,1:k) */
242 
243  akk = a_ref(k, k);
244  bkk = b_ref(k, k);
245  i__2 = k - 1;
246  template_blas_trmv(uplo, "No transpose", "Non-unit", &i__2, &b[b_offset],
247  ldb, &a_ref(1, k), &c__1);
248  ct = akk * .5;
249  i__2 = k - 1;
250  template_blas_axpy(&i__2, &ct, &b_ref(1, k), &c__1, &a_ref(1, k), &c__1);
251  i__2 = k - 1;
252  template_blas_syr2(uplo, &i__2, &c_b27, &a_ref(1, k), &c__1, &b_ref(1, k),
253  &c__1, &a[a_offset], lda);
254  i__2 = k - 1;
255  template_blas_axpy(&i__2, &ct, &b_ref(1, k), &c__1, &a_ref(1, k), &c__1);
256  i__2 = k - 1;
257  template_blas_scal(&i__2, &bkk, &a_ref(1, k), &c__1);
258 /* Computing 2nd power */
259  d__1 = bkk;
260  a_ref(k, k) = akk * (d__1 * d__1);
261 /* L30: */
262  }
263  } else {
264 
265 /* Compute L'*A*L */
266 
267  i__1 = *n;
268  for (k = 1; k <= i__1; ++k) {
269 
270 /* Update the lower triangle of A(1:k,1:k) */
271 
272  akk = a_ref(k, k);
273  bkk = b_ref(k, k);
274  i__2 = k - 1;
275  template_blas_trmv(uplo, "Transpose", "Non-unit", &i__2, &b[b_offset],
276  ldb, &a_ref(k, 1), lda);
277  ct = akk * .5;
278  i__2 = k - 1;
279  template_blas_axpy(&i__2, &ct, &b_ref(k, 1), ldb, &a_ref(k, 1), lda);
280  i__2 = k - 1;
281  template_blas_syr2(uplo, &i__2, &c_b27, &a_ref(k, 1), lda, &b_ref(k, 1),
282  ldb, &a[a_offset], lda);
283  i__2 = k - 1;
284  template_blas_axpy(&i__2, &ct, &b_ref(k, 1), ldb, &a_ref(k, 1), lda);
285  i__2 = k - 1;
286  template_blas_scal(&i__2, &bkk, &a_ref(k, 1), lda);
287 /* Computing 2nd power */
288  d__1 = bkk;
289  a_ref(k, k) = akk * (d__1 * d__1);
290 /* L40: */
291  }
292  }
293  }
294  return 0;
295 
296 /* End of DSYGS2 */
297 
298 } /* dsygs2_ */
299 
300 #undef b_ref
301 #undef a_ref
302 
303 
304 #endif
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:41
int integer
Definition: template_blas_common.h:38
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:40
#define b_ref(a_1, a_2)
#define a_ref(a_1, a_2)
int template_blas_trmv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition: template_blas_trmv.h:40
#define maxMACRO(a, b)
Definition: template_blas_common.h:43
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
int template_blas_trsv(const char *uplo, const char *trans, const char *diag, const integer *n, const Treal *a, const integer *lda, Treal *x, const integer *incx)
Definition: template_blas_trsv.h:40
bool logical
Definition: template_blas_common.h:39
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:41
int template_lapack_sygs2(const integer *itype, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *info)
Definition: template_lapack_sygs2.h:41
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:44