ergo
template_lapack_steqr.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_STEQR_HEADER
38 #define TEMPLATE_LAPACK_STEQR_HEADER
39 
40 
41 template<class Treal>
42 int 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 
224 L10:
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 
248 L30:
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 
297 L40:
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 
315 L60:
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 
403 L80:
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 
418 L90:
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 
436 L110:
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 
525 L130:
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 
538 L140:
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 
572 L160:
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 
605 L190:
606  return 0;
607 
608 /* End of DSTEQR */
609 
610 } /* dsteqr_ */
611 
612 #undef z___ref
613 
614 
615 #endif
#define absMACRO(x)
Definition: template_blas_common.h:47
#define z___ref(a_1, a_2)
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_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition: template_lapack_laev2.h:42
int integer
Definition: template_blas_common.h:40
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
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_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
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
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_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition: template_lapack_lartg.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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
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)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46