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_LAEV2_HEADER 00038 #define TEMPLATE_LAPACK_LAEV2_HEADER 00039 00040 00041 template<class Treal> 00042 int template_lapack_laev2(Treal *a, Treal *b, Treal *c__, 00043 Treal *rt1, Treal *rt2, Treal *cs1, Treal *sn1) 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 DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix 00055 [ A B ] 00056 [ B C ]. 00057 On return, RT1 is the eigenvalue of larger absolute value, RT2 is the 00058 eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right 00059 eigenvector for RT1, giving the decomposition 00060 00061 [ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ] 00062 [-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ]. 00063 00064 Arguments 00065 ========= 00066 00067 A (input) DOUBLE PRECISION 00068 The (1,1) element of the 2-by-2 matrix. 00069 00070 B (input) DOUBLE PRECISION 00071 The (1,2) element and the conjugate of the (2,1) element of 00072 the 2-by-2 matrix. 00073 00074 C (input) DOUBLE PRECISION 00075 The (2,2) element of the 2-by-2 matrix. 00076 00077 RT1 (output) DOUBLE PRECISION 00078 The eigenvalue of larger absolute value. 00079 00080 RT2 (output) DOUBLE PRECISION 00081 The eigenvalue of smaller absolute value. 00082 00083 CS1 (output) DOUBLE PRECISION 00084 SN1 (output) DOUBLE PRECISION 00085 The vector (CS1, SN1) is a unit right eigenvector for RT1. 00086 00087 Further Details 00088 =============== 00089 00090 RT1 is accurate to a few ulps barring over/underflow. 00091 00092 RT2 may be inaccurate if there is massive cancellation in the 00093 determinant A*C-B*B; higher precision or correctly rounded or 00094 correctly truncated arithmetic would be needed to compute RT2 00095 accurately in all cases. 00096 00097 CS1 and SN1 are accurate to a few ulps barring over/underflow. 00098 00099 Overflow is possible only if RT1 is within a factor of 5 of overflow. 00100 Underflow is harmless if the input data is 0 or exceeds 00101 underflow_threshold / macheps. 00102 00103 ===================================================================== 00104 00105 00106 Compute the eigenvalues */ 00107 /* System generated locals */ 00108 Treal d__1; 00109 /* Local variables */ 00110 Treal acmn, acmx, ab, df, cs, ct, tb, sm, tn, rt, adf, acs; 00111 integer sgn1, sgn2; 00112 00113 00114 sm = *a + *c__; 00115 df = *a - *c__; 00116 adf = absMACRO(df); 00117 tb = *b + *b; 00118 ab = absMACRO(tb); 00119 if (absMACRO(*a) > absMACRO(*c__)) { 00120 acmx = *a; 00121 acmn = *c__; 00122 } else { 00123 acmx = *c__; 00124 acmn = *a; 00125 } 00126 if (adf > ab) { 00127 /* Computing 2nd power */ 00128 d__1 = ab / adf; 00129 rt = adf * template_blas_sqrt(d__1 * d__1 + 1.); 00130 } else if (adf < ab) { 00131 /* Computing 2nd power */ 00132 d__1 = adf / ab; 00133 rt = ab * template_blas_sqrt(d__1 * d__1 + 1.); 00134 } else { 00135 00136 /* Includes case AB=ADF=0 */ 00137 00138 rt = ab * template_blas_sqrt(2.); 00139 } 00140 if (sm < 0.) { 00141 *rt1 = (sm - rt) * .5; 00142 sgn1 = -1; 00143 00144 /* Order of execution important. 00145 To get fully accurate smaller eigenvalue, 00146 next line needs to be executed in higher precision. */ 00147 00148 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b; 00149 } else if (sm > 0.) { 00150 *rt1 = (sm + rt) * .5; 00151 sgn1 = 1; 00152 00153 /* Order of execution important. 00154 To get fully accurate smaller eigenvalue, 00155 next line needs to be executed in higher precision. */ 00156 00157 *rt2 = acmx / *rt1 * acmn - *b / *rt1 * *b; 00158 } else { 00159 00160 /* Includes case RT1 = RT2 = 0 */ 00161 00162 *rt1 = rt * .5; 00163 *rt2 = rt * -.5; 00164 sgn1 = 1; 00165 } 00166 00167 /* Compute the eigenvector */ 00168 00169 if (df >= 0.) { 00170 cs = df + rt; 00171 sgn2 = 1; 00172 } else { 00173 cs = df - rt; 00174 sgn2 = -1; 00175 } 00176 acs = absMACRO(cs); 00177 if (acs > ab) { 00178 ct = -tb / cs; 00179 *sn1 = 1. / template_blas_sqrt(ct * ct + 1.); 00180 *cs1 = ct * *sn1; 00181 } else { 00182 if (ab == 0.) { 00183 *cs1 = 1.; 00184 *sn1 = 0.; 00185 } else { 00186 tn = -cs / tb; 00187 *cs1 = 1. / template_blas_sqrt(tn * tn + 1.); 00188 *sn1 = tn * *cs1; 00189 } 00190 } 00191 if (sgn1 == sgn2) { 00192 tn = *cs1; 00193 *cs1 = -(*sn1); 00194 *sn1 = tn; 00195 } 00196 return 0; 00197 00198 /* End of DLAEV2 */ 00199 00200 } /* dlaev2_ */ 00201 00202 #endif