1 #define BOOST_TEST_MODULE TestQR
2 #include <boost/test/unit_test.hpp>
3
4 #include <vector>
5 #include <random>
6 #include <boost/multi_array.hpp>
7
8 #include <amgcl/detail/qr.hpp>
9 #include <amgcl/value_type/interface.hpp>
10 #include <amgcl/value_type/complex.hpp>
11 #include <amgcl/value_type/static_matrix.hpp>
12
13 template <class T>
14 struct make_random {
getmake_random15 static T get() {
16 static std::mt19937 gen;
17 static std::uniform_real_distribution<T> rnd;
18
19 return rnd(gen);
20 }
21 };
22
23 template <class T>
random()24 T random() {
25 return make_random<T>::get();
26 }
27
28 template <class T>
29 struct make_random< std::complex<T> > {
getmake_random30 static std::complex<T> get() {
31 return std::complex<T>( random<T>(), random<T>() );
32 }
33 };
34
35 template <class T, int N, int M>
36 struct make_random< amgcl::static_matrix<T,N,M> > {
37 typedef amgcl::static_matrix<T,N,M> matrix;
getmake_random38 static matrix get() {
39 matrix A = amgcl::math::zero<matrix>();
40 for(int i = 0; i < N; ++i)
41 for(int j = 0; j < M; ++j)
42 A(i,j) = make_random<T>::get();
43 return A;
44 }
45 };
46
47 template <class value_type, amgcl::detail::storage_order order>
qr_factorize(int n,int m)48 void qr_factorize(int n, int m) {
49 std::cout << "factorize " << n << " " << m << std::endl;
50 typedef typename std::conditional<order == amgcl::detail::row_major,
51 boost::c_storage_order,
52 boost::fortran_storage_order
53 >::type ma_storage_order;
54
55 boost::multi_array<value_type, 2> A0(boost::extents[n][m], ma_storage_order());
56
57 for(int i = 0; i < n; ++i)
58 for(int j = 0; j < m; ++j)
59 A0[i][j] = random<value_type>();
60
61 boost::multi_array<value_type, 2> A = A0;
62
63 amgcl::detail::QR<value_type> qr;
64
65 qr.factorize(n, m, A.data(), order);
66
67 // Check that A = QR
68 int p = std::min(n, m);
69 for(int i = 0; i < n; ++i) {
70 for(int j = 0; j < m; ++j) {
71 value_type sum = amgcl::math::zero<value_type>();
72
73 for(int k = 0; k < p; ++k)
74 sum += qr.Q(i,k) * qr.R(k,j);
75
76 sum -= A0[i][j];
77
78 BOOST_CHECK_SMALL(amgcl::math::norm(sum), 1e-8);
79 }
80 }
81 }
82
83 template <class value_type, amgcl::detail::storage_order order>
qr_solve(int n,int m)84 void qr_solve(int n, int m) {
85 std::cout << "solve " << n << " " << m << std::endl;
86 typedef typename std::conditional<order == amgcl::detail::row_major,
87 boost::c_storage_order,
88 boost::fortran_storage_order
89 >::type ma_storage_order;
90
91 typedef typename amgcl::math::rhs_of<value_type>::type rhs_type;
92
93 boost::multi_array<value_type, 2> A0(boost::extents[n][m], ma_storage_order());
94
95 for(int i = 0; i < n; ++i)
96 for(int j = 0; j < m; ++j)
97 A0[i][j] = random<value_type>();
98
99 boost::multi_array<value_type, 2> A = A0;
100
101 amgcl::detail::QR<value_type> qr;
102
103 std::vector<rhs_type> f0(n, amgcl::math::constant<rhs_type>(1));
104 std::vector<rhs_type> f = f0;
105
106 std::vector<rhs_type> x(m);
107
108 qr.solve(n, m, A.data(), f.data(), x.data(), order);
109
110 std::vector<rhs_type> Ax(n);
111 for(int i = 0; i < n; ++i) {
112 rhs_type sum = amgcl::math::zero<rhs_type>();
113 for(int j = 0; j < m; ++j)
114 sum += A0[i][j] * x[j];
115
116 Ax[i] = sum;
117
118 if (n < m) {
119 BOOST_CHECK_SMALL(amgcl::math::norm(sum - f0[i]), 1e-8);
120 }
121 }
122
123 if (n >= m) {
124 for(int i = 0; i < m; ++i) {
125 rhs_type sumx = amgcl::math::zero<rhs_type>();
126 rhs_type sumf = amgcl::math::zero<rhs_type>();
127
128 for(int j = 0; j < n; ++j) {
129 sumx += amgcl::math::adjoint(A0[j][i]) * Ax[j];
130 sumf += amgcl::math::adjoint(A0[j][i]) * f0[j];
131 }
132
133 rhs_type delta = sumx - sumf;
134
135 BOOST_CHECK_SMALL(amgcl::math::norm(delta), 1e-8);
136 }
137 }
138 }
139
140 BOOST_AUTO_TEST_SUITE( test_qr )
141
BOOST_AUTO_TEST_CASE(test_qr_factorize)142 BOOST_AUTO_TEST_CASE( test_qr_factorize ) {
143 const int shape[][2] = {
144 {3, 3},
145 {3, 5},
146 {5, 3},
147 {5, 5}
148 };
149
150 const int n = sizeof(shape) / sizeof(shape[0]);
151
152 for(int i = 0; i < n; ++i) {
153 qr_factorize<double, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
154 qr_factorize<double, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
155 qr_factorize<std::complex<double>, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
156 qr_factorize<std::complex<double>, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
157 qr_factorize<amgcl::static_matrix<double, 2, 2>, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
158 qr_factorize<amgcl::static_matrix<double, 2, 2>, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
159 }
160 }
161
BOOST_AUTO_TEST_CASE(test_qr_solve)162 BOOST_AUTO_TEST_CASE( test_qr_solve ) {
163 const int shape[][2] = {
164 {3, 3},
165 {3, 5},
166 {5, 3},
167 {5, 5}
168 };
169
170 const int n = sizeof(shape) / sizeof(shape[0]);
171
172 for(int i = 0; i < n; ++i) {
173 qr_solve<double, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
174 qr_solve<double, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
175 qr_solve<std::complex<double>, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
176 qr_solve<std::complex<double>, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
177 qr_solve<amgcl::static_matrix<double, 2, 2>, amgcl::detail::row_major>(shape[i][0], shape[i][1]);
178 qr_solve<amgcl::static_matrix<double, 2, 2>, amgcl::detail::col_major>(shape[i][0], shape[i][1]);
179 }
180 }
181
BOOST_AUTO_TEST_CASE(qr_issue_39)182 BOOST_AUTO_TEST_CASE( qr_issue_39 ) {
183 boost::multi_array<double, 2> A0(boost::extents[2][2]);
184 A0[0][0] = 1e+0;
185 A0[0][1] = 1e+0;
186 A0[1][0] = 1e-8;
187 A0[1][1] = 1e+0;
188
189 boost::multi_array<double, 2> A = A0;
190
191 amgcl::detail::QR<double> qr;
192
193 qr.factorize(2, 2, A.data());
194
195 // Check that A = QR
196 for(int i = 0; i < 2; ++i) {
197 for(int j = 0; j < 2; ++j) {
198 double sum = 0;
199 for(int k = 0; k < 2; ++k)
200 sum += qr.Q(i,k) * qr.R(k,j);
201
202 sum -= A0[i][j];
203
204 BOOST_CHECK_SMALL(sum, 1e-8);
205 }
206 }
207 }
208
209 BOOST_AUTO_TEST_SUITE_END()
210