ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
scheduler_matrix.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2015, 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 
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"*/
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 
55 
56 using namespace boost::numeric;
57 
58 template<typename MatrixType, typename VCLMatrixType>
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>
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 
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);
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);
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 
364  GENERATE_UNARY_OP_TEST(floor);
367  GENERATE_UNARY_OP_TEST(log10);
371  //GENERATE_UNARY_OP_TEST(abs); //OpenCL allows abs on integers only
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>
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 
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 
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 
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 
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 
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 
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_div > > element_div(vector_base< T > const &v1, vector_base< T > const &v2)
int main(int, const char **)
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
Class for representing strided submatrices of a bigger matrix A.
Definition: forwards.h:443
void trans(matrix_expression< const matrix_base< NumericT, SizeT, DistanceT >, const matrix_base< NumericT, SizeT, DistanceT >, op_trans > const &proxy, matrix_base< NumericT > &temp_trans)
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Implementation of the dense matrix class.
A tag class representing assignment.
Definition: forwards.h:81
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
A dense matrix class.
Definition: forwards.h:375
void execute(statement const &s)
Definition: execute.hpp:279
Represents a vector consisting of 1 at a given index and zeros otherwise. To be used as an initialize...
Definition: matrix_def.hpp:69
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
A tag class representing inplace addition.
Definition: forwards.h:83
bool check_for_equality(MatrixType const &ublas_A, VCLMatrixType const &vcl_A, double epsilon)
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: matrix_def.hpp:93
#define GENERATE_UNARY_OP_TEST(OPNAME)
Represents a vector consisting of zeros only. To be used as an initializer for viennacl::vector, vector_range, or vector_slize only.
Definition: matrix_def.hpp:81
T::ERROR_CANNOT_DEDUCE_CPU_SCALAR_TYPE_FOR_T type
Definition: result_of.hpp:271
int run_test(double epsilon, UBLASMatrixType &ublas_A, UBLASMatrixType &ublas_B, UBLASMatrixType &ublas_C, ViennaCLMatrixType1 &vcl_A, ViennaCLMatrixType2 &vcl_B, ViennaCLMatrixType3 vcl_C)
Proxy classes for vectors.
Proxy classes for matrices.
A tag class representing inplace subtraction.
Definition: forwards.h:85
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
A range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:424
float ScalarType
Definition: fft_1d.cpp:42
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_prod > > element_prod(vector_base< T > const &v1, vector_base< T > const &v2)
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
Definition: forwards.h:502
A slice class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:429
Implementation of the ViennaCL scalar class.