ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
matrix_vector_int.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 
23 //
24 // *** System
25 //
26 #include <iostream>
27 
28 //
29 // *** Boost
30 //
31 #include <boost/numeric/ublas/io.hpp>
32 #include <boost/numeric/ublas/triangular.hpp>
33 #include <boost/numeric/ublas/matrix_sparse.hpp>
34 #include <boost/numeric/ublas/matrix.hpp>
35 #include <boost/numeric/ublas/matrix_proxy.hpp>
36 #include <boost/numeric/ublas/lu.hpp>
37 #include <boost/numeric/ublas/io.hpp>
38 
39 //
40 // *** ViennaCL
41 //
42 //#define VIENNACL_DEBUG_ALL
43 #define VIENNACL_WITH_UBLAS 1
44 #include "viennacl/scalar.hpp"
45 #include "viennacl/matrix.hpp"
46 #include "viennacl/vector.hpp"
47 #include "viennacl/linalg/prod.hpp"
48 #include "viennacl/linalg/norm_2.hpp"
49 #include "viennacl/linalg/direct_solve.hpp"
50 #include "viennacl/linalg/lu.hpp"
51 #include "viennacl/linalg/sum.hpp"
52 
53 //
54 // -------------------------------------------------------------
55 //
56 using namespace boost::numeric;
57 //
58 // -------------------------------------------------------------
59 //
60 template<typename ScalarType>
61 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
62 {
63  viennacl::backend::finish();
64  if (s1 != s2)
65  return 1;
66  return 0;
67 }
68 
69 template<typename ScalarType, typename VCLVectorType>
70 ScalarType diff(ublas::vector<ScalarType> const & v1, VCLVectorType const & v2)
71 {
72  ublas::vector<ScalarType> v2_cpu(v2.size());
73  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
74  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
75 
76  for (unsigned int i=0;i<v1.size(); ++i)
77  {
78  if (v2_cpu[i] != v1[i])
79  return 1;
80  }
81 
82  return 0;
83 }
84 
85 template<typename ScalarType, typename VCLMatrixType>
86 ScalarType diff(ublas::matrix<ScalarType> const & mat1, VCLMatrixType const & mat2)
87 {
88  ublas::matrix<ScalarType> mat2_cpu(mat2.size1(), mat2.size2());
89  viennacl::backend::finish(); //workaround for a bug in APP SDK 2.7 on Trinity APUs (with Catalyst 12.8)
90  viennacl::copy(mat2, mat2_cpu);
91 
92  for (unsigned int i = 0; i < mat2_cpu.size1(); ++i)
93  {
94  for (unsigned int j = 0; j < mat2_cpu.size2(); ++j)
95  {
96  if (mat2_cpu(i,j) != mat1(i,j))
97  return 1;
98  }
99  }
100  //std::cout << ret << std::endl;
101  return 0;
102 }
103 //
104 // -------------------------------------------------------------
105 //
106 
107 template<typename NumericT,
108  typename UblasMatrixType, typename UblasVectorType,
109  typename VCLMatrixType, typename VCLVectorType1, typename VCLVectorType2>
110 int test_prod_rank1(UblasMatrixType & ublas_m1, UblasVectorType & ublas_v1, UblasVectorType & ublas_v2,
111  VCLMatrixType & vcl_m1, VCLVectorType1 & vcl_v1, VCLVectorType2 & vcl_v2)
112 {
113  int retval = EXIT_SUCCESS;
114 
115  // sync data:
116  ublas_v1 = ublas::scalar_vector<NumericT>(ublas_v1.size(), NumericT(2));
117  ublas_v2 = ublas::scalar_vector<NumericT>(ublas_v2.size(), NumericT(3));
118  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
119  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
120  viennacl::copy(ublas_m1, vcl_m1);
121 
122  // --------------------------------------------------------------------------
123  std::cout << "Rank 1 update" << std::endl;
124 
125  ublas_m1 += ublas::outer_prod(ublas_v1, ublas_v2);
126  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
127  if ( diff(ublas_m1, vcl_m1) != 0 )
128  {
129  std::cout << "# Error at operation: rank 1 update" << std::endl;
130  std::cout << " diff: " << diff(ublas_m1, vcl_m1) << std::endl;
131  return EXIT_FAILURE;
132  }
133 
134 
135 
136  // --------------------------------------------------------------------------
137  std::cout << "Scaled rank 1 update - CPU Scalar" << std::endl;
138  ublas_m1 += NumericT(4) * ublas::outer_prod(ublas_v1, ublas_v2);
139  vcl_m1 += NumericT(2) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
140  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * NumericT(2); //check proper compilation
141  if ( diff(ublas_m1, vcl_m1) != 0 )
142  {
143  std::cout << "# Error at operation: scaled rank 1 update - CPU Scalar" << std::endl;
144  std::cout << " diff: " << diff(ublas_m1, vcl_m1) << std::endl;
145  return EXIT_FAILURE;
146  }
147 
148  // --------------------------------------------------------------------------
149  std::cout << "Scaled rank 1 update - GPU Scalar" << std::endl;
150  ublas_m1 += NumericT(4) * ublas::outer_prod(ublas_v1, ublas_v2);
151  vcl_m1 += viennacl::scalar<NumericT>(2) * viennacl::linalg::outer_prod(vcl_v1, vcl_v2);
152  vcl_m1 += viennacl::linalg::outer_prod(vcl_v1, vcl_v2) * viennacl::scalar<NumericT>(2); //check proper compilation
153  if ( diff(ublas_m1, vcl_m1) != 0 )
154  {
155  std::cout << "# Error at operation: scaled rank 1 update - GPU Scalar" << std::endl;
156  std::cout << " diff: " << diff(ublas_m1, vcl_m1) << std::endl;
157  return EXIT_FAILURE;
158  }
159 
160  //reset vcl_matrix:
161  viennacl::copy(ublas_m1, vcl_m1);
162 
163  // --------------------------------------------------------------------------
164  std::cout << "Matrix-Vector product" << std::endl;
165  ublas_v1 = viennacl::linalg::prod(ublas_m1, ublas_v2);
166  vcl_v1 = viennacl::linalg::prod(vcl_m1, vcl_v2);
167 
168  if ( diff(ublas_v1, vcl_v1) != 0 )
169  {
170  std::cout << "# Error at operation: matrix-vector product" << std::endl;
171  std::cout << " diff: " << diff(ublas_v1, vcl_v1) << std::endl;
172  retval = EXIT_FAILURE;
173  }
174  // --------------------------------------------------------------------------
175  std::cout << "Matrix-Vector product with scaled add" << std::endl;
176  NumericT alpha = static_cast<NumericT>(2);
177  NumericT beta = static_cast<NumericT>(3);
178  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
179  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
180 
181  ublas_v1 = alpha * viennacl::linalg::prod(ublas_m1, ublas_v2) + beta * ublas_v1;
182  vcl_v1 = alpha * viennacl::linalg::prod(vcl_m1, vcl_v2) + beta * vcl_v1;
183 
184  if ( diff(ublas_v1, vcl_v1) != 0 )
185  {
186  std::cout << "# Error at operation: matrix-vector product with scaled additions" << std::endl;
187  std::cout << " diff: " << diff(ublas_v1, vcl_v1) << std::endl;
188  retval = EXIT_FAILURE;
189  }
190  // --------------------------------------------------------------------------
191 
192  viennacl::copy(ublas_v1.begin(), ublas_v1.end(), vcl_v1.begin());
193  viennacl::copy(ublas_v2.begin(), ublas_v2.end(), vcl_v2.begin());
194 
195  std::cout << "Transposed Matrix-Vector product" << std::endl;
196  ublas_v2 = alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1);
197  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1);
198 
199  if ( diff(ublas_v2, vcl_v2) != 0 )
200  {
201  std::cout << "# Error at operation: transposed matrix-vector product" << std::endl;
202  std::cout << " diff: " << diff(ublas_v2, vcl_v2) << std::endl;
203  retval = EXIT_FAILURE;
204  }
205 
206  std::cout << "Transposed Matrix-Vector product with scaled add" << std::endl;
207  ublas_v2 = alpha * viennacl::linalg::prod(trans(ublas_m1), ublas_v1) + beta * ublas_v2;
208  vcl_v2 = alpha * viennacl::linalg::prod(trans(vcl_m1), vcl_v1) + beta * vcl_v2;
209 
210  if ( diff(ublas_v2, vcl_v2) != 0 )
211  {
212  std::cout << "# Error at operation: transposed matrix-vector product with scaled additions" << std::endl;
213  std::cout << " diff: " << diff(ublas_v2, vcl_v2) << std::endl;
214  retval = EXIT_FAILURE;
215  }
216  // --------------------------------------------------------------------------
217 
218 
219  std::cout << "Row sum with matrix" << std::endl;
220  ublas_v1 = ublas::prod(ublas_m1, ublas::scalar_vector<NumericT>(ublas_m1.size2(), NumericT(1)));
221  vcl_v1 = viennacl::linalg::row_sum(vcl_m1);
222 
223  if ( diff(ublas_v1, vcl_v1) != 0 )
224  {
225  std::cout << "# Error at operation: row sum" << std::endl;
226  std::cout << " diff: " << diff(ublas_v1, vcl_v1) << std::endl;
227  retval = EXIT_FAILURE;
228  }
229  // --------------------------------------------------------------------------
230 
231  std::cout << "Row sum with matrix expression" << std::endl;
232  ublas_v1 = ublas::prod(ublas_m1 + ublas_m1, ublas::scalar_vector<NumericT>(ublas_m1.size2(), NumericT(1)));
233  vcl_v1 = viennacl::linalg::row_sum(vcl_m1 + vcl_m1);
234 
235  if ( diff(ublas_v1, vcl_v1) != 0 )
236  {
237  std::cout << "# Error at operation: row sum (with expression)" << std::endl;
238  std::cout << " diff: " << diff(ublas_v1, vcl_v1) << std::endl;
239  retval = EXIT_FAILURE;
240  }
241  // --------------------------------------------------------------------------
242 
243  std::cout << "Column sum with matrix" << std::endl;
244  ublas_v2 = ublas::prod(trans(ublas_m1), ublas::scalar_vector<NumericT>(ublas_m1.size1(), NumericT(1)));
245  vcl_v2 = viennacl::linalg::column_sum(vcl_m1);
246 
247  if ( diff(ublas_v2, vcl_v2) != 0 )
248  {
249  std::cout << "# Error at operation: column sum" << std::endl;
250  std::cout << " diff: " << diff(ublas_v2, vcl_v2) << std::endl;
251  retval = EXIT_FAILURE;
252  }
253  // --------------------------------------------------------------------------
254 
255  std::cout << "Column sum with matrix expression" << std::endl;
256  ublas_v2 = ublas::prod(trans(ublas_m1 + ublas_m1), ublas::scalar_vector<NumericT>(ublas_m1.size1(), NumericT(1)));
257  vcl_v2 = viennacl::linalg::column_sum(vcl_m1 + vcl_m1);
258 
259  if ( diff(ublas_v2, vcl_v2) != 0 )
260  {
261  std::cout << "# Error at operation: column sum (with expression)" << std::endl;
262  std::cout << " diff: " << diff(ublas_v2, vcl_v2) << std::endl;
263  retval = EXIT_FAILURE;
264  }
265  // --------------------------------------------------------------------------
266 
267  return retval;
268 }
269 
270 
271 //
272 // -------------------------------------------------------------
273 //
274 template< typename NumericT, typename F>
275 int test()
276 {
277  int retval = EXIT_SUCCESS;
278 
279  std::size_t num_rows = 141;
280  std::size_t num_cols = 103;
281 
282  // --------------------------------------------------------------------------
283  ublas::vector<NumericT> ublas_v1(num_rows);
284  for (std::size_t i = 0; i < ublas_v1.size(); ++i)
285  ublas_v1(i) = NumericT(i);
286  ublas::vector<NumericT> ublas_v2 = ublas::scalar_vector<NumericT>(num_cols, NumericT(3));
287 
288 
289  ublas::matrix<NumericT> ublas_m1(ublas_v1.size(), ublas_v2.size());
290  ublas::matrix<NumericT> ublas_m2(ublas_v1.size(), ublas_v1.size());
291 
292 
293  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
294  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
295  ublas_m1(i,j) = NumericT(i+j);
296 
297 
298  for (std::size_t i = 0; i < ublas_m2.size1(); ++i)
299  for (std::size_t j = 0; j < ublas_m2.size2(); ++j)
300  ublas_m2(i,j) = NumericT(j - i*j + i);
301 
302 
303  viennacl::vector<NumericT> vcl_v1_native(ublas_v1.size());
304  viennacl::vector<NumericT> vcl_v1_large(4 * ublas_v1.size());
305  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v1_range(vcl_v1_large, viennacl::range(3, ublas_v1.size() + 3));
306  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v1_slice(vcl_v1_large, viennacl::slice(2, 3, ublas_v1.size()));
307 
308  viennacl::vector<NumericT> vcl_v2_native(ublas_v2.size());
309  viennacl::vector<NumericT> vcl_v2_large(4 * ublas_v2.size());
310  viennacl::vector_range< viennacl::vector<NumericT> > vcl_v2_range(vcl_v2_large, viennacl::range(8, ublas_v2.size() + 8));
311  viennacl::vector_slice< viennacl::vector<NumericT> > vcl_v2_slice(vcl_v2_large, viennacl::slice(6, 2, ublas_v2.size()));
312 
313  viennacl::matrix<NumericT, F> vcl_m1_native(ublas_m1.size1(), ublas_m1.size2());
314  viennacl::matrix<NumericT, F> vcl_m1_large(4 * ublas_m1.size1(), 4 * ublas_m1.size2());
315  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m1_range(vcl_m1_large,
316  viennacl::range(8, ublas_m1.size1() + 8),
317  viennacl::range(ublas_m1.size2(), 2 * ublas_m1.size2()) );
318  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m1_slice(vcl_m1_large,
319  viennacl::slice(6, 2, ublas_m1.size1()),
320  viennacl::slice(ublas_m1.size2(), 2, ublas_m1.size2()) );
321 
322  viennacl::matrix<NumericT, F> vcl_m2_native(ublas_m2.size1(), ublas_m2.size2());
323  viennacl::matrix<NumericT, F> vcl_m2_large(4 * ublas_m2.size1(), 4 * ublas_m2.size2());
324  viennacl::matrix_range< viennacl::matrix<NumericT, F> > vcl_m2_range(vcl_m2_large,
325  viennacl::range(8, ublas_m2.size1() + 8),
326  viennacl::range(ublas_m2.size2(), 2 * ublas_m2.size2()) );
327  viennacl::matrix_slice< viennacl::matrix<NumericT, F> > vcl_m2_slice(vcl_m2_large,
328  viennacl::slice(6, 2, ublas_m2.size1()),
329  viennacl::slice(ublas_m2.size2(), 2, ublas_m2.size2()) );
330 
331 
332  //
333  // Run a bunch of tests for rank-1-updates, matrix-vector products
334  //
335  std::cout << "------------ Testing rank-1-updates and matrix-vector products ------------------" << std::endl;
336 
337  std::cout << "* m = full, v1 = full, v2 = full" << std::endl;
338  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
339  vcl_m1_native, vcl_v1_native, vcl_v2_native);
340  if (retval == EXIT_FAILURE)
341  {
342  std::cout << " --- FAILED! ---" << std::endl;
343  return retval;
344  }
345  else
346  std::cout << " --- PASSED ---" << std::endl;
347 
348  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
349  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
350  ublas_m1(i,j) = NumericT(i+j);
351 
352  std::cout << "* m = full, v1 = full, v2 = range" << std::endl;
353  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
354  vcl_m1_native, vcl_v1_native, vcl_v2_range);
355  if (retval == EXIT_FAILURE)
356  {
357  std::cout << " --- FAILED! ---" << std::endl;
358  return retval;
359  }
360  else
361  std::cout << " --- PASSED ---" << std::endl;
362 
363 
364  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
365  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
366  ublas_m1(i,j) = NumericT(i+j);
367 
368  std::cout << "* m = full, v1 = full, v2 = slice" << std::endl;
369  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
370  vcl_m1_native, vcl_v1_native, vcl_v2_slice);
371  if (retval == EXIT_FAILURE)
372  {
373  std::cout << " --- FAILED! ---" << std::endl;
374  return retval;
375  }
376  else
377  std::cout << " --- PASSED ---" << std::endl;
378 
379 
380  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
381  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
382  ublas_m1(i,j) = NumericT(i+j);
383 
384  // v1 = range
385 
386 
387  std::cout << "* m = full, v1 = range, v2 = full" << std::endl;
388  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
389  vcl_m1_native, vcl_v1_range, vcl_v2_native);
390  if (retval == EXIT_FAILURE)
391  {
392  std::cout << " --- FAILED! ---" << std::endl;
393  return retval;
394  }
395  else
396  std::cout << " --- PASSED ---" << std::endl;
397 
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) = NumericT(i+j);
402 
403  std::cout << "* m = full, v1 = range, v2 = range" << std::endl;
404  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
405  vcl_m1_native, vcl_v1_range, vcl_v2_range);
406  if (retval == EXIT_FAILURE)
407  {
408  std::cout << " --- FAILED! ---" << std::endl;
409  return retval;
410  }
411  else
412  std::cout << " --- PASSED ---" << std::endl;
413 
414 
415  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
416  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
417  ublas_m1(i,j) = NumericT(i+j);
418 
419  std::cout << "* m = full, v1 = range, v2 = slice" << std::endl;
420  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
421  vcl_m1_native, vcl_v1_range, vcl_v2_slice);
422  if (retval == EXIT_FAILURE)
423  {
424  std::cout << " --- FAILED! ---" << std::endl;
425  return retval;
426  }
427  else
428  std::cout << " --- PASSED ---" << std::endl;
429 
430 
431  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
432  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
433  ublas_m1(i,j) = NumericT(i+j);
434 
435 
436  // v1 = slice
437 
438  std::cout << "* m = full, v1 = slice, v2 = full" << std::endl;
439  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
440  vcl_m1_native, vcl_v1_slice, vcl_v2_native);
441  if (retval == EXIT_FAILURE)
442  {
443  std::cout << " --- FAILED! ---" << std::endl;
444  return retval;
445  }
446  else
447  std::cout << " --- PASSED ---" << std::endl;
448 
449 
450  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
451  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
452  ublas_m1(i,j) = NumericT(i+j);
453 
454  std::cout << "* m = full, v1 = slice, v2 = range" << std::endl;
455  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
456  vcl_m1_native, vcl_v1_slice, vcl_v2_range);
457  if (retval == EXIT_FAILURE)
458  {
459  std::cout << " --- FAILED! ---" << std::endl;
460  return retval;
461  }
462  else
463  std::cout << " --- PASSED ---" << std::endl;
464 
465 
466  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
467  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
468  ublas_m1(i,j) = NumericT(i+j);
469 
470  std::cout << "* m = full, v1 = slice, v2 = slice" << std::endl;
471  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
472  vcl_m1_native, vcl_v1_slice, vcl_v2_slice);
473  if (retval == EXIT_FAILURE)
474  {
475  std::cout << " --- FAILED! ---" << std::endl;
476  return retval;
477  }
478  else
479  std::cout << " --- PASSED ---" << std::endl;
480 
481 
482  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
483  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
484  ublas_m1(i,j) = NumericT(i+j);
485 
487 
488  std::cout << "* m = range, v1 = full, v2 = full" << std::endl;
489  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
490  vcl_m1_range, vcl_v1_native, vcl_v2_native);
491  if (retval == EXIT_FAILURE)
492  {
493  std::cout << " --- FAILED! ---" << std::endl;
494  return retval;
495  }
496  else
497  std::cout << " --- PASSED ---" << std::endl;
498 
499 
500  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
501  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
502  ublas_m1(i,j) = NumericT(i+j);
503 
504  std::cout << "* m = range, v1 = full, v2 = range" << std::endl;
505  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
506  vcl_m1_range, vcl_v1_native, vcl_v2_range);
507  if (retval == EXIT_FAILURE)
508  {
509  std::cout << " --- FAILED! ---" << std::endl;
510  return retval;
511  }
512  else
513  std::cout << " --- PASSED ---" << std::endl;
514 
515 
516  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
517  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
518  ublas_m1(i,j) = NumericT(i+j);
519 
520  std::cout << "* m = range, v1 = full, v2 = slice" << std::endl;
521  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
522  vcl_m1_range, vcl_v1_native, vcl_v2_slice);
523  if (retval == EXIT_FAILURE)
524  {
525  std::cout << " --- FAILED! ---" << std::endl;
526  return retval;
527  }
528  else
529  std::cout << " --- PASSED ---" << std::endl;
530 
531 
532  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
533  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
534  ublas_m1(i,j) = NumericT(i+j);
535 
536  // v1 = range
537 
538 
539  std::cout << "* m = range, v1 = range, v2 = full" << std::endl;
540  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
541  vcl_m1_range, vcl_v1_range, vcl_v2_native);
542  if (retval == EXIT_FAILURE)
543  {
544  std::cout << " --- FAILED! ---" << std::endl;
545  return retval;
546  }
547  else
548  std::cout << " --- PASSED ---" << std::endl;
549 
550 
551  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
552  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
553  ublas_m1(i,j) = NumericT(i+j);
554 
555  std::cout << "* m = range, v1 = range, v2 = range" << std::endl;
556  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
557  vcl_m1_range, vcl_v1_range, vcl_v2_range);
558  if (retval == EXIT_FAILURE)
559  {
560  std::cout << " --- FAILED! ---" << std::endl;
561  return retval;
562  }
563  else
564  std::cout << " --- PASSED ---" << std::endl;
565 
566 
567  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
568  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
569  ublas_m1(i,j) = NumericT(i+j);
570 
571  std::cout << "* m = range, v1 = range, v2 = slice" << std::endl;
572  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
573  vcl_m1_range, vcl_v1_range, vcl_v2_slice);
574  if (retval == EXIT_FAILURE)
575  {
576  std::cout << " --- FAILED! ---" << std::endl;
577  return retval;
578  }
579  else
580  std::cout << " --- PASSED ---" << std::endl;
581 
582 
583  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
584  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
585  ublas_m1(i,j) = NumericT(i+j);
586 
587 
588  // v1 = slice
589 
590  std::cout << "* m = range, v1 = slice, v2 = full" << std::endl;
591  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
592  vcl_m1_range, vcl_v1_slice, vcl_v2_native);
593  if (retval == EXIT_FAILURE)
594  {
595  std::cout << " --- FAILED! ---" << std::endl;
596  return retval;
597  }
598  else
599  std::cout << " --- PASSED ---" << std::endl;
600 
601 
602  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
603  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
604  ublas_m1(i,j) = NumericT(i+j);
605 
606  std::cout << "* m = range, v1 = slice, v2 = range" << std::endl;
607  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
608  vcl_m1_range, vcl_v1_slice, vcl_v2_range);
609  if (retval == EXIT_FAILURE)
610  {
611  std::cout << " --- FAILED! ---" << std::endl;
612  return retval;
613  }
614  else
615  std::cout << " --- PASSED ---" << std::endl;
616 
617 
618  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
619  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
620  ublas_m1(i,j) = NumericT(i+j);
621 
622  std::cout << "* m = range, v1 = slice, v2 = slice" << std::endl;
623  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
624  vcl_m1_range, vcl_v1_slice, vcl_v2_slice);
625  if (retval == EXIT_FAILURE)
626  {
627  std::cout << " --- FAILED! ---" << std::endl;
628  return retval;
629  }
630  else
631  std::cout << " --- PASSED ---" << std::endl;
632 
633 
634  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
635  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
636  ublas_m1(i,j) = NumericT(i+j);
637 
639 
640  std::cout << "* m = slice, v1 = full, v2 = full" << std::endl;
641  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
642  vcl_m1_slice, vcl_v1_native, vcl_v2_native);
643  if (retval == EXIT_FAILURE)
644  {
645  std::cout << " --- FAILED! ---" << std::endl;
646  return retval;
647  }
648  else
649  std::cout << " --- PASSED ---" << std::endl;
650 
651 
652  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
653  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
654  ublas_m1(i,j) = NumericT(i+j);
655 
656  std::cout << "* m = slice, v1 = full, v2 = range" << std::endl;
657  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
658  vcl_m1_slice, vcl_v1_native, vcl_v2_range);
659  if (retval == EXIT_FAILURE)
660  {
661  std::cout << " --- FAILED! ---" << std::endl;
662  return retval;
663  }
664  else
665  std::cout << " --- PASSED ---" << std::endl;
666 
667 
668  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
669  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
670  ublas_m1(i,j) = NumericT(i+j);
671 
672  std::cout << "* m = slice, v1 = full, v2 = slice" << std::endl;
673  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
674  vcl_m1_slice, vcl_v1_native, vcl_v2_slice);
675  if (retval == EXIT_FAILURE)
676  {
677  std::cout << " --- FAILED! ---" << std::endl;
678  return retval;
679  }
680  else
681  std::cout << " --- PASSED ---" << std::endl;
682 
683 
684  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
685  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
686  ublas_m1(i,j) = NumericT(i+j);
687 
688  // v1 = range
689 
690 
691  std::cout << "* m = slice, v1 = range, v2 = full" << std::endl;
692  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
693  vcl_m1_slice, vcl_v1_range, vcl_v2_native);
694  if (retval == EXIT_FAILURE)
695  {
696  std::cout << " --- FAILED! ---" << std::endl;
697  return retval;
698  }
699  else
700  std::cout << " --- PASSED ---" << std::endl;
701 
702 
703  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
704  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
705  ublas_m1(i,j) = NumericT(i+j);
706 
707  std::cout << "* m = slice, v1 = range, v2 = range" << std::endl;
708  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
709  vcl_m1_slice, vcl_v1_range, vcl_v2_range);
710  if (retval == EXIT_FAILURE)
711  {
712  std::cout << " --- FAILED! ---" << std::endl;
713  return retval;
714  }
715  else
716  std::cout << " --- PASSED ---" << std::endl;
717 
718 
719  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
720  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
721  ublas_m1(i,j) = NumericT(i+j);
722 
723  std::cout << "* m = slice, v1 = range, v2 = slice" << std::endl;
724  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
725  vcl_m1_slice, vcl_v1_range, vcl_v2_slice);
726  if (retval == EXIT_FAILURE)
727  {
728  std::cout << " --- FAILED! ---" << std::endl;
729  return retval;
730  }
731  else
732  std::cout << " --- PASSED ---" << std::endl;
733 
734 
735 
736  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
737  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
738  ublas_m1(i,j) = NumericT(i+j);
739 
740  // v1 = slice
741 
742  std::cout << "* m = slice, v1 = slice, v2 = full" << std::endl;
743  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
744  vcl_m1_slice, 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  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
754  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
755  ublas_m1(i,j) = NumericT(i+j);
756 
757 
758  std::cout << "* m = slice, v1 = slice, v2 = range" << std::endl;
759  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
760  vcl_m1_slice, vcl_v1_slice, vcl_v2_range);
761  if (retval == EXIT_FAILURE)
762  {
763  std::cout << " --- FAILED! ---" << std::endl;
764  return retval;
765  }
766  else
767  std::cout << " --- PASSED ---" << std::endl;
768 
769 
770  for (std::size_t i = 0; i < ublas_m1.size1(); ++i)
771  for (std::size_t j = 0; j < ublas_m1.size2(); ++j)
772  ublas_m1(i,j) = NumericT(i+j);
773 
774  std::cout << "* m = slice, v1 = slice, v2 = slice" << std::endl;
775  retval = test_prod_rank1<NumericT>(ublas_m1, ublas_v1, ublas_v2,
776  vcl_m1_slice, vcl_v1_slice, vcl_v2_slice);
777  if (retval == EXIT_FAILURE)
778  {
779  std::cout << " --- FAILED! ---" << std::endl;
780  return retval;
781  }
782  else
783  std::cout << " --- PASSED ---" << std::endl;
784 
785 
786  return retval;
787 }
788 //
789 // -------------------------------------------------------------
790 //
791 int main()
792 {
793  std::cout << std::endl;
794  std::cout << "----------------------------------------------" << std::endl;
795  std::cout << "----------------------------------------------" << std::endl;
796  std::cout << "## Test :: Matrix" << std::endl;
797  std::cout << "----------------------------------------------" << std::endl;
798  std::cout << "----------------------------------------------" << std::endl;
799  std::cout << std::endl;
800 
801  int retval = EXIT_SUCCESS;
802 
803  std::cout << std::endl;
804  std::cout << "----------------------------------------------" << std::endl;
805  std::cout << std::endl;
806  {
807  typedef int NumericT;
808  std::cout << "# Testing setup:" << std::endl;
809  std::cout << " numeric: int" << std::endl;
810  std::cout << " layout: row-major" << std::endl;
811  retval = test<NumericT, viennacl::row_major>();
812  if ( retval == EXIT_SUCCESS )
813  std::cout << "# Test passed" << std::endl;
814  else
815  return retval;
816  }
817  std::cout << std::endl;
818  std::cout << "----------------------------------------------" << std::endl;
819  std::cout << std::endl;
820  {
821  typedef int NumericT;
822  std::cout << "# Testing setup:" << std::endl;
823  std::cout << " numeric: int" << std::endl;
824  std::cout << " layout: column-major" << std::endl;
825  retval = test<NumericT, viennacl::column_major>();
826  if ( retval == EXIT_SUCCESS )
827  std::cout << "# Test passed" << std::endl;
828  else
829  return retval;
830  }
831  std::cout << std::endl;
832  std::cout << "----------------------------------------------" << std::endl;
833  std::cout << std::endl;
834 
835 
836 #ifdef VIENNACL_WITH_OPENCL
837  if ( viennacl::ocl::current_device().double_support() )
838 #endif
839  {
840  {
841  typedef long NumericT;
842  std::cout << "# Testing setup:" << std::endl;
843  std::cout << " numeric: double" << std::endl;
844  std::cout << " layout: row-major" << std::endl;
845  retval = test<NumericT, viennacl::row_major>();
846  if ( retval == EXIT_SUCCESS )
847  std::cout << "# Test passed" << std::endl;
848  else
849  return retval;
850  }
851  std::cout << std::endl;
852  std::cout << "----------------------------------------------" << std::endl;
853  std::cout << std::endl;
854  {
855  typedef long NumericT;
856  std::cout << "# Testing setup:" << std::endl;
857  std::cout << " numeric: double" << std::endl;
858  std::cout << " layout: column-major" << std::endl;
859  retval = test<NumericT, viennacl::column_major>();
860  if ( retval == EXIT_SUCCESS )
861  std::cout << "# Test passed" << std::endl;
862  else
863  return retval;
864  }
865  std::cout << std::endl;
866  std::cout << "----------------------------------------------" << std::endl;
867  std::cout << std::endl;
868  }
869 
870  std::cout << std::endl;
871  std::cout << "------- Test completed --------" << std::endl;
872  std::cout << std::endl;
873 
874 
875  return retval;
876 }
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...
test_prod_rank1
int test_prod_rank1(UblasMatrixType &ublas_m1, UblasVectorType &ublas_v1, UblasVectorType &ublas_v2, VCLMatrixType &vcl_m1, VCLVectorType1 &vcl_v1, VCLVectorType2 &vcl_v2)
Definition: matrix_vector_int.cpp:110
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.
viennacl::linalg::row_sum
viennacl::vector_expression< const viennacl::matrix_base< NumericT >, const viennacl::matrix_base< NumericT >, viennacl::op_row_sum > row_sum(viennacl::matrix_base< NumericT > const &A)
User interface function for computing the sum of all elements of each row of a matrix.
Definition: sum.hpp:77
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
s2
viennacl::scalar< int > s2
Definition: global_variables.cpp:58
s1
viennacl::scalar< float > s1
Definition: global_variables.cpp:57
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
test
int test()
Definition: matrix_vector_int.cpp:275
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::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
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: matrix_vector_int.cpp:61
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::linalg::outer_prod
viennacl::matrix_expression< const vector_base< NumericT >, const vector_base< NumericT >, op_prod > outer_prod(const vector_base< NumericT > &vec1, const vector_base< NumericT > &vec2)
Returns a proxy class for the operation mat += vec1 * vec2^T, i.e. a rank 1 update.
Definition: matrix_operations.hpp:840
sum.hpp
Stub routines for the summation of elements in a vector, or all elements in either a row or column of...
v2
viennacl::vector< int > v2
Definition: global_variables.cpp:61
prod
void prod(std::vector< std::map< IndexT, NumericT > > const &stl_A, std::vector< std::map< IndexT, NumericT > > const &stl_B, std::vector< std::map< IndexT, NumericT > > &stl_C)
Definition: sparse_prod.cpp:114
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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
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
main
int main()
Definition: matrix_vector_int.cpp:791
viennacl::matrix_range
Class for representing non-strided submatrices of a bigger matrix A.
Definition: forwards.h:440
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
viennacl::linalg::column_sum
viennacl::vector_expression< const viennacl::matrix_base< NumericT >, const viennacl::matrix_base< NumericT >, viennacl::op_col_sum > column_sum(viennacl::matrix_base< NumericT > const &A)
User interface function for computing the sum of all elements of each column of a matrix...
Definition: sum.hpp:109
scalar.hpp
Implementation of the ViennaCL scalar class.