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