ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_LAE2_HEADER 00036 #define TEMPLATE_LAPACK_LAE2_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_lae2(const Treal *a, const Treal *b, const Treal *c__, 00041 Treal *rt1, Treal *rt2) 00042 { 00043 /* -- LAPACK auxiliary routine (version 3.0) -- 00044 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00045 Courant Institute, Argonne National Lab, and Rice University 00046 October 31, 1992 00047 00048 00049 Purpose 00050 ======= 00051 00052 DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix 00053 [ A B ] 00054 [ B C ]. 00055 On return, RT1 is the eigenvalue of larger absolute value, and RT2 00056 is the eigenvalue of smaller absolute value. 00057 00058 Arguments 00059 ========= 00060 00061 A (input) DOUBLE PRECISION 00062 The (1,1) element of the 2-by-2 matrix. 00063 00064 B (input) DOUBLE PRECISION 00065 The (1,2) and (2,1) elements of the 2-by-2 matrix. 00066 00067 C (input) DOUBLE PRECISION 00068 The (2,2) element of the 2-by-2 matrix. 00069 00070 RT1 (output) DOUBLE PRECISION 00071 The eigenvalue of larger absolute value. 00072 00073 RT2 (output) DOUBLE PRECISION 00074 The eigenvalue of smaller absolute value. 00075 00076 Further Details 00077 =============== 00078 00079 RT1 is accurate to a few ulps barring over/underflow. 00080 00081 RT2 may be inaccurate if there is massive cancellation in the 00082 determinant A*C-B*B; higher precision or correctly rounded or 00083 correctly truncated arithmetic would be needed to compute RT2 00084 accurately in all cases. 00085 00086 Overflow is possible only if RT1 is within a factor of 5 of overflow. 00087 Underflow is harmless if the input data is 0 or exceeds 00088 underflow_threshold / macheps. 00089 00090 ===================================================================== 00091 00092 00093 Compute the eigenvalues */ 00094 /* System generated locals */ 00095 Treal d__1; 00096 /* Local variables */ 00097 Treal acmn, acmx, ab, df, tb, sm, rt, adf; 00098 00099 00100 sm = *a + *c__; 00101 df = *a - *c__; 00102 adf = absMACRO(df); 00103 tb = *b + *b; 00104 ab = absMACRO(tb); 00105 if (absMACRO(*a) > absMACRO(*c__)) { 00106 acmx = *a; 00107 acmn = *c__; 00108 } else { 00109 acmx = *c__; 00110 acmn = *a; 00111 } 00112 if (adf > ab) { 00113 /* Computing 2nd power */ 00114 d__1 = ab / adf; 00115 rt = adf * template_blas_sqrt(d__1 * d__1 + 1.); 00116 } else if (adf < ab) { 00117 /* Computing 2nd power */ 00118 d__1 = adf / ab; 00119 rt = ab * template_blas_sqrt(d__1 * d__1 + 1.); 00120 } else { 00121 00122 /* Includes case AB=ADF=0 */ 00123 00124 rt = ab * template_blas_sqrt(2.); 00125 } 00126 if (sm < 0.) { 00127 *rt1 = (sm - rt) * .5; 00128 00129 /* Order of execution important. 00130 To get fully accurate smaller eigenvalue, 00131 next line needs to be executed in higher precision. */ 00132 00133 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b; 00134 } else if (sm > 0.) { 00135 *rt1 = (sm + rt) * .5; 00136 00137 /* Order of execution important. 00138 To get fully accurate smaller eigenvalue, 00139 next line needs to be executed in higher precision. */ 00140 00141 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b; 00142 } else { 00143 00144 /* Includes case RT1 = RT2 = 0 */ 00145 00146 *rt1 = rt * .5; 00147 *rt2 = rt * -.5; 00148 } 00149 return 0; 00150 00151 /* End of DLAE2 */ 00152 00153 } /* dlae2_ */ 00154 00155 #endif