template_lapack_lartg.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_LARTG_HEADER
00036 #define TEMPLATE_LAPACK_LARTG_HEADER
00037 
00038 
00039 template<class Treal>
00040 int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, 
00041         Treal *sn, Treal *r__)
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        September 30, 1994   
00047 
00048 
00049     Purpose   
00050     =======   
00051 
00052     DLARTG generate a plane rotation so that   
00053 
00054        [  CS  SN  ]  .  [ F ]  =  [ R ]   where CS**2 + SN**2 = 1.   
00055        [ -SN  CS  ]     [ G ]     [ 0 ]   
00056 
00057     This is a slower, more accurate version of the BLAS1 routine DROTG,   
00058     with the following other differences:   
00059        F and G are unchanged on return.   
00060        If G=0, then CS=1 and SN=0.   
00061        If F=0 and (G .ne. 0), then CS=0 and SN=1 without doing any   
00062           floating point operations (saves work in DBDSQR when   
00063           there are zeros on the diagonal).   
00064 
00065     If F exceeds G in magnitude, CS will be positive.   
00066 
00067     Arguments   
00068     =========   
00069 
00070     F       (input) DOUBLE PRECISION   
00071             The first component of vector to be rotated.   
00072 
00073     G       (input) DOUBLE PRECISION   
00074             The second component of vector to be rotated.   
00075 
00076     CS      (output) DOUBLE PRECISION   
00077             The cosine of the rotation.   
00078 
00079     SN      (output) DOUBLE PRECISION   
00080             The sine of the rotation.   
00081 
00082     R       (output) DOUBLE PRECISION   
00083             The nonzero component of the rotated vector.   
00084 
00085     ===================================================================== */
00086     /* Initialized data */
00087      logical first = TRUE_;
00088     /* System generated locals */
00089     integer i__1;
00090     Treal d__1, d__2;
00091     /* Local variables */
00092      integer i__;
00093      Treal scale;
00094      integer count;
00095      Treal f1, g1, safmn2, safmx2;
00096      Treal safmin, eps;
00097 
00098 
00099 
00100     if (first) {
00101         first = FALSE_;
00102         safmin = template_lapack_lamch("S", (Treal)0);
00103         eps = template_lapack_lamch("E", (Treal)0);
00104         d__1 = template_lapack_lamch("B", (Treal)0);
00105         i__1 = (integer) (template_blas_log(safmin / eps) / template_blas_log(template_lapack_lamch("B", (Treal)0)) / 
00106                 2.);
00107         safmn2 = template_lapack_pow_di(&d__1, &i__1);
00108         safmx2 = 1. / safmn2;
00109     }
00110     if (*g == 0.) {
00111         *cs = 1.;
00112         *sn = 0.;
00113         *r__ = *f;
00114     } else if (*f == 0.) {
00115         *cs = 0.;
00116         *sn = 1.;
00117         *r__ = *g;
00118     } else {
00119         f1 = *f;
00120         g1 = *g;
00121 /* Computing MAX */
00122         d__1 = absMACRO(f1), d__2 = absMACRO(g1);
00123         scale = maxMACRO(d__1,d__2);
00124         if (scale >= safmx2) {
00125             count = 0;
00126 L10:
00127             ++count;
00128             f1 *= safmn2;
00129             g1 *= safmn2;
00130 /* Computing MAX */
00131             d__1 = absMACRO(f1), d__2 = absMACRO(g1);
00132             scale = maxMACRO(d__1,d__2);
00133             if (scale >= safmx2) {
00134                 goto L10;
00135             }
00136 /* Computing 2nd power */
00137             d__1 = f1;
00138 /* Computing 2nd power */
00139             d__2 = g1;
00140             *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
00141             *cs = f1 / *r__;
00142             *sn = g1 / *r__;
00143             i__1 = count;
00144             for (i__ = 1; i__ <= i__1; ++i__) {
00145                 *r__ *= safmx2;
00146 /* L20: */
00147             }
00148         } else if (scale <= safmn2) {
00149             count = 0;
00150 L30:
00151             ++count;
00152             f1 *= safmx2;
00153             g1 *= safmx2;
00154 /* Computing MAX */
00155             d__1 = absMACRO(f1), d__2 = absMACRO(g1);
00156             scale = maxMACRO(d__1,d__2);
00157             if (scale <= safmn2) {
00158                 goto L30;
00159             }
00160 /* Computing 2nd power */
00161             d__1 = f1;
00162 /* Computing 2nd power */
00163             d__2 = g1;
00164             *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
00165             *cs = f1 / *r__;
00166             *sn = g1 / *r__;
00167             i__1 = count;
00168             for (i__ = 1; i__ <= i__1; ++i__) {
00169                 *r__ *= safmn2;
00170 /* L40: */
00171             }
00172         } else {
00173 /* Computing 2nd power */
00174             d__1 = f1;
00175 /* Computing 2nd power */
00176             d__2 = g1;
00177             *r__ = template_blas_sqrt(d__1 * d__1 + d__2 * d__2);
00178             *cs = f1 / *r__;
00179             *sn = g1 / *r__;
00180         }
00181         if (absMACRO(*f) > absMACRO(*g) && *cs < 0.) {
00182             *cs = -(*cs);
00183             *sn = -(*sn);
00184             *r__ = -(*r__);
00185         }
00186     }
00187     return 0;
00188 
00189 /*     End of DLARTG */
00190 
00191 } /* dlartg_ */
00192 
00193 #endif

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