ergo
template_lapack_sterf.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_STERF_HEADER
36 #define TEMPLATE_LAPACK_STERF_HEADER
37 
38 #include "template_lapack_common.h"
39 
40 template<class Treal>
41 int template_lapack_sterf(const integer *n, Treal *d__, Treal *e,
42  integer *info)
43 {
44 /* -- LAPACK routine (version 3.0) --
45  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
46  Courant Institute, Argonne National Lab, and Rice University
47  June 30, 1999
48 
49 
50  Purpose
51  =======
52 
53  DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
54  using the Pal-Walker-Kahan variant of the QL or QR algorithm.
55 
56  Arguments
57  =========
58 
59  N (input) INTEGER
60  The order of the matrix. N >= 0.
61 
62  D (input/output) DOUBLE PRECISION array, dimension (N)
63  On entry, the n diagonal elements of the tridiagonal matrix.
64  On exit, if INFO = 0, the eigenvalues in ascending order.
65 
66  E (input/output) DOUBLE PRECISION array, dimension (N-1)
67  On entry, the (n-1) subdiagonal elements of the tridiagonal
68  matrix.
69  On exit, E has been destroyed.
70 
71  INFO (output) INTEGER
72  = 0: successful exit
73  < 0: if INFO = -i, the i-th argument had an illegal value
74  > 0: the algorithm failed to find all of the eigenvalues in
75  a total of 30*N iterations; if INFO = i, then i
76  elements of E have not converged to zero.
77 
78  =====================================================================
79 
80 
81  Test the input parameters.
82 
83  Parameter adjustments */
84  /* Table of constant values */
85  integer c__0 = 0;
86  integer c__1 = 1;
87  Treal c_b32 = 1.;
88 
89  /* System generated locals */
90  integer i__1;
91  Treal d__1, d__2, d__3;
92  /* Local variables */
93  Treal oldc;
94  integer lend, jtot;
95  Treal c__;
96  integer i__, l, m;
97  Treal p, gamma, r__, s, alpha, sigma, anorm;
98  integer l1;
99  Treal bb;
100  integer iscale;
101  Treal oldgam, safmin;
102  Treal safmax;
103  integer lendsv;
104  Treal ssfmin;
105  integer nmaxit;
106  Treal ssfmax, rt1, rt2, eps, rte;
107  integer lsv;
108  Treal eps2;
109 
110 
111  --e;
112  --d__;
113 
114  /* Function Body */
115  *info = 0;
116 
117 /* Quick return if possible */
118 
119  if (*n < 0) {
120  *info = -1;
121  i__1 = -(*info);
122  template_blas_erbla("STERF ", &i__1);
123  return 0;
124  }
125  if (*n <= 1) {
126  return 0;
127  }
128 
129 /* Determine the unit roundoff for this environment. */
130 
131  eps = template_lapack_lamch("E", (Treal)0);
132 /* Computing 2nd power */
133  d__1 = eps;
134  eps2 = d__1 * d__1;
135  safmin = template_lapack_lamch("S", (Treal)0);
136  safmax = 1. / safmin;
137  ssfmax = template_blas_sqrt(safmax) / 3.;
138  ssfmin = template_blas_sqrt(safmin) / eps2;
139 
140 /* Compute the eigenvalues of the tridiagonal matrix. */
141 
142  nmaxit = *n * 30;
143  sigma = 0.;
144  jtot = 0;
145 
146 /* Determine where the matrix splits and choose QL or QR iteration
147  for each block, according to whether top or bottom diagonal
148  element is smaller. */
149 
150  l1 = 1;
151 
152 L10:
153  if (l1 > *n) {
154  goto L170;
155  }
156  if (l1 > 1) {
157  e[l1 - 1] = 0.;
158  }
159  i__1 = *n - 1;
160  for (m = l1; m <= i__1; ++m) {
161  if ((d__3 = e[m], absMACRO(d__3)) <= template_blas_sqrt((d__1 = d__[m], absMACRO(d__1))) *
162  template_blas_sqrt((d__2 = d__[m + 1], absMACRO(d__2))) * eps) {
163  e[m] = 0.;
164  goto L30;
165  }
166 /* L20: */
167  }
168  m = *n;
169 
170 L30:
171  l = l1;
172  lsv = l;
173  lend = m;
174  lendsv = lend;
175  l1 = m + 1;
176  if (lend == l) {
177  goto L10;
178  }
179 
180 /* Scale submatrix in rows and columns L to LEND */
181 
182  i__1 = lend - l + 1;
183  anorm = template_lapack_lanst("I", &i__1, &d__[l], &e[l]);
184  iscale = 0;
185  if (anorm > ssfmax) {
186  iscale = 1;
187  i__1 = lend - l + 1;
188  template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
189  info);
190  i__1 = lend - l;
191  template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
192  info);
193  } else if (anorm < ssfmin) {
194  iscale = 2;
195  i__1 = lend - l + 1;
196  template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
197  info);
198  i__1 = lend - l;
199  template_lapack_lascl("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
200  info);
201  }
202 
203  i__1 = lend - 1;
204  for (i__ = l; i__ <= i__1; ++i__) {
205 /* Computing 2nd power */
206  d__1 = e[i__];
207  e[i__] = d__1 * d__1;
208 /* L40: */
209  }
210 
211 /* Choose between QL and QR iteration */
212 
213  if ((d__1 = d__[lend], absMACRO(d__1)) < (d__2 = d__[l], absMACRO(d__2))) {
214  lend = lsv;
215  l = lendsv;
216  }
217 
218  if (lend >= l) {
219 
220 /* QL Iteration
221 
222  Look for small subdiagonal element. */
223 
224 L50:
225  if (l != lend) {
226  i__1 = lend - 1;
227  for (m = l; m <= i__1; ++m) {
228  if ((d__2 = e[m], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
229  + 1], absMACRO(d__1))) {
230  goto L70;
231  }
232 /* L60: */
233  }
234  }
235  m = lend;
236 
237 L70:
238  if (m < lend) {
239  e[m] = 0.;
240  }
241  p = d__[l];
242  if (m == l) {
243  goto L90;
244  }
245 
246 /* If remaining matrix is 2 by 2, use DLAE2 to compute its
247  eigenvalues. */
248 
249  if (m == l + 1) {
250  rte = template_blas_sqrt(e[l]);
251  template_lapack_lae2(&d__[l], &rte, &d__[l + 1], &rt1, &rt2);
252  d__[l] = rt1;
253  d__[l + 1] = rt2;
254  e[l] = 0.;
255  l += 2;
256  if (l <= lend) {
257  goto L50;
258  }
259  goto L150;
260  }
261 
262  if (jtot == nmaxit) {
263  goto L150;
264  }
265  ++jtot;
266 
267 /* Form shift. */
268 
269  rte = template_blas_sqrt(e[l]);
270  sigma = (d__[l + 1] - p) / (rte * 2.);
271  r__ = template_lapack_lapy2(&sigma, &c_b32);
272  sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
273 
274  c__ = 1.;
275  s = 0.;
276  gamma = d__[m] - sigma;
277  p = gamma * gamma;
278 
279 /* Inner loop */
280 
281  i__1 = l;
282  for (i__ = m - 1; i__ >= i__1; --i__) {
283  bb = e[i__];
284  r__ = p + bb;
285  if (i__ != m - 1) {
286  e[i__ + 1] = s * r__;
287  }
288  oldc = c__;
289  c__ = p / r__;
290  s = bb / r__;
291  oldgam = gamma;
292  alpha = d__[i__];
293  gamma = c__ * (alpha - sigma) - s * oldgam;
294  d__[i__ + 1] = oldgam + (alpha - gamma);
295  if (c__ != 0.) {
296  p = gamma * gamma / c__;
297  } else {
298  p = oldc * bb;
299  }
300 /* L80: */
301  }
302 
303  e[l] = s * p;
304  d__[l] = sigma + gamma;
305  goto L50;
306 
307 /* Eigenvalue found. */
308 
309 L90:
310  d__[l] = p;
311 
312  ++l;
313  if (l <= lend) {
314  goto L50;
315  }
316  goto L150;
317 
318  } else {
319 
320 /* QR Iteration
321 
322  Look for small superdiagonal element. */
323 
324 L100:
325  i__1 = lend + 1;
326  for (m = l; m >= i__1; --m) {
327  if ((d__2 = e[m - 1], absMACRO(d__2)) <= eps2 * (d__1 = d__[m] * d__[m
328  - 1], absMACRO(d__1))) {
329  goto L120;
330  }
331 /* L110: */
332  }
333  m = lend;
334 
335 L120:
336  if (m > lend) {
337  e[m - 1] = 0.;
338  }
339  p = d__[l];
340  if (m == l) {
341  goto L140;
342  }
343 
344 /* If remaining matrix is 2 by 2, use DLAE2 to compute its
345  eigenvalues. */
346 
347  if (m == l - 1) {
348  rte = template_blas_sqrt(e[l - 1]);
349  template_lapack_lae2(&d__[l], &rte, &d__[l - 1], &rt1, &rt2);
350  d__[l] = rt1;
351  d__[l - 1] = rt2;
352  e[l - 1] = 0.;
353  l += -2;
354  if (l >= lend) {
355  goto L100;
356  }
357  goto L150;
358  }
359 
360  if (jtot == nmaxit) {
361  goto L150;
362  }
363  ++jtot;
364 
365 /* Form shift. */
366 
367  rte = template_blas_sqrt(e[l - 1]);
368  sigma = (d__[l - 1] - p) / (rte * 2.);
369  r__ = template_lapack_lapy2(&sigma, &c_b32);
370  sigma = p - rte / (sigma + template_lapack_d_sign(&r__, &sigma));
371 
372  c__ = 1.;
373  s = 0.;
374  gamma = d__[m] - sigma;
375  p = gamma * gamma;
376 
377 /* Inner loop */
378 
379  i__1 = l - 1;
380  for (i__ = m; i__ <= i__1; ++i__) {
381  bb = e[i__];
382  r__ = p + bb;
383  if (i__ != m) {
384  e[i__ - 1] = s * r__;
385  }
386  oldc = c__;
387  c__ = p / r__;
388  s = bb / r__;
389  oldgam = gamma;
390  alpha = d__[i__ + 1];
391  gamma = c__ * (alpha - sigma) - s * oldgam;
392  d__[i__] = oldgam + (alpha - gamma);
393  if (c__ != 0.) {
394  p = gamma * gamma / c__;
395  } else {
396  p = oldc * bb;
397  }
398 /* L130: */
399  }
400 
401  e[l - 1] = s * p;
402  d__[l] = sigma + gamma;
403  goto L100;
404 
405 /* Eigenvalue found. */
406 
407 L140:
408  d__[l] = p;
409 
410  --l;
411  if (l >= lend) {
412  goto L100;
413  }
414  goto L150;
415 
416  }
417 
418 /* Undo scaling if necessary */
419 
420 L150:
421  if (iscale == 1) {
422  i__1 = lendsv - lsv + 1;
423  template_lapack_lascl("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
424  n, info);
425  }
426  if (iscale == 2) {
427  i__1 = lendsv - lsv + 1;
428  template_lapack_lascl("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
429  n, info);
430  }
431 
432 /* Check for no convergence to an eigenvalue after a total
433  of N*MAXIT iterations. */
434 
435  if (jtot < nmaxit) {
436  goto L10;
437  }
438  i__1 = *n - 1;
439  for (i__ = 1; i__ <= i__1; ++i__) {
440  if (e[i__] != 0.) {
441  ++(*info);
442  }
443 /* L160: */
444  }
445  goto L180;
446 
447 /* Sort eigenvalues in increasing order. */
448 
449 L170:
450  template_lapack_lasrt("I", n, &d__[1], info);
451 
452 L180:
453  return 0;
454 
455 /* End of DSTERF */
456 
457 } /* dsterf_ */
458 
459 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
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:40
int integer
Definition: template_blas_common.h:38
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:40
int template_lapack_lasrt(const char *id, const integer *n, Treal *d__, integer *info)
Definition: template_lapack_lasrt.h:40
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:46
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:144
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:199
int template_lapack_sterf(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_sterf.h:41
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition: template_lapack_lae2.h:40
Treal template_lapack_lanst(const char *norm, const integer *n, const Treal *d__, const Treal *e)
Definition: template_lapack_lanst.h:40
Treal template_blas_sqrt(Treal x)