ergo
template_lapack_steqr.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_STEQR_HEADER
38#define TEMPLATE_LAPACK_STEQR_HEADER
39
40
41template<class Treal>
42int template_lapack_steqr(const char *compz, const integer *n, Treal *d__,
43 Treal *e, Treal *z__, const integer *ldz, Treal *work,
44 integer *info)
45{
46/* -- LAPACK routine (version 3.0) --
47 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
48 Courant Institute, Argonne National Lab, and Rice University
49 September 30, 1994
50
51
52 Purpose
53 =======
54
55 DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
56 symmetric tridiagonal matrix using the implicit QL or QR method.
57 The eigenvectors of a full or band symmetric matrix can also be found
58 if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
59 tridiagonal form.
60
61 Arguments
62 =========
63
64 COMPZ (input) CHARACTER*1
65 = 'N': Compute eigenvalues only.
66 = 'V': Compute eigenvalues and eigenvectors of the original
67 symmetric matrix. On entry, Z must contain the
68 orthogonal matrix used to reduce the original matrix
69 to tridiagonal form.
70 = 'I': Compute eigenvalues and eigenvectors of the
71 tridiagonal matrix. Z is initialized to the identity
72 matrix.
73
74 N (input) INTEGER
75 The order of the matrix. N >= 0.
76
77 D (input/output) DOUBLE PRECISION array, dimension (N)
78 On entry, the diagonal elements of the tridiagonal matrix.
79 On exit, if INFO = 0, the eigenvalues in ascending order.
80
81 E (input/output) DOUBLE PRECISION array, dimension (N-1)
82 On entry, the (n-1) subdiagonal elements of the tridiagonal
83 matrix.
84 On exit, E has been destroyed.
85
86 Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
87 On entry, if COMPZ = 'V', then Z contains the orthogonal
88 matrix used in the reduction to tridiagonal form.
89 On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
90 orthonormal eigenvectors of the original symmetric matrix,
91 and if COMPZ = 'I', Z contains the orthonormal eigenvectors
92 of the symmetric tridiagonal matrix.
93 If COMPZ = 'N', then Z is not referenced.
94
95 LDZ (input) INTEGER
96 The leading dimension of the array Z. LDZ >= 1, and if
97 eigenvectors are desired, then LDZ >= max(1,N).
98
99 WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
100 If COMPZ = 'N', then WORK is not referenced.
101
102 INFO (output) INTEGER
103 = 0: successful exit
104 < 0: if INFO = -i, the i-th argument had an illegal value
105 > 0: the algorithm has failed to find all the eigenvalues in
106 a total of 30*N iterations; if INFO = i, then i
107 elements of E have not converged to zero; on exit, D
108 and E contain the elements of a symmetric tridiagonal
109 matrix which is orthogonally similar to the original
110 matrix.
111
112 =====================================================================
113
114
115 Test the input parameters.
116
117 Parameter adjustments */
118 /* Table of constant values */
119 Treal c_b9 = 0.;
120 Treal c_b10 = 1.;
121 integer c__0 = 0;
122 integer c__1 = 1;
123 integer c__2 = 2;
124
125 /* System generated locals */
126 integer z_dim1, z_offset, i__1, i__2;
127 Treal d__1, d__2;
128 /* Local variables */
129 integer lend, jtot;
130 Treal b, c__, f, g;
131 integer i__, j, k, l, m;
132 Treal p, r__, s;
133 Treal anorm;
134 integer l1;
135 integer lendm1, lendp1;
136 integer ii;
137 integer mm, iscale;
138 Treal safmin;
139 Treal safmax;
140 integer lendsv;
141 Treal ssfmin;
142 integer nmaxit, icompz;
143 Treal ssfmax;
144 integer lm1, mm1, nm1;
145 Treal rt1, rt2, eps;
146 integer lsv;
147 Treal tst, eps2;
148#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
149
150
151 --d__;
152 --e;
153 z_dim1 = *ldz;
154 z_offset = 1 + z_dim1 * 1;
155 z__ -= z_offset;
156 --work;
157
158 /* Function Body */
159 *info = 0;
160
161 if (template_blas_lsame(compz, "N")) {
162 icompz = 0;
163 } else if (template_blas_lsame(compz, "V")) {
164 icompz = 1;
165 } else if (template_blas_lsame(compz, "I")) {
166 icompz = 2;
167 } else {
168 icompz = -1;
169 }
170 if (icompz < 0) {
171 *info = -1;
172 } else if (*n < 0) {
173 *info = -2;
174 } else if (*ldz < 1 || (icompz > 0 && *ldz < maxMACRO(1,*n) ) ) {
175 *info = -6;
176 }
177 if (*info != 0) {
178 i__1 = -(*info);
179 template_blas_erbla("STEQR ", &i__1);
180 return 0;
181 }
182
183/* Quick return if possible */
184
185 if (*n == 0) {
186 return 0;
187 }
188
189 if (*n == 1) {
190 if (icompz == 2) {
191 z___ref(1, 1) = 1.;
192 }
193 return 0;
194 }
195
196/* Determine the unit roundoff and over/underflow thresholds. */
197
198 eps = template_lapack_lamch("E", (Treal)0);
199/* Computing 2nd power */
200 d__1 = eps;
201 eps2 = d__1 * d__1;
202 safmin = template_lapack_lamch("S", (Treal)0);
203 safmax = 1. / safmin;
204 ssfmax = template_blas_sqrt(safmax) / 3.;
205 ssfmin = template_blas_sqrt(safmin) / eps2;
206
207/* Compute the eigenvalues and eigenvectors of the tridiagonal
208 matrix. */
209
210 if (icompz == 2) {
211 template_lapack_laset("Full", n, n, &c_b9, &c_b10, &z__[z_offset], ldz);
212 }
213
214 nmaxit = *n * 30;
215 jtot = 0;
216
217/* Determine where the matrix splits and choose QL or QR iteration
218 for each block, according to whether top or bottom diagonal
219 element is smaller. */
220
221 l1 = 1;
222 nm1 = *n - 1;
223
224L10:
225 if (l1 > *n) {
226 goto L160;
227 }
228 if (l1 > 1) {
229 e[l1 - 1] = 0.;
230 }
231 if (l1 <= nm1) {
232 i__1 = nm1;
233 for (m = l1; m <= i__1; ++m) {
234 tst = (d__1 = e[m], absMACRO(d__1));
235 if (tst == 0.) {
236 goto L30;
237 }
238 if (tst <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) * template_blas_sqrt((d__2 = d__[m
239 + 1], absMACRO(d__2))) * eps) {
240 e[m] = 0.;
241 goto L30;
242 }
243/* L20: */
244 }
245 }
246 m = *n;
247
248L30:
249 l = l1;
250 lsv = l;
251 lend = m;
252 lendsv = lend;
253 l1 = m + 1;
254 if (lend == l) {
255 goto L10;
256 }
257
258/* Scale submatrix in rows and columns L to LEND */
259
260 i__1 = lend - l + 1;
261 anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
262 iscale = 0;
263 if (anorm == 0.) {
264 goto L10;
265 }
266 if (anorm > ssfmax) {
267 iscale = 1;
268 i__1 = lend - l + 1;
269 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
270 info);
271 i__1 = lend - l;
272 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
273 info);
274 } else if (anorm < ssfmin) {
275 iscale = 2;
276 i__1 = lend - l + 1;
277 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
278 info);
279 i__1 = lend - l;
280 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
281 info);
282 }
283
284/* Choose between QL and QR iteration */
285
286 if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
287 lend = lsv;
288 l = lendsv;
289 }
290
291 if (lend > l) {
292
293/* QL Iteration
294
295 Look for small subdiagonal element. */
296
297L40:
298 if (l != lend) {
299 lendm1 = lend - 1;
300 i__1 = lendm1;
301 for (m = l; m <= i__1; ++m) {
302/* Computing 2nd power */
303 d__2 = (d__1 = e[m], absMACRO(d__1));
304 tst = d__2 * d__2;
305 if (tst <= eps2 * (d__1 = d__[m], absMACRO(d__1)) * (d__2 = d__[m
306 + 1], absMACRO(d__2)) + safmin) {
307 goto L60;
308 }
309/* L50: */
310 }
311 }
312
313 m = lend;
314
315L60:
316 if (m < lend) {
317 e[m] = 0.;
318 }
319 p = d__[l];
320 if (m == l) {
321 goto L80;
322 }
323
324/* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
325 to compute its eigensystem. */
326
327 if (m == l + 1) {
328 if (icompz > 0) {
329 template_lapack_laev2(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
330 work[l] = c__;
331 work[*n - 1 + l] = s;
332 template_lapack_lasr("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
333 z___ref(1, l), ldz);
334 } else {
335 template_lapack_lae2(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
336 }
337 d__[l] = rt1;
338 d__[l + 1] = rt2;
339 e[l] = 0.;
340 l += 2;
341 if (l <= lend) {
342 goto L40;
343 }
344 goto L140;
345 }
346
347 if (jtot == nmaxit) {
348 goto L140;
349 }
350 ++jtot;
351
352/* Form shift. */
353
354 g = (d__[l + 1] - p) / (e[l] * 2.);
355 r__ = template_lapack_lapy2(&g, &c_b10);
356 g = d__[m] - p + e[l] / (g + template_lapack_d_sign(&r__, &g));
357
358 s = 1.;
359 c__ = 1.;
360 p = 0.;
361
362/* Inner loop */
363
364 mm1 = m - 1;
365 i__1 = l;
366 for (i__ = mm1; i__ >= i__1; --i__) {
367 f = s * e[i__];
368 b = c__ * e[i__];
369 template_lapack_lartg(&g, &f, &c__, &s, &r__);
370 if (i__ != m - 1) {
371 e[i__ + 1] = r__;
372 }
373 g = d__[i__ + 1] - p;
374 r__ = (d__[i__] - g) * s + c__ * 2. * b;
375 p = s * r__;
376 d__[i__ + 1] = g + p;
377 g = c__ * r__ - b;
378
379/* If eigenvectors are desired, then save rotations. */
380
381 if (icompz > 0) {
382 work[i__] = c__;
383 work[*n - 1 + i__] = -s;
384 }
385
386/* L70: */
387 }
388
389/* If eigenvectors are desired, then apply saved rotations. */
390
391 if (icompz > 0) {
392 mm = m - l + 1;
393 template_lapack_lasr("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &
394 z___ref(1, l), ldz);
395 }
396
397 d__[l] -= p;
398 e[l] = g;
399 goto L40;
400
401/* Eigenvalue found. */
402
403L80:
404 d__[l] = p;
405
406 ++l;
407 if (l <= lend) {
408 goto L40;
409 }
410 goto L140;
411
412 } else {
413
414/* QR Iteration
415
416 Look for small superdiagonal element. */
417
418L90:
419 if (l != lend) {
420 lendp1 = lend + 1;
421 i__1 = lendp1;
422 for (m = l; m >= i__1; --m) {
423/* Computing 2nd power */
424 d__2 = (d__1 = e[m - 1], absMACRO(d__1));
425 tst = d__2 * d__2;
426 if (tst <= eps2 * (d__1 = d__[m], absMACRO(d__1)) * (d__2 = d__[m
427 - 1], absMACRO(d__2)) + safmin) {
428 goto L110;
429 }
430/* L100: */
431 }
432 }
433
434 m = lend;
435
436L110:
437 if (m > lend) {
438 e[m - 1] = 0.;
439 }
440 p = d__[l];
441 if (m == l) {
442 goto L130;
443 }
444
445/* If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
446 to compute its eigensystem. */
447
448 if (m == l - 1) {
449 if (icompz > 0) {
450 template_lapack_laev2(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
451 ;
452 work[m] = c__;
453 work[*n - 1 + m] = s;
454 template_lapack_lasr("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], &
455 z___ref(1, l - 1), ldz);
456 } else {
457 template_lapack_lae2(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
458 }
459 d__[l - 1] = rt1;
460 d__[l] = rt2;
461 e[l - 1] = 0.;
462 l += -2;
463 if (l >= lend) {
464 goto L90;
465 }
466 goto L140;
467 }
468
469 if (jtot == nmaxit) {
470 goto L140;
471 }
472 ++jtot;
473
474/* Form shift. */
475
476 g = (d__[l - 1] - p) / (e[l - 1] * 2.);
477 r__ = template_lapack_lapy2(&g, &c_b10);
478 g = d__[m] - p + e[l - 1] / (g + template_lapack_d_sign(&r__, &g));
479
480 s = 1.;
481 c__ = 1.;
482 p = 0.;
483
484/* Inner loop */
485
486 lm1 = l - 1;
487 i__1 = lm1;
488 for (i__ = m; i__ <= i__1; ++i__) {
489 f = s * e[i__];
490 b = c__ * e[i__];
491 template_lapack_lartg(&g, &f, &c__, &s, &r__);
492 if (i__ != m) {
493 e[i__ - 1] = r__;
494 }
495 g = d__[i__] - p;
496 r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b;
497 p = s * r__;
498 d__[i__] = g + p;
499 g = c__ * r__ - b;
500
501/* If eigenvectors are desired, then save rotations. */
502
503 if (icompz > 0) {
504 work[i__] = c__;
505 work[*n - 1 + i__] = s;
506 }
507
508/* L120: */
509 }
510
511/* If eigenvectors are desired, then apply saved rotations. */
512
513 if (icompz > 0) {
514 mm = l - m + 1;
515 template_lapack_lasr("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &
516 z___ref(1, m), ldz);
517 }
518
519 d__[l] -= p;
520 e[lm1] = g;
521 goto L90;
522
523/* Eigenvalue found. */
524
525L130:
526 d__[l] = p;
527
528 --l;
529 if (l >= lend) {
530 goto L90;
531 }
532 goto L140;
533
534 }
535
536/* Undo scaling if necessary */
537
538L140:
539 if (iscale == 1) {
540 i__1 = lendsv - lsv + 1;
541 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
542 n, info);
543 i__1 = lendsv - lsv;
544 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
545 info);
546 } else if (iscale == 2) {
547 i__1 = lendsv - lsv + 1;
548 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
549 n, info);
550 i__1 = lendsv - lsv;
551 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
552 info);
553 }
554
555/* Check for no convergence to an eigenvalue after a total
556 of N*MAXIT iterations. */
557
558 if (jtot < nmaxit) {
559 goto L10;
560 }
561 i__1 = *n - 1;
562 for (i__ = 1; i__ <= i__1; ++i__) {
563 if (e[i__] != 0.) {
564 ++(*info);
565 }
566/* L150: */
567 }
568 goto L190;
569
570/* Order eigenvalues and eigenvectors. */
571
572L160:
573 if (icompz == 0) {
574
575/* Use Quick Sort */
576
577 template_lapack_lasrt("I", n, &d__[1], info);
578
579 } else {
580
581/* Use Selection Sort to minimize swaps of eigenvectors */
582
583 i__1 = *n;
584 for (ii = 2; ii <= i__1; ++ii) {
585 i__ = ii - 1;
586 k = i__;
587 p = d__[i__];
588 i__2 = *n;
589 for (j = ii; j <= i__2; ++j) {
590 if (d__[j] < p) {
591 k = j;
592 p = d__[j];
593 }
594/* L170: */
595 }
596 if (k != i__) {
597 d__[k] = d__[i__];
598 d__[i__] = p;
599 template_blas_swap(n, &z___ref(1, i__), &c__1, &z___ref(1, k), &c__1);
600 }
601/* L180: */
602 }
603 }
604
605L190:
606 return 0;
607
608/* End of DSTEQR */
609
610} /* dsteqr_ */
611
612#undef z___ref
613
614
615#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 absMACRO(x)
Definition: template_blas_common.h:47
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
int template_blas_swap(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy)
Definition: template_blas_swap.h:42
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition: template_lapack_lae2.h:42
int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition: template_lapack_laev2.h:42
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:42
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.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_laset(const char *uplo, const integer *m, const integer *n, const Treal *alpha, const Treal *beta, Treal *a, const integer *lda)
Definition: template_lapack_laset.h:42
int template_lapack_lasr(const char *side, const char *pivot, const char *direct, const integer *m, const integer *n, const Treal *c__, const Treal *s, Treal *a, const integer *lda)
Definition: template_lapack_lasr.h:42
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:42
#define z___ref(a_1, a_2)
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