ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
scheduler_matrix_vector.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 //
25 // *** System
26 //
27 #include <iostream>
28 
29 //
30 // *** Boost
31 //
32 #include <boost/numeric/ublas/io.hpp>
33 #include <boost/numeric/ublas/triangular.hpp>
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/matrix.hpp>
36 #include <boost/numeric/ublas/matrix_proxy.hpp>
37 #include <boost/numeric/ublas/lu.hpp>
38 #include <boost/numeric/ublas/io.hpp>
39 
40 //
41 // *** ViennaCL
42 //
43 //#define VIENNACL_DEBUG_ALL
44 #define VIENNACL_WITH_UBLAS 1
45 #include "viennacl/scalar.hpp"
46 #include "viennacl/matrix.hpp"
47 #include "viennacl/vector.hpp"
48 #include "viennacl/linalg/prod.hpp"
49 #include "viennacl/linalg/norm_2.hpp"
50 #include "viennacl/linalg/direct_solve.hpp"
51 #include "viennacl/linalg/lu.hpp"
52 #include "viennacl/tools/random.hpp"
53 
54 #include "viennacl/scheduler/execute.hpp"
55 #include "viennacl/scheduler/io.hpp"
56 
57 //
58 // -------------------------------------------------------------
59 //
60 using namespace boost::numeric;
61 //
62 // -------------------------------------------------------------
63 //
64 template<typename ScalarType>
65 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
66 {
67  viennacl::backend::finish();
68  if (s1 != s2)
69  return (s1 - s2) / std::max(std::fabs(s1), std::fabs(s2));
70  return 0;
71 }
72 
73 template<typename ScalarType, typename VCLVectorType>
74 ScalarType diff(ublas::vector<ScalarType> const & v1, VCLVectorType const & v2)
75 {
76  ublas::vector<ScalarType> v2_cpu(v2.size());
77  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
78  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
79 
80  for (unsigned int i=0;i<v1.size(); ++i)
81  {
82  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
83  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
84  else
85  v2_cpu[i] = 0.0;
86  }
87 
88  return norm_inf(v2_cpu);
89 }
90 
91 template<typename ScalarType, typename VCLMatrixType>
92 ScalarType diff(ublas::matrix<ScalarType> const & mat1, VCLMatrixType const & mat2)
93 {
94  ublas::matrix<ScalarType> mat2_cpu(mat2.size1(), mat2.size2());
95  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
96  viennacl::copy(mat2, mat2_cpu);
97  ScalarType ret = 0;
98  ScalarType act = 0;
99 
100  for (unsigned int i = 0; i < mat2_cpu.size1(); ++i)
101  {
102  for (unsigned int j = 0; j < mat2_cpu.size2(); ++j)
103  {
104  act = std::fabs(mat2_cpu(i,j) - mat1(i,j)) / std::max( std::fabs(mat2_cpu(i, j)), std::fabs(mat1(i,j)) );
105  if (act > ret)
106  ret = act;
107  }
108  }
109  //std::cout << ret << std::endl;
110  return ret;
111 }
112 //
113 // -------------------------------------------------------------
114 //
115 
116 template<typename NumericT, typename Epsilon,
117  typename UblasMatrixType, typename UblasVectorType,
118  typename VCLMatrixType, typename VCLVectorType1, typename VCLVectorType2>
119 int test_prod_rank1(Epsilon const & epsilon,
120  UblasMatrixType & ublas_m1, UblasVectorType & ublas_v1, UblasVectorType & ublas_v2,
121  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1, VCLVectorType2 & vcl_v2)
122 {
123  int retval = EXIT_SUCCESS;
124 
125  // sync data:
126  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
127  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
128  viennacl::copy(ublas_m1, vcl_m1);
129 
130  /* TODO: Add rank-1 operations here */
131 
132  //reset vcl_matrix:
133  viennacl::copy(ublas_m1, vcl_m1);
134 
135  // --------------------------------------------------------------------------
136  std::cout << "Matrix-Vector product" << std::endl;
137  ublas_v1 = viennacl::linalg::prod(ublas_m1, ublas_v2);
138  {
139  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(vcl_m1, vcl_v2));
140  viennacl::scheduler::execute(my_statement);
141  }
142 
143  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
144  {
145  std::cout << "# Error at operation: matrix-vector product" << std::endl;
146  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
147  retval = EXIT_FAILURE;
148  }
149 
150  std::cout << "Matrix-Vector product with inplace-add" << std::endl;
151  ublas_v1 += viennacl::linalg::prod(ublas_m1, ublas_v2);
152  {
153  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), viennacl::linalg::prod(vcl_m1, vcl_v2));
154  viennacl::scheduler::execute(my_statement);
155  }
156 
157  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
158  {
159  std::cout << "# Error at operation: matrix-vector product" << std::endl;
160  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
161  retval = EXIT_FAILURE;
162  }
163 
164  std::cout << "Matrix-Vector product with inplace-sub" << std::endl;
165  ublas_v1 -= viennacl::linalg::prod(ublas_m1, ublas_v2);
166  {
167  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), viennacl::linalg::prod(vcl_m1, vcl_v2));
168  viennacl::scheduler::execute(my_statement);
169  }
170 
171  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
172  {
173  std::cout << "# Error at operation: matrix-vector product" << std::endl;
174  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
175  retval = EXIT_FAILURE;
176  }
177 
178  // --------------------------------------------------------------------------
179  /*
180  std::cout << "Matrix-Vector product with scaled matrix" << std::endl;
181  ublas_v1 = viennacl::linalg::prod(NumericT(2.0) * ublas_m1, ublas_v2);
182  {
183  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(NumericT(2.0) * vcl_m1, vcl_v2));
184  viennacl::scheduler::execute(my_statement);
185  }
186 
187  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
188  {
189  std::cout << "# Error at operation: matrix-vector product" << std::endl;
190  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
191  retval = EXIT_FAILURE;
192  }*/
193 
194  // --------------------------------------------------------------------------
195  std::cout << "Matrix-Vector product with scaled vector" << std::endl;
196  /*
197  ublas_v1 = viennacl::linalg::prod(ublas_m1, NumericT(2.0) * ublas_v2);
198  {
199  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(vcl_m1, NumericT(2.0) * vcl_v2));
200  viennacl::scheduler::execute(my_statement);
201  }
202 
203  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
204  {
205  std::cout << "# Error at operation: matrix-vector product" << std::endl;
206  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
207  retval = EXIT_FAILURE;
208  }*/
209 
210  // --------------------------------------------------------------------------
211  std::cout << "Matrix-Vector product with scaled matrix and scaled vector" << std::endl;
212  /*
213  ublas_v1 = viennacl::linalg::prod(NumericT(2.0) * ublas_m1, NumericT(2.0) * ublas_v2);
214  {
215  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), viennacl::linalg::prod(NumericT(2.0) * vcl_m1, NumericT(2.0) * vcl_v2));
216  viennacl::scheduler::execute(my_statement);
217  }
218 
219  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
220  {
221  std::cout << "# Error at operation: matrix-vector product" << std::endl;
222  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
223  retval = EXIT_FAILURE;
224  }*/
225 
226 
227  // --------------------------------------------------------------------------
228  std::cout << "Matrix-Vector product with scaled add" << std::endl;
229  NumericT alpha = static_cast<NumericT>(2.786);
230  NumericT beta = static_cast<NumericT>(3.1415);
231  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
232  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
233 
234  ublas_v1 = alpha * viennacl::linalg::prod(ublas_m1, ublas_v2) - beta * ublas_v1;
235  {
236  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
237  viennacl::scheduler::execute(my_statement);
238  }
239 
240  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
241  {
242  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
243  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
244  retval = EXIT_FAILURE;
245  }
246 
247  std::cout << "Matrix-Vector product with scaled add, inplace-add" << std::endl;
248  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
249  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
250 
251  ublas_v1 += alpha * viennacl::linalg::prod(ublas_m1, ublas_v2) - beta * ublas_v1;
252  {
253  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_add(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
254  viennacl::scheduler::execute(my_statement);
255  }
256 
257  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
258  {
259  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
260  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
261  retval = EXIT_FAILURE;
262  }
263 
264  std::cout << "Matrix-Vector product with scaled add, inplace-sub" << std::endl;
265  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
266  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
267 
268  ublas_v1 -= alpha * viennacl::linalg::prod(ublas_m1, ublas_v2) - beta * ublas_v1;
269  {
270  viennacl::scheduler::statement my_statement(vcl_v1, viennacl::op_inplace_sub(), alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) - beta * vcl_v1);
271  viennacl::scheduler::execute(my_statement);
272  }
273 
274  if ( std::fabs(diff(ublas_v1, vcl_v1)) > epsilon )
275  {
276  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
277  std::cout << " diff: " << std::fabs(diff(ublas_v1, vcl_v1)) << std::endl;
278  retval = EXIT_FAILURE;
279  }
280 
281  // --------------------------------------------------------------------------
282 
283  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
284  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
285 
286  std::cout << "Transposed Matrix-Vector product" << std::endl;
287  ublas_v2 = viennacl::linalg::prod(trans(ublas_m1), ublas_v1);
288  {
289  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_assign(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
290  viennacl::scheduler::execute(my_statement);
291  }
292 
293  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
294  {
295  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
296  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
297  retval = EXIT_FAILURE;
298  }
299 
300  std::cout << "Transposed Matrix-Vector product, inplace-add" << std::endl;
301  ublas_v2 += viennacl::linalg::prod(trans(ublas_m1), ublas_v1);
302  {
303  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_add(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
304  viennacl::scheduler::execute(my_statement);
305  }
306 
307  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
308  {
309  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
310  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
311  retval = EXIT_FAILURE;
312  }
313 
314  std::cout << "Transposed Matrix-Vector product, inplace-sub" << std::endl;
315  ublas_v2 -= viennacl::linalg::prod(trans(ublas_m1), ublas_v1);
316  {
317  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_sub(), viennacl::linalg::prod(trans(vcl_m1), vcl_v1));
318  viennacl::scheduler::execute(my_statement);
319  }
320 
321  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
322  {
323  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
324  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
325  retval = EXIT_FAILURE;
326  }
327 
328  // --------------------------------------------------------------------------
329  std::cout << "Transposed Matrix-Vector product with scaled add" << std::endl;
330  ublas_v2 = alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1) + beta * ublas_v2;
331  {
332  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_assign(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
333  viennacl::scheduler::execute(my_statement);
334  }
335 
336  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
337  {
338  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
339  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
340  retval = EXIT_FAILURE;
341  }
342 
343  std::cout << "Transposed Matrix-Vector product with scaled add, inplace-add" << std::endl;
344  ublas_v2 += alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1) + beta * ublas_v2;
345  {
346  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_add(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
347  viennacl::scheduler::execute(my_statement);
348  }
349 
350  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
351  {
352  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
353  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
354  retval = EXIT_FAILURE;
355  }
356 
357  std::cout << "Transposed Matrix-Vector product with scaled add, inplace-sub" << std::endl;
358  ublas_v2 -= alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1) + beta * ublas_v2;
359  {
360  viennacl::scheduler::statement my_statement(vcl_v2, viennacl::op_inplace_sub(), alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2);
361  viennacl::scheduler::execute(my_statement);
362  }
363 
364  if ( std::fabs(diff(ublas_v2, vcl_v2)) > epsilon )
365  {
366  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
367  std::cout << " diff: " << std::fabs(diff(ublas_v2, vcl_v2)) << std::endl;
368  retval = EXIT_FAILURE;
369  }
370 
371  // --------------------------------------------------------------------------
372 
373  return retval;
374 }
375 
376 
377 //
378 // -------------------------------------------------------------
379 //
380 template< typename NumericT, typename F, typename Epsilon >
381 int test(Epsilon const& epsilon)
382 {
383  int retval = EXIT_SUCCESS;
384 
385  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
386 
387  std::size_t num_rows = 141;
388  std::size_t num_cols = 79;
389 
390  // --------------------------------------------------------------------------
391  ublas::vector<NumericT> ublas_v1(num_rows);
392  for (std::size_t i = 0; i < ublas_v1.size(); ++i)
393  ublas_v1(i) = randomNumber();
394  ublas::vector<NumericT> ublas_v2 = ublas::scalar_vector<NumericT>(num_cols, NumericT(3.1415));
395 
396 
397  ublas::matrix<NumericT> ublas_m1(ublas_v1.size(), ublas_v2.size());
398 
399  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
400  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
401  ublas_m1(i,j) = static_cast<NumericT>(0.1) * randomNumber();
402 
403 
404  ublas::matrix<NumericT> ublas_m2(ublas_v1.size(), ublas_v1.size());
405 
406  for (std::size_t i = 0; i < ublas_m2.size1(); ++i)
407  {
408  for (std::size_t j = 0; j < ublas_m2.size2(); ++j)
409  ublas_m2(i,j) = static_cast<NumericT>(-0.1) * randomNumber();
410  ublas_m2(i, i) = static_cast<NumericT>(2) + randomNumber();
411  }
412 
413 
414  viennacl::vector<NumericT> vcl_v1_native(ublas_v1.size());
415  viennacl::vector<NumericT> vcl_v1_large(4 * ublas_v1.size());
416  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v1_range(vcl_v1_large, viennacl::range(3, ublas_v1.size() + 3));
417  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v1_slice(vcl_v1_large, viennacl::slice(2, 3, ublas_v1.size()));
418 
419  viennacl::vector<NumericT> vcl_v2_native(ublas_v2.size());
420  viennacl::vector<NumericT> vcl_v2_large(4 * ublas_v2.size());
421  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v2_range(vcl_v2_large, viennacl::range(8, ublas_v2.size() + 8));
422  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v2_slice(vcl_v2_large, viennacl::slice(6, 2, ublas_v2.size()));
423 
424  viennacl::matrix<NumericT, F> vcl_m1_native(ublas_m1.size1(), ublas_m1.size2());
425  viennacl::matrix<NumericT, F> vcl_m1_large(4 * ublas_m1.size1(), 4 * ublas_m1.size2());
426  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m1_range(vcl_m1_large,
427  viennacl::range(8, ublas_m1.size1() + 8),
428  viennacl::range(ublas_m1.size2(), 2 * ublas_m1.size2()) );
429  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m1_slice(vcl_m1_large,
430  viennacl::slice(6, 2, ublas_m1.size1()),
431  viennacl::slice(ublas_m1.size2(), 2, ublas_m1.size2()) );
432 
433  viennacl::matrix<NumericT, F> vcl_m2_native(ublas_m2.size1(), ublas_m2.size2());
434  viennacl::matrix<NumericT, F> vcl_m2_large(4 * ublas_m2.size1(), 4 * ublas_m2.size2());
435  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m2_range(vcl_m2_large,
436  viennacl::range(8, ublas_m2.size1() + 8),
437  viennacl::range(ublas_m2.size2(), 2 * ublas_m2.size2()) );
438  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m2_slice(vcl_m2_large,
439  viennacl::slice(6, 2, ublas_m2.size1()),
440  viennacl::slice(ublas_m2.size2(), 2, ublas_m2.size2()) );
441 
442 
443 /* std::cout << "Matrix resizing (to larger)" << std::endl;
444  matrix.resize(2*num_rows, 2*num_cols, true);
445  for (unsigned int i = 0; i < matrix.size1(); ++i)
446  {
447  for (unsigned int j = (i<result.size() ? rhs.size() : 0); j < matrix.size2(); ++j)
448  matrix(i,j) = 0;
449  }
450  vcl_matrix.resize(2*num_rows, 2*num_cols, true);
451  viennacl::copy(vcl_matrix, matrix);
452  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
453  {
454  std::cout << "# Error at operation: matrix resize (to larger)" << std::endl;
455  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
456  return EXIT_FAILURE;
457  }
458 
459  matrix(12, 14) = NumericT(1.9);
460  matrix(19, 16) = NumericT(1.0);
461  matrix (13, 15) = NumericT(-9);
462  vcl_matrix(12, 14) = NumericT(1.9);
463  vcl_matrix(19, 16) = NumericT(1.0);
464  vcl_matrix (13, 15) = NumericT(-9);
465 
466  std::cout << "Matrix resizing (to smaller)" << std::endl;
467  matrix.resize(result.size(), rhs.size(), true);
468  vcl_matrix.resize(result.size(), rhs.size(), true);
469  if ( std::fabs(diff(matrix, vcl_matrix)) > epsilon )
470  {
471  std::cout << "# Error at operation: matrix resize (to smaller)" << std::endl;
472  std::cout << " diff: " << std::fabs(diff(matrix, vcl_matrix)) << std::endl;
473  return EXIT_FAILURE;
474  }
475  */
476 
477  //
478  // Run a bunch of tests for rank-1-updates, matrix-vector products
479  //
480  std::cout << "------------ Testing rank-1-updates and matrix-vector products ------------------" << std::endl;
481 
482  std::cout << "* m = full, v1 = full, v2 = full" << std::endl;
483  retval = test_prod_rank1<NumericT>(epsilon,
484  ublas_m1, ublas_v1, ublas_v2,
485  vcl_m1_native, vcl_v1_native, vcl_v2_native);
486  if (retval == EXIT_FAILURE)
487  {
488  std::cout << " --- FAILED! ---" << std::endl;
489  return retval;
490  }
491  else
492  std::cout << " --- PASSED ---" << std::endl;
493 
494 
495  std::cout << "* m = full, v1 = full, v2 = range" << std::endl;
496  retval = test_prod_rank1<NumericT>(epsilon,
497  ublas_m1, ublas_v1, ublas_v2,
498  vcl_m1_native, vcl_v1_native, vcl_v2_range);
499  if (retval == EXIT_FAILURE)
500  {
501  std::cout << " --- FAILED! ---" << std::endl;
502  return retval;
503  }
504  else
505  std::cout << " --- PASSED ---" << std::endl;
506 
507 
508  std::cout << "* m = full, v1 = full, v2 = slice" << std::endl;
509  retval = test_prod_rank1<NumericT>(epsilon,
510  ublas_m1, ublas_v1, ublas_v2,
511  vcl_m1_native, vcl_v1_native, vcl_v2_slice);
512  if (retval == EXIT_FAILURE)
513  {
514  std::cout << " --- FAILED! ---" << std::endl;
515  return retval;
516  }
517  else
518  std::cout << " --- PASSED ---" << std::endl;
519 
520 
521  // v1 = range
522 
523 
524  std::cout << "* m = full, v1 = range, v2 = full" << std::endl;
525  retval = test_prod_rank1<NumericT>(epsilon,
526  ublas_m1, ublas_v1, ublas_v2,
527  vcl_m1_native, vcl_v1_range, vcl_v2_native);
528  if (retval == EXIT_FAILURE)
529  {
530  std::cout << " --- FAILED! ---" << std::endl;
531  return retval;
532  }
533  else
534  std::cout << " --- PASSED ---" << std::endl;
535 
536 
537  std::cout << "* m = full, v1 = range, v2 = range" << std::endl;
538  retval = test_prod_rank1<NumericT>(epsilon,
539  ublas_m1, ublas_v1, ublas_v2,
540  vcl_m1_native, vcl_v1_range, vcl_v2_range);
541  if (retval == EXIT_FAILURE)
542  {
543  std::cout << " --- FAILED! ---" << std::endl;
544  return retval;
545  }
546  else
547  std::cout << " --- PASSED ---" << std::endl;
548 
549 
550  std::cout << "* m = full, v1 = range, v2 = slice" << std::endl;
551  retval = test_prod_rank1<NumericT>(epsilon,
552  ublas_m1, ublas_v1, ublas_v2,
553  vcl_m1_native, vcl_v1_range, vcl_v2_slice);
554  if (retval == EXIT_FAILURE)
555  {
556  std::cout << " --- FAILED! ---" << std::endl;
557  return retval;
558  }
559  else
560  std::cout << " --- PASSED ---" << std::endl;
561 
562 
563 
564  // v1 = slice
565 
566  std::cout << "* m = full, v1 = slice, v2 = full" << std::endl;
567  retval = test_prod_rank1<NumericT>(epsilon,
568  ublas_m1, ublas_v1, ublas_v2,
569  vcl_m1_native, vcl_v1_slice, vcl_v2_native);
570  if (retval == EXIT_FAILURE)
571  {
572  std::cout << " --- FAILED! ---" << std::endl;
573  return retval;
574  }
575  else
576  std::cout << " --- PASSED ---" << std::endl;
577 
578 
579  std::cout << "* m = full, v1 = slice, v2 = range" << std::endl;
580  retval = test_prod_rank1<NumericT>(epsilon,
581  ublas_m1, ublas_v1, ublas_v2,
582  vcl_m1_native, vcl_v1_slice, vcl_v2_range);
583  if (retval == EXIT_FAILURE)
584  {
585  std::cout << " --- FAILED! ---" << std::endl;
586  return retval;
587  }
588  else
589  std::cout << " --- PASSED ---" << std::endl;
590 
591 
592  std::cout << "* m = full, v1 = slice, v2 = slice" << std::endl;
593  retval = test_prod_rank1<NumericT>(epsilon,
594  ublas_m1, ublas_v1, ublas_v2,
595  vcl_m1_native, vcl_v1_slice, vcl_v2_slice);
596  if (retval == EXIT_FAILURE)
597  {
598  std::cout << " --- FAILED! ---" << std::endl;
599  return retval;
600  }
601  else
602  std::cout << " --- PASSED ---" << std::endl;
603 
604 
606 
607  std::cout << "* m = range, v1 = full, v2 = full" << std::endl;
608  retval = test_prod_rank1<NumericT>(epsilon,
609  ublas_m1, ublas_v1, ublas_v2,
610  vcl_m1_range, vcl_v1_native, vcl_v2_native);
611  if (retval == EXIT_FAILURE)
612  {
613  std::cout << " --- FAILED! ---" << std::endl;
614  return retval;
615  }
616  else
617  std::cout << " --- PASSED ---" << std::endl;
618 
619 
620  std::cout << "* m = range, v1 = full, v2 = range" << std::endl;
621  retval = test_prod_rank1<NumericT>(epsilon,
622  ublas_m1, ublas_v1, ublas_v2,
623  vcl_m1_range, vcl_v1_native, vcl_v2_range);
624  if (retval == EXIT_FAILURE)
625  {
626  std::cout << " --- FAILED! ---" << std::endl;
627  return retval;
628  }
629  else
630  std::cout << " --- PASSED ---" << std::endl;
631 
632 
633  std::cout << "* m = range, v1 = full, v2 = slice" << std::endl;
634  retval = test_prod_rank1<NumericT>(epsilon,
635  ublas_m1, ublas_v1, ublas_v2,
636  vcl_m1_range, vcl_v1_native, vcl_v2_slice);
637  if (retval == EXIT_FAILURE)
638  {
639  std::cout << " --- FAILED! ---" << std::endl;
640  return retval;
641  }
642  else
643  std::cout << " --- PASSED ---" << std::endl;
644 
645 
646  // v1 = range
647 
648 
649  std::cout << "* m = range, v1 = range, v2 = full" << std::endl;
650  retval = test_prod_rank1<NumericT>(epsilon,
651  ublas_m1, ublas_v1, ublas_v2,
652  vcl_m1_range, vcl_v1_range, vcl_v2_native);
653  if (retval == EXIT_FAILURE)
654  {
655  std::cout << " --- FAILED! ---" << std::endl;
656  return retval;
657  }
658  else
659  std::cout << " --- PASSED ---" << std::endl;
660 
661 
662  std::cout << "* m = range, v1 = range, v2 = range" << std::endl;
663  retval = test_prod_rank1<NumericT>(epsilon,
664  ublas_m1, ublas_v1, ublas_v2,
665  vcl_m1_range, vcl_v1_range, vcl_v2_range);
666  if (retval == EXIT_FAILURE)
667  {
668  std::cout << " --- FAILED! ---" << std::endl;
669  return retval;
670  }
671  else
672  std::cout << " --- PASSED ---" << std::endl;
673 
674 
675  std::cout << "* m = range, v1 = range, v2 = slice" << std::endl;
676  retval = test_prod_rank1<NumericT>(epsilon,
677  ublas_m1, ublas_v1, ublas_v2,
678  vcl_m1_range, vcl_v1_range, vcl_v2_slice);
679  if (retval == EXIT_FAILURE)
680  {
681  std::cout << " --- FAILED! ---" << std::endl;
682  return retval;
683  }
684  else
685  std::cout << " --- PASSED ---" << std::endl;
686 
687 
688 
689  // v1 = slice
690 
691  std::cout << "* m = range, v1 = slice, v2 = full" << std::endl;
692  retval = test_prod_rank1<NumericT>(epsilon,
693  ublas_m1, ublas_v1, ublas_v2,
694  vcl_m1_range, vcl_v1_slice, vcl_v2_native);
695  if (retval == EXIT_FAILURE)
696  {
697  std::cout << " --- FAILED! ---" << std::endl;
698  return retval;
699  }
700  else
701  std::cout << " --- PASSED ---" << std::endl;
702 
703 
704  std::cout << "* m = range, v1 = slice, v2 = range" << std::endl;
705  retval = test_prod_rank1<NumericT>(epsilon,
706  ublas_m1, ublas_v1, ublas_v2,
707  vcl_m1_range, vcl_v1_slice, vcl_v2_range);
708  if (retval == EXIT_FAILURE)
709  {
710  std::cout << " --- FAILED! ---" << std::endl;
711  return retval;
712  }
713  else
714  std::cout << " --- PASSED ---" << std::endl;
715 
716 
717  std::cout << "* m = range, v1 = slice, v2 = slice" << std::endl;
718  retval = test_prod_rank1<NumericT>(epsilon,
719  ublas_m1, ublas_v1, ublas_v2,
720  vcl_m1_range, vcl_v1_slice, vcl_v2_slice);
721  if (retval == EXIT_FAILURE)
722  {
723  std::cout << " --- FAILED! ---" << std::endl;
724  return retval;
725  }
726  else
727  std::cout << " --- PASSED ---" << std::endl;
728 
729 
731 
732  std::cout << "* m = slice, v1 = full, v2 = full" << std::endl;
733  retval = test_prod_rank1<NumericT>(epsilon,
734  ublas_m1, ublas_v1, ublas_v2,
735  vcl_m1_slice, vcl_v1_native, vcl_v2_native);
736  if (retval == EXIT_FAILURE)
737  {
738  std::cout << " --- FAILED! ---" << std::endl;
739  return retval;
740  }
741  else
742  std::cout << " --- PASSED ---" << std::endl;
743 
744 
745  std::cout << "* m = slice, v1 = full, v2 = range" << std::endl;
746  retval = test_prod_rank1<NumericT>(epsilon,
747  ublas_m1, ublas_v1, ublas_v2,
748  vcl_m1_slice, vcl_v1_native, vcl_v2_range);
749  if (retval == EXIT_FAILURE)
750  {
751  std::cout << " --- FAILED! ---" << std::endl;
752  return retval;
753  }
754  else
755  std::cout << " --- PASSED ---" << std::endl;
756 
757 
758  std::cout << "* m = slice, v1 = full, v2 = slice" << std::endl;
759  retval = test_prod_rank1<NumericT>(epsilon,
760  ublas_m1, ublas_v1, ublas_v2,
761  vcl_m1_slice, vcl_v1_native, vcl_v2_slice);
762  if (retval == EXIT_FAILURE)
763  {
764  std::cout << " --- FAILED! ---" << std::endl;
765  return retval;
766  }
767  else
768  std::cout << " --- PASSED ---" << std::endl;
769 
770 
771  // v1 = range
772 
773 
774  std::cout << "* m = slice, v1 = range, v2 = full" << std::endl;
775  retval = test_prod_rank1<NumericT>(epsilon,
776  ublas_m1, ublas_v1, ublas_v2,
777  vcl_m1_slice, vcl_v1_range, vcl_v2_native);
778  if (retval == EXIT_FAILURE)
779  {
780  std::cout << " --- FAILED! ---" << std::endl;
781  return retval;
782  }
783  else
784  std::cout << " --- PASSED ---" << std::endl;
785 
786 
787  std::cout << "* m = slice, v1 = range, v2 = range" << std::endl;
788  retval = test_prod_rank1<NumericT>(epsilon,
789  ublas_m1, ublas_v1, ublas_v2,
790  vcl_m1_slice, vcl_v1_range, vcl_v2_range);
791  if (retval == EXIT_FAILURE)
792  {
793  std::cout << " --- FAILED! ---" << std::endl;
794  return retval;
795  }
796  else
797  std::cout << " --- PASSED ---" << std::endl;
798 
799 
800  std::cout << "* m = slice, v1 = range, v2 = slice" << std::endl;
801  retval = test_prod_rank1<NumericT>(epsilon,
802  ublas_m1, ublas_v1, ublas_v2,
803  vcl_m1_slice, vcl_v1_range, vcl_v2_slice);
804  if (retval == EXIT_FAILURE)
805  {
806  std::cout << " --- FAILED! ---" << std::endl;
807  return retval;
808  }
809  else
810  std::cout << " --- PASSED ---" << std::endl;
811 
812 
813 
814  // v1 = slice
815 
816  std::cout << "* m = slice, v1 = slice, v2 = full" << std::endl;
817  retval = test_prod_rank1<NumericT>(epsilon,
818  ublas_m1, ublas_v1, ublas_v2,
819  vcl_m1_slice, vcl_v1_slice, vcl_v2_native);
820  if (retval == EXIT_FAILURE)
821  {
822  std::cout << " --- FAILED! ---" << std::endl;
823  return retval;
824  }
825  else
826  std::cout << " --- PASSED ---" << std::endl;
827 
828 
829  std::cout << "* m = slice, v1 = slice, v2 = range" << std::endl;
830  retval = test_prod_rank1<NumericT>(epsilon,
831  ublas_m1, ublas_v1, ublas_v2,
832  vcl_m1_slice, vcl_v1_slice, vcl_v2_range);
833  if (retval == EXIT_FAILURE)
834  {
835  std::cout << " --- FAILED! ---" << std::endl;
836  return retval;
837  }
838  else
839  std::cout << " --- PASSED ---" << std::endl;
840 
841 
842  std::cout << "* m = slice, v1 = slice, v2 = slice" << std::endl;
843  retval = test_prod_rank1<NumericT>(epsilon,
844  ublas_m1, ublas_v1, ublas_v2,
845  vcl_m1_slice, vcl_v1_slice, vcl_v2_slice);
846  if (retval == EXIT_FAILURE)
847  {
848  std::cout << " --- FAILED! ---" << std::endl;
849  return retval;
850  }
851  else
852  std::cout << " --- PASSED ---" << std::endl;
853 
854  return retval;
855 }
856 //
857 // -------------------------------------------------------------
858 //
859 int main()
860 {
861  std::cout << std::endl;
862  std::cout << "----------------------------------------------" << std::endl;
863  std::cout << "----------------------------------------------" << std::endl;
864  std::cout << "## Test :: Matrix" << std::endl;
865  std::cout << "----------------------------------------------" << std::endl;
866  std::cout << "----------------------------------------------" << std::endl;
867  std::cout << std::endl;
868 
869  int retval = EXIT_SUCCESS;
870 
871  std::cout << std::endl;
872  std::cout << "----------------------------------------------" << std::endl;
873  std::cout << std::endl;
874  {
875  typedef float NumericT;
876  NumericT epsilon = NumericT(1.0E-3);
877  std::cout << "# Testing setup:" << std::endl;
878  std::cout << " eps: " << epsilon << std::endl;
879  std::cout << " numeric: float" << std::endl;
880  std::cout << " layout: row-major" << std::endl;
881  retval = test<NumericT, viennacl::row_major>(epsilon);
882  if ( retval == EXIT_SUCCESS )
883  std::cout << "# Test passed" << std::endl;
884  else
885  return retval;
886  }
887  std::cout << std::endl;
888  std::cout << "----------------------------------------------" << std::endl;
889  std::cout << std::endl;
890  {
891  typedef float NumericT;
892  NumericT epsilon = NumericT(1.0E-3);
893  std::cout << "# Testing setup:" << std::endl;
894  std::cout << " eps: " << epsilon << std::endl;
895  std::cout << " numeric: float" << std::endl;
896  std::cout << " layout: column-major" << std::endl;
897  retval = test<NumericT, viennacl::column_major>(epsilon);
898  if ( retval == EXIT_SUCCESS )
899  std::cout << "# Test passed" << std::endl;
900  else
901  return retval;
902  }
903  std::cout << std::endl;
904  std::cout << "----------------------------------------------" << std::endl;
905  std::cout << std::endl;
906 
907 
908 #ifdef VIENNACL_WITH_OPENCL
909  if ( viennacl::ocl::current_device().double_support() )
910 #endif
911  {
912  {
913  typedef double NumericT;
914  NumericT epsilon = 1.0E-11;
915  std::cout << "# Testing setup:" << std::endl;
916  std::cout << " eps: " << epsilon << std::endl;
917  std::cout << " numeric: double" << std::endl;
918  std::cout << " layout: row-major" << std::endl;
919  retval = test<NumericT, viennacl::row_major>(epsilon);
920  if ( retval == EXIT_SUCCESS )
921  std::cout << "# Test passed" << std::endl;
922  else
923  return retval;
924  }
925  std::cout << std::endl;
926  std::cout << "----------------------------------------------" << std::endl;
927  std::cout << std::endl;
928  {
929  typedef double NumericT;
930  NumericT epsilon = 1.0E-11;
931  std::cout << "# Testing setup:" << std::endl;
932  std::cout << " eps: " << epsilon << std::endl;
933  std::cout << " numeric: double" << std::endl;
934  std::cout << " layout: column-major" << std::endl;
935  retval = test<NumericT, viennacl::column_major>(epsilon);
936  if ( retval == EXIT_SUCCESS )
937  std::cout << "# Test passed" << std::endl;
938  else
939  return retval;
940  }
941  std::cout << std::endl;
942  std::cout << "----------------------------------------------" << std::endl;
943  std::cout << std::endl;
944  }
945 
946  std::cout << std::endl;
947  std::cout << "------- Test completed --------" << std::endl;
948  std::cout << std::endl;
949 
950 
951  return retval;
952 }
viennacl::scalar
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
norm_2.hpp
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: scheduler_matrix_vector.cpp:65
viennacl::matrix_slice
Class for representing strided submatrices of a bigger matrix A.
Definition: forwards.h:443
viennacl::linalg::cuda::trans
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)
Definition: matrix_operations.hpp:94
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
matrix.hpp
Implementation of the dense matrix class.
io.hpp
Some helper routines for reading/writing/printing scheduler expressions.
viennacl::op_assign
A tag class representing assignment.
Definition: forwards.h:81
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::scheduler::execute
void execute(statement const &s)
Definition: execute.hpp:279
test
int test(Epsilon const &epsilon)
Definition: scheduler_matrix_vector.cpp:381
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
s1
viennacl::scalar< float > s1
Definition: global_variables.cpp:57
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::current_device
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
viennacl::op_inplace_add
A tag class representing inplace addition.
Definition: forwards.h:83
NumericT
float NumericT
Definition: bisect.cpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
main
int main()
Definition: scheduler_matrix_vector.cpp:859
viennacl::vector_range
Class for representing non-strided subvectors of a bigger vector x.
Definition: forwards.h:434
viennacl::vector_slice
Class for representing strided subvectors of a bigger vector x.
Definition: forwards.h:437
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
lu.hpp
Implementations of LU factorization for row-major and column-major dense matrices.
direct_solve.hpp
Implementations of dense direct solvers are found here.
viennacl::op_inplace_sub
A tag class representing inplace subtraction.
Definition: forwards.h:85
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::linalg::norm_inf
T norm_inf(std::vector< T, A > const &v1)
Definition: norm_inf.hpp:60
viennacl::copy
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) ...
Definition: circulant_matrix.hpp:150
random.hpp
A small collection of sequential random number generators.
test_prod_rank1
int test_prod_rank1(Epsilon const &epsilon, UblasMatrixType &ublas_m1, UblasVectorType &ublas_v1, UblasVectorType &ublas_v2, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1, VCLVectorType2 &vcl_v2)
Definition: scheduler_matrix_vector.cpp:119
viennacl::basic_range
A range class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:424
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
execute.hpp
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
viennacl::scheduler::statement
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
Definition: forwards.h:502
viennacl::basic_slice
A slice class that refers to an interval [start, stop), where 'start' is included, and 'stop' is excluded.
Definition: forwards.h:429
scalar.hpp
Implementation of the ViennaCL scalar class.