template_lapack_lamch.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_LAMCH_HEADER
00036 #define TEMPLATE_LAPACK_LAMCH_HEADER
00037 
00038 
00039 #include <stdio.h>
00040 #include <iostream>
00041 #include <limits>
00042 
00043 
00044 
00045 template<class Treal>
00046 Treal template_lapack_d_sign(const Treal *a, const Treal *b)
00047 {
00048   Treal x;
00049   x = (*a >= 0 ? *a : - *a);
00050   return( *b >= 0 ? x : -x);
00051 }
00052 
00053 
00054 
00055 #define log10e 0.43429448190325182765
00056 template<class Treal>
00057 Treal template_blas_lg10(Treal *x)
00058 {
00059   return( log10e * template_blas_log(*x) );
00060 }
00061 
00062 
00063 
00064 
00065 
00066 
00067 
00068 
00069 template<class Treal>
00070 int template_lapack_lassq(const integer *n, const Treal *x, const integer *incx, 
00071         Treal *scale, Treal *sumsq)
00072 {
00073 /*  -- LAPACK auxiliary routine (version 3.0) --   
00074        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00075        Courant Institute, Argonne National Lab, and Rice University   
00076        June 30, 1999   
00077 
00078 
00079     Purpose   
00080     =======   
00081 
00082     DLASSQ  returns the values  scl  and  smsq  such that   
00083 
00084        ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,   
00085 
00086     where  x( i ) = X( 1 + ( i - 1 )*INCX ). The value of  sumsq  is   
00087     assumed to be non-negative and  scl  returns the value   
00088 
00089        scl = max( scale, abs( x( i ) ) ).   
00090 
00091     scale and sumsq must be supplied in SCALE and SUMSQ and   
00092     scl and smsq are overwritten on SCALE and SUMSQ respectively.   
00093 
00094     The routine makes only one pass through the vector x.   
00095 
00096     Arguments   
00097     =========   
00098 
00099     N       (input) INTEGER   
00100             The number of elements to be used from the vector X.   
00101 
00102     X       (input) DOUBLE PRECISION array, dimension (N)   
00103             The vector for which a scaled sum of squares is computed.   
00104                x( i )  = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.   
00105 
00106     INCX    (input) INTEGER   
00107             The increment between successive values of the vector X.   
00108             INCX > 0.   
00109 
00110     SCALE   (input/output) DOUBLE PRECISION   
00111             On entry, the value  scale  in the equation above.   
00112             On exit, SCALE is overwritten with  scl , the scaling factor   
00113             for the sum of squares.   
00114 
00115     SUMSQ   (input/output) DOUBLE PRECISION   
00116             On entry, the value  sumsq  in the equation above.   
00117             On exit, SUMSQ is overwritten with  smsq , the basic sum of   
00118             squares from which  scl  has been factored out.   
00119 
00120    =====================================================================   
00121 
00122 
00123        Parameter adjustments */
00124     /* System generated locals */
00125     integer i__1, i__2;
00126     Treal d__1;
00127     /* Local variables */
00128      Treal absxi;
00129      integer ix;
00130 
00131     --x;
00132 
00133     /* Function Body */
00134     if (*n > 0) {
00135         i__1 = (*n - 1) * *incx + 1;
00136         i__2 = *incx;
00137         for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
00138             if (x[ix] != 0.) {
00139                 absxi = (d__1 = x[ix], absMACRO(d__1));
00140                 if (*scale < absxi) {
00141 /* Computing 2nd power */
00142                     d__1 = *scale / absxi;
00143                     *sumsq = *sumsq * (d__1 * d__1) + 1;
00144                     *scale = absxi;
00145                 } else {
00146 /* Computing 2nd power */
00147                     d__1 = absxi / *scale;
00148                     *sumsq += d__1 * d__1;
00149                 }
00150             }
00151 /* L10: */
00152         }
00153     }
00154     return 0;
00155 
00156 /*     End of DLASSQ */
00157 
00158 } /* dlassq_ */
00159 
00160 
00161 
00162 
00163 
00164 template<class Treal>
00165 double template_lapack_pow_di(Treal *ap, integer *bp)
00166 {
00167   Treal pow, x;
00168   integer n;
00169   unsigned long u;
00170 
00171   pow = 1;
00172   x = *ap;
00173   n = *bp;
00174 
00175   if(n != 0)
00176     {
00177       if(n < 0)
00178         {
00179           n = -n;
00180           x = 1/x;
00181         }
00182       for(u = n; ; )
00183         {
00184           if(u & 01)
00185             pow *= x;
00186           if(u >>= 1)
00187             x *= x;
00188           else
00189             break;
00190         }
00191     }
00192   return(pow);
00193 }
00194 
00195 
00196 
00197 
00198 template<class Treal>
00199 Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
00200 {
00201 
00202 /*  -- LAPACK auxiliary routine (version 3.0) --   
00203        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
00204        Courant Institute, Argonne National Lab, and Rice University   
00205        October 31, 1992   
00206 
00207 
00208     Purpose   
00209     =======   
00210 
00211     DLAMCH determines double precision machine parameters.   
00212 
00213     Arguments   
00214     =========   
00215 
00216     CMACH   (input) CHARACTER*1   
00217             Specifies the value to be returned by DLAMCH:   
00218             = 'E' or 'e',   DLAMCH := eps   
00219             = 'S' or 's ,   DLAMCH := sfmin   
00220             = 'B' or 'b',   DLAMCH := base   
00221             = 'P' or 'p',   DLAMCH := eps*base   
00222             = 'N' or 'n',   DLAMCH := t   
00223             = 'R' or 'r',   DLAMCH := rnd   
00224             = 'M' or 'm',   DLAMCH := emin   
00225             = 'U' or 'u',   DLAMCH := rmin   
00226             = 'L' or 'l',   DLAMCH := emax   
00227             = 'O' or 'o',   DLAMCH := rmax   
00228 
00229             where   
00230 
00231             eps   = relative machine precision   
00232             sfmin = safe minimum, such that 1/sfmin does not overflow   
00233             base  = base of the machine   
00234             prec  = eps*base   
00235             t     = number of (base) digits in the mantissa   
00236             rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise   
00237             emin  = minimum exponent before (gradual) underflow   
00238             rmin  = underflow threshold - base**(emin-1)   
00239             emax  = largest exponent before overflow   
00240             rmax  = overflow threshold  - (base**emax)*(1-eps)   
00241 
00242    ===================================================================== 
00243 */
00244 
00245    Treal rmach, ret_val;
00246 
00247    /* Initialization added by Elias to get rid of compiler warnings. */
00248    rmach = 0;
00249   if (template_blas_lsame(cmach, "E")) { /* Epsilon */
00250     rmach = std::numeric_limits<Treal>::epsilon();
00251   } else if (template_blas_lsame(cmach, "S")) { /* Safe minimum */
00252     rmach = std::numeric_limits<Treal>::min();
00253   } else if (template_blas_lsame(cmach, "B")) { /* Base */
00254     /* Assume "base" is 2 */
00255     rmach = 2.0;
00256   } else if (template_blas_lsame(cmach, "P")) { /* Precision */
00257     /* Assume "base" is 2 */
00258     rmach = 2.0 * std::numeric_limits<Treal>::epsilon();
00259   } else if (template_blas_lsame(cmach, "N")) {
00260     std::cout << "ERROR in template_lapack_lamch: case N not implemented." << std::endl;
00261     throw "ERROR in template_lapack_lamch: case N not implemented.";
00262   } else if (template_blas_lsame(cmach, "R")) {
00263     std::cout << "ERROR in template_lapack_lamch: case R not implemented." << std::endl;
00264     throw "ERROR in template_lapack_lamch: case R not implemented.";
00265   } else if (template_blas_lsame(cmach, "M")) {
00266     std::cout << "ERROR in template_lapack_lamch: case M not implemented." << std::endl;
00267     throw "ERROR in template_lapack_lamch: case M not implemented.";
00268   } else if (template_blas_lsame(cmach, "U")) {
00269     std::cout << "ERROR in template_lapack_lamch: case U not implemented." << std::endl;
00270     throw "ERROR in template_lapack_lamch: case U not implemented.";
00271   } else if (template_blas_lsame(cmach, "L")) {
00272     std::cout << "ERROR in template_lapack_lamch: case L not implemented." << std::endl;
00273     throw "ERROR in template_lapack_lamch: case L not implemented.";
00274   } else if (template_blas_lsame(cmach, "O")) {
00275     std::cout << "ERROR in template_lapack_lamch: case O not implemented." << std::endl;
00276     throw "ERROR in template_lapack_lamch: case O not implemented.";
00277   }
00278 
00279   ret_val = rmach;
00280   return ret_val;
00281 
00282   /*     End of DLAMCH */
00283 
00284 } /* dlamch_ */
00285 
00286 
00287 
00288 #endif

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