ergo
template_lapack_larrr.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_LARRR_HEADER
38#define TEMPLATE_LAPACK_LARRR_HEADER
39
40template<class Treal>
41int 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 }
172L11:
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
Treal template_blas_sqrt(Treal x)
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
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_lapack_lamch(const char *cmach, Treal dummyReal)
Definition template_lapack_lamch.h:202
int template_lapack_larrr(const integer *n, Treal *d__, Treal *e, integer *info)
Definition template_lapack_larrr.h:41