template_lapack_laev2.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_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

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