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