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