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 *
35 SUBROUTINE DLARTG( F, G, CS, SN, R )
36 *
37 * -- LAPACK auxiliary routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LARTG_TCC
44 #define FLENS_LAPACK_EIG_LARTG_TCC 1
45
46 #include <flens/aux/aux.h>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename T>
55 void
56 lartg_generic(const T &f, const T &g, T &cs, T &sn, T &r)
57 {
58 using std::abs;
59 using std::max;
60
61 const T Zero(0), One(1), Two(2);
62 static bool first = true;
63
64 static T eps, safeMin, safeMin2, safeMax2;
65
66 if (first) {
67 safeMin = lamch<T>(SafeMin);
68 eps = lamch<T>(Eps);
69
70 const int exp = explicit_cast<T,int>(log(safeMin/eps)
71 / log(lamch<T>(Base))
72 / Two);
73 safeMin2 = pow(lamch<T>(Base), exp);
74 safeMax2 = One / safeMin2;
75 first = false;
76 }
77
78 if (g==Zero) {
79 cs = One;
80 sn = Zero;
81 r = f;
82 } else if (f==Zero) {
83 cs = Zero;
84 sn = One;
85 r = g;
86 } else {
87 T f1 = f;
88 T g1 = g;
89 T scale = max(abs(f1), abs(g1));
90 if (scale>=safeMax2) {
91 int count = 0;
92 do {
93 ++count;
94 f1 *= safeMin2;
95 g1 *= safeMin2;
96 scale = max(abs(f1), abs(g1));
97 } while (scale>=safeMax2);
98 r = sqrt(pow(f1,2) + pow(g1,2));
99 cs = f1/r;
100 sn = g1/r;
101 for (int i=1; i<=count; ++i) {
102 r *= safeMax2;
103 }
104 } else if (scale<=safeMin2) {
105 int count = 0;
106 do {
107 ++count;
108 f1 *= safeMax2;
109 g1 *= safeMax2;
110 scale = max(abs(f1), abs(g1));
111 } while (scale<=safeMin2);
112 r = sqrt(pow(f1,2) + pow(g1,2));
113 cs = f1/r;
114 sn = g1/r;
115 for (int i=1; i<=count; ++i) {
116 r *= safeMin2;
117 }
118 } else {
119 r = sqrt(pow(f1,2) + pow(g1,2));
120 cs = f1/r;
121 sn = g1/r;
122 }
123 if (abs(f)>abs(g) && cs<Zero) {
124 cs = -cs;
125 sn = -sn;
126 r = -r;
127 }
128 }
129 }
130
131 //== interface for native lapack ===============================================
132
133 #ifdef CHECK_CXXLAPACK
134
135 template <typename T>
136 void
137 lartg_native(const T &f, const T &g, T &cs, T &sn, T &r)
138 {
139 if (IsSame<T, DOUBLE>::value) {
140 LAPACK_IMPL(dlartg)(&f, &g, &cs, &sn, &r);
141 } else {
142 ASSERT(0);
143 }
144 }
145
146 #endif // CHECK_CXXLAPACK
147
148 //== public interface ==========================================================
149
150 template <typename T>
151 void
152 lartg(const T &f, const T &g, T &cs, T &sn, T &r)
153 {
154 LAPACK_DEBUG_OUT("lartg");
155
156 # ifdef CHECK_CXXLAPACK
157 //
158 // Make copies of output arguments
159 //
160 T _cs = cs;
161 T _sn = sn;
162 T _r = r;
163 # endif
164
165 //
166 // Call implementation
167 //
168 lartg_generic(f, g, cs, sn, r);
169
170 # ifdef CHECK_CXXLAPACK
171 //
172 // Compare results
173 //
174 lartg_native(f, g, _cs, _sn, _r);
175
176 bool failed = false;
177 if (! isIdentical(cs, _cs, " cs", "_cs")) {
178 std::cerr << "CXXLAPACK: cs = " << cs << std::endl;
179 std::cerr << "F77LAPACK: _cs = " << _cs << std::endl;
180 failed = true;
181 }
182
183 if (! isIdentical(sn, _sn, " sn", "_sn")) {
184 std::cerr << "CXXLAPACK: sn = " << sn << std::endl;
185 std::cerr << "F77LAPACK: _sn = " << _sn << std::endl;
186 failed = true;
187 }
188
189 if (! isIdentical(r, _r, " r", "_r")) {
190 std::cerr << "CXXLAPACK: r = " << r << std::endl;
191 std::cerr << "F77LAPACK: _r = " << _r << std::endl;
192 failed = true;
193 }
194
195 if (failed) {
196 ASSERT(0);
197 }
198 # endif
199 }
200
201 } } // namespace lapack, flens
202
203 #endif // FLENS_LAPACK_EIG_LARTG_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 *
35 SUBROUTINE DLARTG( F, G, CS, SN, R )
36 *
37 * -- LAPACK auxiliary routine (version 3.2) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * November 2006
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LARTG_TCC
44 #define FLENS_LAPACK_EIG_LARTG_TCC 1
45
46 #include <flens/aux/aux.h>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename T>
55 void
56 lartg_generic(const T &f, const T &g, T &cs, T &sn, T &r)
57 {
58 using std::abs;
59 using std::max;
60
61 const T Zero(0), One(1), Two(2);
62 static bool first = true;
63
64 static T eps, safeMin, safeMin2, safeMax2;
65
66 if (first) {
67 safeMin = lamch<T>(SafeMin);
68 eps = lamch<T>(Eps);
69
70 const int exp = explicit_cast<T,int>(log(safeMin/eps)
71 / log(lamch<T>(Base))
72 / Two);
73 safeMin2 = pow(lamch<T>(Base), exp);
74 safeMax2 = One / safeMin2;
75 first = false;
76 }
77
78 if (g==Zero) {
79 cs = One;
80 sn = Zero;
81 r = f;
82 } else if (f==Zero) {
83 cs = Zero;
84 sn = One;
85 r = g;
86 } else {
87 T f1 = f;
88 T g1 = g;
89 T scale = max(abs(f1), abs(g1));
90 if (scale>=safeMax2) {
91 int count = 0;
92 do {
93 ++count;
94 f1 *= safeMin2;
95 g1 *= safeMin2;
96 scale = max(abs(f1), abs(g1));
97 } while (scale>=safeMax2);
98 r = sqrt(pow(f1,2) + pow(g1,2));
99 cs = f1/r;
100 sn = g1/r;
101 for (int i=1; i<=count; ++i) {
102 r *= safeMax2;
103 }
104 } else if (scale<=safeMin2) {
105 int count = 0;
106 do {
107 ++count;
108 f1 *= safeMax2;
109 g1 *= safeMax2;
110 scale = max(abs(f1), abs(g1));
111 } while (scale<=safeMin2);
112 r = sqrt(pow(f1,2) + pow(g1,2));
113 cs = f1/r;
114 sn = g1/r;
115 for (int i=1; i<=count; ++i) {
116 r *= safeMin2;
117 }
118 } else {
119 r = sqrt(pow(f1,2) + pow(g1,2));
120 cs = f1/r;
121 sn = g1/r;
122 }
123 if (abs(f)>abs(g) && cs<Zero) {
124 cs = -cs;
125 sn = -sn;
126 r = -r;
127 }
128 }
129 }
130
131 //== interface for native lapack ===============================================
132
133 #ifdef CHECK_CXXLAPACK
134
135 template <typename T>
136 void
137 lartg_native(const T &f, const T &g, T &cs, T &sn, T &r)
138 {
139 if (IsSame<T, DOUBLE>::value) {
140 LAPACK_IMPL(dlartg)(&f, &g, &cs, &sn, &r);
141 } else {
142 ASSERT(0);
143 }
144 }
145
146 #endif // CHECK_CXXLAPACK
147
148 //== public interface ==========================================================
149
150 template <typename T>
151 void
152 lartg(const T &f, const T &g, T &cs, T &sn, T &r)
153 {
154 LAPACK_DEBUG_OUT("lartg");
155
156 # ifdef CHECK_CXXLAPACK
157 //
158 // Make copies of output arguments
159 //
160 T _cs = cs;
161 T _sn = sn;
162 T _r = r;
163 # endif
164
165 //
166 // Call implementation
167 //
168 lartg_generic(f, g, cs, sn, r);
169
170 # ifdef CHECK_CXXLAPACK
171 //
172 // Compare results
173 //
174 lartg_native(f, g, _cs, _sn, _r);
175
176 bool failed = false;
177 if (! isIdentical(cs, _cs, " cs", "_cs")) {
178 std::cerr << "CXXLAPACK: cs = " << cs << std::endl;
179 std::cerr << "F77LAPACK: _cs = " << _cs << std::endl;
180 failed = true;
181 }
182
183 if (! isIdentical(sn, _sn, " sn", "_sn")) {
184 std::cerr << "CXXLAPACK: sn = " << sn << std::endl;
185 std::cerr << "F77LAPACK: _sn = " << _sn << std::endl;
186 failed = true;
187 }
188
189 if (! isIdentical(r, _r, " r", "_r")) {
190 std::cerr << "CXXLAPACK: r = " << r << std::endl;
191 std::cerr << "F77LAPACK: _r = " << _r << std::endl;
192 failed = true;
193 }
194
195 if (failed) {
196 ASSERT(0);
197 }
198 # endif
199 }
200
201 } } // namespace lapack, flens
202
203 #endif // FLENS_LAPACK_EIG_LARTG_TCC