ergo
Main Page
Namespaces
Classes
Files
File List
File Members
template_lapack_lae2.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_LAE2_HEADER
36
#define TEMPLATE_LAPACK_LAE2_HEADER
37
38
39
template
<
class
Treal>
40
int
template_lapack_lae2
(
const
Treal *a,
const
Treal *b,
const
Treal *c__,
41
Treal *rt1, Treal *rt2)
42
{
43
/* -- LAPACK auxiliary routine (version 3.0) --
44
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
45
Courant Institute, Argonne National Lab, and Rice University
46
October 31, 1992
47
48
49
Purpose
50
=======
51
52
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
53
[ A B ]
54
[ B C ].
55
On return, RT1 is the eigenvalue of larger absolute value, and RT2
56
is the eigenvalue of smaller absolute value.
57
58
Arguments
59
=========
60
61
A (input) DOUBLE PRECISION
62
The (1,1) element of the 2-by-2 matrix.
63
64
B (input) DOUBLE PRECISION
65
The (1,2) and (2,1) elements of the 2-by-2 matrix.
66
67
C (input) DOUBLE PRECISION
68
The (2,2) element of the 2-by-2 matrix.
69
70
RT1 (output) DOUBLE PRECISION
71
The eigenvalue of larger absolute value.
72
73
RT2 (output) DOUBLE PRECISION
74
The eigenvalue of smaller absolute value.
75
76
Further Details
77
===============
78
79
RT1 is accurate to a few ulps barring over/underflow.
80
81
RT2 may be inaccurate if there is massive cancellation in the
82
determinant A*C-B*B; higher precision or correctly rounded or
83
correctly truncated arithmetic would be needed to compute RT2
84
accurately in all cases.
85
86
Overflow is possible only if RT1 is within a factor of 5 of overflow.
87
Underflow is harmless if the input data is 0 or exceeds
88
underflow_threshold / macheps.
89
90
=====================================================================
91
92
93
Compute the eigenvalues */
94
/* System generated locals */
95
Treal d__1;
96
/* Local variables */
97
Treal acmn, acmx, ab, df, tb, sm, rt, adf;
98
99
100
sm = *a + *c__;
101
df = *a - *c__;
102
adf =
absMACRO
(df);
103
tb = *b + *b;
104
ab =
absMACRO
(tb);
105
if
(
absMACRO
(*a) >
absMACRO
(*c__)) {
106
acmx = *a;
107
acmn = *c__;
108
}
else
{
109
acmx = *c__;
110
acmn = *a;
111
}
112
if
(adf > ab) {
113
/* Computing 2nd power */
114
d__1 = ab / adf;
115
rt = adf *
template_blas_sqrt
(d__1 * d__1 + 1.);
116
}
else
if
(adf < ab) {
117
/* Computing 2nd power */
118
d__1 = adf / ab;
119
rt = ab *
template_blas_sqrt
(d__1 * d__1 + 1.);
120
}
else
{
121
122
/* Includes case AB=ADF=0 */
123
124
rt = ab *
template_blas_sqrt
(2.);
125
}
126
if
(sm < 0.) {
127
*rt1 = (sm - rt) * .5;
128
129
/* Order of execution important.
130
To get fully accurate smaller eigenvalue,
131
next line needs to be executed in higher precision. */
132
133
*rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
134
}
else
if
(sm > 0.) {
135
*rt1 = (sm + rt) * .5;
136
137
/* Order of execution important.
138
To get fully accurate smaller eigenvalue,
139
next line needs to be executed in higher precision. */
140
141
*rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
142
}
else
{
143
144
/* Includes case RT1 = RT2 = 0 */
145
146
*rt1 = rt * .5;
147
*rt2 = rt * -.5;
148
}
149
return
0;
150
151
/* End of DLAE2 */
152
153
}
/* dlae2_ */
154
155
#endif
source
matrix
template_lapack
lapack
template_lapack_lae2.h
Generated on Wed Nov 21 2012 20:40:23 for ergo by
1.8.1.2