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