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