template_lapack_lae2.h

Go to the documentation of this file.
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

Generated on Wed Nov 21 09:32:11 2012 for ergo by  doxygen 1.4.7