ergo
template_lapack_lamch.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_LAMCH_HEADER
36 #define TEMPLATE_LAPACK_LAMCH_HEADER
37 
38 
39 #include <stdio.h>
40 #include <iostream>
41 #include <limits>
42 
43 
44 
45 template<class Treal>
46 Treal template_lapack_d_sign(const Treal *a, const Treal *b)
47 {
48  Treal x;
49  x = (*a >= 0 ? *a : - *a);
50  return( *b >= 0 ? x : -x);
51 }
52 
53 
54 
55 #define log10e 0.43429448190325182765
56 template<class Treal>
57 Treal template_blas_lg10(Treal *x)
58 {
59  return( log10e * template_blas_log(*x) );
60 }
61 
62 
63 
64 
65 
66 
67 
68 
69 template<class Treal>
70 int template_lapack_lassq(const integer *n, const Treal *x, const integer *incx,
71  Treal *scale, Treal *sumsq)
72 {
73 /* -- LAPACK auxiliary routine (version 3.0) --
74  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
75  Courant Institute, Argonne National Lab, and Rice University
76  June 30, 1999
77 
78 
79  Purpose
80  =======
81 
82  DLASSQ returns the values scl and smsq such that
83 
84  ( scl**2 )*smsq = x( 1 )**2 +...+ x( n )**2 + ( scale**2 )*sumsq,
85 
86  where x( i ) = X( 1 + ( i - 1 )*INCX ). The value of sumsq is
87  assumed to be non-negative and scl returns the value
88 
89  scl = max( scale, abs( x( i ) ) ).
90 
91  scale and sumsq must be supplied in SCALE and SUMSQ and
92  scl and smsq are overwritten on SCALE and SUMSQ respectively.
93 
94  The routine makes only one pass through the vector x.
95 
96  Arguments
97  =========
98 
99  N (input) INTEGER
100  The number of elements to be used from the vector X.
101 
102  X (input) DOUBLE PRECISION array, dimension (N)
103  The vector for which a scaled sum of squares is computed.
104  x( i ) = X( 1 + ( i - 1 )*INCX ), 1 <= i <= n.
105 
106  INCX (input) INTEGER
107  The increment between successive values of the vector X.
108  INCX > 0.
109 
110  SCALE (input/output) DOUBLE PRECISION
111  On entry, the value scale in the equation above.
112  On exit, SCALE is overwritten with scl , the scaling factor
113  for the sum of squares.
114 
115  SUMSQ (input/output) DOUBLE PRECISION
116  On entry, the value sumsq in the equation above.
117  On exit, SUMSQ is overwritten with smsq , the basic sum of
118  squares from which scl has been factored out.
119 
120  =====================================================================
121 
122 
123  Parameter adjustments */
124  /* System generated locals */
125  integer i__1, i__2;
126  Treal d__1;
127  /* Local variables */
128  Treal absxi;
129  integer ix;
130 
131  --x;
132 
133  /* Function Body */
134  if (*n > 0) {
135  i__1 = (*n - 1) * *incx + 1;
136  i__2 = *incx;
137  for (ix = 1; i__2 < 0 ? ix >= i__1 : ix <= i__1; ix += i__2) {
138  if (x[ix] != 0.) {
139  absxi = (d__1 = x[ix], absMACRO(d__1));
140  if (*scale < absxi) {
141 /* Computing 2nd power */
142  d__1 = *scale / absxi;
143  *sumsq = *sumsq * (d__1 * d__1) + 1;
144  *scale = absxi;
145  } else {
146 /* Computing 2nd power */
147  d__1 = absxi / *scale;
148  *sumsq += d__1 * d__1;
149  }
150  }
151 /* L10: */
152  }
153  }
154  return 0;
155 
156 /* End of DLASSQ */
157 
158 } /* dlassq_ */
159 
160 
161 
162 
163 
164 template<class Treal>
165 double template_lapack_pow_di(Treal *ap, integer *bp)
166 {
167  Treal pow, x;
168  integer n;
169  unsigned long u;
170 
171  pow = 1;
172  x = *ap;
173  n = *bp;
174 
175  if(n != 0)
176  {
177  if(n < 0)
178  {
179  n = -n;
180  x = 1/x;
181  }
182  for(u = n; ; )
183  {
184  if(u & 01)
185  pow *= x;
186  if(u >>= 1)
187  x *= x;
188  else
189  break;
190  }
191  }
192  return(pow);
193 }
194 
195 
196 
197 
198 template<class Treal>
199 Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
200 {
201 
202 /* -- LAPACK auxiliary routine (version 3.0) --
203  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
204  Courant Institute, Argonne National Lab, and Rice University
205  October 31, 1992
206 
207 
208  Purpose
209  =======
210 
211  DLAMCH determines double precision machine parameters.
212 
213  Arguments
214  =========
215 
216  CMACH (input) CHARACTER*1
217  Specifies the value to be returned by DLAMCH:
218  = 'E' or 'e', DLAMCH := eps
219  = 'S' or 's , DLAMCH := sfmin
220  = 'B' or 'b', DLAMCH := base
221  = 'P' or 'p', DLAMCH := eps*base
222  = 'N' or 'n', DLAMCH := t
223  = 'R' or 'r', DLAMCH := rnd
224  = 'M' or 'm', DLAMCH := emin
225  = 'U' or 'u', DLAMCH := rmin
226  = 'L' or 'l', DLAMCH := emax
227  = 'O' or 'o', DLAMCH := rmax
228 
229  where
230 
231  eps = relative machine precision
232  sfmin = safe minimum, such that 1/sfmin does not overflow
233  base = base of the machine
234  prec = eps*base
235  t = number of (base) digits in the mantissa
236  rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
237  emin = minimum exponent before (gradual) underflow
238  rmin = underflow threshold - base**(emin-1)
239  emax = largest exponent before overflow
240  rmax = overflow threshold - (base**emax)*(1-eps)
241 
242  =====================================================================
243 */
244 
245  Treal rmach, ret_val;
246 
247  /* Initialization added by Elias to get rid of compiler warnings. */
248  rmach = 0;
249  if (template_blas_lsame(cmach, "E")) { /* Epsilon */
250  rmach = std::numeric_limits<Treal>::epsilon();
251  } else if (template_blas_lsame(cmach, "S")) { /* Safe minimum */
253  } else if (template_blas_lsame(cmach, "B")) { /* Base */
254  /* Assume "base" is 2 */
255  rmach = 2.0;
256  } else if (template_blas_lsame(cmach, "P")) { /* Precision */
257  /* Assume "base" is 2 */
258  rmach = 2.0 * std::numeric_limits<Treal>::epsilon();
259  } else if (template_blas_lsame(cmach, "N")) {
260  std::cout << "ERROR in template_lapack_lamch: case N not implemented." << std::endl;
261  throw "ERROR in template_lapack_lamch: case N not implemented.";
262  } else if (template_blas_lsame(cmach, "R")) {
263  std::cout << "ERROR in template_lapack_lamch: case R not implemented." << std::endl;
264  throw "ERROR in template_lapack_lamch: case R not implemented.";
265  } else if (template_blas_lsame(cmach, "M")) {
266  std::cout << "ERROR in template_lapack_lamch: case M not implemented." << std::endl;
267  throw "ERROR in template_lapack_lamch: case M not implemented.";
268  } else if (template_blas_lsame(cmach, "U")) {
269  std::cout << "ERROR in template_lapack_lamch: case U not implemented." << std::endl;
270  throw "ERROR in template_lapack_lamch: case U not implemented.";
271  } else if (template_blas_lsame(cmach, "L")) {
272  std::cout << "ERROR in template_lapack_lamch: case L not implemented." << std::endl;
273  throw "ERROR in template_lapack_lamch: case L not implemented.";
274  } else if (template_blas_lsame(cmach, "O")) {
275  std::cout << "ERROR in template_lapack_lamch: case O not implemented." << std::endl;
276  throw "ERROR in template_lapack_lamch: case O not implemented.";
277  }
278 
279  ret_val = rmach;
280  return ret_val;
281 
282  /* End of DLAMCH */
283 
284 } /* dlamch_ */
285 
286 
287 
288 #endif