ergo
template_lapack_stevr.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_STEVR_HEADER
38#define TEMPLATE_LAPACK_STEVR_HEADER
39
40template<class Treal>
41int template_lapack_stevr(const char *jobz, const char *range, const integer *n,
42 Treal * d__, Treal *e, const Treal *vl,
43 const Treal *vu, const integer *il,
44 const integer *iu, const Treal *abstol,
45 integer *m, Treal *w,
46 Treal *z__, const integer *ldz, integer *isuppz,
47 Treal *work,
48 integer *lwork, integer *iwork, integer *liwork,
49 integer *info)
50{
51 /* System generated locals */
52 integer z_dim1, z_offset, i__1, i__2;
53 Treal d__1, d__2;
54
55 /* Builtin functions */
56
57 /* Local variables */
58 integer i__, j, jj;
59 Treal eps, vll, vuu, tmp1;
60 integer imax;
61 Treal rmin, rmax;
62 logical test;
63 Treal tnrm;
64 integer itmp1;
65 Treal sigma;
66 char order[1];
67 integer lwmin;
68 logical wantz;
69 logical alleig, indeig;
70 integer iscale, ieeeok, indibl, indifl;
71 logical valeig;
72 Treal safmin;
73 Treal bignum;
74 integer indisp;
75 integer indiwo;
76 integer liwmin;
77 logical tryrac;
78 integer nsplit;
79 Treal smlnum;
80 logical lquery;
81
82
83/* -- LAPACK driver routine (version 3.2) -- */
84/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
85/* November 2006 */
86
87/* .. Scalar Arguments .. */
88/* .. */
89/* .. Array Arguments .. */
90/* .. */
91
92/* Purpose */
93/* ======= */
94
95/* DSTEVR computes selected eigenvalues and, optionally, eigenvectors */
96/* of a real symmetric tridiagonal matrix T. Eigenvalues and */
97/* eigenvectors can be selected by specifying either a range of values */
98/* or a range of indices for the desired eigenvalues. */
99
100/* Whenever possible, DSTEVR calls DSTEMR to compute the */
101/* eigenspectrum using Relatively Robust Representations. DSTEMR */
102/* computes eigenvalues by the dqds algorithm, while orthogonal */
103/* eigenvectors are computed from various "good" L D L^T representations */
104/* (also known as Relatively Robust Representations). Gram-Schmidt */
105/* orthogonalization is avoided as far as possible. More specifically, */
106/* the various steps of the algorithm are as follows. For the i-th */
107/* unreduced block of T, */
108/* (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T */
109/* is a relatively robust representation, */
110/* (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high */
111/* relative accuracy by the dqds algorithm, */
112/* (c) If there is a cluster of close eigenvalues, "choose" sigma_i */
113/* close to the cluster, and go to step (a), */
114/* (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, */
115/* compute the corresponding eigenvector by forming a */
116/* rank-revealing twisted factorization. */
117/* The desired accuracy of the output can be specified by the input */
118/* parameter ABSTOL. */
119
120/* For more details, see "A new O(n^2) algorithm for the symmetric */
121/* tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, */
122/* Computer Science Division Technical Report No. UCB//CSD-97-971, */
123/* UC Berkeley, May 1997. */
124
125
126/* Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested */
127/* on machines which conform to the ieee-754 floating point standard. */
128/* DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and */
129/* when partial spectrum requests are made. */
130
131/* Normal execution of DSTEMR may create NaNs and infinities and */
132/* hence may abort due to a floating point exception in environments */
133/* which do not handle NaNs and infinities in the ieee standard default */
134/* manner. */
135
136/* Arguments */
137/* ========= */
138
139/* JOBZ (input) CHARACTER*1 */
140/* = 'N': Compute eigenvalues only; */
141/* = 'V': Compute eigenvalues and eigenvectors. */
142
143/* RANGE (input) CHARACTER*1 */
144/* = 'A': all eigenvalues will be found. */
145/* = 'V': all eigenvalues in the half-open interval (VL,VU] */
146/* will be found. */
147/* = 'I': the IL-th through IU-th eigenvalues will be found. */
148/* ********* For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and */
149/* ********* DSTEIN are called */
150
151/* N (input) INTEGER */
152/* The order of the matrix. N >= 0. */
153
154/* D (input/output) DOUBLE PRECISION array, dimension (N) */
155/* On entry, the n diagonal elements of the tridiagonal matrix */
156/* A. */
157/* On exit, D may be multiplied by a constant factor chosen */
158/* to avoid over/underflow in computing the eigenvalues. */
159
160/* E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1)) */
161/* On entry, the (n-1) subdiagonal elements of the tridiagonal */
162/* matrix A in elements 1 to N-1 of E. */
163/* On exit, E may be multiplied by a constant factor chosen */
164/* to avoid over/underflow in computing the eigenvalues. */
165
166/* VL (input) DOUBLE PRECISION */
167/* VU (input) DOUBLE PRECISION */
168/* If RANGE='V', the lower and upper bounds of the interval to */
169/* be searched for eigenvalues. VL < VU. */
170/* Not referenced if RANGE = 'A' or 'I'. */
171
172/* IL (input) INTEGER */
173/* IU (input) INTEGER */
174/* If RANGE='I', the indices (in ascending order) of the */
175/* smallest and largest eigenvalues to be returned. */
176/* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
177/* Not referenced if RANGE = 'A' or 'V'. */
178
179/* ABSTOL (input) DOUBLE PRECISION */
180/* The absolute error tolerance for the eigenvalues. */
181/* An approximate eigenvalue is accepted as converged */
182/* when it is determined to lie in an interval [a,b] */
183/* of width less than or equal to */
184
185/* ABSTOL + EPS * max( |a|,|b| ) , */
186
187/* where EPS is the machine precision. If ABSTOL is less than */
188/* or equal to zero, then EPS*|T| will be used in its place, */
189/* where |T| is the 1-norm of the tridiagonal matrix obtained */
190/* by reducing A to tridiagonal form. */
191
192/* See "Computing Small Singular Values of Bidiagonal Matrices */
193/* with Guaranteed High Relative Accuracy," by Demmel and */
194/* Kahan, LAPACK Working Note #3. */
195
196/* If high relative accuracy is important, set ABSTOL to */
197/* DLAMCH( 'Safe minimum' ). Doing so will guarantee that */
198/* eigenvalues are computed to high relative accuracy when */
199/* possible in future releases. The current code does not */
200/* make any guarantees about high relative accuracy, but */
201/* future releases will. See J. Barlow and J. Demmel, */
202/* "Computing Accurate Eigensystems of Scaled Diagonally */
203/* Dominant Matrices", LAPACK Working Note #7, for a discussion */
204/* of which matrices define their eigenvalues to high relative */
205/* accuracy. */
206
207/* M (output) INTEGER */
208/* The total number of eigenvalues found. 0 <= M <= N. */
209/* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
210
211/* W (output) DOUBLE PRECISION array, dimension (N) */
212/* The first M elements contain the selected eigenvalues in */
213/* ascending order. */
214
215/* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) */
216/* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
217/* contain the orthonormal eigenvectors of the matrix A */
218/* corresponding to the selected eigenvalues, with the i-th */
219/* column of Z holding the eigenvector associated with W(i). */
220/* Note: the user must ensure that at least max(1,M) columns are */
221/* supplied in the array Z; if RANGE = 'V', the exact value of M */
222/* is not known in advance and an upper bound must be used. */
223
224/* LDZ (input) INTEGER */
225/* The leading dimension of the array Z. LDZ >= 1, and if */
226/* JOBZ = 'V', LDZ >= max(1,N). */
227
228/* ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) */
229/* The support of the eigenvectors in Z, i.e., the indices */
230/* indicating the nonzero elements in Z. The i-th eigenvector */
231/* is nonzero only in elements ISUPPZ( 2*i-1 ) through */
232/* ISUPPZ( 2*i ). */
233/* ********* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 */
234
235/* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
236/* On exit, if INFO = 0, WORK(1) returns the optimal (and */
237/* minimal) LWORK. */
238
239/* LWORK (input) INTEGER */
240/* The dimension of the array WORK. LWORK >= max(1,20*N). */
241
242/* If LWORK = -1, then a workspace query is assumed; the routine */
243/* only calculates the optimal sizes of the WORK and IWORK */
244/* arrays, returns these values as the first entries of the WORK */
245/* and IWORK arrays, and no error message related to LWORK or */
246/* LIWORK is issued by XERBLA. */
247
248/* IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) */
249/* On exit, if INFO = 0, IWORK(1) returns the optimal (and */
250/* minimal) LIWORK. */
251
252/* LIWORK (input) INTEGER */
253/* The dimension of the array IWORK. LIWORK >= max(1,10*N). */
254
255/* If LIWORK = -1, then a workspace query is assumed; the */
256/* routine only calculates the optimal sizes of the WORK and */
257/* IWORK arrays, returns these values as the first entries of */
258/* the WORK and IWORK arrays, and no error message related to */
259/* LWORK or LIWORK is issued by XERBLA. */
260
261/* INFO (output) INTEGER */
262/* = 0: successful exit */
263/* < 0: if INFO = -i, the i-th argument had an illegal value */
264/* > 0: Internal error */
265
266/* Further Details */
267/* =============== */
268
269/* Based on contributions by */
270/* Inderjit Dhillon, IBM Almaden, USA */
271/* Osni Marques, LBNL/NERSC, USA */
272/* Ken Stanley, Computer Science Division, University of */
273/* California at Berkeley, USA */
274
275/* ===================================================================== */
276
277/* .. Parameters .. */
278/* .. */
279/* .. Local Scalars .. */
280/* .. */
281/* .. External Functions .. */
282/* .. */
283/* .. External Subroutines .. */
284/* .. */
285/* .. Intrinsic Functions .. */
286/* .. */
287/* .. Executable Statements .. */
288
289
290/* Test the input parameters. */
291
292 /* Parameter adjustments */
293 /* Table of constant values */
294
295 integer c__10 = 10;
296 integer c__1 = 1;
297 integer c__2 = 2;
298 integer c__3 = 3;
299 integer c__4 = 4;
300
301 --d__;
302 --e;
303 --w;
304 z_dim1 = *ldz;
305 z_offset = 1 + z_dim1;
306 z__ -= z_offset;
307 --isuppz;
308 --work;
309 --iwork;
310
311 /* Function Body */
312 ieeeok = template_lapack_ilaenv(&c__10, "DSTEVR", "N", &c__1, &c__2, &c__3, &c__4, (ftnlen)6, (ftnlen)1);
313
314 wantz = template_blas_lsame(jobz, "V");
315 alleig = template_blas_lsame(range, "A");
316 valeig = template_blas_lsame(range, "V");
317 indeig = template_blas_lsame(range, "I");
318
319 lquery = *lwork == -1 || *liwork == -1;
320/* Computing MAX */
321 i__1 = 1, i__2 = *n * 20;
322 lwmin = maxMACRO(i__1,i__2);
323/* Computing MAX */
324 i__1 = 1, i__2 = *n * 10;
325 liwmin = maxMACRO(i__1,i__2);
326
327
328 *info = 0;
329 if (! (wantz || template_blas_lsame(jobz, "N"))) {
330 *info = -1;
331 } else if (! (alleig || valeig || indeig)) {
332 *info = -2;
333 } else if (*n < 0) {
334 *info = -3;
335 } else {
336 if (valeig) {
337 if (*n > 0 && *vu <= *vl) {
338 *info = -7;
339 }
340 } else if (indeig) {
341 if (*il < 1 || *il > maxMACRO(1,*n)) {
342 *info = -8;
343 } else if (*iu < minMACRO(*n,*il) || *iu > *n) {
344 *info = -9;
345 }
346 }
347 }
348 if (*info == 0) {
349 if (*ldz < 1 || ( wantz && *ldz < *n ) ) {
350 *info = -14;
351 }
352 }
353
354 if (*info == 0) {
355 work[1] = (Treal) lwmin;
356 iwork[1] = liwmin;
357
358 if (*lwork < lwmin && ! lquery) {
359 *info = -17;
360 } else if (*liwork < liwmin && ! lquery) {
361 *info = -19;
362 }
363 }
364
365 if (*info != 0) {
366 i__1 = -(*info);
367 template_blas_erbla("STEVR", &i__1);
368 return 0;
369 } else if (lquery) {
370 return 0;
371 }
372
373/* Quick return if possible */
374
375 *m = 0;
376 if (*n == 0) {
377 return 0;
378 }
379
380 if (*n == 1) {
381 if (alleig || indeig) {
382 *m = 1;
383 w[1] = d__[1];
384 } else {
385 if (*vl < d__[1] && *vu >= d__[1]) {
386 *m = 1;
387 w[1] = d__[1];
388 }
389 }
390 if (wantz) {
391 z__[z_dim1 + 1] = 1.;
392 }
393 return 0;
394 }
395
396/* Get machine constants. */
397
398 safmin = template_lapack_lamch("Safe minimum", (Treal)0);
399 eps = template_lapack_lamch("Precision", (Treal)0);
400 smlnum = safmin / eps;
401 bignum = 1. / smlnum;
402 rmin = template_blas_sqrt(smlnum);
403/* Computing MIN */
404 d__1 = template_blas_sqrt(bignum), d__2 = 1. / template_blas_sqrt(template_blas_sqrt(safmin));
405 rmax = minMACRO(d__1,d__2);
406
407
408/* Scale matrix to allowable range, if necessary. */
409
410 iscale = 0;
411 vll = *vl;
412 vuu = *vu;
413
414 tnrm = template_lapack_lanst("M", n, &d__[1], &e[1]);
415 if (tnrm > 0. && tnrm < rmin) {
416 iscale = 1;
417 sigma = rmin / tnrm;
418 } else if (tnrm > rmax) {
419 iscale = 1;
420 sigma = rmax / tnrm;
421 }
422 if (iscale == 1) {
423 template_blas_scal(n, &sigma, &d__[1], &c__1);
424 i__1 = *n - 1;
425 template_blas_scal(&i__1, &sigma, &e[1], &c__1);
426 if (valeig) {
427 vll = *vl * sigma;
428 vuu = *vu * sigma;
429 }
430 }
431/* Initialize indices into workspaces. Note: These indices are used only */
432/* if DSTERF or DSTEMR fail. */
433/* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and */
434/* stores the block indices of each of the M<=N eigenvalues. */
435 indibl = 1;
436/* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and */
437/* stores the starting and finishing indices of each block. */
438 indisp = indibl + *n;
439/* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors */
440/* that corresponding to eigenvectors that fail to converge in */
441/* DSTEIN. This information is discarded; if any fail, the driver */
442/* returns INFO > 0. */
443 indifl = indisp + *n;
444/* INDIWO is the offset of the remaining integer workspace. */
445 indiwo = indisp + *n;
446
447/* If all eigenvalues are desired, then */
448/* call DSTERF or DSTEMR. If this fails for some eigenvalue, then */
449/* try DSTEBZ. */
450
451
452 test = FALSE_;
453 if (indeig) {
454 if (*il == 1 && *iu == *n) {
455 test = TRUE_;
456 }
457 }
458 if ((alleig || test) && ieeeok == 1) {
459 i__1 = *n - 1;
460 template_blas_copy(&i__1, &e[1], &c__1, &work[1], &c__1);
461 if (! wantz) {
462 template_blas_copy(n, &d__[1], &c__1, &w[1], &c__1);
463 template_lapack_sterf(n, &w[1], &work[1], info);
464 } else {
465 template_blas_copy(n, &d__[1], &c__1, &work[*n + 1], &c__1);
466 if (*abstol <= *n * 2. * eps) {
467 tryrac = TRUE_;
468 } else {
469 tryrac = FALSE_;
470 }
471 i__1 = *lwork - (*n << 1);
472 template_lapack_stemr(jobz, "A", n, &work[*n + 1], &work[1], vl, vu, il, iu, m,
473 &w[1], &z__[z_offset], ldz, n, &isuppz[1], &tryrac, &work[
474 (*n << 1) + 1], &i__1, &iwork[1], liwork, info);
475
476 }
477 if (*info == 0) {
478 *m = *n;
479 goto L10;
480 }
481 *info = 0;
482 }
483
484/* Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN. */
485
486 if (wantz) {
487 *(unsigned char *)order = 'B';
488 } else {
489 *(unsigned char *)order = 'E';
490 }
491 template_lapack_stebz(range, order, n, &vll, &vuu, il, iu, abstol, &d__[1], &e[1], m, &
492 nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[1], &iwork[
493 indiwo], info);
494
495 if (wantz) {
496 template_lapack_stein(n, &d__[1], &e[1], m, &w[1], &iwork[indibl], &iwork[indisp], &
497 z__[z_offset], ldz, &work[1], &iwork[indiwo], &iwork[indifl],
498 info);
499 }
500
501/* If matrix was scaled, then rescale eigenvalues appropriately. */
502
503L10:
504 if (iscale == 1) {
505 if (*info == 0) {
506 imax = *m;
507 } else {
508 imax = *info - 1;
509 }
510 d__1 = 1. / sigma;
511 template_blas_scal(&imax, &d__1, &w[1], &c__1);
512 }
513
514/* If eigenvalues are not in order, then sort them, along with */
515/* eigenvectors. */
516
517 if (wantz) {
518 i__1 = *m - 1;
519 for (j = 1; j <= i__1; ++j) {
520 i__ = 0;
521 tmp1 = w[j];
522 i__2 = *m;
523 for (jj = j + 1; jj <= i__2; ++jj) {
524 if (w[jj] < tmp1) {
525 i__ = jj;
526 tmp1 = w[jj];
527 }
528/* L20: */
529 }
530
531 if (i__ != 0) {
532 itmp1 = iwork[i__];
533 w[i__] = w[j];
534 iwork[i__] = iwork[j];
535 w[j] = tmp1;
536 iwork[j] = itmp1;
537 template_blas_swap(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
538 &c__1);
539 }
540/* L30: */
541 }
542 }
543
544/* Causes problems with tests 19 & 20: */
545/* IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002 */
546
547
548 work[1] = (Treal) lwmin;
549 iwork[1] = liwmin;
550 return 0;
551
552/* End of DSTEVR */
553
554} /* dstevr_ */
555
556#endif
557
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 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
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 TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
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_stemr(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, integer *m, Treal *w, Treal *z__, const integer *ldz, const integer *nzc, integer *isuppz, logical *tryrac, Treal *work, integer *lwork, integer *iwork, integer *liwork, integer *info)
Definition: template_lapack_stemr.h:41
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_sterf.h:43
int template_lapack_stevr(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, integer *isuppz, Treal *work, integer *lwork, integer *iwork, integer *liwork, integer *info)
Definition: template_lapack_stevr.h:41