ergo
template_lapack_sygs2.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_SYGS2_HEADER
38#define TEMPLATE_LAPACK_SYGS2_HEADER
39
41
42template<class Treal>
43int 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
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 maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition template_blas_scal.h:43
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_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
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
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)
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