ergo
template_lapack_syev.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_SYEV_HEADER
38#define TEMPLATE_LAPACK_SYEV_HEADER
39
40
41template<class Treal>
42int 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
Treal template_blas_sqrt(Treal x)
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_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
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_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_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
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
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_sterf.h:43
#define a_ref(a_1, a_2)
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_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