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
19
20 /** \file tests/src/scheduler_matrix.cpp Tests the scheduler for matrix-operations (no matrix-matrix).
21 * \test Tests the scheduler for matrix-operations (no matrix-matrix).
22 **/
23
24 #define VIENNACL_WITH_UBLAS
25 //#define NDEBUG
26 //#define VIENNACL_BUILD_INFO
27
28 #include <utility>
29 #include <iostream>
30 #include <fstream>
31 #include <string>
32 #include <cmath>
33 #include <algorithm>
34 #include <cstdio>
35 #include <ctime>
36
37 #include "viennacl/scalar.hpp"
38 #include "viennacl/matrix.hpp"
39 #include "viennacl/linalg/prod.hpp"
40 /*#include "viennacl/compressed_matrix.hpp"
41 #include "viennacl/linalg/cg.hpp"
42 #include "viennacl/linalg/inner_prod.hpp"
43 #include "viennacl/linalg/ilu.hpp"
44 #include "viennacl/linalg/norm_2.hpp"
45 #include "viennacl/io/matrix_market.hpp"*/
46 #include "viennacl/matrix_proxy.hpp"
47 #include "viennacl/vector_proxy.hpp"
48 #include "boost/numeric/ublas/vector.hpp"
49 #include "boost/numeric/ublas/matrix.hpp"
50 #include "boost/numeric/ublas/matrix_proxy.hpp"
51 #include "boost/numeric/ublas/vector_proxy.hpp"
52 #include "boost/numeric/ublas/io.hpp"
53
54 #include "viennacl/scheduler/execute.hpp"
55
56 using namespace boost::numeric;
57
58 template<typename MatrixType, typename VCLMatrixType>
check_for_equality(MatrixType const & ublas_A,VCLMatrixType const & vcl_A,double epsilon)59 bool check_for_equality(MatrixType const & ublas_A, VCLMatrixType const & vcl_A, double epsilon)
60 {
61 typedef typename MatrixType::value_type value_type;
62
63 boost::numeric::ublas::matrix<value_type> vcl_A_cpu(vcl_A.size1(), vcl_A.size2());
64 viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
65 viennacl::copy(vcl_A, vcl_A_cpu);
66
67 for (std::size_t i=0; i<ublas_A.size1(); ++i)
68 {
69 for (std::size_t j=0; j<ublas_A.size2(); ++j)
70 {
71 if (std::fabs(ublas_A(i,j) - vcl_A_cpu(i,j)) > 0)
72 {
73 if ( (std::fabs(ublas_A(i,j) - vcl_A_cpu(i,j)) / std::max(std::fabs(ublas_A(i,j)), std::fabs(vcl_A_cpu(i,j))) > epsilon) || std::fabs(vcl_A_cpu(i,j) - vcl_A_cpu(i,j)) > 0 )
74 {
75 std::cout << "Error at index (" << i << ", " << j << "): " << ublas_A(i,j) << " vs " << vcl_A_cpu(i,j) << std::endl;
76 std::cout << std::endl << "TEST failed!" << std::endl;
77 return false;
78 }
79 }
80 }
81 }
82
83 std::cout << "PASSED!" << std::endl;
84 return true;
85 }
86
87
88
89
90 template<typename UBLASMatrixType,
91 typename ViennaCLMatrixType1, typename ViennaCLMatrixType2, typename ViennaCLMatrixType3>
run_test(double epsilon,UBLASMatrixType & ublas_A,UBLASMatrixType & ublas_B,UBLASMatrixType & ublas_C,ViennaCLMatrixType1 & vcl_A,ViennaCLMatrixType2 & vcl_B,ViennaCLMatrixType3 vcl_C)92 int run_test(double epsilon,
93 UBLASMatrixType & ublas_A, UBLASMatrixType & ublas_B, UBLASMatrixType & ublas_C,
94 ViennaCLMatrixType1 & vcl_A, ViennaCLMatrixType2 & vcl_B, ViennaCLMatrixType3 vcl_C)
95 {
96
97 typedef typename viennacl::result_of::cpu_value_type<typename ViennaCLMatrixType1::value_type>::type cpu_value_type;
98
99 cpu_value_type alpha = cpu_value_type(3.1415);
100 viennacl::scalar<cpu_value_type> gpu_alpha = alpha;
101
102 cpu_value_type beta = cpu_value_type(2.7182);
103 viennacl::scalar<cpu_value_type> gpu_beta = beta;
104
105
106 //
107 // Initializer:
108 //
109 std::cout << "Checking for zero_matrix initializer..." << std::endl;
110 ublas_A = ublas::zero_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2());
111 vcl_A = viennacl::zero_matrix<cpu_value_type>(vcl_A.size1(), vcl_A.size2());
112 if (!check_for_equality(ublas_A, vcl_A, epsilon))
113 return EXIT_FAILURE;
114
115 std::cout << "Checking for scalar_matrix initializer..." << std::endl;
116 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), alpha);
117 vcl_A = viennacl::scalar_matrix<cpu_value_type>(vcl_A.size1(), vcl_A.size2(), alpha);
118 if (!check_for_equality(ublas_A, vcl_A, epsilon))
119 return EXIT_FAILURE;
120
121 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), gpu_beta);
122 vcl_A = viennacl::scalar_matrix<cpu_value_type>( vcl_A.size1(), vcl_A.size2(), gpu_beta);
123 if (!check_for_equality(ublas_A, vcl_A, epsilon))
124 return EXIT_FAILURE;
125
126 /*std::cout << "Checking for identity initializer..." << std::endl;
127 ublas_A = ublas::identity_matrix<cpu_value_type>(ublas_A.size1());
128 vcl_A = viennacl::identity_matrix<cpu_value_type>(vcl_A.size1());
129 if (!check_for_equality(ublas_A, vcl_A, epsilon))
130 return EXIT_FAILURE;*/
131
132
133 std::cout << std::endl;
134 //std::cout << "//" << std::endl;
135 //std::cout << "////////// Test: Assignments //////////" << std::endl;
136 //std::cout << "//" << std::endl;
137
138 if (!check_for_equality(ublas_B, vcl_B, epsilon))
139 return EXIT_FAILURE;
140
141 std::cout << "Testing matrix assignment... ";
142 //std::cout << ublas_B(0,0) << " vs. " << vcl_B(0,0) << std::endl;
143 ublas_A = ublas_B;
144 vcl_A = vcl_B;
145 if (!check_for_equality(ublas_A, vcl_A, epsilon))
146 return EXIT_FAILURE;
147
148
149
150 //std::cout << std::endl;
151 //std::cout << "//" << std::endl;
152 //std::cout << "////////// Test 1: Copy to GPU //////////" << std::endl;
153 //std::cout << "//" << std::endl;
154
155 ublas_A = ublas_B;
156 viennacl::copy(ublas_B, vcl_A);
157 std::cout << "Testing upper left copy to GPU... ";
158 if (!check_for_equality(ublas_A, vcl_A, epsilon))
159 return EXIT_FAILURE;
160
161
162 ublas_C = ublas_B;
163 viennacl::copy(ublas_B, vcl_C);
164 std::cout << "Testing lower right copy to GPU... ";
165 if (!check_for_equality(ublas_C, vcl_C, epsilon))
166 return EXIT_FAILURE;
167
168
169 //std::cout << std::endl;
170 //std::cout << "//" << std::endl;
171 //std::cout << "////////// Test 2: Copy from GPU //////////" << std::endl;
172 //std::cout << "//" << std::endl;
173
174 std::cout << "Testing upper left copy to A... ";
175 if (!check_for_equality(ublas_A, vcl_A, epsilon))
176 return EXIT_FAILURE;
177
178 std::cout << "Testing lower right copy to C... ";
179 if (!check_for_equality(ublas_C, vcl_C, epsilon))
180 return EXIT_FAILURE;
181
182
183
184 //std::cout << "//" << std::endl;
185 //std::cout << "////////// Test 3: Addition //////////" << std::endl;
186 //std::cout << "//" << std::endl;
187 viennacl::copy(ublas_C, vcl_C);
188
189 std::cout << "Assignment: ";
190 {
191 ublas_C = ublas_B;
192 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), vcl_B); // same as vcl_C = vcl_B;
193 viennacl::scheduler::execute(my_statement);
194
195 if (!check_for_equality(ublas_C, vcl_C, epsilon))
196 return EXIT_FAILURE;
197 }
198
199
200 std::cout << "Inplace add: ";
201 {
202 ublas_C += ublas_C;
203 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_inplace_add(), vcl_C); // same as vcl_C += vcl_C;
204 viennacl::scheduler::execute(my_statement);
205
206 if (!check_for_equality(ublas_C, vcl_C, epsilon))
207 return EXIT_FAILURE;
208 }
209
210 std::cout << "Inplace sub: ";
211 {
212 ublas_C -= ublas_C;
213 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_inplace_sub(), vcl_C); // same as vcl_C -= vcl_C;
214 viennacl::scheduler::execute(my_statement);
215
216 if (!check_for_equality(ublas_C, vcl_C, epsilon))
217 return EXIT_FAILURE;
218 }
219
220
221 std::cout << "Add: ";
222 {
223 ublas_C = ublas_A + ublas_B;
224 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), vcl_A + vcl_B); // same as vcl_C = vcl_A + vcl_B;
225 viennacl::scheduler::execute(my_statement);
226
227 if (!check_for_equality(ublas_C, vcl_C, epsilon))
228 return EXIT_FAILURE;
229 }
230
231 std::cout << "Sub: ";
232 {
233 ublas_C = ublas_A - ublas_B;
234 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), vcl_A - vcl_B); // same as vcl_C = vcl_A - vcl_B;
235 viennacl::scheduler::execute(my_statement);
236
237 if (!check_for_equality(ublas_C, vcl_C, epsilon))
238 return EXIT_FAILURE;
239 }
240
241 std::cout << "Composite assignments: ";
242 {
243 ublas_C += alpha * ublas_A - beta * ublas_B + ublas_A / beta - ublas_B / alpha;
244 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_inplace_add(), alpha * vcl_A - beta * vcl_B + vcl_A / beta - vcl_B / alpha); // same as vcl_C += alpha * vcl_A - beta * vcl_B + vcl_A / beta - vcl_B / alpha;
245 viennacl::scheduler::execute(my_statement);
246
247 if (!check_for_equality(ublas_C, vcl_C, epsilon))
248 return EXIT_FAILURE;
249 }
250
251
252 std::cout << "--- Testing elementwise operations (binary) ---" << std::endl;
253 std::cout << "x = element_prod(x, y)... ";
254 {
255 ublas_C = element_prod(ublas_A, ublas_B);
256 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_A, vcl_B));
257 viennacl::scheduler::execute(my_statement);
258
259 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
260 return EXIT_FAILURE;
261 }
262
263 std::cout << "x = element_prod(x + y, y)... ";
264 {
265 ublas_C = element_prod(ublas_A + ublas_B, ublas_B);
266 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_A + vcl_B, vcl_B));
267 viennacl::scheduler::execute(my_statement);
268
269 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
270 return EXIT_FAILURE;
271 }
272
273 std::cout << "x = element_prod(x, x + y)... ";
274 {
275 ublas_C = element_prod(ublas_A, ublas_A + ublas_B);
276 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_A, vcl_B + vcl_A));
277 viennacl::scheduler::execute(my_statement);
278
279 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
280 return EXIT_FAILURE;
281 }
282
283 std::cout << "x = element_prod(x - y, y + x)... ";
284 {
285 ublas_C = element_prod(ublas_A - ublas_B, ublas_B + ublas_A);
286 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_prod(vcl_A - vcl_B, vcl_B + vcl_A));
287 viennacl::scheduler::execute(my_statement);
288
289 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
290 return EXIT_FAILURE;
291 }
292
293
294
295 std::cout << "x = element_div(x, y)... ";
296 {
297 ublas_C = element_div(ublas_A, ublas_B);
298 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_div(vcl_A, vcl_B));
299 viennacl::scheduler::execute(my_statement);
300
301 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
302 return EXIT_FAILURE;
303 }
304
305 std::cout << "x = element_div(x + y, y)... ";
306 {
307 ublas_C = element_div(ublas_A + ublas_B, ublas_B);
308 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_div(vcl_A + vcl_B, vcl_B));
309 viennacl::scheduler::execute(my_statement);
310
311 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
312 return EXIT_FAILURE;
313 }
314
315 std::cout << "x = element_div(x, x + y)... ";
316 {
317 ublas_C = element_div(ublas_A, ublas_A + ublas_B);
318 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_div(vcl_A, vcl_B + vcl_A));
319 viennacl::scheduler::execute(my_statement);
320
321 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
322 return EXIT_FAILURE;
323 }
324
325 std::cout << "x = element_div(x - y, y + x)... ";
326 {
327 ublas_C = element_div(ublas_A - ublas_B, ublas_B + ublas_A);
328 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_div(vcl_A - vcl_B, vcl_B + vcl_A));
329 viennacl::scheduler::execute(my_statement);
330
331 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
332 return EXIT_FAILURE;
333 }
334
335
336 std::cout << "--- Testing elementwise operations (unary) ---" << std::endl;
337 #define GENERATE_UNARY_OP_TEST(OPNAME) \
338 ublas_A = ublas::scalar_matrix<cpu_value_type>(ublas_A.size1(), ublas_A.size2(), cpu_value_type(0.21)); \
339 ublas_B = cpu_value_type(3.1415) * ublas_A; \
340 viennacl::copy(ublas_A, vcl_A); \
341 viennacl::copy(ublas_B, vcl_B); \
342 { \
343 for (std::size_t i=0; i<ublas_C.size1(); ++i) \
344 for (std::size_t j=0; j<ublas_C.size2(); ++j) \
345 ublas_C(i,j) = static_cast<cpu_value_type>(OPNAME(ublas_A(i,j))); \
346 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_A)); \
347 viennacl::scheduler::execute(my_statement); \
348 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS) \
349 return EXIT_FAILURE; \
350 } \
351 { \
352 for (std::size_t i=0; i<ublas_C.size1(); ++i) \
353 for (std::size_t j=0; j<ublas_C.size2(); ++j) \
354 ublas_C(i,j) = static_cast<cpu_value_type>(OPNAME(ublas_A(i,j) / cpu_value_type(2))); \
355 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_##OPNAME(vcl_A / cpu_value_type(2))); \
356 viennacl::scheduler::execute(my_statement); \
357 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS) \
358 return EXIT_FAILURE; \
359 }
360
361 GENERATE_UNARY_OP_TEST(cos);
362 GENERATE_UNARY_OP_TEST(cosh);
363 GENERATE_UNARY_OP_TEST(exp);
364 GENERATE_UNARY_OP_TEST(floor);
365 GENERATE_UNARY_OP_TEST(fabs);
366 GENERATE_UNARY_OP_TEST(log);
367 GENERATE_UNARY_OP_TEST(log10);
368 GENERATE_UNARY_OP_TEST(sin);
369 GENERATE_UNARY_OP_TEST(sinh);
370 GENERATE_UNARY_OP_TEST(fabs);
371 //GENERATE_UNARY_OP_TEST(abs); //OpenCL allows abs on integers only
372 GENERATE_UNARY_OP_TEST(sqrt);
373 GENERATE_UNARY_OP_TEST(tan);
374 GENERATE_UNARY_OP_TEST(tanh);
375
376 #undef GENERATE_UNARY_OP_TEST
377
378 if (ublas_C.size1() == ublas_C.size2()) // transposition tests
379 {
380 std::cout << "z = element_fabs(x - trans(y))... ";
381 {
382 for (std::size_t i=0; i<ublas_C.size1(); ++i)
383 for (std::size_t j=0; j<ublas_C.size2(); ++j)
384 ublas_C(i,j) = std::fabs(ublas_A(i,j) - ublas_B(j,i));
385 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_fabs((vcl_A) - trans(vcl_B)));
386 viennacl::scheduler::execute(my_statement);
387
388 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
389 return EXIT_FAILURE;
390 }
391
392 std::cout << "z = element_fabs(trans(x) - (y))... ";
393 {
394 for (std::size_t i=0; i<ublas_C.size1(); ++i)
395 for (std::size_t j=0; j<ublas_C.size2(); ++j)
396 ublas_C(i,j) = std::fabs(ublas_A(j,i) - ublas_B(i,j));
397 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_fabs(trans(vcl_A) - (vcl_B)));
398 viennacl::scheduler::execute(my_statement);
399
400 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
401 return EXIT_FAILURE;
402 }
403
404 std::cout << "z = element_fabs(trans(x) - trans(y))... ";
405 {
406 for (std::size_t i=0; i<ublas_C.size1(); ++i)
407 for (std::size_t j=0; j<ublas_C.size2(); ++j)
408 ublas_C(i,j) = std::fabs(ublas_A(j,i) - ublas_B(j,i));
409 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_assign(), viennacl::linalg::element_fabs(trans(vcl_A) - trans(vcl_B)));
410 viennacl::scheduler::execute(my_statement);
411
412 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
413 return EXIT_FAILURE;
414 }
415
416 std::cout << "z += trans(x)... ";
417 {
418 for (std::size_t i=0; i<ublas_C.size1(); ++i)
419 for (std::size_t j=0; j<ublas_C.size2(); ++j)
420 ublas_C(i,j) += ublas_A(j,i);
421 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_inplace_add(), trans(vcl_A));
422 viennacl::scheduler::execute(my_statement);
423
424 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
425 return EXIT_FAILURE;
426 }
427
428 std::cout << "z -= trans(x)... ";
429 {
430 for (std::size_t i=0; i<ublas_C.size1(); ++i)
431 for (std::size_t j=0; j<ublas_C.size2(); ++j)
432 ublas_C(i,j) -= ublas_A(j,i);
433 viennacl::scheduler::statement my_statement(vcl_C, viennacl::op_inplace_sub(), trans(vcl_A));
434 viennacl::scheduler::execute(my_statement);
435
436 if (!check_for_equality(ublas_C, vcl_C, epsilon) != EXIT_SUCCESS)
437 return EXIT_FAILURE;
438 }
439
440 }
441
442 std::cout << std::endl;
443 std::cout << "----------------------------------------------" << std::endl;
444 std::cout << std::endl;
445
446
447 return EXIT_SUCCESS;
448 }
449
450
451
452
453 template<typename T, typename ScalarType>
run_test(double epsilon)454 int run_test(double epsilon)
455 {
456 //typedef float ScalarType;
457 typedef boost::numeric::ublas::matrix<ScalarType> MatrixType;
458
459 typedef viennacl::matrix<ScalarType, T> VCLMatrixType;
460
461 std::size_t dim_rows = 131;
462 std::size_t dim_cols = 33;
463 //std::size_t dim_rows = 4;
464 //std::size_t dim_cols = 4;
465
466 //setup ublas objects:
467 MatrixType ublas_A(dim_rows, dim_cols);
468 MatrixType ublas_B(dim_rows, dim_cols);
469 MatrixType ublas_C(dim_rows, dim_cols);
470
471 for (std::size_t i=0; i<ublas_A.size1(); ++i)
472 for (std::size_t j=0; j<ublas_A.size2(); ++j)
473 {
474 ublas_A(i,j) = ScalarType((i+2) + (j+1)*(i+2));
475 ublas_B(i,j) = ScalarType((j+2) + (j+1)*(j+2));
476 ublas_C(i,j) = ScalarType((i+1) + (i+1)*(i+2));
477 }
478
479 MatrixType ublas_A_large(4 * dim_rows, 4 * dim_cols);
480 for (std::size_t i=0; i<ublas_A_large.size1(); ++i)
481 for (std::size_t j=0; j<ublas_A_large.size2(); ++j)
482 ublas_A_large(i,j) = ScalarType(i * ublas_A_large.size2() + j);
483
484 //Setup ViennaCL objects
485 VCLMatrixType vcl_A_full(4 * dim_rows, 4 * dim_cols);
486 VCLMatrixType vcl_B_full(4 * dim_rows, 4 * dim_cols);
487 VCLMatrixType vcl_C_full(4 * dim_rows, 4 * dim_cols);
488
489 viennacl::copy(ublas_A_large, vcl_A_full);
490 viennacl::copy(ublas_A_large, vcl_B_full);
491 viennacl::copy(ublas_A_large, vcl_C_full);
492
493 //
494 // Create A
495 //
496 VCLMatrixType vcl_A(dim_rows, dim_cols);
497
498 viennacl::range vcl_A_r1(2 * dim_rows, 3 * dim_rows);
499 viennacl::range vcl_A_r2(dim_cols, 2 * dim_cols);
500 viennacl::matrix_range<VCLMatrixType> vcl_range_A(vcl_A_full, vcl_A_r1, vcl_A_r2);
501
502 viennacl::slice vcl_A_s1(2, 3, dim_rows);
503 viennacl::slice vcl_A_s2(2 * dim_cols, 2, dim_cols);
504 viennacl::matrix_slice<VCLMatrixType> vcl_slice_A(vcl_A_full, vcl_A_s1, vcl_A_s2);
505
506
507 //
508 // Create B
509 //
510 VCLMatrixType vcl_B(dim_rows, dim_cols);
511
512 viennacl::range vcl_B_r1(dim_rows, 2 * dim_rows);
513 viennacl::range vcl_B_r2(2 * dim_cols, 3 * dim_cols);
514 viennacl::matrix_range<VCLMatrixType> vcl_range_B(vcl_B_full, vcl_B_r1, vcl_B_r2);
515
516 viennacl::slice vcl_B_s1(2 * dim_rows, 2, dim_rows);
517 viennacl::slice vcl_B_s2(dim_cols, 3, dim_cols);
518 viennacl::matrix_slice<VCLMatrixType> vcl_slice_B(vcl_B_full, vcl_B_s1, vcl_B_s2);
519
520
521 //
522 // Create C
523 //
524 VCLMatrixType vcl_C(dim_rows, dim_cols);
525
526 viennacl::range vcl_C_r1(2 * dim_rows, 3 * dim_rows);
527 viennacl::range vcl_C_r2(3 * dim_cols, 4 * dim_cols);
528 viennacl::matrix_range<VCLMatrixType> vcl_range_C(vcl_C_full, vcl_C_r1, vcl_C_r2);
529
530 viennacl::slice vcl_C_s1(dim_rows, 2, dim_rows);
531 viennacl::slice vcl_C_s2(0, 3, dim_cols);
532 viennacl::matrix_slice<VCLMatrixType> vcl_slice_C(vcl_C_full, vcl_C_s1, vcl_C_s2);
533
534 viennacl::copy(ublas_A, vcl_A);
535 viennacl::copy(ublas_A, vcl_range_A);
536 viennacl::copy(ublas_A, vcl_slice_A);
537
538 viennacl::copy(ublas_B, vcl_B);
539 viennacl::copy(ublas_B, vcl_range_B);
540 viennacl::copy(ublas_B, vcl_slice_B);
541
542 viennacl::copy(ublas_C, vcl_C);
543 viennacl::copy(ublas_C, vcl_range_C);
544 viennacl::copy(ublas_C, vcl_slice_C);
545
546
547 std::cout << std::endl;
548 std::cout << "//" << std::endl;
549 std::cout << "////////// Test: Copy CTOR //////////" << std::endl;
550 std::cout << "//" << std::endl;
551
552 {
553 std::cout << "Testing matrix created from range... ";
554 VCLMatrixType vcl_temp = vcl_range_A;
555 if (check_for_equality(ublas_A, vcl_temp, epsilon))
556 std::cout << "PASSED!" << std::endl;
557 else
558 {
559 std::cout << "ublas_A: " << ublas_A << std::endl;
560 std::cout << "vcl_temp: " << vcl_temp << std::endl;
561 std::cout << "vcl_range_A: " << vcl_range_A << std::endl;
562 std::cout << "vcl_A: " << vcl_A << std::endl;
563 std::cout << std::endl << "TEST failed!" << std::endl;
564 return EXIT_FAILURE;
565 }
566
567 std::cout << "Testing matrix created from slice... ";
568 VCLMatrixType vcl_temp2 = vcl_range_B;
569 if (check_for_equality(ublas_B, vcl_temp2, epsilon))
570 std::cout << "PASSED!" << std::endl;
571 else
572 {
573 std::cout << std::endl << "TEST failed!" << std::endl;
574 return EXIT_FAILURE;
575 }
576 }
577
578 std::cout << "//" << std::endl;
579 std::cout << "////////// Test: Initializer for matrix type //////////" << std::endl;
580 std::cout << "//" << std::endl;
581
582 {
583 ublas::matrix<ScalarType> ublas_dummy1 = ublas::identity_matrix<ScalarType>(ublas_A.size1());
584 ublas::matrix<ScalarType> ublas_dummy2 = ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
585 ublas::matrix<ScalarType> ublas_dummy3 = ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
586
587 viennacl::matrix<ScalarType> vcl_dummy1 = viennacl::identity_matrix<ScalarType>(ublas_A.size1());
588 viennacl::matrix<ScalarType> vcl_dummy2 = viennacl::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
589 viennacl::matrix<ScalarType> vcl_dummy3 = viennacl::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
590
591 std::cout << "Testing initializer CTOR... ";
592 if ( check_for_equality(ublas_dummy1, vcl_dummy1, epsilon)
593 && check_for_equality(ublas_dummy2, vcl_dummy2, epsilon)
594 && check_for_equality(ublas_dummy3, vcl_dummy3, epsilon)
595 )
596 std::cout << "PASSED!" << std::endl;
597 else
598 {
599 std::cout << std::endl << "TEST failed!" << std::endl;
600 return EXIT_FAILURE;
601 }
602
603 ublas_dummy1 = ublas::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
604 ublas_dummy2 = ublas::identity_matrix<ScalarType>(ublas_A.size1());
605 ublas_dummy3 = ublas::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
606
607 vcl_dummy1 = viennacl::zero_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1());
608 vcl_dummy2 = viennacl::identity_matrix<ScalarType>(ublas_A.size1());
609 vcl_dummy3 = viennacl::scalar_matrix<ScalarType>(ublas_A.size1(), ublas_A.size1(), 3.0);
610
611 std::cout << "Testing initializer assignment... ";
612 if ( check_for_equality(ublas_dummy1, vcl_dummy1, epsilon)
613 && check_for_equality(ublas_dummy2, vcl_dummy2, epsilon)
614 && check_for_equality(ublas_dummy3, vcl_dummy3, epsilon)
615 )
616 std::cout << "PASSED!" << std::endl;
617 else
618 {
619 std::cout << std::endl << "TEST failed!" << std::endl;
620 return EXIT_FAILURE;
621 }
622 }
623
624
625 //
626 // run operation tests:
627 //
628
629 /////// A=matrix:
630 std::cout << "Testing A=matrix, B=matrix, C=matrix ..." << std::endl;
631 viennacl::copy(ublas_A, vcl_A);
632 viennacl::copy(ublas_B, vcl_B);
633 viennacl::copy(ublas_C, vcl_C);
634 if (run_test(epsilon,
635 ublas_A, ublas_B, ublas_C,
636 vcl_A, vcl_B, vcl_C) != EXIT_SUCCESS)
637 {
638 return EXIT_FAILURE;
639 }
640
641 std::cout << "Testing A=matrix, B=matrix, C=range ..." << std::endl;
642 viennacl::copy(ublas_A, vcl_A);
643 viennacl::copy(ublas_B, vcl_B);
644 viennacl::copy(ublas_C, vcl_range_C);
645 if (run_test(epsilon,
646 ublas_A, ublas_B, ublas_C,
647 vcl_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
648 {
649 return EXIT_FAILURE;
650 }
651
652 std::cout << "Testing A=matrix, B=matrix, C=slice ..." << std::endl;
653 viennacl::copy(ublas_A, vcl_A);
654 viennacl::copy(ublas_B, vcl_B);
655 viennacl::copy(ublas_C, vcl_slice_C);
656 if (run_test(epsilon,
657 ublas_A, ublas_B, ublas_C,
658 vcl_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
659 {
660 return EXIT_FAILURE;
661 }
662
663 std::cout << "Testing A=matrix, B=range, C=matrix ..." << std::endl;
664 viennacl::copy(ublas_A, vcl_A);
665 viennacl::copy(ublas_B, vcl_range_B);
666 viennacl::copy(ublas_C, vcl_C);
667 if (run_test(epsilon,
668 ublas_A, ublas_B, ublas_C,
669 vcl_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
670 {
671 return EXIT_FAILURE;
672 }
673
674 std::cout << "Testing A=matrix, B=range, C=range ..." << std::endl;
675 viennacl::copy(ublas_A, vcl_A);
676 viennacl::copy(ublas_B, vcl_range_B);
677 viennacl::copy(ublas_C, vcl_range_C);
678 if (run_test(epsilon,
679 ublas_A, ublas_B, ublas_C,
680 vcl_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
681 {
682 return EXIT_FAILURE;
683 }
684
685 std::cout << "Testing A=matrix, B=range, C=slice ..." << std::endl;
686 viennacl::copy(ublas_A, vcl_A);
687 viennacl::copy(ublas_B, vcl_range_B);
688 viennacl::copy(ublas_C, vcl_slice_C);
689 if (run_test(epsilon,
690 ublas_A, ublas_B, ublas_C,
691 vcl_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
692 {
693 return EXIT_FAILURE;
694 }
695
696
697 std::cout << "Testing A=matrix, B=slice, C=matrix ..." << std::endl;
698 viennacl::copy(ublas_A, vcl_A);
699 viennacl::copy(ublas_B, vcl_slice_B);
700 viennacl::copy(ublas_C, vcl_C);
701 if (run_test(epsilon,
702 ublas_A, ublas_B, ublas_C,
703 vcl_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
704 {
705 return EXIT_FAILURE;
706 }
707
708 std::cout << "Testing A=matrix, B=slice, C=range ..." << std::endl;
709 viennacl::copy(ublas_A, vcl_A);
710 viennacl::copy(ublas_B, vcl_slice_B);
711 viennacl::copy(ublas_C, vcl_range_C);
712 if (run_test(epsilon,
713 ublas_A, ublas_B, ublas_C,
714 vcl_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
715 {
716 return EXIT_FAILURE;
717 }
718
719 std::cout << "Testing A=matrix, B=slice, C=slice ..." << std::endl;
720 viennacl::copy(ublas_A, vcl_A);
721 viennacl::copy(ublas_B, vcl_slice_B);
722 viennacl::copy(ublas_C, vcl_slice_C);
723 if (run_test(epsilon,
724 ublas_A, ublas_B, ublas_C,
725 vcl_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
726 {
727 return EXIT_FAILURE;
728 }
729
730
731
732 /////// A=range:
733 std::cout << "Testing A=range, B=matrix, C=matrix ..." << std::endl;
734 viennacl::copy(ublas_A, vcl_range_A);
735 viennacl::copy(ublas_B, vcl_B);
736 viennacl::copy(ublas_C, vcl_C);
737 if (run_test(epsilon,
738 ublas_A, ublas_B, ublas_C,
739 vcl_range_A, vcl_B, vcl_C) != EXIT_SUCCESS)
740 {
741 return EXIT_FAILURE;
742 }
743
744 std::cout << "Testing A=range, B=matrix, C=range ..." << std::endl;
745 viennacl::copy(ublas_A, vcl_range_A);
746 viennacl::copy(ublas_B, vcl_B);
747 viennacl::copy(ublas_C, vcl_range_C);
748 if (run_test(epsilon,
749 ublas_A, ublas_B, ublas_C,
750 vcl_range_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
751 {
752 return EXIT_FAILURE;
753 }
754
755 std::cout << "Testing A=range, B=matrix, C=slice ..." << std::endl;
756 viennacl::copy(ublas_A, vcl_range_A);
757 viennacl::copy(ublas_B, vcl_B);
758 viennacl::copy(ublas_C, vcl_slice_C);
759 if (run_test(epsilon,
760 ublas_A, ublas_B, ublas_C,
761 vcl_range_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
762 {
763 return EXIT_FAILURE;
764 }
765
766
767
768 std::cout << "Testing A=range, B=range, C=matrix ..." << std::endl;
769 viennacl::copy(ublas_A, vcl_range_A);
770 viennacl::copy(ublas_B, vcl_range_B);
771 viennacl::copy(ublas_C, vcl_C);
772 if (run_test(epsilon,
773 ublas_A, ublas_B, ublas_C,
774 vcl_range_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
775 {
776 return EXIT_FAILURE;
777 }
778
779 std::cout << "Testing A=range, B=range, C=range ..." << std::endl;
780 viennacl::copy(ublas_A, vcl_range_A);
781 viennacl::copy(ublas_B, vcl_range_B);
782 viennacl::copy(ublas_C, vcl_range_C);
783 if (run_test(epsilon,
784 ublas_A, ublas_B, ublas_C,
785 vcl_range_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
786 {
787 return EXIT_FAILURE;
788 }
789
790 std::cout << "Testing A=range, B=range, C=slice ..." << std::endl;
791 viennacl::copy(ublas_A, vcl_range_A);
792 viennacl::copy(ublas_B, vcl_range_B);
793 viennacl::copy(ublas_C, vcl_slice_C);
794 if (run_test(epsilon,
795 ublas_A, ublas_B, ublas_C,
796 vcl_range_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
797 {
798 return EXIT_FAILURE;
799 }
800
801
802
803 std::cout << "Testing A=range, B=slice, C=matrix ..." << std::endl;
804 viennacl::copy(ublas_A, vcl_range_A);
805 viennacl::copy(ublas_B, vcl_slice_B);
806 viennacl::copy(ublas_C, vcl_C);
807 if (run_test(epsilon,
808 ublas_A, ublas_B, ublas_C,
809 vcl_range_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
810 {
811 return EXIT_FAILURE;
812 }
813
814 std::cout << "Testing A=range, B=slice, C=range ..." << std::endl;
815 viennacl::copy(ublas_A, vcl_range_A);
816 viennacl::copy(ublas_B, vcl_slice_B);
817 viennacl::copy(ublas_C, vcl_range_C);
818 if (run_test(epsilon,
819 ublas_A, ublas_B, ublas_C,
820 vcl_range_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
821 {
822 return EXIT_FAILURE;
823 }
824
825 std::cout << "Testing A=range, B=slice, C=slice ..." << std::endl;
826 viennacl::copy(ublas_A, vcl_range_A);
827 viennacl::copy(ublas_B, vcl_slice_B);
828 viennacl::copy(ublas_C, vcl_slice_C);
829 if (run_test(epsilon,
830 ublas_A, ublas_B, ublas_C,
831 vcl_range_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
832 {
833 return EXIT_FAILURE;
834 }
835
836
837 /////// A=slice:
838 std::cout << "Testing A=slice, B=matrix, C=matrix ..." << std::endl;
839 viennacl::copy(ublas_A, vcl_slice_A);
840 viennacl::copy(ublas_B, vcl_B);
841 viennacl::copy(ublas_C, vcl_C);
842 if (run_test(epsilon,
843 ublas_A, ublas_B, ublas_C,
844 vcl_slice_A, vcl_B, vcl_C) != EXIT_SUCCESS)
845 {
846 return EXIT_FAILURE;
847 }
848
849 std::cout << "Testing A=slice, B=matrix, C=range ..." << std::endl;
850 viennacl::copy(ublas_A, vcl_slice_A);
851 viennacl::copy(ublas_B, vcl_B);
852 viennacl::copy(ublas_C, vcl_range_C);
853 if (run_test(epsilon,
854 ublas_A, ublas_B, ublas_C,
855 vcl_slice_A, vcl_B, vcl_range_C) != EXIT_SUCCESS)
856 {
857 return EXIT_FAILURE;
858 }
859
860 std::cout << "Testing A=slice, B=matrix, C=slice ..." << std::endl;
861 viennacl::copy(ublas_A, vcl_slice_A);
862 viennacl::copy(ublas_B, vcl_B);
863 viennacl::copy(ublas_C, vcl_slice_C);
864 if (run_test(epsilon,
865 ublas_A, ublas_B, ublas_C,
866 vcl_slice_A, vcl_B, vcl_slice_C) != EXIT_SUCCESS)
867 {
868 return EXIT_FAILURE;
869 }
870
871
872
873 std::cout << "Testing A=slice, B=range, C=matrix ..." << std::endl;
874 viennacl::copy(ublas_A, vcl_slice_A);
875 viennacl::copy(ublas_B, vcl_range_B);
876 viennacl::copy(ublas_C, vcl_C);
877 if (run_test(epsilon,
878 ublas_A, ublas_B, ublas_C,
879 vcl_slice_A, vcl_range_B, vcl_C) != EXIT_SUCCESS)
880 {
881 return EXIT_FAILURE;
882 }
883
884 std::cout << "Testing A=slice, B=range, C=range ..." << std::endl;
885 viennacl::copy(ublas_A, vcl_slice_A);
886 viennacl::copy(ublas_B, vcl_range_B);
887 viennacl::copy(ublas_C, vcl_range_C);
888 if (run_test(epsilon,
889 ublas_A, ublas_B, ublas_C,
890 vcl_slice_A, vcl_range_B, vcl_range_C) != EXIT_SUCCESS)
891 {
892 return EXIT_FAILURE;
893 }
894
895 std::cout << "Testing A=slice, B=range, C=slice ..." << std::endl;
896 viennacl::copy(ublas_A, vcl_slice_A);
897 viennacl::copy(ublas_B, vcl_range_B);
898 viennacl::copy(ublas_C, vcl_slice_C);
899 if (run_test(epsilon,
900 ublas_A, ublas_B, ublas_C,
901 vcl_slice_A, vcl_range_B, vcl_slice_C) != EXIT_SUCCESS)
902 {
903 return EXIT_FAILURE;
904 }
905
906
907
908 std::cout << "Testing A=slice, B=slice, C=matrix ..." << std::endl;
909 viennacl::copy(ublas_A, vcl_slice_A);
910 viennacl::copy(ublas_B, vcl_slice_B);
911 viennacl::copy(ublas_C, vcl_C);
912 if (run_test(epsilon,
913 ublas_A, ublas_B, ublas_C,
914 vcl_slice_A, vcl_slice_B, vcl_C) != EXIT_SUCCESS)
915 {
916 return EXIT_FAILURE;
917 }
918
919 std::cout << "Testing A=slice, B=slice, C=range ..." << std::endl;
920 viennacl::copy(ublas_A, vcl_slice_A);
921 viennacl::copy(ublas_B, vcl_slice_B);
922 viennacl::copy(ublas_C, vcl_range_C);
923 if (run_test(epsilon,
924 ublas_A, ublas_B, ublas_C,
925 vcl_slice_A, vcl_slice_B, vcl_range_C) != EXIT_SUCCESS)
926 {
927 return EXIT_FAILURE;
928 }
929
930 std::cout << "Testing A=slice, B=slice, C=slice ..." << std::endl;
931 viennacl::copy(ublas_A, vcl_slice_A);
932 viennacl::copy(ublas_B, vcl_slice_B);
933 viennacl::copy(ublas_C, vcl_slice_C);
934 if (run_test(epsilon,
935 ublas_A, ublas_B, ublas_C,
936 vcl_slice_A, vcl_slice_B, vcl_slice_C) != EXIT_SUCCESS)
937 {
938 return EXIT_FAILURE;
939 }
940
941
942 return EXIT_SUCCESS;
943 }
944
main(int,const char **)945 int main (int, const char **)
946 {
947 std::cout << std::endl;
948 std::cout << "----------------------------------------------" << std::endl;
949 std::cout << "----------------------------------------------" << std::endl;
950 std::cout << "## Test :: Matrix Range" << std::endl;
951 std::cout << "----------------------------------------------" << std::endl;
952 std::cout << "----------------------------------------------" << std::endl;
953 std::cout << std::endl;
954
955 double epsilon = 1e-4;
956 std::cout << "# Testing setup:" << std::endl;
957 std::cout << " eps: " << epsilon << std::endl;
958 std::cout << " numeric: float" << std::endl;
959 std::cout << " --- row-major ---" << std::endl;
960 if (run_test<viennacl::row_major, float>(epsilon) != EXIT_SUCCESS)
961 return EXIT_FAILURE;
962 std::cout << " --- column-major ---" << std::endl;
963 if (run_test<viennacl::column_major, float>(epsilon) != EXIT_SUCCESS)
964 return EXIT_FAILURE;
965
966
967 #ifdef VIENNACL_WITH_OPENCL
968 if ( viennacl::ocl::current_device().double_support() )
969 #endif
970 {
971 epsilon = 1e-12;
972 std::cout << "# Testing setup:" << std::endl;
973 std::cout << " eps: " << epsilon << std::endl;
974 std::cout << " numeric: double" << std::endl;
975
976 if (run_test<viennacl::row_major, double>(epsilon) != EXIT_SUCCESS)
977 return EXIT_FAILURE;
978 if (run_test<viennacl::column_major, double>(epsilon) != EXIT_SUCCESS)
979 return EXIT_FAILURE;
980 }
981
982 std::cout << std::endl;
983 std::cout << "------- Test completed --------" << std::endl;
984 std::cout << std::endl;
985
986
987 return EXIT_SUCCESS;
988 }
989
990