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 <complex>
51 #include <limits>
52
53 #include <flens/blas/blas.h>
54 #include <flens/lapack/lapack.h>
55
56 namespace flens { namespace lapack {
57
58 //== generic lapack implementation =============================================
59
60 template <typename T>
61 T
62 lamch_generic(MachineParameter machineParameter)
63 {
64 //
65 // Assume rounding, not chopping. Always.
66 //
67 const T eps = std::numeric_limits<T>::epsilon() * T(0.5);
68
69 if (machineParameter==Eps) {
70 return eps;
71 } else if (machineParameter==SafeMin) {
72 T safeMin = std::numeric_limits<T>::min();
73 const T small = T(1) / std::numeric_limits<T>::max();
74 if (small>=safeMin) {
75 //
76 // Use SMALL plus a bit, to avoid the possibility of rounding
77 // causing overflow when computing 1/sfmin.
78 //
79 safeMin = small*(T(1)+eps);
80 }
81 return safeMin;
82 } else if (machineParameter==Base) {
83 return std::numeric_limits<T>::radix;
84 } else if (machineParameter==Precision) {
85 return eps*std::numeric_limits<T>::radix;
86 } else if (machineParameter==Mantissa) {
87 return std::numeric_limits<T>::digits;
88 } else if (machineParameter==Rounding) {
89 return T(1);
90 } else if (machineParameter==UnderflowExp) {
91 return std::numeric_limits<T>::min_exponent;
92 } else if (machineParameter==UnderflowThreshold) {
93 return std::numeric_limits<T>::min();
94 } else if (machineParameter==OverflowExp) {
95 return std::numeric_limits<T>::max_exponent;
96 } else if (machineParameter==OverflowThreshold) {
97 return std::numeric_limits<T>::max();
98 }
99 ASSERT(0);
100 return T(0);
101 }
102
103 //== interface for native lapack ===============================================
104
105 #ifdef CHECK_CXXLAPACK
106
107 template <typename T>
108 T
109 lamch_native(MachineParameter machineParameter)
110 {
111 char c = char(machineParameter);
112
113 if (IsSame<T, double>::value) {
114 return LAPACK_IMPL(dlamch)(&c);
115 }
116
117 ASSERT(0);
118 return 0;
119 }
120
121 #endif // CHECK_CXXLAPACK
122
123 //== public interface ==========================================================
124
125 template <typename TT>
126 typename LAMCH::Real<TT>::Type
127 lamch(MachineParameter machineParameter)
128 {
129 typedef typename LAMCH::Real<TT>::Type T;
130 //
131 // Call implementation
132 //
133 const T cxx_result = lamch_generic<T>(machineParameter);
134
135 # ifdef CHECK_CXXLAPACK
136 //
137 // Compare results
138 //
139 const T f77_result = lamch_native<T>(machineParameter);
140
141 if (! isIdentical(cxx_result, f77_result, "cxx_result", "f77_result")) {
142 ASSERT(0);
143 }
144 # endif
145
146 return cxx_result;
147 }
148
149 } } // namespace lapack, flens
150
151 #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 <complex>
51 #include <limits>
52
53 #include <flens/blas/blas.h>
54 #include <flens/lapack/lapack.h>
55
56 namespace flens { namespace lapack {
57
58 //== generic lapack implementation =============================================
59
60 template <typename T>
61 T
62 lamch_generic(MachineParameter machineParameter)
63 {
64 //
65 // Assume rounding, not chopping. Always.
66 //
67 const T eps = std::numeric_limits<T>::epsilon() * T(0.5);
68
69 if (machineParameter==Eps) {
70 return eps;
71 } else if (machineParameter==SafeMin) {
72 T safeMin = std::numeric_limits<T>::min();
73 const T small = T(1) / std::numeric_limits<T>::max();
74 if (small>=safeMin) {
75 //
76 // Use SMALL plus a bit, to avoid the possibility of rounding
77 // causing overflow when computing 1/sfmin.
78 //
79 safeMin = small*(T(1)+eps);
80 }
81 return safeMin;
82 } else if (machineParameter==Base) {
83 return std::numeric_limits<T>::radix;
84 } else if (machineParameter==Precision) {
85 return eps*std::numeric_limits<T>::radix;
86 } else if (machineParameter==Mantissa) {
87 return std::numeric_limits<T>::digits;
88 } else if (machineParameter==Rounding) {
89 return T(1);
90 } else if (machineParameter==UnderflowExp) {
91 return std::numeric_limits<T>::min_exponent;
92 } else if (machineParameter==UnderflowThreshold) {
93 return std::numeric_limits<T>::min();
94 } else if (machineParameter==OverflowExp) {
95 return std::numeric_limits<T>::max_exponent;
96 } else if (machineParameter==OverflowThreshold) {
97 return std::numeric_limits<T>::max();
98 }
99 ASSERT(0);
100 return T(0);
101 }
102
103 //== interface for native lapack ===============================================
104
105 #ifdef CHECK_CXXLAPACK
106
107 template <typename T>
108 T
109 lamch_native(MachineParameter machineParameter)
110 {
111 char c = char(machineParameter);
112
113 if (IsSame<T, double>::value) {
114 return LAPACK_IMPL(dlamch)(&c);
115 }
116
117 ASSERT(0);
118 return 0;
119 }
120
121 #endif // CHECK_CXXLAPACK
122
123 //== public interface ==========================================================
124
125 template <typename TT>
126 typename LAMCH::Real<TT>::Type
127 lamch(MachineParameter machineParameter)
128 {
129 typedef typename LAMCH::Real<TT>::Type T;
130 //
131 // Call implementation
132 //
133 const T cxx_result = lamch_generic<T>(machineParameter);
134
135 # ifdef CHECK_CXXLAPACK
136 //
137 // Compare results
138 //
139 const T f77_result = lamch_native<T>(machineParameter);
140
141 if (! isIdentical(cxx_result, f77_result, "cxx_result", "f77_result")) {
142 ASSERT(0);
143 }
144 # endif
145
146 return cxx_result;
147 }
148
149 } } // namespace lapack, flens
150
151 #endif // FLENS_LAPACK_AUX_LAMCH_TCC