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