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
sparse.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 #include <map>
30 #include <cmath>
31 
32 //
33 // *** ViennaCL
34 //
35 #include "viennacl/scalar.hpp"
36 #include "viennacl/compressed_matrix.hpp"
37 #include "viennacl/compressed_compressed_matrix.hpp"
38 #include "viennacl/coordinate_matrix.hpp"
39 #include "viennacl/ell_matrix.hpp"
40 #include "viennacl/sliced_ell_matrix.hpp"
41 #include "viennacl/hyb_matrix.hpp"
42 #include "viennacl/vector.hpp"
43 #include "viennacl/vector_proxy.hpp"
44 #include "viennacl/linalg/prod.hpp"
45 #include "viennacl/linalg/norm_2.hpp"
46 #include "viennacl/linalg/ilu.hpp"
47 #include "viennacl/linalg/detail/ilu/common.hpp"
48 #include "viennacl/io/matrix_market.hpp"
49 #include "viennacl/tools/random.hpp"
50 
51 
52 
53 //
54 // -------------------------------------------------------------
55 //
56 template<typename ScalarType>
57 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
58 {
59  if (s1 != s2)
60  return (s1 - s2) / std::max(fabs(s1), std::fabs(s2));
61  return 0;
62 }
63 
64 template<typename ScalarType>
65 ScalarType diff(std::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
66 {
67  std::vector<ScalarType> v2_cpu(v2.size());
68  viennacl::backend::finish();
69  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
70 
71  for (unsigned int i=0;i<v1.size(); ++i)
72  {
73  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
74  {
75  //if (std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) < 1e-10 ) //absolute tolerance (avoid round-off issues)
76  // v2_cpu[i] = 0;
77  //else
78  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
79  }
80  else
81  v2_cpu[i] = 0.0;
82 
83  if (v2_cpu[i] > 0.0001)
84  {
85  //std::cout << "Neighbor: " << i-1 << ": " << v1[i-1] << " vs. " << v2_cpu[i-1] << std::endl;
86  std::cout << "Error at entry " << i << ": Should: " << v1[i] << " vs. Is: " << v2[i] << std::endl;
87  //std::cout << "Neighbor: " << i+1 << ": " << v1[i+1] << " vs. " << v2_cpu[i+1] << std::endl;
88  exit(EXIT_FAILURE);
89  }
90  }
91 
92  ScalarType norm_inf = 0;
93  for (std::size_t i=0; i<v2_cpu.size(); ++i)
94  norm_inf = std::max<ScalarType>(norm_inf, std::fabs(v2_cpu[i]));
95 
96  return norm_inf;
97 }
98 
99 
100 template<typename IndexT, typename NumericT, typename SparseMatrixT>
101 NumericT diff(std::vector<std::map<IndexT, NumericT> > & cpu_A, SparseMatrixT & vcl_A)
102 {
103  typedef typename std::map<IndexT, NumericT>::const_iterator RowIterator;
104 
105  std::vector<std::map<IndexT, NumericT> > from_gpu(vcl_A.size1());
106 
107  viennacl::backend::finish();
108  viennacl::copy(vcl_A, from_gpu);
109 
110  NumericT error = 0;
111 
112  //step 1: compare all entries from cpu_A with vcl_A:
113  for (std::size_t i=0; i<cpu_A.size(); ++i)
114  {
115  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
116  for (RowIterator it = cpu_A[i].begin(); it != cpu_A[i].end(); ++it)
117  {
118  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
119  NumericT current_error = 0;
120  NumericT val_cpu_A = it->second;
121  NumericT val_gpu_A = from_gpu[i][it->first];
122 
123  NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
124  if (max_val > 0)
125  current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
126  if (current_error > error)
127  error = current_error;
128  }
129  }
130 
131  //step 2: compare all entries from gpu_matrix with cpu_matrix (sparsity pattern might differ):
132  //std::cout << "ViennaCL matrix: " << std::endl;
133  for (std::size_t i=0; i<from_gpu.size(); ++i)
134  {
135  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
136  for (RowIterator it = from_gpu[i].begin(); it != from_gpu[i].end(); ++it)
137  {
138  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
139  NumericT current_error = 0;
140  NumericT val_gpu_A = it->second;
141  NumericT val_cpu_A = cpu_A[i][it->first];
142 
143  NumericT max_val = std::max(std::fabs(val_cpu_A), std::fabs(val_gpu_A));
144  if (max_val > 0)
145  current_error = std::fabs(val_cpu_A - val_gpu_A) / max_val;
146  if (current_error > error)
147  error = current_error;
148  }
149  }
150 
151  return error;
152 }
153 
154 
155 template<typename NumericT, typename VCL_MatrixT, typename Epsilon, typename STLVectorT, typename VCLVectorT>
156 int strided_matrix_vector_product_test(Epsilon epsilon,
157  STLVectorT & result, STLVectorT const & rhs,
158  VCLVectorT & vcl_result, VCLVectorT & vcl_rhs)
159 {
160  typedef typename std::map<unsigned int, NumericT>::const_iterator RowIterator;
161  int retval = EXIT_SUCCESS;
162 
163  std::vector<std::map<unsigned int, NumericT> > std_A(5);
164  std_A[0][0] = NumericT(2.0); std_A[0][2] = NumericT(-1.0);
165  std_A[1][0] = NumericT(3.0); std_A[1][2] = NumericT(-5.0);
166  std_A[2][1] = NumericT(5.0); std_A[2][2] = NumericT(-2.0);
167  std_A[3][2] = NumericT(1.0); std_A[3][3] = NumericT(-6.0);
168  std_A[4][1] = NumericT(7.0); std_A[4][2] = NumericT(-5.0);
169  //the following computes project(result, slice(1, 3, 5)) = prod(std_A, project(rhs, slice(3, 2, 4)));
170  for (std::size_t i=0; i<5; ++i)
171  {
172  NumericT val = 0;
173  for (RowIterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
174  val += it->second * rhs[3 + 2*it->first];
175  result[1 + 3*i] = val;
176  }
177 
178  VCL_MatrixT vcl_sparse_matrix2;
179  viennacl::copy(std_A, vcl_sparse_matrix2);
180  viennacl::vector<NumericT> vec(4);
181  vec(0) = rhs[3];
182  vec(1) = rhs[5];
183  vec(2) = rhs[7];
184  vec(3) = rhs[9];
185  viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, viennacl::project(vcl_rhs, viennacl::slice(3, 2, 4)));
186 
187  if ( std::fabs(diff(result, vcl_result)) > epsilon )
188  {
189  std::cout << "# Error at operation: matrix-vector product with strided vectors, part 1" << std::endl;
190  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
191  retval = EXIT_FAILURE;
192  }
193  vcl_result(1) = NumericT(1.0);
194  vcl_result(4) = NumericT(1.0);
195  vcl_result(7) = NumericT(1.0);
196  vcl_result(10) = NumericT(1.0);
197  vcl_result(13) = NumericT(1.0);
198 
199  viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, vec);
200 
201  if ( std::fabs(diff(result, vcl_result)) > epsilon )
202  {
203  std::cout << "# Error at operation: matrix-vector product with strided vectors, part 2" << std::endl;
204  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
205  retval = EXIT_FAILURE;
206  }
207 
208  return retval;
209 }
210 
211 
212 template< typename NumericT, typename VCL_MATRIX, typename Epsilon >
213 int resize_test(Epsilon const& epsilon)
214 {
215  int retval = EXIT_SUCCESS;
216 
217  std::vector<std::map<unsigned int, NumericT> > std_A(5);
218  VCL_MATRIX vcl_matrix;
219 
220  std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
221  std_A[1][0] = NumericT(1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
222  std_A[2][0] = NumericT(2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
223  std_A[3][0] = NumericT(3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
224  std_A[4][0] = NumericT(4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
225 
226  viennacl::copy(std_A, vcl_matrix);
227  std::vector<std::map<unsigned int, NumericT> > std_B(std_A.size());
228  viennacl::copy(vcl_matrix, std_B);
229 
230  std::cout << "Checking for equality after copy..." << std::endl;
231  if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
232  {
233  std::cout << "# Error at operation: equality after copy with sparse matrix" << std::endl;
234  std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
235  return EXIT_FAILURE;
236  }
237 
238  std::cout << "Testing resize to larger..." << std::endl;
239  std_A.resize(10);
240  std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
241  std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
242  std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
243  std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
244  std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
245 
246  vcl_matrix.resize(10, 10, true);
247 
248  if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
249  {
250  std::cout << "# Error at operation: resize (to larger) with sparse matrix" << std::endl;
251  std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
252  return EXIT_FAILURE;
253  }
254 
255  std_A[5][5] = NumericT(5.5); std_A[5][6] = NumericT(5.6); std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
256  std_A[6][5] = NumericT(6.5); std_A[6][6] = NumericT(6.6); std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
257  std_A[7][5] = NumericT(7.5); std_A[7][6] = NumericT(7.6); std_A[7][7] = NumericT(7.7); std_A[7][8] = NumericT(7.8); std_A[7][9] = NumericT(7.9);
258  std_A[8][5] = NumericT(8.5); std_A[8][6] = NumericT(8.6); std_A[8][7] = NumericT(8.7); std_A[8][8] = NumericT(8.8); std_A[8][9] = NumericT(8.9);
259  std_A[9][5] = NumericT(9.5); std_A[9][6] = NumericT(9.6); std_A[9][7] = NumericT(9.7); std_A[9][8] = NumericT(9.8); std_A[9][9] = NumericT(9.9);
260  viennacl::copy(std_A, vcl_matrix);
261 
262  std::cout << "Testing resize to smaller..." << std::endl;
263  std_A.clear();
264  std_A.resize(7);
265  std_A[0][0] = NumericT(10.0); std_A[0][1] = NumericT(0.1); std_A[0][2] = NumericT(0.2); std_A[0][3] = NumericT(0.3); std_A[0][4] = NumericT(0.4);
266  std_A[1][0] = NumericT( 1.0); std_A[1][1] = NumericT(1.1); std_A[1][2] = NumericT(1.2); std_A[1][3] = NumericT(1.3); std_A[1][4] = NumericT(1.4);
267  std_A[2][0] = NumericT( 2.0); std_A[2][1] = NumericT(2.1); std_A[2][2] = NumericT(2.2); std_A[2][3] = NumericT(2.3); std_A[2][4] = NumericT(2.4);
268  std_A[3][0] = NumericT( 3.0); std_A[3][1] = NumericT(3.1); std_A[3][2] = NumericT(3.2); std_A[3][3] = NumericT(3.3); std_A[3][4] = NumericT(3.4);
269  std_A[4][0] = NumericT( 4.0); std_A[4][1] = NumericT(4.1); std_A[4][2] = NumericT(4.2); std_A[4][3] = NumericT(4.3); std_A[4][4] = NumericT(4.4);
270  std_A[5][5] = NumericT( 5.5); std_A[5][6] = NumericT(5.6); //std_A[5][7] = NumericT(5.7); std_A[5][8] = NumericT(5.8); std_A[5][9] = NumericT(5.9);
271  std_A[6][5] = NumericT( 6.5); std_A[6][6] = NumericT(6.6); //std_A[6][7] = NumericT(6.7); std_A[6][8] = NumericT(6.8); std_A[6][9] = NumericT(6.9);
272 
273  vcl_matrix.resize(7, 7);
274 
275  //std::cout << std_A << std::endl;
276  if ( std::fabs(diff(std_A, vcl_matrix)) > epsilon )
277  {
278  std::cout << "# Error at operation: resize (to smaller) with sparse matrix" << std::endl;
279  std::cout << " diff: " << std::fabs(diff(std_A, vcl_matrix)) << std::endl;
280  retval = EXIT_FAILURE;
281  }
282 
283  std::vector<NumericT> std_vec(std_A.size(), NumericT(3.1415));
284  viennacl::vector<NumericT> vcl_vec(std_A.size());
285 
286 
287  std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
288  viennacl::copy(std_vec, vcl_vec);
289 
290  std::cout << "STL..." << std::endl;
291  //boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::unit_lower_tag());
292  for (std::size_t i=1; i<std_A.size(); ++i)
293  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_A[i].begin(); it != std_A[i].end(); ++it)
294  {
295  if (it->first < static_cast<unsigned int>(i))
296  std_vec[i] -= it->second * std_vec[it->first];
297  else
298  continue;
299  }
300 
301  std::cout << "ViennaCL..." << std::endl;
302  viennacl::linalg::inplace_solve((vcl_matrix), vcl_vec, viennacl::linalg::unit_lower_tag());
303 
304  if ( std::fabs(diff(std_vec, vcl_vec)) > epsilon )
305  {
306  std::cout << "# Error at operation: unit lower triangular solve" << std::endl;
307  std::cout << " diff: " << std::fabs(diff(std_vec, vcl_vec)) << std::endl;
308  retval = EXIT_FAILURE;
309  }
310  return retval;
311 }
312 
313 
314 //
315 // -------------------------------------------------------------
316 //
317 template< typename NumericT, typename Epsilon >
318 int test(Epsilon const& epsilon)
319 {
320  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
321 
322  std::cout << "Testing resizing of compressed_matrix..." << std::endl;
323  int retval = resize_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon);
324  if (retval != EXIT_SUCCESS)
325  return retval;
326 
327  // --------------------------------------------------------------------------
328  std::vector<NumericT> rhs;
329  std::vector<NumericT> result;
330  std::vector<std::map<unsigned int, NumericT> > std_matrix;
331 
332  if (viennacl::io::read_matrix_market_file(std_matrix, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
333  {
334  std::cout << "Error reading Matrix file" << std::endl;
335  return EXIT_FAILURE;
336  }
337 
338  //unsigned int cg_mat_size = cg_mat.size();
339  std::cout << "done reading matrix" << std::endl;
340 
341 
342  rhs.resize(std_matrix.size());
343  for (std::size_t i=0; i<rhs.size(); ++i)
344  {
345  std_matrix[i][static_cast<unsigned int>(i)] = NumericT(0.5); // Get rid of round-off errors by making row-sums unequal to zero:
346  rhs[i] = NumericT(1) + randomNumber();
347  }
348 
349  // add some random numbers to the double-compressed matrix:
350  std::vector<std::map<unsigned int, NumericT> > std_cc_matrix(std_matrix.size());
351  std_cc_matrix[42][199] = NumericT(3.1415);
352  std_cc_matrix[31][69] = NumericT(2.71);
353  std_cc_matrix[23][32] = NumericT(6);
354  std_cc_matrix[177][57] = NumericT(4);
355  std_cc_matrix[21][97] = NumericT(-4);
356  std_cc_matrix[92][25] = NumericT(2);
357  std_cc_matrix[89][62] = NumericT(11);
358  std_cc_matrix[ 1][ 7] = NumericT(8);
359  std_cc_matrix[85][41] = NumericT(13);
360  std_cc_matrix[66][28] = NumericT(8);
361  std_cc_matrix[21][74] = NumericT(-2);
362  viennacl::tools::sparse_matrix_adapter<NumericT> adapted_std_cc_matrix(std_cc_matrix, std_matrix.size(), std_matrix.size());
363 
364 
365  result = rhs;
366 
367 
368  viennacl::vector<NumericT> vcl_rhs(rhs.size());
369  viennacl::vector<NumericT> vcl_result(result.size());
370  viennacl::vector<NumericT> vcl_result2(result.size());
371  viennacl::compressed_matrix<NumericT> vcl_compressed_matrix(rhs.size(), rhs.size());
372  viennacl::compressed_compressed_matrix<NumericT> vcl_compressed_compressed_matrix(rhs.size(), rhs.size());
373  viennacl::coordinate_matrix<NumericT> vcl_coordinate_matrix(rhs.size(), rhs.size());
374  viennacl::ell_matrix<NumericT> vcl_ell_matrix;
375  viennacl::sliced_ell_matrix<NumericT> vcl_sliced_ell_matrix;
376  viennacl::hyb_matrix<NumericT> vcl_hyb_matrix;
377 
378  viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
379  viennacl::copy(std_matrix, vcl_compressed_matrix);
380  viennacl::copy(adapted_std_cc_matrix, vcl_compressed_compressed_matrix);
381  viennacl::copy(std_matrix, vcl_coordinate_matrix);
382 
383  // --------------------------------------------------------------------------
384  std::cout << "Testing products: STL" << std::endl;
385  result = viennacl::linalg::prod(std_matrix, rhs);
386 
387  std::cout << "Testing products: compressed_matrix" << std::endl;
388  vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
389 
390  if ( std::fabs(diff(result, vcl_result)) > epsilon )
391  {
392  std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
393  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
394  return EXIT_FAILURE;
395  }
396 
397  std::cout << "Testing products: compressed_matrix, strided vectors" << std::endl;
398  retval = strided_matrix_vector_product_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
399  if (retval != EXIT_SUCCESS)
400  return retval;
401 
402  result = rhs;
403  result = viennacl::linalg::prod(std_matrix, rhs);
404  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
405  vcl_result = vcl_rhs;
406  vcl_result += viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
407 
408  if ( std::fabs(diff(result, vcl_result)) > epsilon )
409  {
410  std::cout << "# Error at operation: matrix-vector product with compressed_matrix (+=)" << std::endl;
411  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
412  return EXIT_FAILURE;
413  }
414 
415  result = rhs;
416  result = viennacl::linalg::prod(std_matrix, rhs);
417  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
418  vcl_result = vcl_rhs;
419  vcl_result -= viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
420 
421  if ( std::fabs(diff(result, vcl_result)) > epsilon )
422  {
423  std::cout << "# Error at operation: matrix-vector product with compressed_matrix (-=)" << std::endl;
424  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
425  return EXIT_FAILURE;
426  }
427 
428  //
429  // Triangular solvers for A \ b:
430  //
431 
432  std::cout << "Testing unit upper triangular solve: compressed_matrix" << std::endl;
433  result = rhs;
434  viennacl::copy(result, vcl_result);
435  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_upper_tag());
436  for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
437  {
438  std::size_t row = std_matrix.size() - i2 - 1;
439  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
440  {
441  if (it->first > static_cast<unsigned int>(row))
442  result[row] -= it->second * result[it->first];
443  else
444  continue;
445  }
446  }
447 
448  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_upper_tag());
449 
450  if ( std::fabs(diff(result, vcl_result)) > epsilon )
451  {
452  std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
453  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
454  return EXIT_FAILURE;
455  }
456 
458 
459  std::cout << "Testing upper triangular solve: compressed_matrix" << std::endl;
460  result = rhs;
461  viennacl::copy(result, vcl_result);
462  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::upper_tag());
463  for (std::size_t i2=0; i2<std_matrix.size(); ++i2)
464  {
465  std::size_t row = std_matrix.size() - i2 - 1;
466  NumericT diag = 0;
467  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[row].begin(); it != std_matrix[row].end(); ++it)
468  {
469  if (it->first > static_cast<unsigned int>(row))
470  result[row] -= it->second * result[it->first];
471  else if (it->first == static_cast<unsigned int>(row))
472  diag = it->second;
473  else
474  continue;
475  }
476  result[row] /= diag;
477  }
478 
479  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::upper_tag());
480 
481  if ( std::fabs(diff(result, vcl_result)) > epsilon )
482  {
483  std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
484  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
485  return EXIT_FAILURE;
486  }
487 
489 
490  std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
491  result = rhs;
492  viennacl::copy(result, vcl_result);
493  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_lower_tag());
494  for (std::size_t i=1; i<std_matrix.size(); ++i)
495  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
496  {
497  if (it->first < static_cast<unsigned int>(i))
498  result[i] -= it->second * result[it->first];
499  else
500  continue;
501  }
502  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_lower_tag());
503 
504 
505  if ( std::fabs(diff(result, vcl_result)) > epsilon )
506  {
507  std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
508  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
509  return EXIT_FAILURE;
510  }
511 
512 
513  std::cout << "Testing lower triangular solve: compressed_matrix" << std::endl;
514  result = rhs;
515  viennacl::copy(result, vcl_result);
516  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::lower_tag());
517  for (std::size_t i=0; i<std_matrix.size(); ++i)
518  {
519  NumericT diag = 0;
520  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
521  {
522  if (it->first < static_cast<unsigned int>(i))
523  result[i] -= it->second * result[it->first];
524  else if (it->first == static_cast<unsigned int>(i))
525  diag = it->second;
526  else
527  continue;
528  }
529  result[i] /= diag;
530  }
531  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::lower_tag());
532 
533 
534  if ( std::fabs(diff(result, vcl_result)) > epsilon )
535  {
536  std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
537  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
538  return EXIT_FAILURE;
539  }
540 
541 
542 
543  //
544  // Triangular solvers for A^T \ b
545  //
546  std::vector<std::map<unsigned int, NumericT> > std_matrix_trans(std_matrix.size());
547 
548  // compute transpose:
549  for (std::size_t i=0; i<std_matrix.size(); ++i)
550  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix[i].begin(); it != std_matrix[i].end(); ++it)
551  std_matrix_trans[i][it->first] = it->second;
552 
553  std::cout << "Testing transposed unit upper triangular solve: compressed_matrix" << std::endl;
554  result = rhs;
555  viennacl::copy(result, vcl_result);
556  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_upper_tag());
557  for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
558  {
559  std::size_t row = std_matrix_trans.size() - i2 - 1;
560  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
561  {
562  if (it->first > static_cast<unsigned int>(row))
563  result[row] -= it->second * result[it->first];
564  else
565  continue;
566  }
567  }
568  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_upper_tag());
569 
570  if ( std::fabs(diff(result, vcl_result)) > epsilon )
571  {
572  std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
573  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
574  return EXIT_FAILURE;
575  }
576 
578 
579  std::cout << "Testing transposed upper triangular solve: compressed_matrix" << std::endl;
580  result = rhs;
581  viennacl::copy(result, vcl_result);
582  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::upper_tag());
583  for (std::size_t i2=0; i2<std_matrix_trans.size(); ++i2)
584  {
585  std::size_t row = std_matrix_trans.size() - i2 - 1;
586  NumericT diag = 0;
587  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[row].begin(); it != std_matrix_trans[row].end(); ++it)
588  {
589  if (it->first > static_cast<unsigned int>(row))
590  result[row] -= it->second * result[it->first];
591  else if (it->first == static_cast<unsigned int>(row))
592  diag = it->second;
593  else
594  continue;
595  }
596  result[row] /= diag;
597  }
598  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::upper_tag());
599 
600  if ( std::fabs(diff(result, vcl_result)) > epsilon )
601  {
602  std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
603  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
604  return EXIT_FAILURE;
605  }
606 
608 
609  std::cout << "Testing transposed unit lower triangular solve: compressed_matrix" << std::endl;
610  result = rhs;
611  viennacl::copy(result, vcl_result);
612  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_lower_tag());
613  for (std::size_t i=1; i<std_matrix_trans.size(); ++i)
614  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
615  {
616  if (it->first < static_cast<unsigned int>(i))
617  result[i] -= it->second * result[it->first];
618  else
619  continue;
620  }
621  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_lower_tag());
622 
623  if ( std::fabs(diff(result, vcl_result)) > epsilon )
624  {
625  std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
626  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
627  return EXIT_FAILURE;
628  }
629 
631 
632  std::cout << "Testing transposed lower triangular solve: compressed_matrix" << std::endl;
633  result = rhs;
634  viennacl::copy(result, vcl_result);
635  //boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::lower_tag());
636  for (std::size_t i=0; i<std_matrix_trans.size(); ++i)
637  {
638  NumericT diag = 0;
639  for (typename std::map<unsigned int, NumericT>::const_iterator it = std_matrix_trans[i].begin(); it != std_matrix_trans[i].end(); ++it)
640  {
641  if (it->first < static_cast<unsigned int>(i))
642  result[i] -= it->second * result[it->first];
643  else if (it->first == static_cast<unsigned int>(i))
644  diag = it->second;
645  else
646  continue;
647  }
648  result[i] /= diag;
649  }
650  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::lower_tag());
651 
652  if ( std::fabs(diff(result, vcl_result)) > epsilon )
653  {
654  std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
655  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
656  return EXIT_FAILURE;
657  }
658 
659 
660  //
662  //
663 
664 
665  std::cout << "Testing products: compressed_compressed_matrix" << std::endl;
666  result = viennacl::linalg::prod(std_cc_matrix, rhs);
667  vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
668 
669  if ( std::fabs(diff(result, vcl_result)) > epsilon )
670  {
671  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (=)" << std::endl;
672  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
673  return EXIT_FAILURE;
674  }
675 
676  {
677  std::vector<std::map<unsigned int, NumericT> > temp(vcl_compressed_compressed_matrix.size1());
678  viennacl::copy(vcl_compressed_compressed_matrix, temp);
679 
680  // check that entries are correct by computing the product again:
681  result = viennacl::linalg::prod(temp, rhs);
682 
683  if ( std::fabs(diff(result, vcl_result)) > epsilon )
684  {
685  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (after copy back)" << std::endl;
686  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
687  return EXIT_FAILURE;
688  }
689 
690  }
691 
692  result = rhs;
693  result = viennacl::linalg::prod(std_cc_matrix, rhs);
694  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
695  vcl_result = vcl_rhs;
696  vcl_result += viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
697 
698  if ( std::fabs(diff(result, vcl_result)) > epsilon )
699  {
700  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (+=)" << std::endl;
701  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
702  return EXIT_FAILURE;
703  }
704 
705  result = rhs;
706  result = viennacl::linalg::prod(std_cc_matrix, rhs);
707  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
708  vcl_result = vcl_rhs;
709  vcl_result -= viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
710 
711  if ( std::fabs(diff(result, vcl_result)) > epsilon )
712  {
713  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (-=)" << std::endl;
714  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
715  return EXIT_FAILURE;
716  }
717 
718 
719  //
721  //
722 
723 
724  std::cout << "Testing products: coordinate_matrix" << std::endl;
725  result = viennacl::linalg::prod(std_matrix, rhs);
726  vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
727 
728  if ( std::fabs(diff(result, vcl_result)) > epsilon )
729  {
730  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
731  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
732  return EXIT_FAILURE;
733  }
734 
735  std::cout << "Testing products: coordinate_matrix, strided vectors" << std::endl;
736  //std::cout << " --> SKIPPING <--" << std::endl;
737  retval = strided_matrix_vector_product_test<NumericT, viennacl::coordinate_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
738  if (retval != EXIT_SUCCESS)
739  return retval;
740 
741  result = rhs;
742  result = viennacl::linalg::prod(std_matrix, rhs);
743  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
744  vcl_result = vcl_rhs;
745  vcl_result += viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
746 
747  if ( std::fabs(diff(result, vcl_result)) > epsilon )
748  {
749  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (+=)" << std::endl;
750  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
751  return EXIT_FAILURE;
752  }
753 
754  result = rhs;
755  result = viennacl::linalg::prod(std_matrix, rhs);
756  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
757  vcl_result = vcl_rhs;
758  vcl_result -= viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
759 
760  if ( std::fabs(diff(result, vcl_result)) > epsilon )
761  {
762  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix (-=)" << std::endl;
763  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
764  return EXIT_FAILURE;
765  }
766 
767  //
769  //
770 
771 
772  //std::cout << "Copying ell_matrix" << std::endl;
773  viennacl::copy(std_matrix, vcl_ell_matrix);
774  std_matrix.clear();
775  viennacl::copy(vcl_ell_matrix, std_matrix);// just to check that it works
776 
777 
778  std::cout << "Testing products: ell_matrix" << std::endl;
779  result = viennacl::linalg::prod(std_matrix, rhs);
780  vcl_result.clear();
781  vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
782  //viennacl::linalg::prod_impl(vcl_ell_matrix, vcl_rhs, vcl_result);
783  //std::cout << vcl_result << "\n";
784  //std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
785  //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
786 
787  if ( std::fabs(diff(result, vcl_result)) > epsilon )
788  {
789  std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
790  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
791  return EXIT_FAILURE;
792  }
793 
794  std::cout << "Testing products: ell_matrix, strided vectors" << std::endl;
795  retval = strided_matrix_vector_product_test<NumericT, viennacl::ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
796  if (retval != EXIT_SUCCESS)
797  return retval;
798 
799  result = rhs;
800  result = viennacl::linalg::prod(std_matrix, rhs);
801  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
802  vcl_result = vcl_rhs;
803  vcl_result += viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
804 
805  if ( std::fabs(diff(result, vcl_result)) > epsilon )
806  {
807  std::cout << "# Error at operation: matrix-vector product with ell_matrix (+=)" << std::endl;
808  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
809  return EXIT_FAILURE;
810  }
811 
812  result = rhs;
813  result = viennacl::linalg::prod(std_matrix, rhs);
814  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
815  vcl_result = vcl_rhs;
816  vcl_result -= viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
817 
818  if ( std::fabs(diff(result, vcl_result)) > epsilon )
819  {
820  std::cout << "# Error at operation: matrix-vector product with ell_matrix (-=)" << std::endl;
821  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
822  return EXIT_FAILURE;
823  }
824 
825  //
827  //
828 
829 
830  //std::cout << "Copying sliced_ell_matrix" << std::endl;
831  viennacl::copy(std_matrix, vcl_sliced_ell_matrix);
832 
833  std::cout << "Testing products: sliced_ell_matrix" << std::endl;
834  result = viennacl::linalg::prod(std_matrix, rhs);
835  vcl_result.clear();
836  vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
837 
838  if ( std::fabs(diff(result, vcl_result)) > epsilon )
839  {
840  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
841  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
842  return EXIT_FAILURE;
843  }
844 
845  std::cout << "Testing products: sliced_ell_matrix, strided vectors" << std::endl;
846  retval = strided_matrix_vector_product_test<NumericT, viennacl::sliced_ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
847  if (retval != EXIT_SUCCESS)
848  return retval;
849 
850  result = rhs;
851  result = viennacl::linalg::prod(std_matrix, rhs);
852  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
853  vcl_result = vcl_rhs;
854  vcl_result += viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
855 
856  if ( std::fabs(diff(result, vcl_result)) > epsilon )
857  {
858  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (+=)" << std::endl;
859  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
860  return EXIT_FAILURE;
861  }
862 
863  result = rhs;
864  result = viennacl::linalg::prod(std_matrix, rhs);
865  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
866  vcl_result = vcl_rhs;
867  vcl_result -= viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
868 
869  if ( std::fabs(diff(result, vcl_result)) > epsilon )
870  {
871  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix (-=)" << std::endl;
872  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
873  return EXIT_FAILURE;
874  }
875 
876  //
878  //
879 
880 
881  //std::cout << "Copying hyb_matrix" << std::endl;
882  viennacl::copy(std_matrix, vcl_hyb_matrix);
883  std_matrix.clear();
884  viennacl::copy(vcl_hyb_matrix, std_matrix);// just to check that it works
885  viennacl::copy(std_matrix, vcl_hyb_matrix);
886 
887  std::cout << "Testing products: hyb_matrix" << std::endl;
888  result = viennacl::linalg::prod(std_matrix, rhs);
889  vcl_result.clear();
890  vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
891  //viennacl::linalg::prod_impl(vcl_hyb_matrix, vcl_rhs, vcl_result);
892  //std::cout << vcl_result << "\n";
893  //std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
894  //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
895 
896  if ( std::fabs(diff(result, vcl_result)) > epsilon )
897  {
898  std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
899  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
900  return EXIT_FAILURE;
901  }
902 
903  std::cout << "Testing products: hyb_matrix, strided vectors" << std::endl;
904  retval = strided_matrix_vector_product_test<NumericT, viennacl::hyb_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
905  if (retval != EXIT_SUCCESS)
906  return retval;
907 
908  result = rhs;
909  result = viennacl::linalg::prod(std_matrix, rhs);
910  for (std::size_t i=0; i<result.size(); ++i) result[i] += rhs[i];
911  vcl_result = vcl_rhs;
912  vcl_result += viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
913 
914  if ( std::fabs(diff(result, vcl_result)) > epsilon )
915  {
916  std::cout << "# Error at operation: matrix-vector product with hyb_matrix (+=)" << std::endl;
917  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
918  return EXIT_FAILURE;
919  }
920 
921  result = rhs;
922  result = viennacl::linalg::prod(std_matrix, rhs);
923  for (std::size_t i=0; i<result.size(); ++i) result[i] = rhs[i] - result[i];
924  vcl_result = vcl_rhs;
925  vcl_result -= viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
926 
927  if ( std::fabs(diff(result, vcl_result)) > epsilon )
928  {
929  std::cout << "# Error at operation: matrix-vector product with hyb_matrix (-=)" << std::endl;
930  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
931  return EXIT_FAILURE;
932  }
933 
934 
935  // --------------------------------------------------------------------------
936  // --------------------------------------------------------------------------
937  NumericT alpha = static_cast<NumericT>(2.786);
938  NumericT beta = static_cast<NumericT>(1.432);
939  copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
940  copy(result.begin(), result.end(), vcl_result.begin());
941  copy(result.begin(), result.end(), vcl_result2.begin());
942 
943  std::cout << "Testing scaled additions of products and vectors" << std::endl;
944  std::vector<NumericT> result2(result);
945  result2 = viennacl::linalg::prod(std_matrix, rhs);
946  for (std::size_t i=0; i<result.size(); ++i)
947  result[i] = alpha * result2[i] + beta * result[i];
948  vcl_result2 = alpha * viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) + beta * vcl_result;
949 
950  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
951  {
952  std::cout << "# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
953  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
954  return EXIT_FAILURE;
955  }
956 
957 
958  vcl_result2.clear();
959  vcl_result2 = alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result;
960 
961  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
962  {
963  std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
964  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
965  return EXIT_FAILURE;
966  }
967 
968  vcl_result2.clear();
969  vcl_result2 = alpha * viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs) + beta * vcl_result;
970 
971  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
972  {
973  std::cout << "# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
974  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
975  return EXIT_FAILURE;
976  }
977 
978  vcl_result2.clear();
979  vcl_result2 = alpha * viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs) + beta * vcl_result;
980 
981  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
982  {
983  std::cout << "# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
984  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
985  return EXIT_FAILURE;
986  }
987 
989  std_matrix.clear();
990 
991  std::cout << "Testing products after clear(): compressed_matrix" << std::endl;
992  vcl_compressed_matrix.clear();
993  result = std::vector<NumericT>(result.size(), NumericT(1));
994  result = viennacl::linalg::prod(std_matrix, rhs);
995  vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
996 
997  if ( std::fabs(diff(result, vcl_result)) > epsilon )
998  {
999  std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
1000  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1001  return EXIT_FAILURE;
1002  }
1003 
1004  std::cout << "Testing products after clear(): compressed_compressed_matrix" << std::endl;
1005  vcl_compressed_compressed_matrix.clear();
1006  result = std::vector<NumericT>(result.size(), NumericT(1));
1007  result = viennacl::linalg::prod(std_matrix, rhs);
1008  vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
1009 
1010  if ( std::fabs(diff(result, vcl_result)) > epsilon )
1011  {
1012  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix" << std::endl;
1013  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1014  return EXIT_FAILURE;
1015  }
1016 
1017  std::cout << "Testing products after clear(): coordinate_matrix" << std::endl;
1018  vcl_coordinate_matrix.clear();
1019  result = std::vector<NumericT>(result.size(), NumericT(1));
1020  result = viennacl::linalg::prod(std_matrix, rhs);
1021  vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
1022 
1023  if ( std::fabs(diff(result, vcl_result)) > epsilon )
1024  {
1025  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
1026  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1027  return EXIT_FAILURE;
1028  }
1029 
1030  std::cout << "Testing products after clear(): ell_matrix" << std::endl;
1031  vcl_ell_matrix.clear();
1032  result = std::vector<NumericT>(result.size(), NumericT(1));
1033  result = viennacl::linalg::prod(std_matrix, rhs);
1034  vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
1035 
1036  if ( std::fabs(diff(result, vcl_result)) > epsilon )
1037  {
1038  std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
1039  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1040  return EXIT_FAILURE;
1041  }
1042 
1043  std::cout << "Testing products after clear(): hyb_matrix" << std::endl;
1044  vcl_hyb_matrix.clear();
1045  result = std::vector<NumericT>(result.size(), NumericT(1));
1046  result = viennacl::linalg::prod(std_matrix, rhs);
1047  vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
1048 
1049  if ( std::fabs(diff(result, vcl_result)) > epsilon )
1050  {
1051  std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
1052  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1053  return EXIT_FAILURE;
1054  }
1055 
1056  std::cout << "Testing products after clear(): sliced_ell_matrix" << std::endl;
1057  vcl_sliced_ell_matrix.clear();
1058  result = std::vector<NumericT>(result.size(), NumericT(1));
1059  result = viennacl::linalg::prod(std_matrix, rhs);
1060  vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
1061 
1062  if ( std::fabs(diff(result, vcl_result)) > epsilon )
1063  {
1064  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
1065  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
1066  return EXIT_FAILURE;
1067  }
1068 
1069 
1070  // --------------------------------------------------------------------------
1071  return retval;
1072 }
1073 //
1074 // -------------------------------------------------------------
1075 //
1076 int main()
1077 {
1078  std::cout << std::endl;
1079  std::cout << "----------------------------------------------" << std::endl;
1080  std::cout << "----------------------------------------------" << std::endl;
1081  std::cout << "## Test :: Sparse Matrices" << std::endl;
1082  std::cout << "----------------------------------------------" << std::endl;
1083  std::cout << "----------------------------------------------" << std::endl;
1084  std::cout << std::endl;
1085 
1086  int retval = EXIT_SUCCESS;
1087 
1088  std::cout << std::endl;
1089  std::cout << "----------------------------------------------" << std::endl;
1090  std::cout << std::endl;
1091  {
1092  typedef float NumericT;
1093  NumericT epsilon = static_cast<NumericT>(1E-4);
1094  std::cout << "# Testing setup:" << std::endl;
1095  std::cout << " eps: " << epsilon << std::endl;
1096  std::cout << " numeric: float" << std::endl;
1097  retval = test<NumericT>(epsilon);
1098  if ( retval == EXIT_SUCCESS )
1099  std::cout << "# Test passed" << std::endl;
1100  else
1101  return retval;
1102  }
1103  std::cout << std::endl;
1104  std::cout << "----------------------------------------------" << std::endl;
1105  std::cout << std::endl;
1106 
1107 #ifdef VIENNACL_WITH_OPENCL
1108  if ( viennacl::ocl::current_device().double_support() )
1109 #endif
1110  {
1111  {
1112  typedef double NumericT;
1113  NumericT epsilon = 1.0E-12;
1114  std::cout << "# Testing setup:" << std::endl;
1115  std::cout << " eps: " << epsilon << std::endl;
1116  std::cout << " numeric: double" << std::endl;
1117  retval = test<NumericT>(epsilon);
1118  if ( retval == EXIT_SUCCESS )
1119  std::cout << "# Test passed" << std::endl;
1120  else
1121  return retval;
1122  }
1123  std::cout << std::endl;
1124  std::cout << "----------------------------------------------" << std::endl;
1125  std::cout << std::endl;
1126  }
1127 #ifdef VIENNACL_WITH_OPENCL
1128  else
1129  std::cout << "No double precision support, skipping test..." << std::endl;
1130 #endif
1131 
1132 
1133  std::cout << std::endl;
1134  std::cout << "------- Test completed --------" << std::endl;
1135  std::cout << std::endl;
1136 
1137  return retval;
1138 }
viennacl::hyb_matrix
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
Definition: forwards.h:406
viennacl::sliced_ell_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: sliced_ell_matrix.hpp:85
viennacl::linalg::inplace_solve
void inplace_solve(const matrix_base< NumericT > &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:217
matrix_market.hpp
A reader and writer for the matrix market format is implemented here.
strided_matrix_vector_product_test
int strided_matrix_vector_product_test(Epsilon epsilon, STLVectorT &result, STLVectorT const &rhs, VCLVectorT &vcl_result, VCLVectorT &vcl_rhs)
Definition: sparse.cpp:156
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: sparse.cpp:57
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
int test(Epsilon const &epsilon)
Definition: sparse.cpp:318
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...
main
int main()
Definition: sparse.cpp:1076
viennacl::hyb_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: hyb_matrix.hpp:69
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
viennacl::backend::finish
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
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
coordinate_matrix.hpp
Implementation of the coordinate_matrix class.
NumericT
float NumericT
Definition: bisect.cpp:40
v1
viennacl::vector< float > v1
Definition: global_variables.cpp:60
hyb_matrix.hpp
Implementation of the hyb_matrix class.
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
viennacl::ell_matrix
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
viennacl::vector_base::begin
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:841
ilu.hpp
Implementations of incomplete factorization preconditioners. Convenience header file.
viennacl::linalg::upper_tag
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
viennacl::sliced_ell_matrix
Sparse matrix class using the sliced ELLPACK with parameters C, .
Definition: forwards.h:403
compressed_matrix.hpp
Implementation of the compressed_matrix class.
viennacl::compressed_compressed_matrix
A sparse square matrix in compressed sparse rows format optimized for the case that only a few rows c...
Definition: compressed_compressed_matrix.hpp:265
sliced_ell_matrix.hpp
Implementation of the sliced_ell_matrix class.
viennacl::ocl::device::double_support
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
viennacl::project
matrix_range< MatrixType > project(MatrixType const &A, viennacl::range const &r1, viennacl::range const &r2)
Definition: matrix_proxy.hpp:326
ell_matrix.hpp
Implementation of the ell_matrix class.
viennacl::diag
vector_expression< const matrix_base< NumericT >, const int, op_matrix_diag > diag(const matrix_base< NumericT > &A, int k=0)
Definition: matrix.hpp:895
vector_proxy.hpp
Proxy classes for vectors.
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
compressed_compressed_matrix.hpp
Implementation of the compressed_compressed_matrix class (CSR format with a relatively small number o...
viennacl::vector_base< NumericT >::clear
void clear()
Resets all entries to zero. Does not change the size of the vector.
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::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
viennacl::linalg::unit_lower_tag
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
viennacl::compressed_matrix< NumericT >
viennacl::tools::sparse_matrix_adapter
Adapts a non-const sparse matrix type made up from std::vector > to basic u...
Definition: adapter.hpp:357
viennacl::io::read_matrix_market_file
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
Definition: matrix_market.hpp:339
viennacl::vector_base::end
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Definition: vector.hpp:848
common.hpp
Common routines used within ILU-type preconditioners.
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.
resize_test
int resize_test(Epsilon const &epsilon)
Definition: sparse.cpp:213
viennacl::linalg::unit_upper_tag
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864
viennacl::coordinate_matrix
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...
Definition: coordinate_matrix.hpp:186