ergo
template_lapack_laev2.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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
#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
int integer
Definition: template_blas_common.h:40
Treal template_blas_sqrt(Treal x)