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