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