1 /* =========================================================================
2 Copyright (c) 2010-2016, Institute for Microelectronics,
3 Institute for Analysis and Scientific Computing,
4 TU Wien.
5 Portions of this software are copyright by UChicago Argonne, LLC.
6
7 -----------------
8 ViennaCL - The Vienna Computing Library
9 -----------------
10
11 Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12
13 (A list of authors and contributors can be found in the PDF manual)
14
15 License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17
18 #ifndef TEST_MATRIX_PRODUCT_FLOAT_DOUBLE_HPP_
19 #define TEST_MATRIX_PRODUCT_FLOAT_DOUBLE_HPP_
20
21 // We don't need debug mode in UBLAS:
22 #define BOOST_UBLAS_NDEBUG
23
24 #include <cstddef>
25
26 #include "viennacl/matrix.hpp"
27 #include "viennacl/matrix_proxy.hpp"
28 #include "viennacl/linalg/prod.hpp"
29
30 #include "boost/numeric/ublas/matrix.hpp"
31 #include "boost/numeric/ublas/matrix_proxy.hpp"
32 #include "boost/numeric/ublas/io.hpp"
33
34 #include "viennacl/tools/random.hpp"
35
36 template<typename ScalarType, typename VCLMatrixType>
37 ScalarType diff(boost::numeric::ublas::matrix<ScalarType> const & mat1, VCLMatrixType const & mat2)
38 {
39 boost::numeric::ublas::matrix<ScalarType> mat2_cpu(mat2.size1(), mat2.size2());
40 viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
41 viennacl::copy(mat2, mat2_cpu);
42 ScalarType ret = 0;
43 ScalarType act = 0;
44
45 for (unsigned int i = 0; i < mat2_cpu.size1(); ++i)
46 {
47 for (unsigned int j = 0; j < mat2_cpu.size2(); ++j)
48 {
49 act = std::fabs(mat2_cpu(i,j) - mat1(i,j)) / std::max( std::fabs(mat2_cpu(i, j)), std::fabs(mat1(i,j)) );
50 if (act > ret)
51 ret = act;
52 }
53 }
54
55 return ret;
56 }
57
58
59 template<class UBlasType, class F>
60 struct matrix_maker;
61
62 template<class T, class F>
63 struct matrix_maker< boost::numeric::ublas::matrix<T>, F>
64 {
65 typedef viennacl::matrix<T, F> result_type;
makematrix_maker66 static result_type make(viennacl::matrix<T, F> const &, boost::numeric::ublas::matrix<T> & base)
67 {
68 viennacl::matrix<T, F> result(base.size1(), base.size2());
69 viennacl::copy(base, result);
70 return result;
71 }
72 };
73
74 template<class MatrixT, class F>
75 struct matrix_maker< boost::numeric::ublas::matrix_range<MatrixT>, F>
76 {
77 typedef typename MatrixT::value_type T;
78 typedef viennacl::matrix_range< viennacl::matrix<T, F> > result_type;
79
makematrix_maker80 static result_type make(viennacl::matrix<T, F> & M, boost::numeric::ublas::matrix_range<MatrixT> & base)
81 {
82 viennacl::range r0(base.start1(), base.start1() + base.size1());
83 viennacl::range r1(base.start2(), base.start2() + base.size2());
84 result_type result(M, r0, r1);
85 viennacl::copy(base, result);
86 return result;
87 }
88 };
89
90 template<class MatrixT, class F>
91 struct matrix_maker< boost::numeric::ublas::matrix_slice<MatrixT>, F>
92 {
93 typedef typename MatrixT::value_type T;
94 typedef viennacl::matrix_slice< viennacl::matrix<T, F> > result_type;
95
makematrix_maker96 static result_type make(viennacl::matrix<T, F> & M, boost::numeric::ublas::matrix_slice<MatrixT> & base)
97 {
98 viennacl::slice s0(base.start1(), std::size_t(base.stride1()), base.size1());
99 viennacl::slice s1(base.start2(), std::size_t(base.stride2()), base.size2());
100 result_type result(M, s0, s1);
101 viennacl::copy(base, result);
102 return result;
103 }
104 };
105
106 template<typename T, typename CType, typename AType, typename BType>
test_layout(CType & C,AType const & A,AType const & AT,BType const & B,BType const & BT,boost::numeric::ublas::matrix<T> const & ground,T epsilon,bool with_composite)107 int test_layout(CType & C, AType const & A, AType const & AT, BType const & B, BType const & BT,
108 boost::numeric::ublas::matrix<T> const & ground, T epsilon, bool with_composite)
109 {
110 using viennacl::linalg::prod;
111 using viennacl::trans;
112
113 std::cout << "C = A.B" << std::endl;
114 C = prod(A, B);
115 if (diff(ground, C)>epsilon)
116 return EXIT_FAILURE;
117
118 std::cout << "C = A'.B" << std::endl;
119 C = prod(trans(AT), B);
120 if (diff(ground, C)>epsilon)
121 return EXIT_FAILURE;
122
123 std::cout << "C = A.B'" << std::endl;
124 C = prod(A, trans(BT));
125 if (diff(ground, C)>epsilon)
126 return EXIT_FAILURE;
127
128 std::cout << "C = A'.B'" << std::endl;
129 C = prod(trans(AT), trans(BT));
130 if (diff(ground, C)>epsilon)
131 return EXIT_FAILURE;
132
133 // composite operations:
134 if (with_composite)
135 {
136 boost::numeric::ublas::matrix<T> ground2 = T(2) * ground;
137
138 std::cout << "C = (A + A).B" << std::endl;
139 C = prod(A + A, B);
140 if (diff(ground2, C)>epsilon)
141 return EXIT_FAILURE;
142
143 std::cout << "C = trans(AT + AT).B" << std::endl;
144 C = prod(viennacl::trans(AT + AT), B);
145 if (diff(ground2, C)>epsilon)
146 return EXIT_FAILURE;
147
148 std::cout << "C = A.(B + B)" << std::endl;
149 C = prod(A, T(2) * B);
150 if (diff(ground2, C)>epsilon)
151 return EXIT_FAILURE;
152
153 std::cout << "C = A.trans(BT + BT)" << std::endl;
154 C = prod(A, trans(BT + BT));
155 if (diff(ground2, C)>epsilon)
156 return EXIT_FAILURE;
157
158 std::cout << "C = (A + A).(B + B)" << std::endl;
159 C = T(0.25) * prod(A + A, B + B);
160 if (diff(ground, C)>epsilon)
161 return EXIT_FAILURE;
162
163 std::cout << "C = trans(AT + AT).trans(BT + BT)" << std::endl;
164 C += prod(trans(AT + AT), trans(BT + BT));
165 C -= prod(trans(AT + AT), trans(BT + BT));
166 if (diff(ground, C)>epsilon)
167 return EXIT_FAILURE;
168
169 }
170
171 return EXIT_SUCCESS;
172 }
173
174 template<typename T, typename RefAType, typename RefBType, typename RefCType>
test_all_layouts(std::size_t CM,std::size_t CN,RefCType & cC,std::size_t AM,std::size_t AK,RefAType & cA,RefAType & cAT,std::size_t BK,std::size_t BN,RefBType & cB,RefBType & cBT,T epsilon)175 int test_all_layouts(std::size_t CM, std::size_t CN, RefCType & cC,
176 std::size_t AM, std::size_t AK, RefAType & cA, RefAType & cAT,
177 std::size_t BK, std::size_t BN, RefBType & cB, RefBType & cBT,
178 T epsilon)
179 {
180 viennacl::matrix<T, viennacl::row_major> ArowTmp(AM, AK);
181 viennacl::matrix<T, viennacl::row_major> ATrowTmp(AK, AM);
182 viennacl::matrix<T, viennacl::row_major> BrowTmp(BK, BN);
183 viennacl::matrix<T, viennacl::row_major> BTrowTmp(BN, BK);
184 viennacl::matrix<T, viennacl::row_major> CrowTmp(CM, CN);
185
186 viennacl::matrix<T, viennacl::column_major> AcolTmp(AM, AK);
187 viennacl::matrix<T, viennacl::column_major> ATcolTmp(AK, AM);
188 viennacl::matrix<T, viennacl::column_major> BcolTmp(BK, BN);
189 viennacl::matrix<T, viennacl::column_major> BTcolTmp(BN, BK);
190 viennacl::matrix<T, viennacl::column_major> CcolTmp(CM, CN);
191
192
193 typename matrix_maker<RefCType, viennacl::row_major>::result_type Crow = matrix_maker<RefCType, viennacl::row_major>::make(CrowTmp, cC);
194 typename matrix_maker<RefAType, viennacl::row_major>::result_type Arow = matrix_maker<RefAType, viennacl::row_major>::make(ArowTmp, cA);
195 typename matrix_maker<RefAType, viennacl::row_major>::result_type ATrow = matrix_maker<RefAType, viennacl::row_major>::make(ATrowTmp, cAT);
196 typename matrix_maker<RefBType, viennacl::row_major>::result_type Brow = matrix_maker<RefBType, viennacl::row_major>::make(BrowTmp, cB);
197 typename matrix_maker<RefBType, viennacl::row_major>::result_type BTrow = matrix_maker<RefBType, viennacl::row_major>::make(BTrowTmp, cBT);
198
199 typename matrix_maker<RefCType, viennacl::column_major>::result_type Ccol = matrix_maker<RefCType, viennacl::column_major>::make(CcolTmp, cC);
200 typename matrix_maker<RefAType, viennacl::column_major>::result_type Acol = matrix_maker<RefAType, viennacl::column_major>::make(AcolTmp, cA);
201 typename matrix_maker<RefAType, viennacl::column_major>::result_type ATcol = matrix_maker<RefAType, viennacl::column_major>::make(ATcolTmp, cAT);
202 typename matrix_maker<RefBType, viennacl::column_major>::result_type Bcol = matrix_maker<RefBType, viennacl::column_major>::make(BcolTmp, cB);
203 typename matrix_maker<RefBType, viennacl::column_major>::result_type BTcol = matrix_maker<RefBType, viennacl::column_major>::make(BTcolTmp, cBT);
204
205
206 boost::numeric::ublas::matrix<T> ground = boost::numeric::ublas::prod(cA, cB);
207
208 #define TEST_LAYOUT(Clayout, Alayout, Blayout, composite) \
209 std::cout << "> " #Clayout " = " #Alayout "." #Blayout << std::endl; \
210 if (test_layout(C ## Clayout, A ## Alayout, AT ## Alayout, B ## Blayout, BT ## Blayout, ground, epsilon, composite) != EXIT_SUCCESS) \
211 return EXIT_FAILURE; \
212
213 TEST_LAYOUT(row, row, row, true);
214 TEST_LAYOUT(row, row, col, false);
215 TEST_LAYOUT(row, col, row, false);
216 TEST_LAYOUT(row, col, col, false);
217 TEST_LAYOUT(col, row, row, false);
218 TEST_LAYOUT(col, row, col, false);
219 TEST_LAYOUT(col, col, row, false);
220 TEST_LAYOUT(col, col, col, true);
221
222 #undef TEST_LAYOUT
223
224 return EXIT_SUCCESS;
225 }
226
227 template<class MatrixType>
init_rand(MatrixType & A)228 void init_rand(MatrixType & A)
229 {
230 typedef typename MatrixType::value_type T;
231
232 viennacl::tools::uniform_random_numbers<T> randomNumber;
233
234 for (unsigned int i = 0; i < A.size1(); ++i)
235 for (unsigned int j = 0; j < A.size2(); ++j)
236 A(i, j) = static_cast<T>(0.1) * randomNumber();
237 }
238
239 template<typename T>
run_test(T epsilon)240 int run_test(T epsilon)
241 {
242 typedef boost::numeric::ublas::range range_type;
243 typedef boost::numeric::ublas::slice slice_type;
244 typedef boost::numeric::ublas::matrix<T> matrix_type;
245 typedef boost::numeric::ublas::matrix_range<matrix_type> matrix_range_type;
246 typedef boost::numeric::ublas::matrix_slice<matrix_type> matrix_slice_type;
247
248 typedef typename matrix_type::difference_type difference_type;
249
250 std::size_t matrix_holder_M = 143;
251 std::size_t matrix_holder_N = 124;
252 std::size_t matrix_holder_K = 184;
253
254 std::size_t start_M = 14;
255 std::size_t start_N = 20;
256 std::size_t start_K = 73;
257
258 std::size_t range_holder_M = start_M + matrix_holder_M;
259 std::size_t range_holder_N = start_N + matrix_holder_N;
260 std::size_t range_holder_K = start_K + matrix_holder_K;
261
262 range_type range_M(start_M, range_holder_M);
263 range_type range_N(start_N, range_holder_N);
264 range_type range_K(start_K, range_holder_K);
265
266 difference_type stride_M = 9;
267 difference_type stride_N = 13;
268 difference_type stride_K = 4;
269
270 std::size_t slice_holder_M = start_M + std::size_t(stride_M)*matrix_holder_M;
271 std::size_t slice_holder_N = start_N + std::size_t(stride_N)*matrix_holder_N;
272 std::size_t slice_holder_K = start_K + std::size_t(stride_K)*matrix_holder_K;
273
274 slice_type slice_M(start_M, stride_M, matrix_holder_M);
275 slice_type slice_N(start_N, stride_N, matrix_holder_N);
276 slice_type slice_K(start_K, stride_K, matrix_holder_K);
277
278 #define DECLARE(NAME, size1, size2) \
279 matrix_type NAME ## _matrix(matrix_holder_ ## size1, matrix_holder_ ## size2);\
280 init_rand(NAME ## _matrix);\
281 matrix_type NAME ## T_matrix = boost::numeric::ublas::trans(NAME ## _matrix);\
282 \
283 matrix_type NAME ## _range_holder(range_holder_ ## size1, range_holder_ ## size2);\
284 init_rand(NAME ## _range_holder);\
285 matrix_range_type NAME ## _range(NAME ## _range_holder, range_ ## size1, range_ ## size2);\
286 matrix_type NAME ## T_range_holder = boost::numeric::ublas::trans(NAME ## _range_holder);\
287 matrix_range_type NAME ## T_range(NAME ## T_range_holder, range_ ## size2, range_ ## size1);\
288 \
289 matrix_type NAME ## _slice_holder(slice_holder_ ## size1, slice_holder_ ## size2);\
290 init_rand(NAME ## _slice_holder);\
291 matrix_slice_type NAME ## _slice(NAME ## _slice_holder, slice_ ## size1, slice_ ## size2);\
292 matrix_type NAME ## T_slice_holder = boost::numeric::ublas::trans(NAME ## _slice_holder);\
293 matrix_slice_type NAME ## T_slice(NAME ## T_slice_holder, slice_ ## size2, slice_ ## size1);\
294
295 DECLARE(A, M, K);
296 DECLARE(B, K, N);
297 DECLARE(C, M, N);
298 #undef DECLARE
299
300 #define TEST_ALL_LAYOUTS(C_TYPE, A_TYPE, B_TYPE)\
301 std::cout << ">> " #C_TYPE " = " #A_TYPE "." #B_TYPE << std::endl;\
302 if (test_all_layouts<T>(C_TYPE ## _holder_M, C_TYPE ## _holder_N, C_ ## C_TYPE,\
303 A_TYPE ## _holder_M, A_TYPE ## _holder_K, A_ ## A_TYPE, AT_ ## A_TYPE,\
304 B_TYPE ## _holder_K, B_TYPE ## _holder_N, B_ ## B_TYPE, BT_ ## B_TYPE, epsilon) != EXIT_SUCCESS)\
305 return EXIT_FAILURE;\
306
307 // //C=matrix
308 TEST_ALL_LAYOUTS(matrix, matrix, matrix)
309 TEST_ALL_LAYOUTS(matrix, matrix, range)
310 TEST_ALL_LAYOUTS(matrix, matrix, slice)
311
312 TEST_ALL_LAYOUTS(matrix, range, matrix)
313 TEST_ALL_LAYOUTS(matrix, range, range)
314 TEST_ALL_LAYOUTS(matrix, range, slice)
315
316 TEST_ALL_LAYOUTS(matrix, slice, matrix)
317 TEST_ALL_LAYOUTS(matrix, slice, range)
318 TEST_ALL_LAYOUTS(matrix, slice, slice)
319
320 // C = range
321 TEST_ALL_LAYOUTS(range, matrix, matrix)
322 TEST_ALL_LAYOUTS(range, matrix, range)
323 TEST_ALL_LAYOUTS(range, matrix, slice)
324
325 TEST_ALL_LAYOUTS(range, range, matrix)
326 TEST_ALL_LAYOUTS(range, range, range)
327 TEST_ALL_LAYOUTS(range, range, slice)
328
329 TEST_ALL_LAYOUTS(range, slice, matrix)
330 TEST_ALL_LAYOUTS(range, slice, range)
331 TEST_ALL_LAYOUTS(range, slice, slice)
332
333 // C = slice
334 TEST_ALL_LAYOUTS(slice, matrix, matrix)
335 TEST_ALL_LAYOUTS(slice, matrix, range)
336 TEST_ALL_LAYOUTS(slice, matrix, slice)
337
338 TEST_ALL_LAYOUTS(slice, range, matrix)
339 TEST_ALL_LAYOUTS(slice, range, range)
340 TEST_ALL_LAYOUTS(slice, range, slice)
341
342 TEST_ALL_LAYOUTS(slice, slice, matrix)
343 TEST_ALL_LAYOUTS(slice, slice, range)
344 TEST_ALL_LAYOUTS(slice, slice, slice)
345
346 #undef TEST_ALL_LAYOUTS
347
348 return EXIT_SUCCESS;
349 }
350
351 #endif
352