1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
/*! \file
Copyright (c) 2003, The Regents of the University of California, through
Lawrence Berkeley National Laboratory (subject to receipt of any required
approvals from U.S. Dept. of Energy)
All rights reserved.
The source code is distributed under BSD license, see the file License.txt
at the top-level directory.
*/
#include <float.h>
#include <math.h>
#include <stdio.h>
double dmach(char *cmach)
{
/* -- SuperLU auxiliary routine (version 5.0) --
This uses C99 standard constants, and is thread safe.
Must be compiled with -std=c99 flag.
Purpose
=======
DMACH returns double precision machine parameters.
Arguments
=========
CMACH (input) CHARACTER*1
Specifies the value to be returned by DMACH:
= 'E' or 'e', DMACH := eps
= 'S' or 's , DMACH := sfmin
= 'B' or 'b', DMACH := base
= 'P' or 'p', DMACH := eps*base
= 'N' or 'n', DMACH := t
= 'R' or 'r', DMACH := rnd
= 'M' or 'm', DMACH := emin
= 'U' or 'u', DMACH := rmin
= 'L' or 'l', DMACH := emax
= 'O' or 'o', DMACH := rmax
where
eps = relative machine precision
sfmin = safe minimum, such that 1/sfmin does not overflow
base = base of the machine
prec = eps*base
t = number of (base) digits in the mantissa
rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
emin = minimum exponent before (gradual) underflow
rmin = underflow threshold - base**(emin-1)
emax = largest exponent before overflow
rmax = overflow threshold - (base**emax)*(1-eps)
=====================================================================
*/
double sfmin, small, rmach;
if (strncmp(cmach, "E", 1)==0) {
rmach = DBL_EPSILON * 0.5;
} else if (strncmp(cmach, "S", 1)==0) {
sfmin = DBL_MIN;
small = 1. / DBL_MAX;
if (small >= sfmin) {
/* Use SMALL plus a bit, to avoid the possibility of rounding
causing overflow when computing 1/sfmin. */
sfmin = small * (DBL_EPSILON*0.5 + 1.);
}
rmach = sfmin;
} else if (strncmp(cmach, "B", 1)==0) {
rmach = FLT_RADIX;
} else if (strncmp(cmach, "P", 1)==0) {
rmach = DBL_EPSILON * 0.5 * FLT_RADIX;
} else if (strncmp(cmach, "N", 1)==0) {
rmach = DBL_MANT_DIG;
} else if (strncmp(cmach, "R", 1)==0) {
rmach = FLT_ROUNDS;
} else if (strncmp(cmach, "M", 1)==0) {
rmach = DBL_MIN_EXP;
} else if (strncmp(cmach, "U", 1)==0) {
rmach = DBL_MIN;
} else if (strncmp(cmach, "L", 1)==0) {
rmach = DBL_MAX_EXP;
} else if (strncmp(cmach, "O", 1)==0) {
rmach = DBL_MAX;
}
return rmach;
} /* end dmach */
|
