ergo
template_lapack_lamch.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 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
47template<class Treal>
48Treal 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
59template<class Treal>
60Treal template_blas_lg10(Treal *x)
61{
62 return( log10e * template_blas_log(*x) );
63}
64
65
66
67
68
69
70
71
72template<class Treal>
73int 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
167template<class Treal>
168double 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
201template<class Treal>
202Treal 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
Treal template_blas_log(Treal x)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
int template_lapack_lassq(const integer *n, const Treal *x, const integer *incx, Treal *scale, Treal *sumsq)
Definition: template_lapack_lamch.h:73
Treal template_blas_lg10(Treal *x)
Definition: template_lapack_lamch.h:60
#define log10e
Definition: template_lapack_lamch.h:58
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
double template_lapack_pow_di(Treal *ap, integer *bp)
Definition: template_lapack_lamch.h:168
Treal template_lapack_d_sign(const Treal *a, const Treal *b)
Definition: template_lapack_lamch.h:48