1 /*
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33
34 /* Based on
35 *
36 DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
37 *
38 * -- LAPACK auxiliary routine (version 3.3.0) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * Based on LAPACK DLAMCH but with Fortran 95 query functions
42 * See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
43 * and http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
44 * July 2010
45 */
46
47 #ifndef FLENS_LAPACK_AUX_LAMCH_TCC
48 #define FLENS_LAPACK_AUX_LAMCH_TCC 1
49
50 #include <limits>
51
52 #include <flens/blas/blas.h>
53 #include <flens/lapack/lapack.h>
54
55 namespace flens { namespace lapack {
56
57 //== generic lapack implementation =============================================
58
59 template <typename T>
60 T
61 lamch_generic(MachineParameter machineParameter)
62 {
63 //
64 // Assume rounding, not chopping. Always.
65 //
66 const T eps = std::numeric_limits<T>::epsilon() * T(0.5);
67
68 if (machineParameter==Eps) {
69 return eps;
70 } else if (machineParameter==SafeMin) {
71 T safeMin = std::numeric_limits<T>::min();
72 const T small = T(1) / std::numeric_limits<T>::max();
73 if (small>=safeMin) {
74 //
75 // Use SMALL plus a bit, to avoid the possibility of rounding
76 // causing overflow when computing 1/sfmin.
77 //
78 safeMin = small*(T(1)+eps);
79 }
80 return safeMin;
81 } else if (machineParameter==Base) {
82 return std::numeric_limits<T>::radix;
83 } else if (machineParameter==Precision) {
84 return eps*std::numeric_limits<T>::radix;
85 } else if (machineParameter==Mantissa) {
86 return std::numeric_limits<T>::digits;
87 } else if (machineParameter==Rounding) {
88 return T(1);
89 } else if (machineParameter==UnderflowExp) {
90 return std::numeric_limits<T>::min_exponent;
91 } else if (machineParameter==UnderflowThreshold) {
92 return std::numeric_limits<T>::min();
93 } else if (machineParameter==OverflowExp) {
94 return std::numeric_limits<T>::max_exponent;
95 } else if (machineParameter==OverflowThreshold) {
96 return std::numeric_limits<T>::max();
97 }
98 ASSERT(0);
99 return T(0);
100 }
101
102 //== interface for native lapack ===============================================
103
104 #ifdef CHECK_CXXLAPACK
105
106 template <typename T>
107 T
108 lamch_native(MachineParameter machineParameter)
109 {
110 char c = char(machineParameter);
111
112 if (IsSame<T, double>::value) {
113 return LAPACK_IMPL(dlamch)(&c);
114 }
115
116 ASSERT(0);
117 return 0;
118 }
119
120 #endif // CHECK_CXXLAPACK
121
122 //== public interface ==========================================================
123
124 template <typename T>
125 T
126 lamch(MachineParameter machineParameter)
127 {
128 //
129 // Call implementation
130 //
131 const T cxx_result = lamch_generic<T>(machineParameter);
132
133 # ifdef CHECK_CXXLAPACK
134 //
135 // Compare results
136 //
137 const T f77_result = lamch_native<T>(machineParameter);
138
139 if (! isIdentical(cxx_result, f77_result, "cxx_result", "f77_result")) {
140 ASSERT(0);
141 }
142 # endif
143
144 return cxx_result;
145 }
146
147 } } // namespace lapack, flens
148
149 #endif // FLENS_LAPACK_AUX_LAMCH_TCC
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33
34 /* Based on
35 *
36 DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
37 *
38 * -- LAPACK auxiliary routine (version 3.3.0) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * Based on LAPACK DLAMCH but with Fortran 95 query functions
42 * See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
43 * and http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
44 * July 2010
45 */
46
47 #ifndef FLENS_LAPACK_AUX_LAMCH_TCC
48 #define FLENS_LAPACK_AUX_LAMCH_TCC 1
49
50 #include <limits>
51
52 #include <flens/blas/blas.h>
53 #include <flens/lapack/lapack.h>
54
55 namespace flens { namespace lapack {
56
57 //== generic lapack implementation =============================================
58
59 template <typename T>
60 T
61 lamch_generic(MachineParameter machineParameter)
62 {
63 //
64 // Assume rounding, not chopping. Always.
65 //
66 const T eps = std::numeric_limits<T>::epsilon() * T(0.5);
67
68 if (machineParameter==Eps) {
69 return eps;
70 } else if (machineParameter==SafeMin) {
71 T safeMin = std::numeric_limits<T>::min();
72 const T small = T(1) / std::numeric_limits<T>::max();
73 if (small>=safeMin) {
74 //
75 // Use SMALL plus a bit, to avoid the possibility of rounding
76 // causing overflow when computing 1/sfmin.
77 //
78 safeMin = small*(T(1)+eps);
79 }
80 return safeMin;
81 } else if (machineParameter==Base) {
82 return std::numeric_limits<T>::radix;
83 } else if (machineParameter==Precision) {
84 return eps*std::numeric_limits<T>::radix;
85 } else if (machineParameter==Mantissa) {
86 return std::numeric_limits<T>::digits;
87 } else if (machineParameter==Rounding) {
88 return T(1);
89 } else if (machineParameter==UnderflowExp) {
90 return std::numeric_limits<T>::min_exponent;
91 } else if (machineParameter==UnderflowThreshold) {
92 return std::numeric_limits<T>::min();
93 } else if (machineParameter==OverflowExp) {
94 return std::numeric_limits<T>::max_exponent;
95 } else if (machineParameter==OverflowThreshold) {
96 return std::numeric_limits<T>::max();
97 }
98 ASSERT(0);
99 return T(0);
100 }
101
102 //== interface for native lapack ===============================================
103
104 #ifdef CHECK_CXXLAPACK
105
106 template <typename T>
107 T
108 lamch_native(MachineParameter machineParameter)
109 {
110 char c = char(machineParameter);
111
112 if (IsSame<T, double>::value) {
113 return LAPACK_IMPL(dlamch)(&c);
114 }
115
116 ASSERT(0);
117 return 0;
118 }
119
120 #endif // CHECK_CXXLAPACK
121
122 //== public interface ==========================================================
123
124 template <typename T>
125 T
126 lamch(MachineParameter machineParameter)
127 {
128 //
129 // Call implementation
130 //
131 const T cxx_result = lamch_generic<T>(machineParameter);
132
133 # ifdef CHECK_CXXLAPACK
134 //
135 // Compare results
136 //
137 const T f77_result = lamch_native<T>(machineParameter);
138
139 if (! isIdentical(cxx_result, f77_result, "cxx_result", "f77_result")) {
140 ASSERT(0);
141 }
142 # endif
143
144 return cxx_result;
145 }
146
147 } } // namespace lapack, flens
148
149 #endif // FLENS_LAPACK_AUX_LAMCH_TCC