template_lapack_larfg.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_LARFG_HEADER
00036 #define TEMPLATE_LAPACK_LARFG_HEADER
00037 
00038 #include "template_lapack_lapy2.h"
00039 
00040 template<class Treal>
00041 int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, 
00042         const integer *incx, Treal *tau)
00043 {
00044 /*  -- LAPACK auxiliary routine (version 3.0) --   
00045        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00046        Courant Institute, Argonne National Lab, and Rice University   
00047        September 30, 1994   
00048 
00049 
00050     Purpose   
00051     =======   
00052 
00053     DLARFG generates a real elementary reflector H of order n, such   
00054     that   
00055 
00056           H * ( alpha ) = ( beta ),   H' * H = I.   
00057               (   x   )   (   0  )   
00058 
00059     where alpha and beta are scalars, and x is an (n-1)-element real   
00060     vector. H is represented in the form   
00061 
00062           H = I - tau * ( 1 ) * ( 1 v' ) ,   
00063                         ( v )   
00064 
00065     where tau is a real scalar and v is a real (n-1)-element   
00066     vector.   
00067 
00068     If the elements of x are all zero, then tau = 0 and H is taken to be   
00069     the unit matrix.   
00070 
00071     Otherwise  1 <= tau <= 2.   
00072 
00073     Arguments   
00074     =========   
00075 
00076     N       (input) INTEGER   
00077             The order of the elementary reflector.   
00078 
00079     ALPHA   (input/output) DOUBLE PRECISION   
00080             On entry, the value alpha.   
00081             On exit, it is overwritten with the value beta.   
00082 
00083     X       (input/output) DOUBLE PRECISION array, dimension   
00084                            (1+(N-2)*abs(INCX))   
00085             On entry, the vector x.   
00086             On exit, it is overwritten with the vector v.   
00087 
00088     INCX    (input) INTEGER   
00089             The increment between elements of X. INCX > 0.   
00090 
00091     TAU     (output) DOUBLE PRECISION   
00092             The value tau.   
00093 
00094     =====================================================================   
00095 
00096 
00097        Parameter adjustments */
00098     /* System generated locals */
00099     integer i__1;
00100     Treal d__1;
00101     /* Local variables */
00102      Treal beta;
00103      integer j;
00104      Treal xnorm;
00105      Treal safmin, rsafmn;
00106      integer knt;
00107 
00108     --x;
00109 
00110     /* Function Body */
00111     if (*n <= 1) {
00112         *tau = 0.;
00113         return 0;
00114     }
00115 
00116     i__1 = *n - 1;
00117     xnorm = template_blas_nrm2(&i__1, &x[1], incx);
00118 
00119     if (xnorm == 0.) {
00120 
00121 /*        H  =  I */
00122 
00123         *tau = 0.;
00124     } else {
00125 
00126 /*        general case */
00127 
00128         d__1 = template_lapack_lapy2(alpha, &xnorm);
00129         beta = -template_lapack_d_sign(&d__1, alpha);
00130         safmin = template_lapack_lamch("S", (Treal)0) / template_lapack_lamch("E", (Treal)0);
00131         if (absMACRO(beta) < safmin) {
00132 
00133 /*           XNORM, BETA may be inaccurate; scale X and recompute them */
00134 
00135             rsafmn = 1. / safmin;
00136             knt = 0;
00137 L10:
00138             ++knt;
00139             i__1 = *n - 1;
00140             template_blas_scal(&i__1, &rsafmn, &x[1], incx);
00141             beta *= rsafmn;
00142             *alpha *= rsafmn;
00143             if (absMACRO(beta) < safmin) {
00144                 goto L10;
00145             }
00146 
00147 /*           New BETA is at most 1, at least SAFMIN */
00148 
00149             i__1 = *n - 1;
00150             xnorm = template_blas_nrm2(&i__1, &x[1], incx);
00151             d__1 = template_lapack_lapy2(alpha, &xnorm);
00152             beta = -template_lapack_d_sign(&d__1, alpha);
00153             *tau = (beta - *alpha) / beta;
00154             i__1 = *n - 1;
00155             d__1 = 1. / (*alpha - beta);
00156             template_blas_scal(&i__1, &d__1, &x[1], incx);
00157 
00158 /*           If ALPHA is subnormal, it may lose relative accuracy */
00159 
00160             *alpha = beta;
00161             i__1 = knt;
00162             for (j = 1; j <= i__1; ++j) {
00163                 *alpha *= safmin;
00164 /* L20: */
00165             }
00166         } else {
00167             *tau = (beta - *alpha) / beta;
00168             i__1 = *n - 1;
00169             d__1 = 1. / (*alpha - beta);
00170             template_blas_scal(&i__1, &d__1, &x[1], incx);
00171             *alpha = beta;
00172         }
00173     }
00174 
00175     return 0;
00176 
00177 /*     End of DLARFG */
00178 
00179 } /* dlarfg_ */
00180 
00181 #endif

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