ergo
template_lapack_lae2.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_LAE2_HEADER
36 #define TEMPLATE_LAPACK_LAE2_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__,
41  Treal *rt1, Treal *rt2)
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  DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
53  [ A B ]
54  [ B C ].
55  On return, RT1 is the eigenvalue of larger absolute value, and RT2
56  is the eigenvalue of smaller absolute value.
57 
58  Arguments
59  =========
60 
61  A (input) DOUBLE PRECISION
62  The (1,1) element of the 2-by-2 matrix.
63 
64  B (input) DOUBLE PRECISION
65  The (1,2) and (2,1) elements of the 2-by-2 matrix.
66 
67  C (input) DOUBLE PRECISION
68  The (2,2) element of the 2-by-2 matrix.
69 
70  RT1 (output) DOUBLE PRECISION
71  The eigenvalue of larger absolute value.
72 
73  RT2 (output) DOUBLE PRECISION
74  The eigenvalue of smaller absolute value.
75 
76  Further Details
77  ===============
78 
79  RT1 is accurate to a few ulps barring over/underflow.
80 
81  RT2 may be inaccurate if there is massive cancellation in the
82  determinant A*C-B*B; higher precision or correctly rounded or
83  correctly truncated arithmetic would be needed to compute RT2
84  accurately in all cases.
85 
86  Overflow is possible only if RT1 is within a factor of 5 of overflow.
87  Underflow is harmless if the input data is 0 or exceeds
88  underflow_threshold / macheps.
89 
90  =====================================================================
91 
92 
93  Compute the eigenvalues */
94  /* System generated locals */
95  Treal d__1;
96  /* Local variables */
97  Treal acmn, acmx, ab, df, tb, sm, rt, adf;
98 
99 
100  sm = *a + *c__;
101  df = *a - *c__;
102  adf = absMACRO(df);
103  tb = *b + *b;
104  ab = absMACRO(tb);
105  if (absMACRO(*a) > absMACRO(*c__)) {
106  acmx = *a;
107  acmn = *c__;
108  } else {
109  acmx = *c__;
110  acmn = *a;
111  }
112  if (adf > ab) {
113 /* Computing 2nd power */
114  d__1 = ab / adf;
115  rt = adf * template_blas_sqrt(d__1 * d__1 + 1.);
116  } else if (adf < ab) {
117 /* Computing 2nd power */
118  d__1 = adf / ab;
119  rt = ab * template_blas_sqrt(d__1 * d__1 + 1.);
120  } else {
121 
122 /* Includes case AB=ADF=0 */
123 
124  rt = ab * template_blas_sqrt(2.);
125  }
126  if (sm < 0.) {
127  *rt1 = (sm - rt) * .5;
128 
129 /* Order of execution important.
130  To get fully accurate smaller eigenvalue,
131  next line needs to be executed in higher precision. */
132 
133  *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
134  } else if (sm > 0.) {
135  *rt1 = (sm + rt) * .5;
136 
137 /* Order of execution important.
138  To get fully accurate smaller eigenvalue,
139  next line needs to be executed in higher precision. */
140 
141  *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b;
142  } else {
143 
144 /* Includes case RT1 = RT2 = 0 */
145 
146  *rt1 = rt * .5;
147  *rt2 = rt * -.5;
148  }
149  return 0;
150 
151 /* End of DLAE2 */
152 
153 } /* dlae2_ */
154 
155 #endif
#define absMACRO(x)
Definition: template_blas_common.h:45
int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, Treal *rt1, Treal *rt2)
Definition: template_lapack_lae2.h:40
Treal template_blas_sqrt(Treal x)