ergo
template_lapack_laev2.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
41template<class Treal>
42int 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
Treal template_blas_sqrt(Treal x)
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1)
Definition template_lapack_laev2.h:42