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