ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
sparse.cpp
Go to the documentation of this file.
1 /* =========================================================================
2  Copyright (c) 2010-2015, Institute for Microelectronics,
3  Institute for Analysis and Scientific Computing,
4  TU Wien.
5  Portions of this software are copyright by UChicago Argonne, LLC.
6 
7  -----------------
8  ViennaCL - The Vienna Computing Library
9  -----------------
10 
11  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12 
13  (A list of authors and contributors can be found in the PDF manual)
14 
15  License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17 
18 
19 
24 #ifndef NDEBUG
25  #define NDEBUG
26 #endif
27 
28 //
29 // *** System
30 //
31 #include <iostream>
32 
33 //
34 // *** Boost
35 //
36 #include <boost/numeric/ublas/io.hpp>
37 #include <boost/numeric/ublas/triangular.hpp>
38 #include <boost/numeric/ublas/matrix_sparse.hpp>
39 #include <boost/numeric/ublas/matrix.hpp>
40 #include <boost/numeric/ublas/matrix_proxy.hpp>
41 #include <boost/numeric/ublas/lu.hpp>
42 #include <boost/numeric/ublas/io.hpp>
43 #include <boost/numeric/ublas/operation_sparse.hpp>
44 #include <boost/numeric/ublas/vector_proxy.hpp>
45 
46 //
47 // *** ViennaCL
48 //
49 //#define VIENNACL_DEBUG_ALL
50 #define VIENNACL_WITH_UBLAS 1
51 #include "viennacl/scalar.hpp"
52 #include "viennacl/compressed_matrix.hpp"
53 #include "viennacl/compressed_compressed_matrix.hpp"
54 #include "viennacl/coordinate_matrix.hpp"
55 #include "viennacl/ell_matrix.hpp"
56 #include "viennacl/sliced_ell_matrix.hpp"
57 #include "viennacl/hyb_matrix.hpp"
58 #include "viennacl/vector.hpp"
59 #include "viennacl/vector_proxy.hpp"
60 #include "viennacl/linalg/prod.hpp"
61 #include "viennacl/linalg/norm_2.hpp"
62 #include "viennacl/linalg/ilu.hpp"
63 #include "viennacl/linalg/detail/ilu/common.hpp"
64 #include "viennacl/io/matrix_market.hpp"
65 #include "viennacl/tools/random.hpp"
66 
67 
68 //
69 // -------------------------------------------------------------
70 //
71 using namespace boost::numeric;
72 //
73 // -------------------------------------------------------------
74 //
75 template<typename ScalarType>
76 ScalarType diff(ScalarType & s1, viennacl::scalar<ScalarType> & s2)
77 {
78  if (s1 != s2)
79  return (s1 - s2) / std::max(fabs(s1), std::fabs(s2));
80  return 0;
81 }
82 
83 template<typename ScalarType>
84 ScalarType diff(ublas::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
85 {
86  ublas::vector<ScalarType> v2_cpu(v2.size());
87  viennacl::backend::finish();
88  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
89 
90  for (unsigned int i=0;i<v1.size(); ++i)
91  {
92  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
93  {
94  //if (std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) < 1e-10 ) //absolute tolerance (avoid round-off issues)
95  // v2_cpu[i] = 0;
96  //else
97  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
98  }
99  else
100  v2_cpu[i] = 0.0;
101 
102  if (v2_cpu[i] > 0.0001)
103  {
104  //std::cout << "Neighbor: " << i-1 << ": " << v1[i-1] << " vs. " << v2_cpu[i-1] << std::endl;
105  std::cout << "Error at entry " << i << ": Should: " << v1[i] << " vs. Is: " << v2[i] << std::endl;
106  //std::cout << "Neighbor: " << i+1 << ": " << v1[i+1] << " vs. " << v2_cpu[i+1] << std::endl;
107  exit(EXIT_FAILURE);
108  }
109  }
110 
111  return norm_inf(v2_cpu);
112 }
113 
114 
115 template<typename ScalarType, typename VCL_MATRIX>
116 ScalarType diff(ublas::compressed_matrix<ScalarType> & cpu_matrix, VCL_MATRIX & gpu_matrix)
117 {
118  typedef ublas::compressed_matrix<ScalarType> CPU_MATRIX;
119  CPU_MATRIX from_gpu(gpu_matrix.size1(), gpu_matrix.size2());
120 
121  viennacl::backend::finish();
122  viennacl::copy(gpu_matrix, from_gpu);
123 
124  ScalarType error = 0;
125 
126  //step 1: compare all entries from cpu_matrix with gpu_matrix:
127  //std::cout << "Ublas matrix: " << std::endl;
128  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
129  row_it != cpu_matrix.end1();
130  ++row_it)
131  {
132  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
133  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
134  col_it != row_it.end();
135  ++col_it)
136  {
137  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
138  ScalarType current_error = 0;
139 
140  if ( std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
141  std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
142  current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
143  / std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
144  std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
145  if (current_error > error)
146  error = current_error;
147  }
148  }
149 
150  //step 2: compare all entries from gpu_matrix with cpu_matrix (sparsity pattern might differ):
151  //std::cout << "ViennaCL matrix: " << std::endl;
152  for (typename CPU_MATRIX::const_iterator1 row_it = from_gpu.begin1();
153  row_it != from_gpu.end1();
154  ++row_it)
155  {
156  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
157  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
158  col_it != row_it.end();
159  ++col_it)
160  {
161  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
162  ScalarType current_error = 0;
163 
164  if ( std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
165  std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
166  current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
167  / std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
168  std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
169  if (current_error > error)
170  error = current_error;
171  }
172  }
173 
174  return error;
175 }
176 
177 
178 template<typename NumericT, typename VCL_MatrixT, typename Epsilon, typename UblasVectorT, typename VCLVectorT>
179 int strided_matrix_vector_product_test(Epsilon epsilon,
180  UblasVectorT & result, UblasVectorT const & rhs,
181  VCLVectorT & vcl_result, VCLVectorT & vcl_rhs)
182 {
183  int retval = EXIT_SUCCESS;
184 
185  ublas::compressed_matrix<NumericT> ublas_matrix2(5, 4);
186  ublas_matrix2(0, 0) = NumericT(2.0); ublas_matrix2(0, 2) = NumericT(-1.0);
187  ublas_matrix2(1, 0) = NumericT(3.0); ublas_matrix2(1, 2) = NumericT(-5.0);
188  ublas_matrix2(2, 1) = NumericT(5.0); ublas_matrix2(2, 2) = NumericT(-2.0);
189  ublas_matrix2(3, 2) = NumericT(1.0); ublas_matrix2(3, 3) = NumericT(-6.0);
190  ublas_matrix2(4, 1) = NumericT(7.0); ublas_matrix2(4, 2) = NumericT(-5.0);
191  project(result, ublas::slice(1, 3, 5)) = ublas::prod(ublas_matrix2, project(rhs, ublas::slice(3, 2, 4)));
192 
193  VCL_MatrixT vcl_sparse_matrix2;
194  viennacl::copy(ublas_matrix2, vcl_sparse_matrix2);
195  viennacl::vector<NumericT> vec(4);
196  vec(0) = rhs(3);
197  vec(1) = rhs(5);
198  vec(2) = rhs(7);
199  vec(3) = rhs(9);
200  viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, viennacl::project(vcl_rhs, viennacl::slice(3, 2, 4)));
201 
202  if ( std::fabs(diff(result, vcl_result)) > epsilon )
203  {
204  std::cout << "# Error at operation: matrix-vector product with stided vectors, part 1" << std::endl;
205  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
206  retval = EXIT_FAILURE;
207  }
208  vcl_result(1) = NumericT(1.0);
209  vcl_result(4) = NumericT(1.0);
210  vcl_result(7) = NumericT(1.0);
211  vcl_result(10) = NumericT(1.0);
212  vcl_result(13) = NumericT(1.0);
213 
214  viennacl::project(vcl_result, viennacl::slice(1, 3, 5)) = viennacl::linalg::prod(vcl_sparse_matrix2, vec);
215 
216  if ( std::fabs(diff(result, vcl_result)) > epsilon )
217  {
218  std::cout << "# Error at operation: matrix-vector product with strided vectors, part 2" << std::endl;
219  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
220  retval = EXIT_FAILURE;
221  }
222 
223  return retval;
224 }
225 
226 
227 template< typename NumericT, typename VCL_MATRIX, typename Epsilon >
228 int resize_test(Epsilon const& epsilon)
229 {
230  int retval = EXIT_SUCCESS;
231 
232  ublas::compressed_matrix<NumericT> ublas_matrix(5,5);
233  VCL_MATRIX vcl_matrix;
234 
235  ublas_matrix(0,0) = NumericT(10.0); ublas_matrix(0, 1) = NumericT(0.1); ublas_matrix(0, 2) = NumericT(0.2); ublas_matrix(0, 3) = NumericT(0.3); ublas_matrix(0, 4) = NumericT(0.4);
236  ublas_matrix(1,0) = NumericT(1.0); ublas_matrix(1, 1) = NumericT(1.1); ublas_matrix(1, 2) = NumericT(1.2); ublas_matrix(1, 3) = NumericT(1.3); ublas_matrix(1, 4) = NumericT(1.4);
237  ublas_matrix(2,0) = NumericT(2.0); ublas_matrix(2, 1) = NumericT(2.1); ublas_matrix(2, 2) = NumericT(2.2); ublas_matrix(2, 3) = NumericT(2.3); ublas_matrix(2, 4) = NumericT(2.4);
238  ublas_matrix(3,0) = NumericT(3.0); ublas_matrix(3, 1) = NumericT(3.1); ublas_matrix(3, 2) = NumericT(3.2); ublas_matrix(3, 3) = NumericT(3.3); ublas_matrix(3, 4) = NumericT(3.4);
239  ublas_matrix(4,0) = NumericT(4.0); ublas_matrix(4, 1) = NumericT(4.1); ublas_matrix(4, 2) = NumericT(4.2); ublas_matrix(4, 3) = NumericT(4.3); ublas_matrix(4, 4) = NumericT(4.4);
240 
241  viennacl::copy(ublas_matrix, vcl_matrix);
242  ublas::compressed_matrix<NumericT> other_matrix(ublas_matrix.size1(), ublas_matrix.size2());
243  viennacl::copy(vcl_matrix, other_matrix);
244 
245  std::cout << "Checking for equality after copy..." << std::endl;
246  if ( std::fabs(diff(ublas_matrix, vcl_matrix)) > epsilon )
247  {
248  std::cout << "# Error at operation: equality after copy with sparse matrix" << std::endl;
249  std::cout << " diff: " << std::fabs(diff(ublas_matrix, vcl_matrix)) << std::endl;
250  return EXIT_FAILURE;
251  }
252 
253  std::cout << "Testing resize to larger..." << std::endl;
254  ublas_matrix.resize(10, 10, false); //ublas does not allow preserve = true here
255  ublas_matrix(0,0) = NumericT(10.0); ublas_matrix(0, 1) = NumericT(0.1); ublas_matrix(0, 2) = NumericT(0.2); ublas_matrix(0, 3) = NumericT(0.3); ublas_matrix(0, 4) = NumericT(0.4);
256  ublas_matrix(1,0) = NumericT( 1.0); ublas_matrix(1, 1) = NumericT(1.1); ublas_matrix(1, 2) = NumericT(1.2); ublas_matrix(1, 3) = NumericT(1.3); ublas_matrix(1, 4) = NumericT(1.4);
257  ublas_matrix(2,0) = NumericT( 2.0); ublas_matrix(2, 1) = NumericT(2.1); ublas_matrix(2, 2) = NumericT(2.2); ublas_matrix(2, 3) = NumericT(2.3); ublas_matrix(2, 4) = NumericT(2.4);
258  ublas_matrix(3,0) = NumericT( 3.0); ublas_matrix(3, 1) = NumericT(3.1); ublas_matrix(3, 2) = NumericT(3.2); ublas_matrix(3, 3) = NumericT(3.3); ublas_matrix(3, 4) = NumericT(3.4);
259  ublas_matrix(4,0) = NumericT( 4.0); ublas_matrix(4, 1) = NumericT(4.1); ublas_matrix(4, 2) = NumericT(4.2); ublas_matrix(4, 3) = NumericT(4.3); ublas_matrix(4, 4) = NumericT(4.4);
260  //std::cout << ublas_matrix << std::endl;
261 
262  vcl_matrix.resize(10, 10, true);
263 
264  if ( std::fabs(diff(ublas_matrix, vcl_matrix)) > epsilon )
265  {
266  std::cout << "# Error at operation: resize (to larger) with sparse matrix" << std::endl;
267  std::cout << " diff: " << std::fabs(diff(ublas_matrix, vcl_matrix)) << std::endl;
268  return EXIT_FAILURE;
269  }
270 
271  ublas_matrix(5,5) = NumericT(5.5); ublas_matrix(5, 6) = NumericT(5.6); ublas_matrix(5, 7) = NumericT(5.7); ublas_matrix(5, 8) = NumericT(5.8); ublas_matrix(5, 9) = NumericT(5.9);
272  ublas_matrix(6,5) = NumericT(6.5); ublas_matrix(6, 6) = NumericT(6.6); ublas_matrix(6, 7) = NumericT(6.7); ublas_matrix(6, 8) = NumericT(6.8); ublas_matrix(6, 9) = NumericT(6.9);
273  ublas_matrix(7,5) = NumericT(7.5); ublas_matrix(7, 6) = NumericT(7.6); ublas_matrix(7, 7) = NumericT(7.7); ublas_matrix(7, 8) = NumericT(7.8); ublas_matrix(7, 9) = NumericT(7.9);
274  ublas_matrix(8,5) = NumericT(8.5); ublas_matrix(8, 6) = NumericT(8.6); ublas_matrix(8, 7) = NumericT(8.7); ublas_matrix(8, 8) = NumericT(8.8); ublas_matrix(8, 9) = NumericT(8.9);
275  ublas_matrix(9,5) = NumericT(9.5); ublas_matrix(9, 6) = NumericT(9.6); ublas_matrix(9, 7) = NumericT(9.7); ublas_matrix(9, 8) = NumericT(9.8); ublas_matrix(9, 9) = NumericT(9.9);
276  viennacl::copy(ublas_matrix, vcl_matrix);
277 
278  std::cout << "Testing resize to smaller..." << std::endl;
279  ublas_matrix.resize(7, 7, false); //ublas does not allow preserve = true here
280  ublas_matrix(0,0) = NumericT(10.0); ublas_matrix(0, 1) = NumericT(0.1); ublas_matrix(0, 2) = NumericT(0.2); ublas_matrix(0, 3) = NumericT(0.3); ublas_matrix(0, 4) = NumericT(0.4);
281  ublas_matrix(1,0) = NumericT( 1.0); ublas_matrix(1, 1) = NumericT(1.1); ublas_matrix(1, 2) = NumericT(1.2); ublas_matrix(1, 3) = NumericT(1.3); ublas_matrix(1, 4) = NumericT(1.4);
282  ublas_matrix(2,0) = NumericT( 2.0); ublas_matrix(2, 1) = NumericT(2.1); ublas_matrix(2, 2) = NumericT(2.2); ublas_matrix(2, 3) = NumericT(2.3); ublas_matrix(2, 4) = NumericT(2.4);
283  ublas_matrix(3,0) = NumericT( 3.0); ublas_matrix(3, 1) = NumericT(3.1); ublas_matrix(3, 2) = NumericT(3.2); ublas_matrix(3, 3) = NumericT(3.3); ublas_matrix(3, 4) = NumericT(3.4);
284  ublas_matrix(4,0) = NumericT( 4.0); ublas_matrix(4, 1) = NumericT(4.1); ublas_matrix(4, 2) = NumericT(4.2); ublas_matrix(4, 3) = NumericT(4.3); ublas_matrix(4, 4) = NumericT(4.4);
285  ublas_matrix(5,5) = NumericT( 5.5); ublas_matrix(5, 6) = NumericT(5.6); ublas_matrix(5, 7) = NumericT(5.7); ublas_matrix(5, 8) = NumericT(5.8); ublas_matrix(5, 9) = NumericT(5.9);
286  ublas_matrix(6,5) = NumericT( 6.5); ublas_matrix(6, 6) = NumericT(6.6); ublas_matrix(6, 7) = NumericT(6.7); ublas_matrix(6, 8) = NumericT(6.8); ublas_matrix(6, 9) = NumericT(6.9);
287 
288  vcl_matrix.resize(7, 7);
289 
290  //std::cout << ublas_matrix << std::endl;
291  if ( std::fabs(diff(ublas_matrix, vcl_matrix)) > epsilon )
292  {
293  std::cout << "# Error at operation: resize (to smaller) with sparse matrix" << std::endl;
294  std::cout << " diff: " << std::fabs(diff(ublas_matrix, vcl_matrix)) << std::endl;
295  retval = EXIT_FAILURE;
296  }
297 
298  ublas::vector<NumericT> ublas_vec = ublas::scalar_vector<NumericT>(ublas_matrix.size1(), NumericT(3.1415));
299  viennacl::vector<NumericT> vcl_vec(ublas_matrix.size1());
300 
301 
302  std::cout << "Testing transposed unit lower triangular solve: compressed_matrix" << std::endl;
303  viennacl::copy(ublas_vec, vcl_vec);
304  std::cout << "matrix: " << ublas_matrix << std::endl;
305  std::cout << "vector: " << ublas_vec << std::endl;
306  std::cout << "ViennaCL matrix size: " << vcl_matrix.size1() << " x " << vcl_matrix.size2() << std::endl;
307 
308  std::cout << "ublas..." << std::endl;
309  boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::unit_lower_tag());
310  std::cout << "ViennaCL..." << std::endl;
311  viennacl::linalg::inplace_solve((vcl_matrix), vcl_vec, viennacl::linalg::unit_lower_tag());
312 
313  /*
314  std::list< viennacl::backend::mem_handle > multifrontal_L_row_index_arrays_;
315  std::list< viennacl::backend::mem_handle > multifrontal_L_row_buffers_;
316  std::list< viennacl::backend::mem_handle > multifrontal_L_col_buffers_;
317  std::list< viennacl::backend::mem_handle > multifrontal_L_element_buffers_;
318  std::list< std::size_t > multifrontal_L_row_elimination_num_list_;
319 
320  viennacl::vector<NumericT> multifrontal_U_diagonal_;
321 
322  viennacl::linalg::detail::multifrontal_setup_L(vcl_matrix,
323  multifrontal_U_diagonal_, //dummy
324  multifrontal_L_row_index_arrays_,
325  multifrontal_L_row_buffers_,
326  multifrontal_L_col_buffers_,
327  multifrontal_L_element_buffers_,
328  multifrontal_L_row_elimination_num_list_);
329 
330  viennacl::linalg::detail::multifrontal_substitute(vcl_vec,
331  multifrontal_L_row_index_arrays_,
332  multifrontal_L_row_buffers_,
333  multifrontal_L_col_buffers_,
334  multifrontal_L_element_buffers_,
335  multifrontal_L_row_elimination_num_list_);
336 
337 
338  std::cout << "ublas..." << std::endl;
339  boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::upper_tag());
340  std::cout << "ViennaCL..." << std::endl;
341  std::list< viennacl::backend::mem_handle > multifrontal_U_row_index_arrays_;
342  std::list< viennacl::backend::mem_handle > multifrontal_U_row_buffers_;
343  std::list< viennacl::backend::mem_handle > multifrontal_U_col_buffers_;
344  std::list< viennacl::backend::mem_handle > multifrontal_U_element_buffers_;
345  std::list< std::size_t > multifrontal_U_row_elimination_num_list_;
346 
347  multifrontal_U_diagonal_.resize(vcl_matrix.size1(), false);
348  viennacl::linalg::single_threaded::detail::row_info(vcl_matrix, multifrontal_U_diagonal_, viennacl::linalg::detail::SPARSE_ROW_DIAGONAL);
349  viennacl::linalg::detail::multifrontal_setup_U(vcl_matrix,
350  multifrontal_U_diagonal_,
351  multifrontal_U_row_index_arrays_,
352  multifrontal_U_row_buffers_,
353  multifrontal_U_col_buffers_,
354  multifrontal_U_element_buffers_,
355  multifrontal_U_row_elimination_num_list_);
356 
357  vcl_vec = viennacl::linalg::element_div(vcl_vec, multifrontal_U_diagonal_);
358  viennacl::linalg::detail::multifrontal_substitute(vcl_vec,
359  multifrontal_U_row_index_arrays_,
360  multifrontal_U_row_buffers_,
361  multifrontal_U_col_buffers_,
362  multifrontal_U_element_buffers_,
363  multifrontal_U_row_elimination_num_list_);
364  */
365  for (std::size_t i=0; i<ublas_vec.size(); ++i)
366  {
367  std::cout << ublas_vec[i] << " vs. " << vcl_vec[i] << std::endl;
368  }
369 
370  /*std::cout << "Testing transposed unit upper triangular solve: compressed_matrix" << std::endl;
371  viennacl::copy(ublas_vec, vcl_vec);
372  std::cout << "matrix: " << ublas_matrix << std::endl;
373  std::cout << "vector: " << ublas_vec << std::endl;
374  std::cout << "ViennaCL matrix size: " << vcl_matrix.size1() << " x " << vcl_matrix.size2() << std::endl;
375 
376  std::cout << "ublas..." << std::endl;
377  boost::numeric::ublas::inplace_solve((ublas_matrix), ublas_vec, boost::numeric::ublas::lower_tag());
378  std::cout << "ViennaCL..." << std::endl;
379  viennacl::linalg::inplace_solve((vcl_matrix), vcl_vec, viennacl::linalg::lower_tag());
380 
381  for (std::size_t i=0; i<ublas_vec.size(); ++i)
382  {
383  std::cout << ublas_vec[i] << " vs. " << vcl_vec[i] << std::endl;
384  }*/
385 
386  return retval;
387 }
388 
389 
390 //
391 // -------------------------------------------------------------
392 //
393 template< typename NumericT, typename Epsilon >
394 int test(Epsilon const& epsilon)
395 {
396  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
397 
398  std::cout << "Testing resizing of compressed_matrix..." << std::endl;
399  int retval = resize_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon);
400  if (retval != EXIT_SUCCESS)
401  return retval;
402  std::cout << "Testing resizing of coordinate_matrix..." << std::endl;
403  //if (retval != EXIT_FAILURE)
404  // retval = resize_test<NumericT, viennacl::coordinate_matrix<NumericT> >(epsilon);
405  //else
406  // return retval;
407 
408  // --------------------------------------------------------------------------
409  ublas::vector<NumericT> rhs;
410  ublas::vector<NumericT> result;
411  ublas::compressed_matrix<NumericT> ublas_matrix;
412 
413  if (viennacl::io::read_matrix_market_file(ublas_matrix, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
414  {
415  std::cout << "Error reading Matrix file" << std::endl;
416  return EXIT_FAILURE;
417  }
418 
419  //unsigned int cg_mat_size = cg_mat.size();
420  std::cout << "done reading matrix" << std::endl;
421 
422 
423  rhs.resize(ublas_matrix.size2());
424  for (std::size_t i=0; i<rhs.size(); ++i)
425  {
426  ublas_matrix(i,i) = NumericT(0.5); // Get rid of round-off errors by making row-sums unequal to zero:
427  rhs[i] = NumericT(1) + randomNumber();
428  }
429 
430  // add some random numbers to the double-compressed matrix:
431  ublas::compressed_matrix<NumericT> ublas_cc_matrix(ublas_matrix.size1(), ublas_matrix.size2());
432  ublas_cc_matrix(42,199) = NumericT(3.1415);
433  ublas_cc_matrix(31, 69) = NumericT(2.71);
434  ublas_cc_matrix(23, 32) = NumericT(6);
435  ublas_cc_matrix(177,57) = NumericT(4);
436  ublas_cc_matrix(21, 97) = NumericT(-4);
437  ublas_cc_matrix(92, 25) = NumericT(2);
438  ublas_cc_matrix(89, 62) = NumericT(11);
439  ublas_cc_matrix(1, 7) = NumericT(8);
440  ublas_cc_matrix(85, 41) = NumericT(13);
441  ublas_cc_matrix(66, 28) = NumericT(8);
442  ublas_cc_matrix(21, 74) = NumericT(-2);
443 
444 
445  result = rhs;
446 
447 
448  viennacl::vector<NumericT> vcl_rhs(rhs.size());
449  viennacl::vector<NumericT> vcl_result(result.size());
450  viennacl::vector<NumericT> vcl_result2(result.size());
451  viennacl::compressed_matrix<NumericT> vcl_compressed_matrix(rhs.size(), rhs.size());
452  viennacl::compressed_compressed_matrix<NumericT> vcl_compressed_compressed_matrix(rhs.size(), rhs.size());
453  viennacl::coordinate_matrix<NumericT> vcl_coordinate_matrix(rhs.size(), rhs.size());
454  viennacl::ell_matrix<NumericT> vcl_ell_matrix;
455  viennacl::sliced_ell_matrix<NumericT> vcl_sliced_ell_matrix;
456  viennacl::hyb_matrix<NumericT> vcl_hyb_matrix;
457 
458  viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
459  viennacl::copy(ublas_matrix, vcl_compressed_matrix);
460  viennacl::copy(ublas_cc_matrix, vcl_compressed_compressed_matrix);
461  viennacl::copy(ublas_matrix, vcl_coordinate_matrix);
462 
463  // --------------------------------------------------------------------------
464  std::cout << "Testing products: ublas" << std::endl;
465  result = viennacl::linalg::prod(ublas_matrix, rhs);
466 
467  std::cout << "Testing products: compressed_matrix" << std::endl;
468  vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
469 
470  if ( std::fabs(diff(result, vcl_result)) > epsilon )
471  {
472  std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
473  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
474  retval = EXIT_FAILURE;
475  }
476 
477  std::cout << "Testing products: compressed_matrix, strided vectors" << std::endl;
478  retval = strided_matrix_vector_product_test<NumericT, viennacl::compressed_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
479  if (retval != EXIT_SUCCESS)
480  return retval;
481 
482  //
483  // Triangular solvers for A \ b:
484  //
485  ublas::compressed_matrix<NumericT> ublas_matrix_trans(ublas_matrix.size2(), ublas_matrix.size1(), ublas_matrix.nnz()); // = trans(ublas_matrix); //note: triangular solvers with uBLAS show atrocious performance, while transposed solvers are quite okay. To keep execution times short, we use a double-transpose-trick in the following.
486 
487  // fast transpose:
488  for (typename ublas::compressed_matrix<NumericT>::iterator1 row_it = ublas_matrix.begin1();
489  row_it != ublas_matrix.end1();
490  ++row_it)
491  {
492  for (typename ublas::compressed_matrix<NumericT>::iterator2 col_it = row_it.begin();
493  col_it != row_it.end();
494  ++col_it)
495  {
496  ublas_matrix_trans(col_it.index1(), col_it.index2()) = *col_it;
497  }
498  }
499 
500 
501  std::cout << "Testing unit upper triangular solve: compressed_matrix" << std::endl;
502  result = rhs;
503  viennacl::copy(result, vcl_result);
504  boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_upper_tag());
505  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_upper_tag());
506 
507  if ( std::fabs(diff(result, vcl_result)) > epsilon )
508  {
509  std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
510  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
511  retval = EXIT_FAILURE;
512  }
513 
514  std::cout << "Testing upper triangular solve: compressed_matrix" << std::endl;
515  result = rhs;
516  viennacl::copy(result, vcl_result);
517  boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::upper_tag());
518  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::upper_tag());
519 
520  if ( std::fabs(diff(result, vcl_result)) > epsilon )
521  {
522  std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
523  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
524  retval = EXIT_FAILURE;
525  }
526 
527  std::cout << "Testing unit lower triangular solve: compressed_matrix" << std::endl;
528  result = rhs;
529  viennacl::copy(result, vcl_result);
530  boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::unit_lower_tag());
531  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::unit_lower_tag());
532 
533  /*std::list< viennacl::backend::mem_handle > multifrontal_L_row_index_arrays_;
534  std::list< viennacl::backend::mem_handle > multifrontal_L_row_buffers_;
535  std::list< viennacl::backend::mem_handle > multifrontal_L_col_buffers_;
536  std::list< viennacl::backend::mem_handle > multifrontal_L_element_buffers_;
537  std::list< std::size_t > multifrontal_L_row_elimination_num_list_;
538 
539  viennacl::vector<NumericT> multifrontal_U_diagonal_;
540 
541  viennacl::switch_memory_domain(multifrontal_U_diagonal_, viennacl::MAIN_MEMORY);
542  multifrontal_U_diagonal_.resize(vcl_compressed_matrix.size1(), false);
543  viennacl::linalg::single_threaded::detail::row_info(vcl_compressed_matrix, multifrontal_U_diagonal_, viennacl::linalg::detail::SPARSE_ROW_DIAGONAL);
544 
545  viennacl::linalg::detail::multifrontal_setup_L(vcl_compressed_matrix,
546  multifrontal_U_diagonal_, //dummy
547  multifrontal_L_row_index_arrays_,
548  multifrontal_L_row_buffers_,
549  multifrontal_L_col_buffers_,
550  multifrontal_L_element_buffers_,
551  multifrontal_L_row_elimination_num_list_);
552 
553  viennacl::linalg::detail::multifrontal_substitute(vcl_result,
554  multifrontal_L_row_index_arrays_,
555  multifrontal_L_row_buffers_,
556  multifrontal_L_col_buffers_,
557  multifrontal_L_element_buffers_,
558  multifrontal_L_row_elimination_num_list_);*/
559 
560 
561  if ( std::fabs(diff(result, vcl_result)) > epsilon )
562  {
563  std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
564  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
565  retval = EXIT_FAILURE;
566  }
567 
568 
569  std::cout << "Testing lower triangular solve: compressed_matrix" << std::endl;
570  result = rhs;
571  viennacl::copy(result, vcl_result);
572  boost::numeric::ublas::inplace_solve(trans(ublas_matrix_trans), result, boost::numeric::ublas::lower_tag());
573  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::lower_tag());
574 
575  /*std::list< viennacl::backend::mem_handle > multifrontal_U_row_index_arrays_;
576  std::list< viennacl::backend::mem_handle > multifrontal_U_row_buffers_;
577  std::list< viennacl::backend::mem_handle > multifrontal_U_col_buffers_;
578  std::list< viennacl::backend::mem_handle > multifrontal_U_element_buffers_;
579  std::list< std::size_t > multifrontal_U_row_elimination_num_list_;
580 
581  multifrontal_U_diagonal_.resize(vcl_compressed_matrix.size1(), false);
582  viennacl::linalg::single_threaded::detail::row_info(vcl_compressed_matrix, multifrontal_U_diagonal_, viennacl::linalg::detail::SPARSE_ROW_DIAGONAL);
583  viennacl::linalg::detail::multifrontal_setup_U(vcl_compressed_matrix,
584  multifrontal_U_diagonal_,
585  multifrontal_U_row_index_arrays_,
586  multifrontal_U_row_buffers_,
587  multifrontal_U_col_buffers_,
588  multifrontal_U_element_buffers_,
589  multifrontal_U_row_elimination_num_list_);
590 
591  vcl_result = viennacl::linalg::element_div(vcl_result, multifrontal_U_diagonal_);
592  viennacl::linalg::detail::multifrontal_substitute(vcl_result,
593  multifrontal_U_row_index_arrays_,
594  multifrontal_U_row_buffers_,
595  multifrontal_U_col_buffers_,
596  multifrontal_U_element_buffers_,
597  multifrontal_U_row_elimination_num_list_);*/
598 
599 
600  if ( std::fabs(diff(result, vcl_result)) > epsilon )
601  {
602  std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
603  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
604  retval = EXIT_FAILURE;
605  }
606 
607 /*
608  std::cout << "Testing lower triangular solve: compressed_matrix" << std::endl;
609  result = rhs;
610  viennacl::copy(result, vcl_result);
611  boost::numeric::ublas::inplace_solve(ublas_matrix, result, boost::numeric::ublas::lower_tag());
612  viennacl::linalg::inplace_solve(vcl_compressed_matrix, vcl_result, viennacl::linalg::lower_tag());
613 
614  if ( std::fabs(diff(result, vcl_result)) > epsilon )
615  {
616  std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
617  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
618  retval = EXIT_FAILURE;
619  }*/
620 
621  //
622  // Triangular solvers for A^T \ b
623  //
624 
625  std::cout << "Testing transposed unit upper triangular solve: compressed_matrix" << std::endl;
626  result = rhs;
627  viennacl::copy(result, vcl_result);
628  boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_upper_tag());
629  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_upper_tag());
630 
631  if ( std::fabs(diff(result, vcl_result)) > epsilon )
632  {
633  std::cout << "# Error at operation: unit upper triangular solve with compressed_matrix" << std::endl;
634  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
635  retval = EXIT_FAILURE;
636  }
637 
638  std::cout << "Testing transposed upper triangular solve: compressed_matrix" << std::endl;
639  result = rhs;
640  viennacl::copy(result, vcl_result);
641  boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::upper_tag());
642  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::upper_tag());
643 
644  if ( std::fabs(diff(result, vcl_result)) > epsilon )
645  {
646  std::cout << "# Error at operation: upper triangular solve with compressed_matrix" << std::endl;
647  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
648  retval = EXIT_FAILURE;
649  }
650 
651 
652  std::cout << "Testing transposed unit lower triangular solve: compressed_matrix" << std::endl;
653  result = rhs;
654  viennacl::copy(result, vcl_result);
655  boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::unit_lower_tag());
656  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::unit_lower_tag());
657 
658  if ( std::fabs(diff(result, vcl_result)) > epsilon )
659  {
660  std::cout << "# Error at operation: unit lower triangular solve with compressed_matrix" << std::endl;
661  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
662  retval = EXIT_FAILURE;
663  }
664 
665  std::cout << "Testing transposed lower triangular solve: compressed_matrix" << std::endl;
666  result = rhs;
667  viennacl::copy(result, vcl_result);
668  boost::numeric::ublas::inplace_solve(trans(ublas_matrix), result, boost::numeric::ublas::lower_tag());
669  viennacl::linalg::inplace_solve(trans(vcl_compressed_matrix), vcl_result, viennacl::linalg::lower_tag());
670 
671  if ( std::fabs(diff(result, vcl_result)) > epsilon )
672  {
673  std::cout << "# Error at operation: lower triangular solve with compressed_matrix" << std::endl;
674  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
675  retval = EXIT_FAILURE;
676  }
677 
678  std::cout << "Testing products: compressed_compressed_matrix" << std::endl;
679  result = viennacl::linalg::prod(ublas_cc_matrix, rhs);
680  vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
681 
682  if ( std::fabs(diff(result, vcl_result)) > epsilon )
683  {
684  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix" << std::endl;
685  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
686  retval = EXIT_FAILURE;
687  }
688 
689  {
690  ublas::compressed_matrix<NumericT> temp(vcl_compressed_compressed_matrix.size1(), vcl_compressed_compressed_matrix.size2());
691  viennacl::copy(vcl_compressed_compressed_matrix, temp);
692 
693  // check that entries are correct by computing the product again:
694  result = viennacl::linalg::prod(temp, rhs);
695 
696  if ( std::fabs(diff(result, vcl_result)) > epsilon )
697  {
698  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix (after copy back)" << std::endl;
699  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
700  retval = EXIT_FAILURE;
701  }
702 
703  }
704 
705 
706 
707 
708  std::cout << "Testing products: coordinate_matrix" << std::endl;
709  result = viennacl::linalg::prod(ublas_matrix, rhs);
710  vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
711 
712  if ( std::fabs(diff(result, vcl_result)) > epsilon )
713  {
714  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
715  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
716  retval = EXIT_FAILURE;
717  }
718 
719  std::cout << "Testing products: coordinate_matrix, strided vectors" << std::endl;
720  //std::cout << " --> SKIPPING <--" << std::endl;
721  retval = strided_matrix_vector_product_test<NumericT, viennacl::coordinate_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
722  if (retval != EXIT_SUCCESS)
723  return retval;
724 
725 
726  //std::cout << "Copying ell_matrix" << std::endl;
727  viennacl::copy(ublas_matrix, vcl_ell_matrix);
728  ublas_matrix.clear();
729  viennacl::copy(vcl_ell_matrix, ublas_matrix);// just to check that it's works
730 
731 
732  std::cout << "Testing products: ell_matrix" << std::endl;
733  result = viennacl::linalg::prod(ublas_matrix, rhs);
734  vcl_result.clear();
735  vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
736  //viennacl::linalg::prod_impl(vcl_ell_matrix, vcl_rhs, vcl_result);
737  //std::cout << vcl_result << "\n";
738  //std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
739  //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
740 
741  if ( std::fabs(diff(result, vcl_result)) > epsilon )
742  {
743  std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
744  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
745  retval = EXIT_FAILURE;
746  }
747 
748  std::cout << "Testing products: ell_matrix, strided vectors" << std::endl;
749  retval = strided_matrix_vector_product_test<NumericT, viennacl::ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
750  if (retval != EXIT_SUCCESS)
751  return retval;
752 
753  //std::cout << "Copying sliced_ell_matrix" << std::endl;
754  viennacl::copy(ublas_matrix, vcl_sliced_ell_matrix);
755  //ublas_matrix.clear();
756  //viennacl::copy(vcl_hyb_matrix, ublas_matrix);// just to check that it's works
757  //viennacl::copy(ublas_matrix, vcl_hyb_matrix);
758 
759  std::cout << "Testing products: sliced_ell_matrix" << std::endl;
760  result = viennacl::linalg::prod(ublas_matrix, rhs);
761  vcl_result.clear();
762  vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
763 
764  if ( std::fabs(diff(result, vcl_result)) > epsilon )
765  {
766  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
767  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
768  retval = EXIT_FAILURE;
769  }
770 
771  std::cout << "Testing products: sliced_ell_matrix, strided vectors" << std::endl;
772  retval = strided_matrix_vector_product_test<NumericT, viennacl::sliced_ell_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
773  if (retval != EXIT_SUCCESS)
774  return retval;
775 
776 
777  //std::cout << "Copying hyb_matrix" << std::endl;
778  viennacl::copy(ublas_matrix, vcl_hyb_matrix);
779  ublas_matrix.clear();
780  viennacl::copy(vcl_hyb_matrix, ublas_matrix);// just to check that it's works
781  viennacl::copy(ublas_matrix, vcl_hyb_matrix);
782 
783  std::cout << "Testing products: hyb_matrix" << std::endl;
784  result = viennacl::linalg::prod(ublas_matrix, rhs);
785  vcl_result.clear();
786  vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
787  //viennacl::linalg::prod_impl(vcl_hyb_matrix, vcl_rhs, vcl_result);
788  //std::cout << vcl_result << "\n";
789  //std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
790  //std::cout << "First entry of result vector: " << vcl_result[0] << std::endl;
791 
792  if ( std::fabs(diff(result, vcl_result)) > epsilon )
793  {
794  std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
795  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
796  retval = EXIT_FAILURE;
797  }
798 
799  std::cout << "Testing products: hyb_matrix, strided vectors" << std::endl;
800  retval = strided_matrix_vector_product_test<NumericT, viennacl::hyb_matrix<NumericT> >(epsilon, result, rhs, vcl_result, vcl_rhs);
801  if (retval != EXIT_SUCCESS)
802  return retval;
803 
804 
805  // --------------------------------------------------------------------------
806  // --------------------------------------------------------------------------
807  NumericT alpha = static_cast<NumericT>(2.786);
808  NumericT beta = static_cast<NumericT>(1.432);
809  copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
810  copy(result.begin(), result.end(), vcl_result.begin());
811  copy(result.begin(), result.end(), vcl_result2.begin());
812 
813  std::cout << "Testing scaled additions of products and vectors" << std::endl;
814  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
815  vcl_result2 = alpha * viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) + beta * vcl_result;
816 
817  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
818  {
819  std::cout << "# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
820  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
821  retval = EXIT_FAILURE;
822  }
823 
824 
825  vcl_result2.clear();
826  vcl_result2 = alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result;
827 
828  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
829  {
830  std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
831  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
832  retval = EXIT_FAILURE;
833  }
834 
835  vcl_result2.clear();
836  vcl_result2 = alpha * viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs) + beta * vcl_result;
837 
838  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
839  {
840  std::cout << "# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
841  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
842  retval = EXIT_FAILURE;
843  }
844 
845  vcl_result2.clear();
846  vcl_result2 = alpha * viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs) + beta * vcl_result;
847 
848  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
849  {
850  std::cout << "# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
851  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
852  retval = EXIT_FAILURE;
853  }
854 
856  ublas_matrix.clear();
857 
858  std::cout << "Testing products after clear(): compressed_matrix" << std::endl;
859  vcl_compressed_matrix.clear();
860  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
861  result = viennacl::linalg::prod(ublas_matrix, rhs);
862  vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
863 
864  if ( std::fabs(diff(result, vcl_result)) > epsilon )
865  {
866  std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
867  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
868  retval = EXIT_FAILURE;
869  }
870 
871  std::cout << "Testing products after clear(): compressed_compressed_matrix" << std::endl;
872  vcl_compressed_compressed_matrix.clear();
873  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
874  result = viennacl::linalg::prod(ublas_matrix, rhs);
875  vcl_result = viennacl::linalg::prod(vcl_compressed_compressed_matrix, vcl_rhs);
876 
877  if ( std::fabs(diff(result, vcl_result)) > epsilon )
878  {
879  std::cout << "# Error at operation: matrix-vector product with compressed_compressed_matrix" << std::endl;
880  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
881  retval = EXIT_FAILURE;
882  }
883 
884  std::cout << "Testing products after clear(): coordinate_matrix" << std::endl;
885  vcl_coordinate_matrix.clear();
886  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
887  result = viennacl::linalg::prod(ublas_matrix, rhs);
888  vcl_result = viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs);
889 
890  if ( std::fabs(diff(result, vcl_result)) > epsilon )
891  {
892  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
893  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
894  retval = EXIT_FAILURE;
895  }
896 
897  std::cout << "Testing products after clear(): ell_matrix" << std::endl;
898  vcl_ell_matrix.clear();
899  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
900  result = viennacl::linalg::prod(ublas_matrix, rhs);
901  vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
902 
903  if ( std::fabs(diff(result, vcl_result)) > epsilon )
904  {
905  std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
906  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
907  retval = EXIT_FAILURE;
908  }
909 
910  std::cout << "Testing products after clear(): hyb_matrix" << std::endl;
911  vcl_hyb_matrix.clear();
912  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
913  result = viennacl::linalg::prod(ublas_matrix, rhs);
914  vcl_result = viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs);
915 
916  if ( std::fabs(diff(result, vcl_result)) > epsilon )
917  {
918  std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
919  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
920  retval = EXIT_FAILURE;
921  }
922 
923  std::cout << "Testing products after clear(): sliced_ell_matrix" << std::endl;
924  vcl_sliced_ell_matrix.clear();
925  result = ublas::scalar_vector<NumericT>(result.size(), NumericT(1));
926  result = viennacl::linalg::prod(ublas_matrix, rhs);
927  vcl_result = viennacl::linalg::prod(vcl_sliced_ell_matrix, vcl_rhs);
928 
929  if ( std::fabs(diff(result, vcl_result)) > epsilon )
930  {
931  std::cout << "# Error at operation: matrix-vector product with sliced_ell_matrix" << std::endl;
932  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
933  retval = EXIT_FAILURE;
934  }
935 
936 
937  // --------------------------------------------------------------------------
938  return retval;
939 }
940 //
941 // -------------------------------------------------------------
942 //
943 int main()
944 {
945  std::cout << std::endl;
946  std::cout << "----------------------------------------------" << std::endl;
947  std::cout << "----------------------------------------------" << std::endl;
948  std::cout << "## Test :: Sparse Matrices" << std::endl;
949  std::cout << "----------------------------------------------" << std::endl;
950  std::cout << "----------------------------------------------" << std::endl;
951  std::cout << std::endl;
952 
953  int retval = EXIT_SUCCESS;
954 
955  std::cout << std::endl;
956  std::cout << "----------------------------------------------" << std::endl;
957  std::cout << std::endl;
958  {
959  typedef float NumericT;
960  NumericT epsilon = static_cast<NumericT>(1E-4);
961  std::cout << "# Testing setup:" << std::endl;
962  std::cout << " eps: " << epsilon << std::endl;
963  std::cout << " numeric: float" << std::endl;
964  retval = test<NumericT>(epsilon);
965  if ( retval == EXIT_SUCCESS )
966  std::cout << "# Test passed" << std::endl;
967  else
968  return retval;
969  }
970  std::cout << std::endl;
971  std::cout << "----------------------------------------------" << std::endl;
972  std::cout << std::endl;
973 
974 #ifdef VIENNACL_WITH_OPENCL
975  if ( viennacl::ocl::current_device().double_support() )
976 #endif
977  {
978  {
979  typedef double NumericT;
980  NumericT epsilon = 1.0E-12;
981  std::cout << "# Testing setup:" << std::endl;
982  std::cout << " eps: " << epsilon << std::endl;
983  std::cout << " numeric: double" << std::endl;
984  retval = test<NumericT>(epsilon);
985  if ( retval == EXIT_SUCCESS )
986  std::cout << "# Test passed" << std::endl;
987  else
988  return retval;
989  }
990  std::cout << std::endl;
991  std::cout << "----------------------------------------------" << std::endl;
992  std::cout << std::endl;
993  }
994 #ifdef VIENNACL_WITH_OPENCL
995  else
996  std::cout << "No double precision support, skipping test..." << std::endl;
997 #endif
998 
999 
1000  std::cout << std::endl;
1001  std::cout << "------- Test completed --------" << std::endl;
1002  std::cout << std::endl;
1003 
1004  return retval;
1005 }
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::cuda::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:253
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.
diff
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
Definition: sparse.cpp:76
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:394
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...
main
int main()
Definition: sparse.cpp:943
viennacl::linalg::lower_tag
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
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::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:293
ell_matrix.hpp
Implementation of the ell_matrix class.
vector_proxy.hpp
Proxy classes for vectors.
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
viennacl::slice
basic_slice slice
Definition: forwards.h:429
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::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.
strided_matrix_vector_product_test
int strided_matrix_vector_product_test(Epsilon epsilon, UblasVectorT &result, UblasVectorT const &rhs, VCLVectorT &vcl_result, VCLVectorT &vcl_rhs)
Definition: sparse.cpp:179
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::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:228
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:174