ergo
template_lapack_larrr.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_LARRR_HEADER
38 #define TEMPLATE_LAPACK_LARRR_HEADER
39 
40 template<class Treal>
41 int template_lapack_larrr(const integer *n, Treal *d__, Treal *e,
42  integer *info)
43 {
44  /* System generated locals */
45  integer i__1;
46  Treal d__1;
47 
48 
49  /* Local variables */
50  integer i__;
51  Treal eps, tmp, tmp2, rmin;
52  Treal offdig, safmin;
53  logical yesrel;
54  Treal smlnum, offdig2;
55 
56 
57 /* -- LAPACK auxiliary routine (version 3.2) -- */
58 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
59 /* November 2006 */
60 
61 /* .. Scalar Arguments .. */
62 /* .. */
63 /* .. Array Arguments .. */
64 /* .. */
65 
66 
67 /* Purpose */
68 /* ======= */
69 
70 /* Perform tests to decide whether the symmetric tridiagonal matrix T */
71 /* warrants expensive computations which guarantee high relative accuracy */
72 /* in the eigenvalues. */
73 
74 /* Arguments */
75 /* ========= */
76 
77 /* N (input) INTEGER */
78 /* The order of the matrix. N > 0. */
79 
80 /* D (input) DOUBLE PRECISION array, dimension (N) */
81 /* The N diagonal elements of the tridiagonal matrix T. */
82 
83 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
84 /* On entry, the first (N-1) entries contain the subdiagonal */
85 /* elements of the tridiagonal matrix T; E(N) is set to ZERO. */
86 
87 /* INFO (output) INTEGER */
88 /* INFO = 0(default) : the matrix warrants computations preserving */
89 /* relative accuracy. */
90 /* INFO = 1 : the matrix warrants computations guaranteeing */
91 /* only absolute accuracy. */
92 
93 /* Further Details */
94 /* =============== */
95 
96 /* Based on contributions by */
97 /* Beresford Parlett, University of California, Berkeley, USA */
98 /* Jim Demmel, University of California, Berkeley, USA */
99 /* Inderjit Dhillon, University of Texas, Austin, USA */
100 /* Osni Marques, LBNL/NERSC, USA */
101 /* Christof Voemel, University of California, Berkeley, USA */
102 
103 /* ===================================================================== */
104 
105 /* .. Parameters .. */
106 /* .. */
107 /* .. Local Scalars .. */
108 /* .. */
109 /* .. External Functions .. */
110 /* .. */
111 /* .. Intrinsic Functions .. */
112 /* .. */
113 /* .. Executable Statements .. */
114 
115 /* As a default, do NOT go for relative-accuracy preserving computations. */
116  /* Parameter adjustments */
117  --e;
118  --d__;
119 
120  /* Function Body */
121  *info = 1;
122  safmin = template_lapack_lamch("Safe minimum", (Treal)0);
123  eps = template_lapack_lamch("Precision", (Treal)0);
124  smlnum = safmin / eps;
125  rmin = template_blas_sqrt(smlnum);
126 /* Tests for relative accuracy */
127 
128 /* Test for scaled diagonal dominance */
129 /* Scale the diagonal entries to one and check whether the sum of the */
130 /* off-diagonals is less than one */
131 
132 /* The sdd relative error bounds have a 1/(1- 2*x) factor in them, */
133 /* x = max(OFFDIG + OFFDIG2), so when x is close to 1/2, no relative */
134 /* accuracy is promised. In the notation of the code fragment below, */
135 /* 1/(1 - (OFFDIG + OFFDIG2)) is the condition number. */
136 /* We don't think it is worth going into "sdd mode" unless the relative */
137 /* condition number is reasonable, not 1/macheps. */
138 /* The threshold should be compatible with other thresholds used in the */
139 /* code. We set OFFDIG + OFFDIG2 <= .999 =: RELCOND, it corresponds */
140 /* to losing at most 3 decimal digits: 1 / (1 - (OFFDIG + OFFDIG2)) <= 1000 */
141 /* instead of the current OFFDIG + OFFDIG2 < 1 */
142 
143  yesrel = TRUE_;
144  offdig = 0.;
145  tmp = template_blas_sqrt((absMACRO(d__[1])));
146  if (tmp < rmin) {
147  yesrel = FALSE_;
148  }
149  if (! yesrel) {
150  goto L11;
151  }
152  i__1 = *n;
153  for (i__ = 2; i__ <= i__1; ++i__) {
154  tmp2 = template_blas_sqrt((d__1 = d__[i__], absMACRO(d__1)));
155  if (tmp2 < rmin) {
156  yesrel = FALSE_;
157  }
158  if (! yesrel) {
159  goto L11;
160  }
161  offdig2 = (d__1 = e[i__ - 1], absMACRO(d__1)) / (tmp * tmp2);
162  if (offdig + offdig2 >= .999) {
163  yesrel = FALSE_;
164  }
165  if (! yesrel) {
166  goto L11;
167  }
168  tmp = tmp2;
169  offdig = offdig2;
170 /* L10: */
171  }
172 L11:
173  if (yesrel) {
174  *info = 0;
175  return 0;
176  } else {
177  }
178 
179 
180 /* *** MORE TO BE IMPLEMENTED *** */
181 
182 
183 /* Test if the lower bidiagonal matrix L from T = L D L^T */
184 /* (zero shift facto) is well conditioned */
185 
186 
187 /* Test if the upper bidiagonal matrix U from T = U D U^T */
188 /* (zero shift facto) is well conditioned. */
189 /* In this case, the matrix needs to be flipped and, at the end */
190 /* of the eigenvector computation, the flip needs to be applied */
191 /* to the computed eigenvectors (and the support) */
192 
193 
194  return 0;
195 
196 /* END OF DLARRR */
197 
198 } /* dlarrr_ */
199 
200 #endif
#define absMACRO(x)
Definition: template_blas_common.h:47
int integer
Definition: template_blas_common.h:40
int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, integer *info)
Definition: template_lapack_larrr.h:41
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
bool logical
Definition: template_blas_common.h:41
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
Treal template_blas_sqrt(Treal x)