ergo
template_lapack_sterf.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
40 #include "template_lapack_common.h"
41 
42 template<class Treal>
43 int 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 
154 L10:
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 
172 L30:
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 
226 L50:
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 
239 L70:
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 
311 L90:
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 
326 L100:
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 
337 L120:
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 
409 L140:
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 
422 L150:
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 
451 L170:
452  template_lapack_lasrt("I", n, &d__[1], info);
453 
454 L180:
455  return 0;
456 
457 /* End of DSTERF */
458 
459 } /* dsterf_ */
460 
461 #endif
#define absMACRO(x)
Definition: template_blas_common.h:47
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 integer
Definition: template_blas_common.h:40
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:42
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_sterf.h:43
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_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:42
Treal template_blas_sqrt(Treal x)