ergo
template_lapack_stevx.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_STEVX_HEADER
38#define TEMPLATE_LAPACK_STEVX_HEADER
39
40
41template<class Treal>
42int template_lapack_stevx(const char *jobz, const char *range, const integer *n, Treal *
43 d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il,
44 const integer *iu, const Treal *abstol, integer *m, Treal *w,
45 Treal *z__, const integer *ldz, Treal *work, integer *iwork,
46 integer *ifail, integer *info)
47{
48/* -- LAPACK driver routine (version 3.0) --
49 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
50 Courant Institute, Argonne National Lab, and Rice University
51 June 30, 1999
52
53
54 Purpose
55 =======
56
57 DSTEVX computes selected eigenvalues and, optionally, eigenvectors
58 of a real symmetric tridiagonal matrix A. Eigenvalues and
59 eigenvectors can be selected by specifying either a range of values
60 or a range of indices for the desired eigenvalues.
61
62 Arguments
63 =========
64
65 JOBZ (input) CHARACTER*1
66 = 'N': Compute eigenvalues only;
67 = 'V': Compute eigenvalues and eigenvectors.
68
69 RANGE (input) CHARACTER*1
70 = 'A': all eigenvalues will be found.
71 = 'V': all eigenvalues in the half-open interval (VL,VU]
72 will be found.
73 = 'I': the IL-th through IU-th eigenvalues will be found.
74
75 N (input) INTEGER
76 The order of the matrix. N >= 0.
77
78 D (input/output) DOUBLE PRECISION array, dimension (N)
79 On entry, the n diagonal elements of the tridiagonal matrix
80 A.
81 On exit, D may be multiplied by a constant factor chosen
82 to avoid over/underflow in computing the eigenvalues.
83
84 E (input/output) DOUBLE PRECISION array, dimension (N)
85 On entry, the (n-1) subdiagonal elements of the tridiagonal
86 matrix A in elements 1 to N-1 of E; E(N) need not be set.
87 On exit, E may be multiplied by a constant factor chosen
88 to avoid over/underflow in computing the eigenvalues.
89
90 VL (input) DOUBLE PRECISION
91 VU (input) DOUBLE PRECISION
92 If RANGE='V', the lower and upper bounds of the interval to
93 be searched for eigenvalues. VL < VU.
94 Not referenced if RANGE = 'A' or 'I'.
95
96 IL (input) INTEGER
97 IU (input) INTEGER
98 If RANGE='I', the indices (in ascending order) of the
99 smallest and largest eigenvalues to be returned.
100 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
101 Not referenced if RANGE = 'A' or 'V'.
102
103 ABSTOL (input) DOUBLE PRECISION
104 The absolute error tolerance for the eigenvalues.
105 An approximate eigenvalue is accepted as converged
106 when it is determined to lie in an interval [a,b]
107 of width less than or equal to
108
109 ABSTOL + EPS * max( |a|,|b| ) ,
110
111 where EPS is the machine precision. If ABSTOL is less
112 than or equal to zero, then EPS*|T| will be used in
113 its place, where |T| is the 1-norm of the tridiagonal
114 matrix.
115
116 Eigenvalues will be computed most accurately when ABSTOL is
117 set to twice the underflow threshold 2*DLAMCH('S'), not zero.
118 If this routine returns with INFO>0, indicating that some
119 eigenvectors did not converge, try setting ABSTOL to
120 2*DLAMCH('S').
121
122 See "Computing Small Singular Values of Bidiagonal Matrices
123 with Guaranteed High Relative Accuracy," by Demmel and
124 Kahan, LAPACK Working Note #3.
125
126 M (output) INTEGER
127 The total number of eigenvalues found. 0 <= M <= N.
128 If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
129
130 W (output) DOUBLE PRECISION array, dimension (N)
131 The first M elements contain the selected eigenvalues in
132 ascending order.
133
134 Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
135 If JOBZ = 'V', then if INFO = 0, the first M columns of Z
136 contain the orthonormal eigenvectors of the matrix A
137 corresponding to the selected eigenvalues, with the i-th
138 column of Z holding the eigenvector associated with W(i).
139 If an eigenvector fails to converge (INFO > 0), then that
140 column of Z contains the latest approximation to the
141 eigenvector, and the index of the eigenvector is returned
142 in IFAIL. If JOBZ = 'N', then Z is not referenced.
143 Note: the user must ensure that at least max(1,M) columns are
144 supplied in the array Z; if RANGE = 'V', the exact value of M
145 is not known in advance and an upper bound must be used.
146
147 LDZ (input) INTEGER
148 The leading dimension of the array Z. LDZ >= 1, and if
149 JOBZ = 'V', LDZ >= max(1,N).
150
151 WORK (workspace) DOUBLE PRECISION array, dimension (5*N)
152
153 IWORK (workspace) INTEGER array, dimension (5*N)
154
155 IFAIL (output) INTEGER array, dimension (N)
156 If JOBZ = 'V', then if INFO = 0, the first M elements of
157 IFAIL are zero. If INFO > 0, then IFAIL contains the
158 indices of the eigenvectors that failed to converge.
159 If JOBZ = 'N', then IFAIL is not referenced.
160
161 INFO (output) INTEGER
162 = 0: successful exit
163 < 0: if INFO = -i, the i-th argument had an illegal value
164 > 0: if INFO = i, then i eigenvectors failed to converge.
165 Their indices are stored in array IFAIL.
166
167 =====================================================================
168
169
170 Test the input parameters.
171
172 Parameter adjustments */
173 /* Table of constant values */
174 integer c__1 = 1;
175
176 /* System generated locals */
177 integer z_dim1, z_offset, i__1, i__2;
178 Treal d__1, d__2;
179 /* Local variables */
180 integer imax;
181 Treal rmin, rmax, tnrm;
182 integer itmp1, i__, j;
183 Treal sigma;
184 char order[1];
185 logical wantz;
186 integer jj;
187 logical alleig, indeig;
188 integer iscale, indibl;
189 logical valeig;
190 Treal safmin;
191 Treal bignum;
192 integer indisp;
193 integer indiwo;
194 integer indwrk;
195 integer nsplit;
196 Treal smlnum, eps, vll, vuu, tmp1;
197#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
198
199
200 --d__;
201 --e;
202 --w;
203 z_dim1 = *ldz;
204 z_offset = 1 + z_dim1 * 1;
205 z__ -= z_offset;
206 --work;
207 --iwork;
208 --ifail;
209
210 /* Function Body */
211 wantz = template_blas_lsame(jobz, "V");
212 alleig = template_blas_lsame(range, "A");
213 valeig = template_blas_lsame(range, "V");
214 indeig = template_blas_lsame(range, "I");
215
216 *info = 0;
217 if (! (wantz || template_blas_lsame(jobz, "N"))) {
218 *info = -1;
219 } else if (! (alleig || valeig || indeig)) {
220 *info = -2;
221 } else if (*n < 0) {
222 *info = -3;
223 } else {
224 if (valeig) {
225 if (*n > 0 && *vu <= *vl) {
226 *info = -7;
227 }
228 } else if (indeig) {
229 if (*il < 1 || *il > maxMACRO(1,*n)) {
230 *info = -8;
231 } else if (*iu < minMACRO(*n,*il) || *iu > *n) {
232 *info = -9;
233 }
234 }
235 }
236 if (*info == 0) {
237 if (*ldz < 1 || (wantz && *ldz < *n) ) {
238 *info = -14;
239 }
240 }
241
242 if (*info != 0) {
243 i__1 = -(*info);
244 template_blas_erbla("STEVX ", &i__1);
245 return 0;
246 }
247
248/* Quick return if possible */
249
250 *m = 0;
251 if (*n == 0) {
252 return 0;
253 }
254
255 if (*n == 1) {
256 if (alleig || indeig) {
257 *m = 1;
258 w[1] = d__[1];
259 } else {
260 if (*vl < d__[1] && *vu >= d__[1]) {
261 *m = 1;
262 w[1] = d__[1];
263 }
264 }
265 if (wantz) {
266 z___ref(1, 1) = 1.;
267 }
268 return 0;
269 }
270
271/* Get machine constants. */
272
273 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
274 eps = template_lapack_lamch("Precision", (Treal)0);
275 smlnum = safmin / eps;
276 bignum = 1. / smlnum;
277 rmin = template_blas_sqrt(smlnum);
278/* Computing MIN */
279 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
280 rmax = minMACRO(d__1,d__2);
281
282/* Scale matrix to allowable range, if necessary. */
283
284 iscale = 0;
285 if (valeig) {
286 vll = *vl;
287 vuu = *vu;
288 } else {
289 vll = 0.;
290 vuu = 0.;
291 }
292 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
293 if (tnrm > 0. && tnrm < rmin) {
294 iscale = 1;
295 sigma = rmin / tnrm;
296 } else if (tnrm > rmax) {
297 iscale = 1;
298 sigma = rmax / tnrm;
299 }
300 if (iscale == 1) {
301 template_blas_scal(n, &sigma, &d__[1], &c__1);
302 i__1 = *n - 1;
303 template_blas_scal(&i__1, &sigma, &e[1], &c__1);
304 if (valeig) {
305 vll = *vl * sigma;
306 vuu = *vu * sigma;
307 }
308 }
309
310/* If all eigenvalues are desired and ABSTOL is less than zero, then
311 call DSTERF or SSTEQR. If this fails for some eigenvalue, then
312 try DSTEBZ. */
313
314 if ((alleig || (indeig && *il == 1 && *iu == *n) ) && *abstol <= 0.) {
315 template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1);
316 i__1 = *n - 1;
317 template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1);
318 indwrk = *n + 1;
319 if (! wantz) {
320 template_lapack_sterf(n, &w[1], &work[1], info);
321 } else {
322 template_lapack_steqr("I", n, &w[1], &work[1], &z__[z_offset], ldz, &work[
323 indwrk], info);
324 if (*info == 0) {
325 i__1 = *n;
326 for (i__ = 1; i__ <= i__1; ++i__) {
327 ifail[i__] = 0;
328/* L10: */
329 }
330 }
331 }
332 if (*info == 0) {
333 *m = *n;
334 goto L20;
335 }
336 *info = 0;
337 }
338
339/* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. */
340
341 if (wantz) {
342 *(unsigned char *)order = 'B';
343 } else {
344 *(unsigned char *)order = 'E';
345 }
346 indwrk = 1;
347 indibl = 1;
348 indisp = indibl + *n;
349 indiwo = indisp + *n;
350 template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, &
351 nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[indwrk], &
352 iwork[indiwo], info);
353
354 if (wantz) {
355 template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], &
356 z__[z_offset], ldz, &work[indwrk], &iwork[indiwo], &ifail[1],
357 info);
358 }
359
360/* If matrix was scaled, then rescale eigenvalues appropriately. */
361
362L20:
363 if (iscale == 1) {
364 if (*info == 0) {
365 imax = *m;
366 } else {
367 imax = *info - 1;
368 }
369 d__1 = 1. / sigma;
370 template_blas_scal(&imax, &d__1, &w[1], &c__1);
371 }
372
373/* If eigenvalues are not in order, then sort them, along with
374 eigenvectors. */
375
376 if (wantz) {
377 i__1 = *m - 1;
378 for (j = 1; j <= i__1; ++j) {
379 i__ = 0;
380 tmp1 = w[j];
381 i__2 = *m;
382 for (jj = j + 1; jj <= i__2; ++jj) {
383 if (w[jj] < tmp1) {
384 i__ = jj;
385 tmp1 = w[jj];
386 }
387/* L30: */
388 }
389
390 if (i__ != 0) {
391 itmp1 = iwork[indibl + i__ - 1];
392 w[i__] = w[j];
393 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
394 w[j] = tmp1;
395 iwork[indibl + j - 1] = itmp1;
396 template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, j), &c__1);
397 if (*info != 0) {
398 itmp1 = ifail[i__];
399 ifail[i__] = ifail[j];
400 ifail[j] = itmp1;
401 }
402 }
403/* L40: */
404 }
405 }
406
407 return 0;
408
409/* End of DSTEVX */
410
411} /* dstevx_ */
412
413#undef z___ref
414
415
416#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
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
int template_blas_copy(const integer *n, const Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_copy.h:42
int template_blas_scal(const integer *n, const Treal *da, Treal *dx, const integer *incx)
Definition: template_blas_scal.h:43
int template_blas_swap(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_swap.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:42
int template_lapack_stebz(const char *range, const char *order, const integer *n, const Treal *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, const Treal *d__, const Treal *e, integer *m, integer *nsplit, Treal *w, integer *iblock, integer *isplit, Treal *work, integer *iwork, integer *info)
Definition: template_lapack_stebz.h:42
int template_lapack_stein(const integer *n, const Treal *d__, const Treal *e, const integer *m, const Treal *w, const integer *iblock, const integer *isplit, Treal *z__, const integer *ldz, Treal *work, integer *iwork, integer *ifail, integer *info)
Definition: template_lapack_stein.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 z___ref(a_1, a_2)
int template_lapack_stevx(const char *jobz, const char *range, const integer *n, Treal *d__, Treal *e, const Treal *vl, const Treal *vu, const integer *il, const integer *iu, const Treal *abstol, integer *m, Treal *w, Treal *z__, const integer *ldz, Treal *work, integer *iwork, integer *ifail, integer *info)
Definition: template_lapack_stevx.h:42