MyraMath
csy_rotate2


Source: tests/dense/rotate2.cpp

1 // ========================================================================= //
2 // This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton //
3 // and distributed by MyraCore, LLC. See LICENSE.txt for license terms. //
4 // ========================================================================= //
5 
11 // Containers.
12 #include <myramath/utility/Number.h>
13 #include <myramath/dense/Matrix22.h>
14 #include <myramath/dense/MatrixRange.h>
15 
16 // Algorithms.
17 #include <myramath/utility/random.h>
18 #include <myramath/dense/hermitian.h>
19 #include <myramath/dense/transpose.h>
20 #include <myramath/dense/gemm.h>
21 #include <myramath/dense/rotate.h>
22 
23 #include <myramath/io/FileStreams.h>
24 
25 // Reporting.
26 #include <tests/myratest.h>
27 
28 using namespace myra;
29 
30 namespace {
31 
32 // Given random symmetric A, compute an sy_rotate2() rotation to diagonalize it.
33 template<class Number> void test_sy(const Matrix22<Number>& A, typename ReflectPrecision<Number>::type tolerance)
34  {
35  typedef typename ReflectPrecision<Number>::type Precision;
36  myra::out() << typestring<Number>() << std::endl;
37  auto R = sy_rotate2(A);
38  auto B = R*A*transpose(R);
39  auto I = Matrix22<Number>::identity();
40  Precision B_error = scalar_norm1( B(1,0) );
41  Precision R_error = frobenius(hermitian(R)*R-I);
42  myra::out() << "|(R*A*Rt)(1,0)| = " << B_error << std::endl;
43  myra::out() << "|Rh*R-I| = " << R_error << std::endl;
44  REQUIRE(B_error < tolerance);
45  REQUIRE(R_error < tolerance);
46  }
47 
48 // Given random hermitian A, compute an he_rotate2() rotation to diagonalize it.
49 template<class Precision> void test_he(const Matrix22<typename ReflectComplex<Precision>::type>& A, Precision tolerance)
50  {
51  typedef typename ReflectComplex<Precision>::type Number;
52  myra::out() << typestring<Number>() << std::endl;
53  auto R = he_rotate2(A);
54  auto B = R*A*hermitian(R);
55  auto I = Matrix22<Number>::identity();
56  Precision B_error = scalar_norm1( B(1,0) );
57  Precision R_error = frobenius(hermitian(R)*R-I);
58  myra::out() << "|(R*A*Rh)(1,0)| = " << B_error << std::endl;
59  myra::out() << "|Rh*R-I| = " << R_error << std::endl;
60  REQUIRE(B_error < tolerance);
61  REQUIRE(R_error < tolerance);
62  }
63 
64 // Given random A, compute its gesvd2()
65 template<class Number> void test_gesvd2(const Matrix22<Number>& A, typename ReflectPrecision<Number>::type tolerance)
66  {
67  typedef typename ReflectPrecision<Number>::type Precision;
68  auto USV = gesvd2(A);
69  const Matrix22<Number>& U = USV.U;
70  const Matrix22<Precision>& S = USV.S;
71  const Matrix22<Number>& V = USV.V;
72  auto I = Matrix22<Number>::identity();
73  Precision A_error = frobenius(A-U*S*V);
74  Precision U_error = frobenius(hermitian(U)*U-I);
75  Precision V_error = frobenius(V*hermitian(V)-I);
76  myra::out() << "|A-U*S*V| = " << A_error << std::endl;
77  myra::out() << "|U'*U-I| = " << U_error << std::endl;
78  myra::out() << "|V*V'-I| = " << V_error << std::endl;
79  REQUIRE(A_error < tolerance);
80  REQUIRE(U_error < tolerance);
81  REQUIRE(V_error < tolerance);
82  }
83 
84 // Returns a random symmetric Matrix22 that is singular.
85 template<class Number> Matrix22<Number> sy_singular()
86  {
87  auto A = Matrix22<Number>::random();
88  auto D = Matrix22<Number>::fill_cmajor(1,0,0,0);
89  return A*D*transpose(A);
90  }
91 
92 // Returns a random hermitian Matrix22 that is singular.
93 template<class Number> Matrix22<Number> he_singular()
94  {
95  auto A = Matrix22<Number>::random();
96  auto D = Matrix22<Number>::fill_cmajor(1,0,0,0);
97  return A*D*hermitian(A);
98  }
99 
100 } // namespace
101 
102 // Symmetric rotations.
103 ADD_TEST("ssy_rotate2","[dense][lapack]")
104  {
105  for (int t = 0; t < 10; ++t)
106  test_sy(Matrix22<NumberS>::random_symmetric(), 1.0e-5f);
107  }
108 
109 ADD_TEST("dsy_rotate2","[dense][lapack]")
110  {
111  for (int t = 0; t < 10; ++t)
112  test_sy(Matrix22<NumberD>::random_symmetric(), 1.0e-10);
113  }
114 ADD_TEST("csy_rotate2","[dense][lapack]")
115  {
116  for (int t = 0; t < 10; ++t)
117  test_sy(Matrix22<NumberC>::random_symmetric(), 1.0e-5f);
118  }
119 ADD_TEST("zsy_rotate2","[dense][lapack]")
120  {
121  for (int t = 0; t < 10; ++t)
122  test_sy(Matrix22<NumberZ>::random_symmetric(), 1.0e-10);
123  }
124 
125 // Special cases.
126 ADD_TEST("sy_rotate2_special","[dense][lapack]")
127  {
128  // Identity matrix.
129  test_sy(Matrix22<NumberS>::identity(), 1.0e-5f);
130  test_sy(Matrix22<NumberD>::identity(), 1.0e-10);
131  test_sy(Matrix22<NumberC>::identity(), 1.0e-5f);
132  test_sy(Matrix22<NumberZ>::identity(), 1.0e-10);
133  // Zero matrix.
134  test_sy(Matrix22<NumberS>::zeros(), 1.0e-5f);
135  test_sy(Matrix22<NumberD>::zeros(), 1.0e-10);
136  test_sy(Matrix22<NumberC>::zeros(), 1.0e-5f);
137  test_sy(Matrix22<NumberZ>::zeros(), 1.0e-10);
138  // One matrix.
139  test_sy(Matrix22<NumberS>::ones(), 1.0e-5f);
140  test_sy(Matrix22<NumberD>::ones(), 1.0e-10);
141  test_sy(Matrix22<NumberC>::ones(), 1.0e-5f);
142  test_sy(Matrix22<NumberZ>::ones(), 1.0e-10);
143  // Skew matrix, [0 1; 1 0]
144  test_sy(Matrix22<NumberS>::fill_rmajor(0,1,1,0), 1.0e-5f);
145  test_sy(Matrix22<NumberD>::fill_rmajor(0,1,1,0), 1.0e-10);
146  test_sy(Matrix22<NumberC>::fill_rmajor(0,1,1,0), 1.0e-5f);
147  test_sy(Matrix22<NumberZ>::fill_rmajor(0,1,1,0), 1.0e-10);
148  // Random singular matrix.
149  for (int t = 0; t < 10; ++t)
150  test_sy(sy_singular<NumberS>(), 1.0e-5f);
151  for (int t = 0; t < 10; ++t)
152  test_sy(sy_singular<NumberD>(), 1.0e-10);
153  for (int t = 0; t < 10; ++t)
154  test_sy(sy_singular<NumberC>(), 1.0e-5f);
155  for (int t = 0; t < 10; ++t)
156  test_sy(sy_singular<NumberZ>(), 1.0e-10);
157  }
158 
159 // Hermitian rotations.
160 ADD_TEST("che_rotate2","[dense][lapack]")
161  {
162  for (int t = 0; t < 1; ++t)
163  test_he<float>(Matrix22<NumberC>::random_hermitian(), 1.0e-5f);
164  }
165 ADD_TEST("zhe_rotate2","[dense][lapack]")
166  {
167  for (int t = 0; t < 1; ++t)
168  test_he<double>(Matrix22<NumberZ>::random_hermitian(), 1.0e-10);
169  }
170 
171 // Special cases.
172 ADD_TEST("he_rotate2_special","[dense][lapack]")
173  {
174  // Identity matrix.
175  test_he<float>(Matrix22<NumberC>::identity(), 1.0e-5f);
176  test_he<double>(Matrix22<NumberZ>::identity(), 1.0e-10);
177  // Zero matrix.
178  test_he<float>(Matrix22<NumberC>::zeros(), 1.0e-5f);
179  test_he<double>(Matrix22<NumberZ>::zeros(), 1.0e-10);
180  // One matrix.
181  test_he<float>(Matrix22<NumberC>::ones(), 1.0e-5f);
182  test_he<double>(Matrix22<NumberZ>::ones(), 1.0e-10);
183  // Skew matrix, [0 1; 1 0]
184  test_he<float>(Matrix22<NumberC>::fill_rmajor(0,1,1,0), 1.0e-5f);
185  test_he<double>(Matrix22<NumberZ>::fill_rmajor(0,1,1,0), 1.0e-10);
186  // Random singular matrix.
187  for (int t = 0; t < 10; ++t)
188  test_he<float>(he_singular<NumberC>(), 1.0e-5f);
189  for (int t = 0; t < 10; ++t)
190  test_he<double>(he_singular<NumberZ>(), 1.0e-10);
191  }
192 
193 // Singular value decompositions.
194 ADD_TEST("gesvd2","[dense][lapack]")
195  {
196  test_gesvd2(Matrix22<NumberS>::random(), 1.0e-5f);
197  test_gesvd2(Matrix22<NumberD>::random(), 1.0e-10);
198  test_gesvd2(Matrix22<NumberC>::random(), 1.0e-5f);
199  test_gesvd2(Matrix22<NumberZ>::random(), 1.0e-10);
200  }
201 
202 
203 // Special cases.
204 ADD_TEST("gesvd2_special","[dense][lapack]")
205  {
206  // Identity matrix.
207  test_gesvd2(Matrix22<NumberS>::identity(), 1.0e-5f);
208  test_gesvd2(Matrix22<NumberD>::identity(), 1.0e-10);
209  test_gesvd2(Matrix22<NumberC>::identity(), 1.0e-5f);
210  test_gesvd2(Matrix22<NumberZ>::identity(), 1.0e-10);
211  // Zero matrix.
212  test_gesvd2(Matrix22<NumberS>::zeros(), 1.0e-5f);
213  test_gesvd2(Matrix22<NumberD>::zeros(), 1.0e-10);
214  test_gesvd2(Matrix22<NumberC>::zeros(), 1.0e-5f);
215  test_gesvd2(Matrix22<NumberZ>::zeros(), 1.0e-10);
216  // One matrix.
217  test_gesvd2(Matrix22<NumberS>::ones(), 1.0e-5f);
218  test_gesvd2(Matrix22<NumberD>::ones(), 1.0e-10);
219  test_gesvd2(Matrix22<NumberC>::ones(), 1.0e-5f);
220  test_gesvd2(Matrix22<NumberZ>::ones(), 1.0e-10);
221  // Skew matrix, [0 1; 1 0]
222  test_gesvd2(Matrix22<NumberS>::fill_rmajor(0,1,1,0), 1.0e-5f);
223  test_gesvd2(Matrix22<NumberD>::fill_rmajor(0,1,1,0), 1.0e-10);
224  test_gesvd2(Matrix22<NumberC>::fill_rmajor(0,1,1,0), 1.0e-5f);
225  test_gesvd2(Matrix22<NumberZ>::fill_rmajor(0,1,1,0), 1.0e-10);
226  // Singular matrix, [1 0; 0 0]
227  test_gesvd2(Matrix22<NumberS>::fill_rmajor(1,0,0,0), 1.0e-5f);
228  test_gesvd2(Matrix22<NumberD>::fill_rmajor(1,0,0,0), 1.0e-10);
229  test_gesvd2(Matrix22<NumberC>::fill_rmajor(1,0,0,0), 1.0e-5f);
230  test_gesvd2(Matrix22<NumberZ>::fill_rmajor(1,0,0,0), 1.0e-10);
231  // Singular matrix, [0 1; 0 0]
232  test_gesvd2(Matrix22<NumberS>::fill_rmajor(0,1,0,0), 1.0e-5f);
233  test_gesvd2(Matrix22<NumberD>::fill_rmajor(0,1,0,0), 1.0e-10);
234  test_gesvd2(Matrix22<NumberC>::fill_rmajor(0,1,0,0), 1.0e-5f);
235  test_gesvd2(Matrix22<NumberZ>::fill_rmajor(0,1,0,0), 1.0e-10);
236  // Singular matrix, [0 0; 1 0]
237  test_gesvd2(Matrix22<NumberS>::fill_rmajor(0,0,1,0), 1.0e-5f);
238  test_gesvd2(Matrix22<NumberD>::fill_rmajor(0,0,1,0), 1.0e-10);
239  test_gesvd2(Matrix22<NumberC>::fill_rmajor(0,0,1,0), 1.0e-5f);
240  test_gesvd2(Matrix22<NumberZ>::fill_rmajor(0,0,1,0), 1.0e-10);
241  // Singular matrix, [0 0; 0 1]
242  test_gesvd2(Matrix22<NumberS>::fill_rmajor(0,0,0,1), 1.0e-5f);
243  test_gesvd2(Matrix22<NumberD>::fill_rmajor(0,0,0,1), 1.0e-10);
244  test_gesvd2(Matrix22<NumberC>::fill_rmajor(0,0,0,1), 1.0e-5f);
245  test_gesvd2(Matrix22<NumberZ>::fill_rmajor(0,0,0,1), 1.0e-10);
246  }
247 
myra::Matrix22::random
static Matrix22< Number > random()
Generates a random Matrix22.
Definition: Matrix22.cpp:150
MatrixRange.h
Interface class for representing subranges of dense Matrix&#39;s.
myra::Matrix22::zeros
static Matrix22< Number > zeros()
Generates a zeros Matrix of specified size.
Definition: Matrix22.cpp:175
myra
Definition: syntax.dox:1
Matrix22.h
Matrix type with fixed size of 2x2, algorithms that operate upon them.
transpose.h
Returns a transposed copy of a Matrix. The inplace version only works on a square operand...
myra::ReflectComplex
Reflects a Precision type into a complex type.
Definition: Number.h:46
myra::Matrix22::identity
static Matrix22< Number > identity()
Generates an identity Matrix22.
Definition: Matrix22.cpp:137
Number.h
Various utility functions/classes related to scalar Number types.
myra::Matrix22::fill_cmajor
static Matrix22< Number > fill_cmajor(Number a00, Number a10, Number a01, Number a11)
Fills from four aij values in column major order, A = [a00 a01; a10 a11].
Definition: Matrix22.cpp:205
myra::Matrix22::ones
static Matrix22< Number > ones()
Generates a ones Matrix of specified size.
Definition: Matrix22.cpp:179
myra::Matrix22
Matrix type with fixed size 2x2.
Definition: Matrix22.h:25
myra::Matrix22::fill_rmajor
static Matrix22< Number > fill_rmajor(Number a00, Number a01, Number a10, Number a11)
Fills from four aij values in row major order, A = [a00 a01; a10 a11].
Definition: Matrix22.cpp:194
myra::ReflectPrecision
Reflects Precision trait for a Number, scalar Number types should specialize it.
Definition: Number.h:33
hermitian.h
Returns a hermitian copy of a Matrix. The inplace version only works on a square operand.
random.h
Simplistic random number functions.
FileStreams.h
Streams that serialize/deserialize to/from files.
gemm.h
Variety of routines all for dense Matrix*Matrix multiplies. Delegates to the BLAS.
rotate.h
Computes a plane rotation, either single sided (rotate1) or double sided (rotate2).


Results: [PASS]

std::complex<float>
|(R*A*Rt)(1,0)| = 3.6365e-07
|Rh*R-I| = 2.07549e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 5.96046e-08
|Rh*R-I| = 1.49384e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 2.45756e-07
|Rh*R-I| = 2.10734e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 1.8128e-07
|Rh*R-I| = 1.68587e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 6.664e-08
|Rh*R-I| = 2.80364e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 8.42937e-08
|Rh*R-I| = 1.88486e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 4.47035e-08
|Rh*R-I| = 1.76627e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 2.68635e-08
|Rh*R-I| = 1.26879e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 3.332e-08
|Rh*R-I| = 1.21058e-07
std::complex<float>
|(R*A*Rt)(1,0)| = 1.63404e-07
|Rh*R-I| = 1.91119e-07


Go back to Summary of /test programs.

Generated on Sun Aug 4 2024 12:26:08 for MyraMath by   doxygen 1.8.13