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