ergo
template_lapack_laev2.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_LAEV2_HEADER
38
#define TEMPLATE_LAPACK_LAEV2_HEADER
39
40
41
template
<
class
Treal>
42
int
template_lapack_laev2
(Treal *a, Treal *b, Treal *c__,
43
Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
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
DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
55
[ A B ]
56
[ B C ].
57
On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
58
eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
59
eigenvector for RT1, giving the decomposition
60
61
[ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
62
[-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
63
64
Arguments
65
=========
66
67
A (input) DOUBLE PRECISION
68
The (1,1) element of the 2-by-2 matrix.
69
70
B (input) DOUBLE PRECISION
71
The (1,2) element and the conjugate of the (2,1) element of
72
the 2-by-2 matrix.
73
74
C (input) DOUBLE PRECISION
75
The (2,2) element of the 2-by-2 matrix.
76
77
RT1 (output) DOUBLE PRECISION
78
The eigenvalue of larger absolute value.
79
80
RT2 (output) DOUBLE PRECISION
81
The eigenvalue of smaller absolute value.
82
83
CS1 (output) DOUBLE PRECISION
84
SN1 (output) DOUBLE PRECISION
85
The vector (CS1, SN1) is a unit right eigenvector for RT1.
86
87
Further Details
88
===============
89
90
RT1 is accurate to a few ulps barring over/underflow.
91
92
RT2 may be inaccurate if there is massive cancellation in the
93
determinant A*C-B*B; higher precision or correctly rounded or
94
correctly truncated arithmetic would be needed to compute RT2
95
accurately in all cases.
96
97
CS1 and SN1 are accurate to a few ulps barring over/underflow.
98
99
Overflow is possible only if RT1 is within a factor of 5 of overflow.
100
Underflow is harmless if the input data is 0 or exceeds
101
underflow_threshold / macheps.
102
103
=====================================================================
104
105
106
Compute the eigenvalues */
107
/* System generated locals */
108
Treal d__1;
109
/* Local variables */
110
Treal acmn, acmx, ab, df, cs, ct, tb, sm, tn, rt, adf, acs;
111
integer
sgn1, sgn2;
112
113
114
sm = *a + *c__;
115
df = *a - *c__;
116
adf =
absMACRO
(df);
117
tb = *b + *b;
118
ab =
absMACRO
(tb);
119
if
(
absMACRO
(*a) >
absMACRO
(*c__)) {
120
acmx = *a;
121
acmn = *c__;
122
}
else
{
123
acmx = *c__;
124
acmn = *a;
125
}
126
if
(adf > ab) {
127
/* Computing 2nd power */
128
d__1 = ab / adf;
129
rt = adf *
template_blas_sqrt
(d__1 * d__1 + 1.);
130
}
else
if
(adf < ab) {
131
/* Computing 2nd power */
132
d__1 = adf / ab;
133
rt = ab *
template_blas_sqrt
(d__1 * d__1 + 1.);
134
}
else
{
135
136
/* Includes case AB=ADF=0 */
137
138
rt = ab *
template_blas_sqrt
(2.);
139
}
140
if
(sm < 0.) {
141
*rt1 = (sm - rt) * .5;
142
sgn1 = -1;
143
144
/* Order of execution important.
145
To get fully accurate smaller eigenvalue,
146
next line needs to be executed in higher precision. */
147
148
*rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
149
}
else
if
(sm > 0.) {
150
*rt1 = (sm + rt) * .5;
151
sgn1 = 1;
152
153
/* Order of execution important.
154
To get fully accurate smaller eigenvalue,
155
next line needs to be executed in higher precision. */
156
157
*rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
158
}
else
{
159
160
/* Includes case RT1 = RT2 = 0 */
161
162
*rt1 = rt * .5;
163
*rt2 = rt * -.5;
164
sgn1 = 1;
165
}
166
167
/* Compute the eigenvector */
168
169
if
(df >= 0.) {
170
cs = df + rt;
171
sgn2 = 1;
172
}
else
{
173
cs = df - rt;
174
sgn2 = -1;
175
}
176
acs =
absMACRO
(cs);
177
if
(acs > ab) {
178
ct = -tb / cs;
179
*sn1 = 1. /
template_blas_sqrt
(ct * ct + 1.);
180
*cs1 = ct * *sn1;
181
}
else
{
182
if
(ab == 0.) {
183
*cs1 = 1.;
184
*sn1 = 0.;
185
}
else
{
186
tn = -cs / tb;
187
*cs1 = 1. /
template_blas_sqrt
(tn * tn + 1.);
188
*sn1 = tn * *cs1;
189
}
190
}
191
if
(sgn1 == sgn2) {
192
tn = *cs1;
193
*cs1 = -(*sn1);
194
*sn1 = tn;
195
}
196
return
0;
197
198
/* End of DLAEV2 */
199
200
}
/* dlaev2_ */
201
202
#endif
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
template_lapack_laev2
int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition:
template_lapack_laev2.h:42
absMACRO
#define absMACRO(x)
Definition:
template_blas_common.h:47
integer
int integer
Definition:
template_blas_common.h:40
source
matrix
template_lapack
lapack
template_lapack_laev2.h
Generated on Sun Aug 16 2020 00:00:00 for ergo by
1.8.18