ergo
template_lapack_sygv.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_SYGV_HEADER
36 #define TEMPLATE_LAPACK_SYGV_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_sygv(const integer *itype, const char *jobz, const char *uplo, const integer *
41  n, Treal *a, const integer *lda, Treal *b, const integer *ldb,
42  Treal *w, Treal *work, const integer *lwork, integer *info)
43 {
44 
45  //printf("entering template_lapack_sygv\n");
46 
47 /* -- LAPACK driver routine (version 3.0) --
48  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
49  Courant Institute, Argonne National Lab, and Rice University
50  June 30, 1999
51 
52 
53  Purpose
54  =======
55 
56  DSYGV computes all the eigenvalues, and optionally, the eigenvectors
57  of a real generalized symmetric-definite eigenproblem, of the form
58  A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
59  Here A and B are assumed to be symmetric and B is also
60  positive definite.
61 
62  Arguments
63  =========
64 
65  ITYPE (input) INTEGER
66  Specifies the problem type to be solved:
67  = 1: A*x = (lambda)*B*x
68  = 2: A*B*x = (lambda)*x
69  = 3: B*A*x = (lambda)*x
70 
71  JOBZ (input) CHARACTER*1
72  = 'N': Compute eigenvalues only;
73  = 'V': Compute eigenvalues and eigenvectors.
74 
75  UPLO (input) CHARACTER*1
76  = 'U': Upper triangles of A and B are stored;
77  = 'L': Lower triangles of A and B are stored.
78 
79  N (input) INTEGER
80  The order of the matrices A and B. N >= 0.
81 
82  A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
83  On entry, the symmetric matrix A. If UPLO = 'U', the
84  leading N-by-N upper triangular part of A contains the
85  upper triangular part of the matrix A. If UPLO = 'L',
86  the leading N-by-N lower triangular part of A contains
87  the lower triangular part of the matrix A.
88 
89  On exit, if JOBZ = 'V', then if INFO = 0, A contains the
90  matrix Z of eigenvectors. The eigenvectors are normalized
91  as follows:
92  if ITYPE = 1 or 2, Z**T*B*Z = I;
93  if ITYPE = 3, Z**T*inv(B)*Z = I.
94  If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
95  or the lower triangle (if UPLO='L') of A, including the
96  diagonal, is destroyed.
97 
98  LDA (input) INTEGER
99  The leading dimension of the array A. LDA >= max(1,N).
100 
101  B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
102  On entry, the symmetric positive definite matrix B.
103  If UPLO = 'U', the leading N-by-N upper triangular part of B
104  contains the upper triangular part of the matrix B.
105  If UPLO = 'L', the leading N-by-N lower triangular part of B
106  contains the lower triangular part of the matrix B.
107 
108  On exit, if INFO <= N, the part of B containing the matrix is
109  overwritten by the triangular factor U or L from the Cholesky
110  factorization B = U**T*U or B = L*L**T.
111 
112  LDB (input) INTEGER
113  The leading dimension of the array B. LDB >= max(1,N).
114 
115  W (output) DOUBLE PRECISION array, dimension (N)
116  If INFO = 0, the eigenvalues in ascending order.
117 
118  WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
119  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
120 
121  LWORK (input) INTEGER
122  The length of the array WORK. LWORK >= max(1,3*N-1).
123  For optimal efficiency, LWORK >= (NB+2)*N,
124  where NB is the blocksize for DSYTRD returned by ILAENV.
125 
126  If LWORK = -1, then a workspace query is assumed; the routine
127  only calculates the optimal size of the WORK array, returns
128  this value as the first entry of the WORK array, and no error
129  message related to LWORK is issued by XERBLA.
130 
131  INFO (output) INTEGER
132  = 0: successful exit
133  < 0: if INFO = -i, the i-th argument had an illegal value
134  > 0: DPOTRF or DSYEV returned an error code:
135  <= N: if INFO = i, DSYEV failed to converge;
136  i off-diagonal elements of an intermediate
137  tridiagonal form did not converge to zero;
138  > N: if INFO = N + i, for 1 <= i <= N, then the leading
139  minor of order i of B is not positive definite.
140  The factorization of B could not be completed and
141  no eigenvalues or eigenvectors were computed.
142 
143  =====================================================================
144 
145 
146  Test the input parameters.
147 
148  Parameter adjustments */
149  /* Table of constant values */
150  integer c__1 = 1;
151  integer c_n1 = -1;
152  Treal c_b16 = 1.;
153 
154  /* System generated locals */
155  integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
156  /* Local variables */
157  integer neig;
158  char trans[1];
159  logical upper;
160  logical wantz;
161  integer nb;
162  integer lwkopt;
163  logical lquery;
164 
165 
166  a_dim1 = *lda;
167  a_offset = 1 + a_dim1 * 1;
168  a -= a_offset;
169  b_dim1 = *ldb;
170  b_offset = 1 + b_dim1 * 1;
171  b -= b_offset;
172  --w;
173  --work;
174 
175  /* Initialization added by Elias to get rid of compiler warnings. */
176  lwkopt = 0;
177  /* Function Body */
178  wantz = template_blas_lsame(jobz, "V");
179  upper = template_blas_lsame(uplo, "U");
180  lquery = *lwork == -1;
181 
182  *info = 0;
183  if (*itype < 1 || *itype > 3) {
184  *info = -1;
185  } else if (! (wantz || template_blas_lsame(jobz, "N"))) {
186  *info = -2;
187  } else if (! (upper || template_blas_lsame(uplo, "L"))) {
188  *info = -3;
189  } else if (*n < 0) {
190  *info = -4;
191  } else if (*lda < maxMACRO(1,*n)) {
192  *info = -6;
193  } else if (*ldb < maxMACRO(1,*n)) {
194  *info = -8;
195  } else /* if(complicated condition) */ {
196 /* Computing MAX */
197  i__1 = 1, i__2 = *n * 3 - 1;
198  if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
199  *info = -11;
200  }
201  }
202 
203  if (*info == 0) {
204  nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
205  (ftnlen)1);
206  lwkopt = (nb + 2) * *n;
207  work[1] = (Treal) lwkopt;
208  }
209 
210  if (*info != 0) {
211  i__1 = -(*info);
212  template_blas_erbla("SYGV ", &i__1);
213  return 0;
214  } else if (lquery) {
215  return 0;
216  }
217 
218 /* Quick return if possible */
219 
220  if (*n == 0) {
221  return 0;
222  }
223 
224 /* Form a Cholesky factorization of B. */
225 
226  //printf("calling template_lapack_potrf\n");
227  template_lapack_potrf(uplo, n, &b[b_offset], ldb, info);
228  //printf("template_lapack_potrf returned\n");
229 
230 
231  if (*info != 0) {
232  *info = *n + *info;
233  return 0;
234  }
235 
236 /* Transform problem to standard eigenvalue problem and solve. */
237 
238  //printf("calling template_lapack_sygst\n");
239  template_lapack_sygst(itype, uplo, n, &a[a_offset], lda, &b[b_offset], ldb, info);
240  //printf("template_lapack_sygst returned\n");
241 
242  //printf("calling template_lapack_syev\n");
243  template_lapack_syev(jobz, uplo, n, &a[a_offset], lda, &w[1], &work[1], lwork, info);
244  //printf("template_lapack_syev returned\n");
245 
246  if (wantz) {
247 
248 /* Backtransform eigenvectors to the original problem. */
249 
250  neig = *n;
251  if (*info > 0) {
252  neig = *info - 1;
253  }
254  if (*itype == 1 || *itype == 2) {
255 
256 /* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
257  backtransform eigenvectors: x = inv(L)'*y or inv(U)*y */
258 
259  if (upper) {
260  *(unsigned char *)trans = 'N';
261  } else {
262  *(unsigned char *)trans = 'T';
263  }
264 
265  template_blas_trsm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[
266  b_offset], ldb, &a[a_offset], lda);
267 
268  } else if (*itype == 3) {
269 
270 /* For B*A*x=(lambda)*x;
271  backtransform eigenvectors: x = L*y or U'*y */
272 
273  if (upper) {
274  *(unsigned char *)trans = 'T';
275  } else {
276  *(unsigned char *)trans = 'N';
277  }
278 
279  template_blas_trmm("Left", uplo, trans, "Non-unit", n, &neig, &c_b16, &b[
280  b_offset], ldb, &a[a_offset], lda);
281  }
282  }
283 
284  work[1] = (Treal) lwkopt;
285  return 0;
286 
287 /* End of DSYGV */
288 
289 } /* dsygv_ */
290 
291 #endif