ergo
template_lapack_syev.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_SYEV_HEADER
36 #define TEMPLATE_LAPACK_SYEV_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_syev(const char *jobz, const char *uplo, const integer *n, Treal *a,
41  const integer *lda, Treal *w, Treal *work, const integer *lwork,
42  integer *info)
43 {
44 
45  //printf("entering template_lapack_syev\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  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
57  real symmetric matrix A.
58 
59  Arguments
60  =========
61 
62  JOBZ (input) CHARACTER*1
63  = 'N': Compute eigenvalues only;
64  = 'V': Compute eigenvalues and eigenvectors.
65 
66  UPLO (input) CHARACTER*1
67  = 'U': Upper triangle of A is stored;
68  = 'L': Lower triangle of A is stored.
69 
70  N (input) INTEGER
71  The order of the matrix A. N >= 0.
72 
73  A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
74  On entry, the symmetric matrix A. If UPLO = 'U', the
75  leading N-by-N upper triangular part of A contains the
76  upper triangular part of the matrix A. If UPLO = 'L',
77  the leading N-by-N lower triangular part of A contains
78  the lower triangular part of the matrix A.
79  On exit, if JOBZ = 'V', then if INFO = 0, A contains the
80  orthonormal eigenvectors of the matrix A.
81  If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
82  or the upper triangle (if UPLO='U') of A, including the
83  diagonal, is destroyed.
84 
85  LDA (input) INTEGER
86  The leading dimension of the array A. LDA >= max(1,N).
87 
88  W (output) DOUBLE PRECISION array, dimension (N)
89  If INFO = 0, the eigenvalues in ascending order.
90 
91  WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
92  On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
93 
94  LWORK (input) INTEGER
95  The length of the array WORK. LWORK >= max(1,3*N-1).
96  For optimal efficiency, LWORK >= (NB+2)*N,
97  where NB is the blocksize for DSYTRD returned by ILAENV.
98 
99  If LWORK = -1, then a workspace query is assumed; the routine
100  only calculates the optimal size of the WORK array, returns
101  this value as the first entry of the WORK array, and no error
102  message related to LWORK is issued by XERBLA.
103 
104  INFO (output) INTEGER
105  = 0: successful exit
106  < 0: if INFO = -i, the i-th argument had an illegal value
107  > 0: if INFO = i, the algorithm failed to converge; i
108  off-diagonal elements of an intermediate tridiagonal
109  form did not converge to zero.
110 
111  =====================================================================
112 
113 
114  Test the input parameters.
115 
116  Parameter adjustments */
117  /* Table of constant values */
118  integer c__1 = 1;
119  integer c_n1 = -1;
120  integer c__0 = 0;
121  Treal c_b17 = 1.;
122 
123  /* System generated locals */
124  integer a_dim1, a_offset, i__1, i__2;
125  Treal d__1;
126  /* Local variables */
127  integer inde;
128  Treal anrm;
129  integer imax;
130  Treal rmin, rmax;
131  Treal sigma;
132  integer iinfo;
133  logical lower, wantz;
134  integer nb;
135  integer iscale;
136  Treal safmin;
137  Treal bignum;
138  integer indtau;
139  integer indwrk;
140  integer llwork;
141  Treal smlnum;
142  integer lwkopt;
143  logical lquery;
144  Treal eps;
145 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
146 
147 
148  a_dim1 = *lda;
149  a_offset = 1 + a_dim1 * 1;
150  a -= a_offset;
151  --w;
152  --work;
153 
154  /* Initialization added by Elias to get rid of compiler warnings. */
155  lwkopt = 0;
156  /* Function Body */
157  wantz = template_blas_lsame(jobz, "V");
158  lower = template_blas_lsame(uplo, "L");
159  lquery = *lwork == -1;
160 
161  *info = 0;
162  if (! (wantz || template_blas_lsame(jobz, "N"))) {
163  *info = -1;
164  } else if (! (lower || template_blas_lsame(uplo, "U"))) {
165  *info = -2;
166  } else if (*n < 0) {
167  *info = -3;
168  } else if (*lda < maxMACRO(1,*n)) {
169  *info = -5;
170  } else /* if(complicated condition) */ {
171 /* Computing MAX */
172  i__1 = 1, i__2 = *n * 3 - 1;
173  if (*lwork < maxMACRO(i__1,i__2) && ! lquery) {
174  *info = -8;
175  }
176  }
177 
178  if (*info == 0) {
179  nb = template_lapack_ilaenv(&c__1, "DSYTRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
180  (ftnlen)1);
181 /* Computing MAX */
182  i__1 = 1, i__2 = (nb + 2) * *n;
183  lwkopt = maxMACRO(i__1,i__2);
184  work[1] = (Treal) lwkopt;
185  }
186 
187  if (*info != 0) {
188  i__1 = -(*info);
189  template_blas_erbla("SYEV ", &i__1);
190  return 0;
191  } else if (lquery) {
192  return 0;
193  }
194 
195 /* Quick return if possible */
196 
197  if (*n == 0) {
198  work[1] = 1.;
199  return 0;
200  }
201 
202  if (*n == 1) {
203  w[1] = a_ref(1, 1);
204  work[1] = 3.;
205  if (wantz) {
206  a_ref(1, 1) = 1.;
207  }
208  return 0;
209  }
210 
211 /* Get machine constants. */
212 
213  //printf("before getting machine constants.\n");
214 
215  //printf("calling template_lapack_lamch for Safe minimum\n");
216  safmin = template_lapack_lamch("Safe minimum", (Treal)0);
217  //printf("template_lapack_lamch for Safe minimum returned\n");
218 
219  eps = template_lapack_lamch("Precision", (Treal)0);
220 
221  //printf("after getting machine constants.\n");
222 
223  smlnum = safmin / eps;
224  bignum = 1. / smlnum;
225  rmin = template_blas_sqrt(smlnum);
226  rmax = template_blas_sqrt(bignum);
227 
228 /* Scale matrix to allowable range, if necessary. */
229 
230  anrm = template_lapack_lansy("M", uplo, n, &a[a_offset], lda, &work[1]);
231  iscale = 0;
232  if (anrm > 0. && anrm < rmin) {
233  iscale = 1;
234  sigma = rmin / anrm;
235  } else if (anrm > rmax) {
236  iscale = 1;
237  sigma = rmax / anrm;
238  }
239  if (iscale == 1) {
240  template_lapack_lascl(uplo, &c__0, &c__0, &c_b17, &sigma, n, n, &a[a_offset], lda,
241  info);
242  }
243 
244 /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */
245 
246  inde = 1;
247  indtau = inde + *n;
248  indwrk = indtau + *n;
249  llwork = *lwork - indwrk + 1;
250  template_lapack_sytrd(uplo, n, &a[a_offset], lda, &w[1], &work[inde], &work[indtau], &
251  work[indwrk], &llwork, &iinfo);
252 
253 /* For eigenvalues only, call DSTERF. For eigenvectors, first call
254  DORGTR to generate the orthogonal matrix, then call DSTEQR. */
255 
256  if (! wantz) {
257  template_lapack_sterf(n, &w[1], &work[inde], info);
258  } else {
259  template_lapack_orgtr(uplo, n, &a[a_offset], lda, &work[indtau], &work[indwrk], &
260  llwork, &iinfo);
261  template_lapack_steqr(jobz, n, &w[1], &work[inde], &a[a_offset], lda, &work[indtau],
262  info);
263  }
264 
265 /* If matrix was scaled, then rescale eigenvalues appropriately. */
266 
267  if (iscale == 1) {
268  if (*info == 0) {
269  imax = *n;
270  } else {
271  imax = *info - 1;
272  }
273  d__1 = 1. / sigma;
274  template_blas_scal(&imax, &d__1, &w[1], &c__1);
275  }
276 
277 /* Set WORK(1) to optimal workspace size. */
278 
279  work[1] = (Treal) lwkopt;
280 
281  return 0;
282 
283 /* End of DSYEV */
284 
285 } /* dsyev_ */
286 
287 #undef a_ref
288 
289 
290 #endif