ergo
|
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