ergo
template_lapack_sterf.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_STERF_HEADER
38#define TEMPLATE_LAPACK_STERF_HEADER
39
41
42template<class Treal>
43int template_lapack_sterf(const integer *n, Treal *d__, Treal *e,
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 June 30, 1999
50
51
52 Purpose
53 =======
54
55 DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
56 using the Pal-Walker-Kahan variant of the QL or QR algorithm.
57
58 Arguments
59 =========
60
61 N (input) INTEGER
62 The order of the matrix. N >= 0.
63
64 D (input/output) DOUBLE PRECISION array, dimension (N)
65 On entry, the n diagonal elements of the tridiagonal matrix.
66 On exit, if INFO = 0, the eigenvalues in ascending order.
67
68 E (input/output) DOUBLE PRECISION array, dimension (N-1)
69 On entry, the (n-1) subdiagonal elements of the tridiagonal
70 matrix.
71 On exit, E has been destroyed.
72
73 INFO (output) INTEGER
74 = 0: successful exit
75 < 0: if INFO = -i, the i-th argument had an illegal value
76 > 0: the algorithm failed to find all of the eigenvalues in
77 a total of 30*N iterations; if INFO = i, then i
78 elements of E have not converged to zero.
79
80 =====================================================================
81
82
83 Test the input parameters.
84
85 Parameter adjustments */
86 /* Table of constant values */
87 integer c__0 = 0;
88 integer c__1 = 1;
89 Treal c_b32 = 1.;
90
91 /* System generated locals */
92 integer i__1;
93 Treal d__1, d__2, d__3;
94 /* Local variables */
95 Treal oldc;
96 integer lend, jtot;
97 Treal c__;
98 integer i__, l, m;
99 Treal p, gamma, r__, s, alpha, sigma, anorm;
100 integer l1;
101 Treal bb;
102 integer iscale;
103 Treal oldgam, safmin;
104 Treal safmax;
105 integer lendsv;
106 Treal ssfmin;
107 integer nmaxit;
108 Treal ssfmax, rt1, rt2, eps, rte;
109 integer lsv;
110 Treal eps2;
111
112
113 --e;
114 --d__;
115
116 /* Function Body */
117 *info = 0;
118
119/* Quick return if possible */
120
121 if (*n < 0) {
122 *info = -1;
123 i__1 = -(*info);
124 template_blas_erbla("STERF ", &i__1);
125 return 0;
126 }
127 if (*n <= 1) {
128 return 0;
129 }
130
131/* Determine the unit roundoff for this environment. */
132
133 eps = template_lapack_lamch("E", (Treal)0);
134/* Computing 2nd power */
135 d__1 = eps;
136 eps2 = d__1 * d__1;
137 safmin = template_lapack_lamch("S", (Treal)0);
138 safmax = 1. / safmin;
139 ssfmax = template_blas_sqrt(safmax) / 3.;
140 ssfmin = template_blas_sqrt(safmin) / eps2;
141
142/* Compute the eigenvalues of the tridiagonal matrix. */
143
144 nmaxit = *n * 30;
145 sigma = 0.;
146 jtot = 0;
147
148/* Determine where the matrix splits and choose QL or QR iteration
149 for each block, according to whether top or bottom diagonal
150 element is smaller. */
151
152 l1 = 1;
153
154L10:
155 if (l1 > *n) {
156 goto L170;
157 }
158 if (l1 > 1) {
159 e[l1 - 1] = 0.;
160 }
161 i__1 = *n - 1;
162 for (m = l1; m <= i__1; ++m) {
163 if ((d__3 = e[m], absMACRO(d__3)) <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) *
164 template_blas_sqrt((d__2 = d__[m + 1], absMACRO(d__2))) * eps) {
165 e[m] = 0.;
166 goto L30;
167 }
168/* L20: */
169 }
170 m = *n;
171
172L30:
173 l = l1;
174 lsv = l;
175 lend = m;
176 lendsv = lend;
177 l1 = m + 1;
178 if (lend == l) {
179 goto L10;
180 }
181
182/* Scale submatrix in rows and columns L to LEND */
183
184 i__1 = lend - l + 1;
185 anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
186 iscale = 0;
187 if (anorm > ssfmax) {
188 iscale = 1;
189 i__1 = lend - l + 1;
190 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
191 info);
192 i__1 = lend - l;
193 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
194 info);
195 } else if (anorm < ssfmin) {
196 iscale = 2;
197 i__1 = lend - l + 1;
198 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
199 info);
200 i__1 = lend - l;
201 template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
202 info);
203 }
204
205 i__1 = lend - 1;
206 for (i__ = l; i__ <= i__1; ++i__) {
207/* Computing 2nd power */
208 d__1 = e[i__];
209 e[i__] = d__1 * d__1;
210/* L40: */
211 }
212
213/* Choose between QL and QR iteration */
214
215 if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
216 lend = lsv;
217 l = lendsv;
218 }
219
220 if (lend >= l) {
221
222/* QL Iteration
223
224 Look for small subdiagonal element. */
225
226L50:
227 if (l != lend) {
228 i__1 = lend - 1;
229 for (m = l; m <= i__1; ++m) {
230 if ((d__2 = e[m], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
231 + 1], absMACRO(d__1))) {
232 goto L70;
233 }
234/* L60: */
235 }
236 }
237 m = lend;
238
239L70:
240 if (m < lend) {
241 e[m] = 0.;
242 }
243 p = d__[l];
244 if (m == l) {
245 goto L90;
246 }
247
248/* If remaining matrix is 2 by 2, use DLAE2 to compute its
249 eigenvalues. */
250
251 if (m == l + 1) {
252 rte = template_blas_sqrt(e[l]);
253 template_lapack_lae2(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
254 d__[l] = rt1;
255 d__[l + 1] = rt2;
256 e[l] = 0.;
257 l += 2;
258 if (l <= lend) {
259 goto L50;
260 }
261 goto L150;
262 }
263
264 if (jtot == nmaxit) {
265 goto L150;
266 }
267 ++jtot;
268
269/* Form shift. */
270
271 rte = template_blas_sqrt(e[l]);
272 sigma = (d__[l + 1] - p) / (rte * 2.);
273 r__ = template_lapack_lapy2(&sigma, &c_b32);
274 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
275
276 c__ = 1.;
277 s = 0.;
278 gamma = d__[m] - sigma;
279 p = gamma * gamma;
280
281/* Inner loop */
282
283 i__1 = l;
284 for (i__ = m - 1; i__ >= i__1; --i__) {
285 bb = e[i__];
286 r__ = p + bb;
287 if (i__ != m - 1) {
288 e[i__ + 1] = s * r__;
289 }
290 oldc = c__;
291 c__ = p / r__;
292 s = bb / r__;
293 oldgam = gamma;
294 alpha = d__[i__];
295 gamma = c__ * (alpha - sigma) - s * oldgam;
296 d__[i__ + 1] = oldgam + (alpha - gamma);
297 if (c__ != 0.) {
298 p = gamma * gamma / c__;
299 } else {
300 p = oldc * bb;
301 }
302/* L80: */
303 }
304
305 e[l] = s * p;
306 d__[l] = sigma + gamma;
307 goto L50;
308
309/* Eigenvalue found. */
310
311L90:
312 d__[l] = p;
313
314 ++l;
315 if (l <= lend) {
316 goto L50;
317 }
318 goto L150;
319
320 } else {
321
322/* QR Iteration
323
324 Look for small superdiagonal element. */
325
326L100:
327 i__1 = lend + 1;
328 for (m = l; m >= i__1; --m) {
329 if ((d__2 = e[m - 1], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
330 - 1], absMACRO(d__1))) {
331 goto L120;
332 }
333/* L110: */
334 }
335 m = lend;
336
337L120:
338 if (m > lend) {
339 e[m - 1] = 0.;
340 }
341 p = d__[l];
342 if (m == l) {
343 goto L140;
344 }
345
346/* If remaining matrix is 2 by 2, use DLAE2 to compute its
347 eigenvalues. */
348
349 if (m == l - 1) {
350 rte = template_blas_sqrt(e[l - 1]);
351 template_lapack_lae2(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
352 d__[l] = rt1;
353 d__[l - 1] = rt2;
354 e[l - 1] = 0.;
355 l += -2;
356 if (l >= lend) {
357 goto L100;
358 }
359 goto L150;
360 }
361
362 if (jtot == nmaxit) {
363 goto L150;
364 }
365 ++jtot;
366
367/* Form shift. */
368
369 rte = template_blas_sqrt(e[l - 1]);
370 sigma = (d__[l - 1] - p) / (rte * 2.);
371 r__ = template_lapack_lapy2(&sigma, &c_b32);
372 sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
373
374 c__ = 1.;
375 s = 0.;
376 gamma = d__[m] - sigma;
377 p = gamma * gamma;
378
379/* Inner loop */
380
381 i__1 = l - 1;
382 for (i__ = m; i__ <= i__1; ++i__) {
383 bb = e[i__];
384 r__ = p + bb;
385 if (i__ != m) {
386 e[i__ - 1] = s * r__;
387 }
388 oldc = c__;
389 c__ = p / r__;
390 s = bb / r__;
391 oldgam = gamma;
392 alpha = d__[i__ + 1];
393 gamma = c__ * (alpha - sigma) - s * oldgam;
394 d__[i__] = oldgam + (alpha - gamma);
395 if (c__ != 0.) {
396 p = gamma * gamma / c__;
397 } else {
398 p = oldc * bb;
399 }
400/* L130: */
401 }
402
403 e[l - 1] = s * p;
404 d__[l] = sigma + gamma;
405 goto L100;
406
407/* Eigenvalue found. */
408
409L140:
410 d__[l] = p;
411
412 --l;
413 if (l >= lend) {
414 goto L100;
415 }
416 goto L150;
417
418 }
419
420/* Undo scaling if necessary */
421
422L150:
423 if (iscale == 1) {
424 i__1 = lendsv - lsv + 1;
425 template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
426 n, info);
427 }
428 if (iscale == 2) {
429 i__1 = lendsv - lsv + 1;
430 template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
431 n, info);
432 }
433
434/* Check for no convergence to an eigenvalue after a total
435 of N*MAXIT iterations. */
436
437 if (jtot < nmaxit) {
438 goto L10;
439 }
440 i__1 = *n - 1;
441 for (i__ = 1; i__ <= i__1; ++i__) {
442 if (e[i__] != 0.) {
443 ++(*info);
444 }
445/* L160: */
446 }
447 goto L180;
448
449/* Sort eigenvalues in increasing order. */
450
451L170:
452 template_lapack_lasrt("I", n, &d__[1], info);
453
454L180:
455 return 0;
456
457/* End of DSTERF */
458
459} /* dsterf_ */
460
461#endif
Treal template_blas_sqrt(Treal x)
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition template_lapack_lae2.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_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_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition template_lapack_lasrt.h:42
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition template_lapack_sterf.h:43