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 DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
36 *
37 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LANV2_TCC
44 #define FLENS_LAPACK_EIG_LANV2_TCC 1
45
46 namespace flens { namespace lapack {
47
48 //== generic lapack implementation =============================================
49
50 template <typename T>
51 void
52 lanv2_generic(T &a, T &b, T &c, T &d,
53 T &rt1r, T &rt1i,
54 T &rt2r, T &rt2i,
55 T &cs, T &sn)
56 {
57 using std::abs;
58 using std::max;
59 using std::min;
60 using std::sqrt;
61 using std::swap;
62
63 const T Zero(0), Half(0.5), One(1), Multpl(4);
64
65 const T eps = lamch<T>(Precision);
66 if (c==Zero) {
67 cs = One;
68 sn = Zero;
69 } else if (b==Zero) {
70 //
71 // Swap rows and columns
72 //
73 cs = Zero;
74 sn = One;
75 swap(a, d);
76 b = -c;
77 c = Zero;
78 } else if ((a-d==Zero) && (sign(One,b)!=sign(One,c))) {
79 cs = One;
80 sn = Zero;
81 } else {
82 T temp = a-d;
83 const T p = Half*temp;
84 const T bcMax = max(abs(b), abs(c));
85 const T bcMis = min(abs(b), abs(c)) * sign(One,b) * sign (One, c);
86 const T scale = max(abs(p), bcMax);
87 T z = (p/scale)*p + (bcMax/scale)*bcMis;
88 //
89 // If Z is of the order of the machine accuracy, postpone the
90 // decision on the nature of eigenvalues
91 //
92 if (z>=Multpl*eps) {
93 //
94 // Real eigenvalues. Compute A and D.
95 //
96 z = p + sign(sqrt(scale)*sqrt(z), p);
97 a = d + z;
98 d = d - (bcMax/z)*bcMis;
99 //
100 // Compute B and the rotation matrix
101 //
102 const T tau = lapy2(c, z);
103 cs = z/tau;
104 sn = c/tau;
105 b = b -c;
106 c = Zero;
107 } else {
108 //
109 // Complex eigenvalues, or real (almost) equal eigenvalues.
110 // Make diagonal elements equal.
111 //
112 const T sigma = b + c;
113 const T tau = lapy2(sigma, temp);
114 cs = sqrt(Half*(One+abs(sigma)/tau));
115 sn = -(p/(tau*cs)) * sign(One, sigma);
116 //
117 // Compute [ AA BB ] = [ A B ] [ CS -SN ]
118 // [ CC DD ] [ C D ] [ SN CS ]
119 //
120 const T aa = a*cs + b*sn;
121 const T bb = -a*sn + b*cs;
122 const T cc = c*cs + d*sn;
123 const T dd = -c*sn + d*cs;
124 //
125 // Compute [ A B ] = [ CS SN ] [ AA BB ]
126 // [ C D ] [-SN CS ] [ CC DD ]
127 //
128 a = aa*cs + cc*sn;
129 b = bb*cs + dd*sn;
130 c = -aa*sn + cc*cs;
131 d = -bb*sn + dd*cs;
132
133 temp = Half * (a+d);
134 a = temp;
135 d = temp;
136
137 if (c!=Zero) {
138 if (b!=Zero) {
139 if (sign(One, b)==sign(One,c)) {
140 //
141 // Real eigenvalues: reduce to upper triangular form
142 //
143 const T sab = sqrt(abs(b));
144 const T sac = sqrt(abs(c));
145 const T p = sign(sab*sac, c);
146 const T tau = One / sqrt(abs(b+c));
147 a = temp + p;
148 d = temp - p;
149 b = b -c;
150 c = Zero;
151 const T cs1 = sab*tau;
152 const T sn1 = sac*tau;
153 temp = cs*cs1 - sn*sn1;
154 sn = cs*sn1 + sn*cs1;
155 cs = temp;
156 }
157 } else {
158 b = -c;
159 c = Zero;
160 temp = cs;
161 cs = -sn;
162 sn = temp;
163 }
164 }
165 }
166 }
167 //
168 // Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
169 //
170 rt1r = a;
171 rt2r = d;
172 if (c==Zero) {
173 rt1i = Zero;
174 rt2i = Zero;
175 } else {
176 rt1i = sqrt(abs(b))*sqrt(abs(c));
177 rt2i = -rt1i;
178 }
179 }
180
181 //== interface for native lapack ===============================================
182
183 #ifdef CHECK_CXXLAPACK
184
185 template <typename T>
186 void
187 lanv2_native(T &a, T &b, T &c, T &d,
188 T &rt1r, T &rt1i,
189 T &rt2r, T &rt2i,
190 T &cs, T &sn)
191 {
192 if (IsSame<T,DOUBLE>::value) {
193 LAPACK_IMPL(dlanv2)(&a, &b, &c, &d,
194 &rt1r, &rt1i,
195 &rt2r, &rt2i,
196 &cs, &sn);
197 } else {
198 ASSERT(0);
199 }
200 }
201
202 #endif // CHECK_CXXLAPACK
203
204 //== public interface ==========================================================
205
206 template <typename T>
207 void
208 lanv2(T &a, T &b, T &c, T &d, T &rt1r, T &rt1i, T &rt2r, T &rt2i, T &cs, T &sn)
209 {
210 LAPACK_DEBUG_OUT("lanv2");
211
212 # ifdef CHECK_CXXLAPACK
213 //
214 // Make copies of output arguments
215 //
216 T _a = a;
217 T _b = b;
218 T _c = c;
219 T _d = d;
220 T _rt1r = rt1r;
221 T _rt1i = rt1i;
222 T _rt2r = rt2r;
223 T _rt2i = rt2i;
224 T _cs = cs;
225 T _sn = sn;
226 # endif
227
228 //
229 // Call implementation
230 //
231 lanv2_generic(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn);
232
233 # ifdef CHECK_CXXLAPACK
234 //
235 // Compare results
236 //
237 lanv2_native(_a, _b, _c, _d, _rt1r, _rt1i, _rt2r, _rt2i, _cs, _sn);
238
239 bool failed = false;
240 if (! isIdentical(a, _a, " a", "_a")) {
241 std::cerr << "CXXLAPACK: a = " << a << std::endl;
242 std::cerr << "F77LAPACK: _a = " << _a << std::endl;
243 failed = true;
244 }
245
246 if (! isIdentical(b, _b, " b", "_b")) {
247 std::cerr << "CXXLAPACK: b = " << b << std::endl;
248 std::cerr << "F77LAPACK: _b = " << _b << std::endl;
249 failed = true;
250 }
251
252 if (! isIdentical(c, _c, " c", "_c")) {
253 std::cerr << "CXXLAPACK: c = " << c << std::endl;
254 std::cerr << "F77LAPACK: _c = " << _c << std::endl;
255 failed = true;
256 }
257
258 if (! isIdentical(d, _d, " d", "_d")) {
259 std::cerr << "CXXLAPACK: d = " << d << std::endl;
260 std::cerr << "F77LAPACK: _d = " << _d << std::endl;
261 failed = true;
262 }
263
264 if (! isIdentical(rt1r, _rt1r, " rt1r", "_rt1r")) {
265 std::cerr << "CXXLAPACK: rt1r = " << rt1r << std::endl;
266 std::cerr << "F77LAPACK: _rt1r = " << _rt1r << std::endl;
267 failed = true;
268 }
269
270 if (! isIdentical(rt1i, _rt1i, " rt1i", "_rt1i")) {
271 std::cerr << "CXXLAPACK: rt1i = " << rt1i << std::endl;
272 std::cerr << "F77LAPACK: _rt1i = " << _rt1i << std::endl;
273 failed = true;
274 }
275
276 if (! isIdentical(rt2r, _rt2r, " rt2r", "_rt2r")) {
277 std::cerr << "CXXLAPACK: rt2r = " << rt2r << std::endl;
278 std::cerr << "F77LAPACK: _rt2r = " << _rt2r << std::endl;
279 failed = true;
280 }
281
282 if (! isIdentical(rt2i, _rt2i, " rt2i", "_rt2i")) {
283 std::cerr << "CXXLAPACK: rt2i = " << rt2i << std::endl;
284 std::cerr << "F77LAPACK: _rt2i = " << _rt2i << std::endl;
285 failed = true;
286 }
287
288 if (! isIdentical(cs, _cs, " cs", "_cs")) {
289 std::cerr << "CXXLAPACK: cs = " << cs << std::endl;
290 std::cerr << "F77LAPACK: _cs = " << _cs << std::endl;
291 failed = true;
292 }
293
294 if (! isIdentical(sn, _sn, " sn", "_sn")) {
295 std::cerr << "CXXLAPACK: sn = " << sn << std::endl;
296 std::cerr << "F77LAPACK: _sn = " << _sn << std::endl;
297 failed = true;
298 }
299
300 if (failed) {
301 ASSERT(0);
302 }
303 # endif
304 }
305
306 } } // namespace lapack, flens
307
308 #endif // FLENS_LAPACK_EIG_LANV2_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 DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
36 *
37 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LANV2_TCC
44 #define FLENS_LAPACK_EIG_LANV2_TCC 1
45
46 namespace flens { namespace lapack {
47
48 //== generic lapack implementation =============================================
49
50 template <typename T>
51 void
52 lanv2_generic(T &a, T &b, T &c, T &d,
53 T &rt1r, T &rt1i,
54 T &rt2r, T &rt2i,
55 T &cs, T &sn)
56 {
57 using std::abs;
58 using std::max;
59 using std::min;
60 using std::sqrt;
61 using std::swap;
62
63 const T Zero(0), Half(0.5), One(1), Multpl(4);
64
65 const T eps = lamch<T>(Precision);
66 if (c==Zero) {
67 cs = One;
68 sn = Zero;
69 } else if (b==Zero) {
70 //
71 // Swap rows and columns
72 //
73 cs = Zero;
74 sn = One;
75 swap(a, d);
76 b = -c;
77 c = Zero;
78 } else if ((a-d==Zero) && (sign(One,b)!=sign(One,c))) {
79 cs = One;
80 sn = Zero;
81 } else {
82 T temp = a-d;
83 const T p = Half*temp;
84 const T bcMax = max(abs(b), abs(c));
85 const T bcMis = min(abs(b), abs(c)) * sign(One,b) * sign (One, c);
86 const T scale = max(abs(p), bcMax);
87 T z = (p/scale)*p + (bcMax/scale)*bcMis;
88 //
89 // If Z is of the order of the machine accuracy, postpone the
90 // decision on the nature of eigenvalues
91 //
92 if (z>=Multpl*eps) {
93 //
94 // Real eigenvalues. Compute A and D.
95 //
96 z = p + sign(sqrt(scale)*sqrt(z), p);
97 a = d + z;
98 d = d - (bcMax/z)*bcMis;
99 //
100 // Compute B and the rotation matrix
101 //
102 const T tau = lapy2(c, z);
103 cs = z/tau;
104 sn = c/tau;
105 b = b -c;
106 c = Zero;
107 } else {
108 //
109 // Complex eigenvalues, or real (almost) equal eigenvalues.
110 // Make diagonal elements equal.
111 //
112 const T sigma = b + c;
113 const T tau = lapy2(sigma, temp);
114 cs = sqrt(Half*(One+abs(sigma)/tau));
115 sn = -(p/(tau*cs)) * sign(One, sigma);
116 //
117 // Compute [ AA BB ] = [ A B ] [ CS -SN ]
118 // [ CC DD ] [ C D ] [ SN CS ]
119 //
120 const T aa = a*cs + b*sn;
121 const T bb = -a*sn + b*cs;
122 const T cc = c*cs + d*sn;
123 const T dd = -c*sn + d*cs;
124 //
125 // Compute [ A B ] = [ CS SN ] [ AA BB ]
126 // [ C D ] [-SN CS ] [ CC DD ]
127 //
128 a = aa*cs + cc*sn;
129 b = bb*cs + dd*sn;
130 c = -aa*sn + cc*cs;
131 d = -bb*sn + dd*cs;
132
133 temp = Half * (a+d);
134 a = temp;
135 d = temp;
136
137 if (c!=Zero) {
138 if (b!=Zero) {
139 if (sign(One, b)==sign(One,c)) {
140 //
141 // Real eigenvalues: reduce to upper triangular form
142 //
143 const T sab = sqrt(abs(b));
144 const T sac = sqrt(abs(c));
145 const T p = sign(sab*sac, c);
146 const T tau = One / sqrt(abs(b+c));
147 a = temp + p;
148 d = temp - p;
149 b = b -c;
150 c = Zero;
151 const T cs1 = sab*tau;
152 const T sn1 = sac*tau;
153 temp = cs*cs1 - sn*sn1;
154 sn = cs*sn1 + sn*cs1;
155 cs = temp;
156 }
157 } else {
158 b = -c;
159 c = Zero;
160 temp = cs;
161 cs = -sn;
162 sn = temp;
163 }
164 }
165 }
166 }
167 //
168 // Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
169 //
170 rt1r = a;
171 rt2r = d;
172 if (c==Zero) {
173 rt1i = Zero;
174 rt2i = Zero;
175 } else {
176 rt1i = sqrt(abs(b))*sqrt(abs(c));
177 rt2i = -rt1i;
178 }
179 }
180
181 //== interface for native lapack ===============================================
182
183 #ifdef CHECK_CXXLAPACK
184
185 template <typename T>
186 void
187 lanv2_native(T &a, T &b, T &c, T &d,
188 T &rt1r, T &rt1i,
189 T &rt2r, T &rt2i,
190 T &cs, T &sn)
191 {
192 if (IsSame<T,DOUBLE>::value) {
193 LAPACK_IMPL(dlanv2)(&a, &b, &c, &d,
194 &rt1r, &rt1i,
195 &rt2r, &rt2i,
196 &cs, &sn);
197 } else {
198 ASSERT(0);
199 }
200 }
201
202 #endif // CHECK_CXXLAPACK
203
204 //== public interface ==========================================================
205
206 template <typename T>
207 void
208 lanv2(T &a, T &b, T &c, T &d, T &rt1r, T &rt1i, T &rt2r, T &rt2i, T &cs, T &sn)
209 {
210 LAPACK_DEBUG_OUT("lanv2");
211
212 # ifdef CHECK_CXXLAPACK
213 //
214 // Make copies of output arguments
215 //
216 T _a = a;
217 T _b = b;
218 T _c = c;
219 T _d = d;
220 T _rt1r = rt1r;
221 T _rt1i = rt1i;
222 T _rt2r = rt2r;
223 T _rt2i = rt2i;
224 T _cs = cs;
225 T _sn = sn;
226 # endif
227
228 //
229 // Call implementation
230 //
231 lanv2_generic(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn);
232
233 # ifdef CHECK_CXXLAPACK
234 //
235 // Compare results
236 //
237 lanv2_native(_a, _b, _c, _d, _rt1r, _rt1i, _rt2r, _rt2i, _cs, _sn);
238
239 bool failed = false;
240 if (! isIdentical(a, _a, " a", "_a")) {
241 std::cerr << "CXXLAPACK: a = " << a << std::endl;
242 std::cerr << "F77LAPACK: _a = " << _a << std::endl;
243 failed = true;
244 }
245
246 if (! isIdentical(b, _b, " b", "_b")) {
247 std::cerr << "CXXLAPACK: b = " << b << std::endl;
248 std::cerr << "F77LAPACK: _b = " << _b << std::endl;
249 failed = true;
250 }
251
252 if (! isIdentical(c, _c, " c", "_c")) {
253 std::cerr << "CXXLAPACK: c = " << c << std::endl;
254 std::cerr << "F77LAPACK: _c = " << _c << std::endl;
255 failed = true;
256 }
257
258 if (! isIdentical(d, _d, " d", "_d")) {
259 std::cerr << "CXXLAPACK: d = " << d << std::endl;
260 std::cerr << "F77LAPACK: _d = " << _d << std::endl;
261 failed = true;
262 }
263
264 if (! isIdentical(rt1r, _rt1r, " rt1r", "_rt1r")) {
265 std::cerr << "CXXLAPACK: rt1r = " << rt1r << std::endl;
266 std::cerr << "F77LAPACK: _rt1r = " << _rt1r << std::endl;
267 failed = true;
268 }
269
270 if (! isIdentical(rt1i, _rt1i, " rt1i", "_rt1i")) {
271 std::cerr << "CXXLAPACK: rt1i = " << rt1i << std::endl;
272 std::cerr << "F77LAPACK: _rt1i = " << _rt1i << std::endl;
273 failed = true;
274 }
275
276 if (! isIdentical(rt2r, _rt2r, " rt2r", "_rt2r")) {
277 std::cerr << "CXXLAPACK: rt2r = " << rt2r << std::endl;
278 std::cerr << "F77LAPACK: _rt2r = " << _rt2r << std::endl;
279 failed = true;
280 }
281
282 if (! isIdentical(rt2i, _rt2i, " rt2i", "_rt2i")) {
283 std::cerr << "CXXLAPACK: rt2i = " << rt2i << std::endl;
284 std::cerr << "F77LAPACK: _rt2i = " << _rt2i << std::endl;
285 failed = true;
286 }
287
288 if (! isIdentical(cs, _cs, " cs", "_cs")) {
289 std::cerr << "CXXLAPACK: cs = " << cs << std::endl;
290 std::cerr << "F77LAPACK: _cs = " << _cs << std::endl;
291 failed = true;
292 }
293
294 if (! isIdentical(sn, _sn, " sn", "_sn")) {
295 std::cerr << "CXXLAPACK: sn = " << sn << std::endl;
296 std::cerr << "F77LAPACK: _sn = " << _sn << std::endl;
297 failed = true;
298 }
299
300 if (failed) {
301 ASSERT(0);
302 }
303 # endif
304 }
305
306 } } // namespace lapack, flens
307
308 #endif // FLENS_LAPACK_EIG_LANV2_TCC