ergo
template_lapack_sygv.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_SYGV_HEADER
38 #define TEMPLATE_LAPACK_SYGV_HEADER
39 
40 
41 template<class Treal>
42 int 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
ftnlen
int ftnlen
Definition: template_blas_common.h:42
template_lapack_sygst
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
template_lapack_ilaenv
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
template_blas_trsm
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
template_lapack_potrf
int template_lapack_potrf(const char *uplo, const integer *n, Treal *a, const integer *lda, integer *info)
Definition: template_lapack_potrf.h:43
logical
bool logical
Definition: template_blas_common.h:41
template_lapack_syev
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
template_blas_trmm
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
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_lapack_sygv
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
integer
int integer
Definition: template_blas_common.h:40
maxMACRO
#define maxMACRO(a, b)
Definition: template_blas_common.h:45