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