ergo
template_lapack_sygst.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_SYGST_HEADER
38#define TEMPLATE_LAPACK_SYGST_HEADER
39
41
42template<class Treal>
43int template_lapack_sygst(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 September 30, 1994
51
52
53 Purpose
54 =======
55
56 DSYGST 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**T)*A*inv(U) or inv(L)*A*inv(L**T)
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**T or L**T*A*L.
64
65 B must have been previously factorized as U**T*U or L*L**T by DPOTRF.
66
67 Arguments
68 =========
69
70 ITYPE (input) INTEGER
71 = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
72 = 2 or 3: compute U*A*U**T or L**T*A*L.
73
74 UPLO (input) CHARACTER
75 = 'U': Upper triangle of A is stored and B is factored as
76 U**T*U;
77 = 'L': Lower triangle of A is stored and B is factored as
78 L*L**T.
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 integer c__1 = 1;
117 integer c_n1 = -1;
118 Treal c_b14 = 1.;
119 Treal c_b16 = -.5;
120 Treal c_b19 = -1.;
121 Treal c_b52 = .5;
122
123 /* System generated locals */
124 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
125 /* Local variables */
126 integer k;
127 logical upper;
128 integer kb;
129 integer nb;
130#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
131#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
132
133
134 a_dim1 = *lda;
135 a_offset = 1 + a_dim1 * 1;
136 a -= a_offset;
137 b_dim1 = *ldb;
138 b_offset = 1 + b_dim1 * 1;
139 b -= b_offset;
140
141 /* Function Body */
142 *info = 0;
143 upper = template_blas_lsame(uplo, "U");
144 if (*itype < 1 || *itype > 3) {
145 *info = -1;
146 } else if (! upper && ! template_blas_lsame(uplo, "L")) {
147 *info = -2;
148 } else if (*n < 0) {
149 *info = -3;
150 } else if (*lda < maxMACRO(1,*n)) {
151 *info = -5;
152 } else if (*ldb < maxMACRO(1,*n)) {
153 *info = -7;
154 }
155 if (*info != 0) {
156 i__1 = -(*info);
157 template_blas_erbla("SYGST ", &i__1);
158 return 0;
159 }
160
161/* Quick return if possible */
162
163 if (*n == 0) {
164 return 0;
165 }
166
167/* Determine the block size for this environment. */
168
169 nb = template_lapack_ilaenv(&c__1, "DSYGST", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
170 ftnlen)1);
171
172 if (nb <= 1 || nb >= *n) {
173
174/* Use unblocked code */
175
176 template_lapack_sygs2(itype, uplo, n, &a[a_offset], lda, &b[b_offset], ldb, info);
177 } else {
178
179/* Use blocked code */
180
181 if (*itype == 1) {
182 if (upper) {
183
184/* Compute inv(U')*A*inv(U) */
185
186 i__1 = *n;
187 i__2 = nb;
188 for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
189/* Computing MIN */
190 i__3 = *n - k + 1;
191 kb = minMACRO(i__3,nb);
192
193/* Update the upper triangle of A(k:n,k:n) */
194
195 template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
196 ldb, info);
197 if (k + kb <= *n) {
198 i__3 = *n - k - kb + 1;
199 template_blas_trsm("Left", uplo, "Transpose", "Non-unit", &kb, &
200 i__3, &c_b14, &b_ref(k, k), ldb, &a_ref(k, k
201 + kb), lda);
202 i__3 = *n - k - kb + 1;
203 template_blas_symm("Left", uplo, &kb, &i__3, &c_b16, &a_ref(k, k),
204 lda, &b_ref(k, k + kb), ldb, &c_b14, &a_ref(
205 k, k + kb), lda);
206 i__3 = *n - k - kb + 1;
207 template_blas_syr2k(uplo, "Transpose", &i__3, &kb, &c_b19, &a_ref(
208 k, k + kb), lda, &b_ref(k, k + kb), ldb, &
209 c_b14, &a_ref(k + kb, k + kb), lda);
210 i__3 = *n - k - kb + 1;
211 template_blas_symm("Left", uplo, &kb, &i__3, &c_b16, &a_ref(k, k),
212 lda, &b_ref(k, k + kb), ldb, &c_b14, &a_ref(
213 k, k + kb), lda);
214 i__3 = *n - k - kb + 1;
215 template_blas_trsm("Right", uplo, "No transpose", "Non-unit", &kb,
216 &i__3, &c_b14, &b_ref(k + kb, k + kb), ldb, &
217 a_ref(k, k + kb), lda);
218 }
219/* L10: */
220 }
221 } else {
222
223/* Compute inv(L)*A*inv(L') */
224
225 i__2 = *n;
226 i__1 = nb;
227 for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
228/* Computing MIN */
229 i__3 = *n - k + 1;
230 kb = minMACRO(i__3,nb);
231
232/* Update the lower triangle of A(k:n,k:n) */
233
234 template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
235 ldb, info);
236 if (k + kb <= *n) {
237 i__3 = *n - k - kb + 1;
238 template_blas_trsm("Right", uplo, "Transpose", "Non-unit", &i__3,
239 &kb, &c_b14, &b_ref(k, k), ldb, &a_ref(k + kb,
240 k), lda);
241 i__3 = *n - k - kb + 1;
242 template_blas_symm("Right", uplo, &i__3, &kb, &c_b16, &a_ref(k, k)
243 , lda, &b_ref(k + kb, k), ldb, &c_b14, &a_ref(
244 k + kb, k), lda);
245 i__3 = *n - k - kb + 1;
246 template_blas_syr2k(uplo, "No transpose", &i__3, &kb, &c_b19, &
247 a_ref(k + kb, k), lda, &b_ref(k + kb, k), ldb,
248 &c_b14, &a_ref(k + kb, k + kb), lda);
249 i__3 = *n - k - kb + 1;
250 template_blas_symm("Right", uplo, &i__3, &kb, &c_b16, &a_ref(k, k)
251 , lda, &b_ref(k + kb, k), ldb, &c_b14, &a_ref(
252 k + kb, k), lda);
253 i__3 = *n - k - kb + 1;
254 template_blas_trsm("Left", uplo, "No transpose", "Non-unit", &
255 i__3, &kb, &c_b14, &b_ref(k + kb, k + kb),
256 ldb, &a_ref(k + kb, k), lda);
257 }
258/* L20: */
259 }
260 }
261 } else {
262 if (upper) {
263
264/* Compute U*A*U' */
265
266 i__1 = *n;
267 i__2 = nb;
268 for (k = 1; i__2 < 0 ? k >= i__1 : k <= i__1; k += i__2) {
269/* Computing MIN */
270 i__3 = *n - k + 1;
271 kb = minMACRO(i__3,nb);
272
273/* Update the upper triangle of A(1:k+kb-1,1:k+kb-1) */
274
275 i__3 = k - 1;
276 template_blas_trmm("Left", uplo, "No transpose", "Non-unit", &i__3, &
277 kb, &c_b14, &b[b_offset], ldb, &a_ref(1, k), lda);
278 i__3 = k - 1;
279 template_blas_symm("Right", uplo, &i__3, &kb, &c_b52, &a_ref(k, k),
280 lda, &b_ref(1, k), ldb, &c_b14, &a_ref(1, k), lda);
281 i__3 = k - 1;
282 template_blas_syr2k(uplo, "No transpose", &i__3, &kb, &c_b14, &a_ref(
283 1, k), lda, &b_ref(1, k), ldb, &c_b14, &a[
284 a_offset], lda);
285 i__3 = k - 1;
286 template_blas_symm("Right", uplo, &i__3, &kb, &c_b52, &a_ref(k, k),
287 lda, &b_ref(1, k), ldb, &c_b14, &a_ref(1, k), lda);
288 i__3 = k - 1;
289 template_blas_trmm("Right", uplo, "Transpose", "Non-unit", &i__3, &kb,
290 &c_b14, &b_ref(k, k), ldb, &a_ref(1, k), lda);
291 template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
292 ldb, info);
293/* L30: */
294 }
295 } else {
296
297/* Compute L'*A*L */
298
299 i__2 = *n;
300 i__1 = nb;
301 for (k = 1; i__1 < 0 ? k >= i__2 : k <= i__2; k += i__1) {
302/* Computing MIN */
303 i__3 = *n - k + 1;
304 kb = minMACRO(i__3,nb);
305
306/* Update the lower triangle of A(1:k+kb-1,1:k+kb-1) */
307
308 i__3 = k - 1;
309 template_blas_trmm("Right", uplo, "No transpose", "Non-unit", &kb, &
310 i__3, &c_b14, &b[b_offset], ldb, &a_ref(k, 1),
311 lda);
312 i__3 = k - 1;
313 template_blas_symm("Left", uplo, &kb, &i__3, &c_b52, &a_ref(k, k),
314 lda, &b_ref(k, 1), ldb, &c_b14, &a_ref(k, 1), lda);
315 i__3 = k - 1;
316 template_blas_syr2k(uplo, "Transpose", &i__3, &kb, &c_b14, &a_ref(k,
317 1), lda, &b_ref(k, 1), ldb, &c_b14, &a[a_offset],
318 lda);
319 i__3 = k - 1;
320 template_blas_symm("Left", uplo, &kb, &i__3, &c_b52, &a_ref(k, k),
321 lda, &b_ref(k, 1), ldb, &c_b14, &a_ref(k, 1), lda);
322 i__3 = k - 1;
323 template_blas_trmm("Left", uplo, "Transpose", "Non-unit", &kb, &i__3,
324 &c_b14, &b_ref(k, k), ldb, &a_ref(k, 1), lda);
325 template_lapack_sygs2(itype, uplo, &kb, &a_ref(k, k), lda, &b_ref(k, k),
326 ldb, info);
327/* L40: */
328 }
329 }
330 }
331 }
332 return 0;
333
334/* End of DSYGST */
335
336} /* dsygst_ */
337
338#undef b_ref
339#undef a_ref
340
341
342#endif
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
int ftnlen
Definition template_blas_common.h:42
#define minMACRO(a, b)
Definition template_blas_common.h:46
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
int template_blas_symm(const char *side, const char *uplo, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition template_blas_symm.h:42
int template_blas_syr2k(const char *uplo, const char *trans, const integer *n, const integer *k, const Treal *alpha, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *beta, Treal *c__, const integer *ldc)
Definition template_blas_syr2k.h:42
int template_blas_trmm(const char *side, const char *uplo, const char *transa, const char *diag, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, Treal *b, const integer *ldb)
Definition template_blas_trmm.h:43
int template_blas_trsm(const char *side, const char *uplo, const char *transa, const char *diag, const integer *m, const integer *n, const Treal *alpha, const Treal *a, const integer *lda, Treal *b, const integer *ldb)
Definition template_blas_trsm.h:43
integer template_lapack_ilaenv(const integer *ispec, const char *name__, const char *opts, const integer *n1, const integer *n2, const integer *n3, const integer *n4, ftnlen name_len, ftnlen opts_len)
Definition template_lapack_common.cc:281
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
int template_lapack_sygst(const integer *itype, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *b, const integer *ldb, integer *info)
Definition template_lapack_sygst.h:43
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)