ergo
template_lapack_sygv.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_SYGV_HEADER
38#define TEMPLATE_LAPACK_SYGV_HEADER
39
40
41template<class Treal>
42int template_lapack_sygv(const integer *itype, const char *jobz, const char *uplo, const integer *
43 n, Treal *a, const integer *lda, Treal *b, const integer *ldb,
44 Treal *w, Treal *work, const integer *lwork, integer *info)
45{
46
47 //printf("entering template_lapack_sygv\n");
48
49/* -- LAPACK driver routine (version 3.0) --
50 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
51 Courant Institute, Argonne National Lab, and Rice University
52 June 30, 1999
53
54
55 Purpose
56 =======
57
58 DSYGV computes all the eigenvalues, and optionally, the eigenvectors
59 of a real generalized symmetric-definite eigenproblem, of the form
60 A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
61 Here A and B are assumed to be symmetric and B is also
62 positive definite.
63
64 Arguments
65 =========
66
67 ITYPE (input) INTEGER
68 Specifies the problem type to be solved:
69 = 1: A*x = (lambda)*B*x
70 = 2: A*B*x = (lambda)*x
71 = 3: B*A*x = (lambda)*x
72
73 JOBZ (input) CHARACTER*1
74 = 'N': Compute eigenvalues only;
75 = 'V': Compute eigenvalues and eigenvectors.
76
77 UPLO (input) CHARACTER*1
78 = 'U': Upper triangles of A and B are stored;
79 = 'L': Lower triangles of A and B are stored.
80
81 N (input) INTEGER
82 The order of the matrices A and B. N >= 0.
83
84 A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
85 On entry, the symmetric matrix A. If UPLO = 'U', the
86 leading N-by-N upper triangular part of A contains the
87 upper triangular part of the matrix A. If UPLO = 'L',
88 the leading N-by-N lower triangular part of A contains
89 the lower triangular part of the matrix A.
90
91 On exit, if JOBZ = 'V', then if INFO = 0, A contains the
92 matrix Z of eigenvectors. The eigenvectors are normalized
93 as follows:
94 if ITYPE = 1 or 2, Z**T*B*Z = I;
95 if ITYPE = 3, Z**T*inv(B)*Z = I.
96 If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
97 or the lower triangle (if UPLO='L') of A, including the
98 diagonal, is destroyed.
99
100 LDA (input) INTEGER
101 The leading dimension of the array A. LDA >= max(1,N).
102
103 B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
104 On entry, the symmetric positive definite matrix B.
105 If UPLO = 'U', the leading N-by-N upper triangular part of B
106 contains the upper triangular part of the matrix B.
107 If UPLO = 'L', the leading N-by-N lower triangular part of B
108 contains the lower triangular part of the matrix B.
109
110 On exit, if INFO <= N, the part of B containing the matrix is
111 overwritten by the triangular factor U or L from the Cholesky
112 factorization B = U**T*U or B = L*L**T.
113
114 LDB (input) INTEGER
115 The leading dimension of the array B. LDB >= max(1,N).
116
117 W (output) DOUBLE PRECISION array, dimension (N)
118 If INFO = 0, the eigenvalues in ascending order.
119
120 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
121 On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
122
123 LWORK (input) INTEGER
124 The length of the array WORK. LWORK >= max(1,3*N-1).
125 For optimal efficiency, LWORK >= (NB+2)*N,
126 where NB is the blocksize for DSYTRD returned by ILAENV.
127
128 If LWORK = -1, then a workspace query is assumed; the routine
129 only calculates the optimal size of the WORK array, returns
130 this value as the first entry of the WORK array, and no error
131 message related to LWORK is issued by XERBLA.
132
133 INFO (output) INTEGER
134 = 0: successful exit
135 < 0: if INFO = -i, the i-th argument had an illegal value
136 > 0: DPOTRF or DSYEV returned an error code:
137 <= N: if INFO = i, DSYEV failed to converge;
138 i off-diagonal elements of an intermediate
139 tridiagonal form did not converge to zero;
140 > N: if INFO = N + i, for 1 <= i <= N, then the leading
141 minor of order i of B is not positive definite.
142 The factorization of B could not be completed and
143 no eigenvalues or eigenvectors were computed.
144
145 =====================================================================
146
147
148 Test the input parameters.
149
150 Parameter adjustments */
151 /* Table of constant values */
152 integer c__1 = 1;
153 integer c_n1 = -1;
154 Treal c_b16 = 1.;
155
156 /* System generated locals */
157 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
158 /* Local variables */
159 integer neig;
160 char trans[1];
161 logical upper;
162 logical wantz;
163 integer nb;
164 integer lwkopt;
165 logical lquery;
166
167
168 a_dim1 = *lda;
169 a_offset = 1 + a_dim1 * 1;
170 a -= a_offset;
171 b_dim1 = *ldb;
172 b_offset = 1 + b_dim1 * 1;
173 b -= b_offset;
174 --w;
175 --work;
176
177 /* Initialization added by Elias to get rid of compiler warnings. */
178 lwkopt = 0;
179 /* Function Body */
180 wantz = template_blas_lsame(jobz, "V");
181 upper = template_blas_lsame(uplo, "U");
182 lquery = *lwork == -1;
183
184 *info = 0;
185 if (*itype < 1 || *itype > 3) {
186 *info = -1;
187 } else if (! (wantz || template_blas_lsame(jobz, "N"))) {
188 *info = -2;
189 } else if (! (upper || template_blas_lsame(uplo, "L"))) {
190 *info = -3;
191 } else if (*n < 0) {
192 *info = -4;
193 } else if (*lda < maxMACRO(1,*n)) {
194 *info = -6;
195 } else if (*ldb < maxMACRO(1,*n)) {
196 *info = -8;
197 } else /* if(complicated condition) */ {
198/* Computing MAX */
199 i__1 = 1, i__2 = *n * 3 - 1;
200 if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
201 *info = -11;
202 }
203 }
204
205 if (*info == 0) {
206 nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
207 (ftnlen)1);
208 lwkopt = (nb + 2) * *n;
209 work[1] = (Treal) lwkopt;
210 }
211
212 if (*info != 0) {
213 i__1 = -(*info);
214 template_blas_erbla("SYGV ", &i__1);
215 return 0;
216 } else if (lquery) {
217 return 0;
218 }
219
220/* Quick return if possible */
221
222 if (*n == 0) {
223 return 0;
224 }
225
226/* Form a Cholesky factorization of B. */
227
228 //printf("calling template_lapack_potrf\n");
229 template_lapack_potrf(uplo, n, &b[b_offset], ldb, info);
230 //printf("template_lapack_potrf returned\n");
231
232
233 if (*info != 0) {
234 *info = *n + *info;
235 return 0;
236 }
237
238/* Transform problem to standard eigenvalue problem and solve. */
239
240 //printf("calling template_lapack_sygst\n");
241 template_lapack_sygst(itype, uplo, n, &a[a_offset], lda, &b[b_offset], ldb, info);
242 //printf("template_lapack_sygst returned\n");
243
244 //printf("calling template_lapack_syev\n");
245 template_lapack_syev(jobz, uplo, n, &a[a_offset], lda, &w[1], &work[1], lwork, info);
246 //printf("template_lapack_syev returned\n");
247
248 if (wantz) {
249
250/* Backtransform eigenvectors to the original problem. */
251
252 neig = *n;
253 if (*info > 0) {
254 neig = *info - 1;
255 }
256 if (*itype == 1 || *itype == 2) {
257
258/* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
259 backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */
260
261 if (upper) {
262 *(unsigned char *)trans = 'N';
263 } else {
264 *(unsigned char *)trans = 'T';
265 }
266
267 template_blas_trsm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[
268 b_offset], ldb, &a[a_offset], lda);
269
270 } else if (*itype == 3) {
271
272/* For B*A*x=(lambda)*x;
273 backtransform eigenvectors: x = L*y or U'*y */
274
275 if (upper) {
276 *(unsigned char *)trans = 'T';
277 } else {
278 *(unsigned char *)trans = 'N';
279 }
280
281 template_blas_trmm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[
282 b_offset], ldb, &a[a_offset], lda);
283 }
284 }
285
286 work[1] = (Treal) lwkopt;
287 return 0;
288
289/* End of DSYGV */
290
291} /* dsygv_ */
292
293#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 maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
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_potrf(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_potrf.h:43
int template_lapack_syev(const char *jobz, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *w, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_syev.h:42
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
int template_lapack_sygv(const integer *itype, const char *jobz, const char *uplo, const integer *n, Treal *a, const integer *lda, Treal *b, const integer *ldb, Treal *w, Treal *work, const integer *lwork, integer *info)
Definition: template_lapack_sygv.h:42