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 /* Based on
34 DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
35 *
36 * -- LAPACK auxiliary routine (version 3.2) --
37 * -- LAPACK is a software package provided by Univ. of Tennessee, --
38 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
39 * November 2006
40 */
41
42 #ifndef FLENS_LAPACK_AUX_LAPY2_TCC
43 #define FLENS_LAPACK_AUX_LAPY2_TCC 1
44
45 #include <cmath>
46
47 namespace flens { namespace lapack {
48
49 //== generic lapack implementation =============================================
50
51 template <typename T>
52 T
53 lapy2_generic(const T &x, const T &y)
54 {
55 using std::abs;
56 using std::max;
57 using std::min;
58 using std::pow;
59 using std::sqrt;
60
61 const T xAbs = abs(x);
62 const T yAbs = abs(y);
63
64 const T w = max(xAbs, yAbs);
65 const T z = min(xAbs, yAbs);
66
67 if (z==T(0)) {
68 return w;
69 }
70
71 return w*sqrt(T(1) + pow(z/w,2));
72 }
73
74 //== interface for native lapack ===============================================
75
76 #ifdef CHECK_CXXLAPACK
77
78 template <typename T>
79 T
80 lapy2_native(const T &x, const T &y)
81 {
82 if (IsSame<T, DOUBLE>::value) {
83 return LAPACK_IMPL(dlapy2)(&x,
84 &y);
85 } else {
86 ASSERT(0);
87 }
88 }
89
90 #endif // CHECK_CXXLAPACK
91
92 //== public interface ==========================================================
93
94 template <typename T>
95 T
96 lapy2(const T &x, const T &y)
97 {
98 LAPACK_DEBUG_OUT("lapy2");
99
100 # ifdef CHECK_CXXLAPACK
101 //
102 // Make copies of output arguments
103 //
104 T _x = x;
105 T _y = y;
106 # endif
107
108 //
109 // Call implementation
110 //
111 const T result = lapy2_generic(x, y);
112
113 # ifdef CHECK_CXXLAPACK
114 //
115 // Compare results
116 //
117 const T _result = lapy2_native(_x, _y);
118
119 bool failed = false;
120 if (! isIdentical(x, _x, " x", "_x")) {
121 std::cerr << "CXXLAPACK: x = " << x << std::endl;
122 std::cerr << "F77LAPACK: _x = " << _x << std::endl;
123 failed = true;
124 }
125 if (! isIdentical(y, _y, " y", "_y")) {
126 std::cerr << "CXXLAPACK: y = " << y << std::endl;
127 std::cerr << "F77LAPACK: _y = " << _y << std::endl;
128 failed = true;
129 }
130 if (! isIdentical(result, _result, " result", "_result")) {
131 std::cerr << "CXXLAPACK: result = " << result << std::endl;
132 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
133 failed = true;
134 }
135
136 if (failed) {
137 ASSERT(0);
138 }
139 # endif
140
141 return result;
142 }
143
144 } } // namespace lapack, flens
145
146 #endif // FLENS_LAPACK_AUX_LAPY2_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 /* Based on
34 DOUBLE PRECISION FUNCTION DLAPY2( X, Y )
35 *
36 * -- LAPACK auxiliary routine (version 3.2) --
37 * -- LAPACK is a software package provided by Univ. of Tennessee, --
38 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
39 * November 2006
40 */
41
42 #ifndef FLENS_LAPACK_AUX_LAPY2_TCC
43 #define FLENS_LAPACK_AUX_LAPY2_TCC 1
44
45 #include <cmath>
46
47 namespace flens { namespace lapack {
48
49 //== generic lapack implementation =============================================
50
51 template <typename T>
52 T
53 lapy2_generic(const T &x, const T &y)
54 {
55 using std::abs;
56 using std::max;
57 using std::min;
58 using std::pow;
59 using std::sqrt;
60
61 const T xAbs = abs(x);
62 const T yAbs = abs(y);
63
64 const T w = max(xAbs, yAbs);
65 const T z = min(xAbs, yAbs);
66
67 if (z==T(0)) {
68 return w;
69 }
70
71 return w*sqrt(T(1) + pow(z/w,2));
72 }
73
74 //== interface for native lapack ===============================================
75
76 #ifdef CHECK_CXXLAPACK
77
78 template <typename T>
79 T
80 lapy2_native(const T &x, const T &y)
81 {
82 if (IsSame<T, DOUBLE>::value) {
83 return LAPACK_IMPL(dlapy2)(&x,
84 &y);
85 } else {
86 ASSERT(0);
87 }
88 }
89
90 #endif // CHECK_CXXLAPACK
91
92 //== public interface ==========================================================
93
94 template <typename T>
95 T
96 lapy2(const T &x, const T &y)
97 {
98 LAPACK_DEBUG_OUT("lapy2");
99
100 # ifdef CHECK_CXXLAPACK
101 //
102 // Make copies of output arguments
103 //
104 T _x = x;
105 T _y = y;
106 # endif
107
108 //
109 // Call implementation
110 //
111 const T result = lapy2_generic(x, y);
112
113 # ifdef CHECK_CXXLAPACK
114 //
115 // Compare results
116 //
117 const T _result = lapy2_native(_x, _y);
118
119 bool failed = false;
120 if (! isIdentical(x, _x, " x", "_x")) {
121 std::cerr << "CXXLAPACK: x = " << x << std::endl;
122 std::cerr << "F77LAPACK: _x = " << _x << std::endl;
123 failed = true;
124 }
125 if (! isIdentical(y, _y, " y", "_y")) {
126 std::cerr << "CXXLAPACK: y = " << y << std::endl;
127 std::cerr << "F77LAPACK: _y = " << _y << std::endl;
128 failed = true;
129 }
130 if (! isIdentical(result, _result, " result", "_result")) {
131 std::cerr << "CXXLAPACK: result = " << result << std::endl;
132 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
133 failed = true;
134 }
135
136 if (failed) {
137 ASSERT(0);
138 }
139 # endif
140
141 return result;
142 }
143
144 } } // namespace lapack, flens
145
146 #endif // FLENS_LAPACK_AUX_LAPY2_TCC