ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
scheduler_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 
23 #ifndef NDEBUG
24  #define NDEBUG
25 #endif
26 
27 //
28 // *** System
29 //
30 #include <iostream>
31 
32 //
33 // *** Boost
34 //
35 #include <boost/numeric/ublas/io.hpp>
36 #include <boost/numeric/ublas/triangular.hpp>
37 #include <boost/numeric/ublas/matrix_sparse.hpp>
38 #include <boost/numeric/ublas/matrix.hpp>
39 #include <boost/numeric/ublas/matrix_proxy.hpp>
40 #include <boost/numeric/ublas/lu.hpp>
41 #include <boost/numeric/ublas/io.hpp>
42 #include <boost/numeric/ublas/operation_sparse.hpp>
43 
44 //
45 // *** ViennaCL
46 //
47 //#define VIENNACL_DEBUG_ALL
48 #define VIENNACL_WITH_UBLAS 1
49 #include "viennacl/scalar.hpp"
52 #include "viennacl/ell_matrix.hpp"
53 #include "viennacl/hyb_matrix.hpp"
54 #include "viennacl/vector.hpp"
56 #include "viennacl/linalg/prod.hpp"
58 #include "viennacl/linalg/ilu.hpp"
62 
65 
66 //
67 // -------------------------------------------------------------
68 //
69 using namespace boost::numeric;
70 //
71 // -------------------------------------------------------------
72 //
73 template<typename ScalarType>
75 {
76  if (s1 != s2)
77  return (s1 - s2) / std::max(fabs(s1), std::fabs(s2));
78  return 0;
79 }
80 
81 template<typename ScalarType>
82 ScalarType diff(ublas::vector<ScalarType> & v1, viennacl::vector<ScalarType> & v2)
83 {
84  ublas::vector<ScalarType> v2_cpu(v2.size());
86  viennacl::copy(v2.begin(), v2.end(), v2_cpu.begin());
87 
88  for (unsigned int i=0;i<v1.size(); ++i)
89  {
90  if ( std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) > 0 )
91  {
92  //if (std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) ) < 1e-10 ) //absolute tolerance (avoid round-off issues)
93  // v2_cpu[i] = 0;
94  //else
95  v2_cpu[i] = std::fabs(v2_cpu[i] - v1[i]) / std::max( std::fabs(v2_cpu[i]), std::fabs(v1[i]) );
96  }
97  else
98  v2_cpu[i] = 0.0;
99 
100  if (v2_cpu[i] > 0.0001)
101  {
102  //std::cout << "Neighbor: " << i-1 << ": " << v1[i-1] << " vs. " << v2_cpu[i-1] << std::endl;
103  std::cout << "Error at entry " << i << ": " << v1[i] << " vs. " << v2_cpu[i] << std::endl;
104  //std::cout << "Neighbor: " << i+1 << ": " << v1[i+1] << " vs. " << v2_cpu[i+1] << std::endl;
105  exit(EXIT_FAILURE);
106  }
107  }
108 
109  return norm_inf(v2_cpu);
110 }
111 
112 
113 template<typename ScalarType, typename VCL_MATRIX>
114 ScalarType diff(ublas::compressed_matrix<ScalarType> & cpu_matrix, VCL_MATRIX & gpu_matrix)
115 {
116  typedef ublas::compressed_matrix<ScalarType> CPU_MATRIX;
117  CPU_MATRIX from_gpu;
118 
120  viennacl::copy(gpu_matrix, from_gpu);
121 
122  ScalarType error = 0;
123 
124  //step 1: compare all entries from cpu_matrix with gpu_matrix:
125  //std::cout << "Ublas matrix: " << std::endl;
126  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
127  row_it != cpu_matrix.end1();
128  ++row_it)
129  {
130  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
131  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
132  col_it != row_it.end();
133  ++col_it)
134  {
135  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
136  ScalarType current_error = 0;
137 
138  if ( std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
139  std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
140  current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
141  / std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
142  std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
143  if (current_error > error)
144  error = current_error;
145  }
146  }
147 
148  //step 2: compare all entries from gpu_matrix with cpu_matrix (sparsity pattern might differ):
149  //std::cout << "ViennaCL matrix: " << std::endl;
150  for (typename CPU_MATRIX::const_iterator1 row_it = from_gpu.begin1();
151  row_it != from_gpu.end1();
152  ++row_it)
153  {
154  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
155  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
156  col_it != row_it.end();
157  ++col_it)
158  {
159  //std::cout << "(" << col_it.index2() << ", " << *col_it << std::endl;
160  ScalarType current_error = 0;
161 
162  if ( std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
163  std::fabs(from_gpu(col_it.index1(), col_it.index2())) ) > 0 )
164  current_error = std::fabs(cpu_matrix(col_it.index1(), col_it.index2()) - from_gpu(col_it.index1(), col_it.index2()))
165  / std::max( std::fabs(cpu_matrix(col_it.index1(), col_it.index2())),
166  std::fabs(from_gpu(col_it.index1(), col_it.index2())) );
167  if (current_error > error)
168  error = current_error;
169  }
170  }
171 
172  return error;
173 }
174 
175 //
176 // -------------------------------------------------------------
177 //
178 template< typename NumericT, typename Epsilon >
179 int test(Epsilon const& epsilon)
180 {
181  int retval = EXIT_SUCCESS;
182 
184 
185  // --------------------------------------------------------------------------
186  NumericT alpha = static_cast<NumericT>(2.786);
187  NumericT beta = static_cast<NumericT>(1.432);
188 
189  ublas::vector<NumericT> rhs;
190  ublas::vector<NumericT> result;
191  ublas::compressed_matrix<NumericT> ublas_matrix;
192 
193  if (viennacl::io::read_matrix_market_file(ublas_matrix, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
194  {
195  std::cout << "Error reading Matrix file" << std::endl;
196  return EXIT_FAILURE;
197  }
198  //unsigned int cg_mat_size = cg_mat.size();
199  std::cout << "done reading matrix" << std::endl;
200 
201 
202  rhs.resize(ublas_matrix.size2());
203  for (std::size_t i=0; i<rhs.size(); ++i)
204  {
205  ublas_matrix(i,i) = NumericT(0.5); // Get rid of round-off errors by making row-sums unequal to zero:
206  rhs[i] = NumericT(1) + randomNumber();
207  }
208 
209  result = rhs;
210 
211 
212  viennacl::vector<NumericT> vcl_rhs(rhs.size());
213  viennacl::vector<NumericT> vcl_result(result.size());
214  viennacl::vector<NumericT> vcl_result2(result.size());
215  viennacl::compressed_matrix<NumericT> vcl_compressed_matrix(rhs.size(), rhs.size());
216  viennacl::coordinate_matrix<NumericT> vcl_coordinate_matrix(rhs.size(), rhs.size());
217  viennacl::ell_matrix<NumericT> vcl_ell_matrix;
218  viennacl::hyb_matrix<NumericT> vcl_hyb_matrix;
219 
220  viennacl::copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
221  viennacl::copy(ublas_matrix, vcl_compressed_matrix);
222  viennacl::copy(ublas_matrix, vcl_coordinate_matrix);
223 
224  // --------------------------------------------------------------------------
225  std::cout << "Testing products: compressed_matrix" << std::endl;
226  result = viennacl::linalg::prod(ublas_matrix, rhs);
227  {
228  viennacl::scheduler::statement my_statement(vcl_result, viennacl::op_assign(), viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs));
229  viennacl::scheduler::execute(my_statement);
230  }
231  vcl_result = viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs);
232 
233  if ( std::fabs(diff(result, vcl_result)) > epsilon )
234  {
235  std::cout << "# Error at operation: matrix-vector product with compressed_matrix" << std::endl;
236  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
237  retval = EXIT_FAILURE;
238  }
239 
240  std::cout << "Testing products: coordinate_matrix" << std::endl;
241  rhs *= NumericT(1.1);
242  vcl_rhs *= NumericT(1.1);
243  result = viennacl::linalg::prod(ublas_matrix, rhs);
244  {
245  viennacl::scheduler::statement my_statement(vcl_result, viennacl::op_assign(), viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs));
246  viennacl::scheduler::execute(my_statement);
247  }
248 
249  if ( std::fabs(diff(result, vcl_result)) > epsilon )
250  {
251  std::cout << "# Error at operation: matrix-vector product with coordinate_matrix" << std::endl;
252  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
253  retval = EXIT_FAILURE;
254  }
255 
256  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
257  {
258  viennacl::scheduler::statement my_statement(vcl_result2, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result);
259  viennacl::scheduler::execute(my_statement);
260  }
261 
262  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
263  {
264  std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
265  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
266  retval = EXIT_FAILURE;
267  }
268 
269  //std::cout << "Copying ell_matrix" << std::endl;
270  viennacl::copy(ublas_matrix, vcl_ell_matrix);
271  ublas_matrix.clear();
272  viennacl::copy(vcl_ell_matrix, ublas_matrix);// just to check that it's works
273 
274 
275  std::cout << "Testing products: ell_matrix" << std::endl;
276  rhs *= NumericT(1.1);
277  vcl_rhs *= NumericT(1.1);
278  result = viennacl::linalg::prod(ublas_matrix, rhs);
279  {
280  //viennacl::scheduler::statement my_statement(vcl_result, viennacl::op_assign(), viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs));
281  //viennacl::scheduler::execute(my_statement);
282  }
283  vcl_result = viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs);
284 
285  if ( std::fabs(diff(result, vcl_result)) > epsilon )
286  {
287  std::cout << "# Error at operation: matrix-vector product with ell_matrix" << std::endl;
288  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
289  retval = EXIT_FAILURE;
290  }
291 
292  //std::cout << "Copying hyb_matrix" << std::endl;
293  viennacl::copy(ublas_matrix, vcl_hyb_matrix);
294  ublas_matrix.clear();
295  viennacl::copy(vcl_hyb_matrix, ublas_matrix);// just to check that it's works
296  viennacl::copy(ublas_matrix, vcl_hyb_matrix);
297 
298  std::cout << "Testing products: hyb_matrix" << std::endl;
299  rhs *= NumericT(1.1);
300  vcl_rhs *= NumericT(1.1);
301  result = viennacl::linalg::prod(ublas_matrix, rhs);
302  {
303  viennacl::scheduler::statement my_statement(vcl_result, viennacl::op_assign(), viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs));
304  viennacl::scheduler::execute(my_statement);
305  }
306 
307  if ( std::fabs(diff(result, vcl_result)) > epsilon )
308  {
309  std::cout << "# Error at operation: matrix-vector product with hyb_matrix" << std::endl;
310  std::cout << " diff: " << std::fabs(diff(result, vcl_result)) << std::endl;
311  retval = EXIT_FAILURE;
312  }
313 
314 
315  // --------------------------------------------------------------------------
316  // --------------------------------------------------------------------------
317  copy(rhs.begin(), rhs.end(), vcl_rhs.begin());
318  copy(result.begin(), result.end(), vcl_result.begin());
319  copy(result.begin(), result.end(), vcl_result2.begin());
320  copy(ublas_matrix, vcl_compressed_matrix);
321  copy(ublas_matrix, vcl_coordinate_matrix);
322  copy(ublas_matrix, vcl_ell_matrix);
323  copy(ublas_matrix, vcl_hyb_matrix);
324 
325  std::cout << "Testing scaled additions of products and vectors: compressed_matrix" << std::endl;
326  rhs *= NumericT(1.1);
327  vcl_rhs *= NumericT(1.1);
328  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
329  {
330  viennacl::scheduler::statement my_statement(vcl_result2, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_compressed_matrix, vcl_rhs) + beta * vcl_result);
331  viennacl::scheduler::execute(my_statement);
332  }
333 
334  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
335  {
336  std::cout << "# Error at operation: matrix-vector product (compressed_matrix) with scaled additions" << std::endl;
337  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
338  retval = EXIT_FAILURE;
339  }
340 
341 
342  std::cout << "Testing scaled additions of products and vectors: coordinate_matrix" << std::endl;
343  copy(result.begin(), result.end(), vcl_result.begin());
344  rhs *= NumericT(1.1);
345  vcl_rhs *= NumericT(1.1);
346  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
347  {
348  viennacl::scheduler::statement my_statement(vcl_result2, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_coordinate_matrix, vcl_rhs) + beta * vcl_result);
349  viennacl::scheduler::execute(my_statement);
350  }
351 
352  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
353  {
354  std::cout << "# Error at operation: matrix-vector product (coordinate_matrix) with scaled additions" << std::endl;
355  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
356  retval = EXIT_FAILURE;
357  }
358 
359  std::cout << "Testing scaled additions of products and vectors: ell_matrix" << std::endl;
360  copy(result.begin(), result.end(), vcl_result.begin());
361  rhs *= NumericT(1.1);
362  vcl_rhs *= NumericT(1.1);
363  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
364  {
365  viennacl::scheduler::statement my_statement(vcl_result2, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_ell_matrix, vcl_rhs) + beta * vcl_result);
366  viennacl::scheduler::execute(my_statement);
367  }
368 
369  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
370  {
371  std::cout << "# Error at operation: matrix-vector product (ell_matrix) with scaled additions" << std::endl;
372  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
373  retval = EXIT_FAILURE;
374  }
375 
376  std::cout << "Testing scaled additions of products and vectors: hyb_matrix" << std::endl;
377  copy(result.begin(), result.end(), vcl_result.begin());
378  rhs *= NumericT(1.1);
379  vcl_rhs *= NumericT(1.1);
380  result = alpha * viennacl::linalg::prod(ublas_matrix, rhs) + beta * result;
381  {
382  viennacl::scheduler::statement my_statement(vcl_result2, viennacl::op_assign(), alpha * viennacl::linalg::prod(vcl_hyb_matrix, vcl_rhs) + beta * vcl_result);
383  viennacl::scheduler::execute(my_statement);
384  }
385 
386  if ( std::fabs(diff(result, vcl_result2)) > epsilon )
387  {
388  std::cout << "# Error at operation: matrix-vector product (hyb_matrix) with scaled additions" << std::endl;
389  std::cout << " diff: " << std::fabs(diff(result, vcl_result2)) << std::endl;
390  retval = EXIT_FAILURE;
391  }
392 
393 
394  // --------------------------------------------------------------------------
395  return retval;
396 }
397 //
398 // -------------------------------------------------------------
399 //
400 int main()
401 {
402  std::cout << std::endl;
403  std::cout << "----------------------------------------------" << std::endl;
404  std::cout << "----------------------------------------------" << std::endl;
405  std::cout << "## Test :: Sparse Matrices" << std::endl;
406  std::cout << "----------------------------------------------" << std::endl;
407  std::cout << "----------------------------------------------" << std::endl;
408  std::cout << std::endl;
409 
410  int retval = EXIT_SUCCESS;
411 
412  std::cout << std::endl;
413  std::cout << "----------------------------------------------" << std::endl;
414  std::cout << std::endl;
415  {
416  typedef float NumericT;
417  NumericT epsilon = static_cast<NumericT>(1E-4);
418  std::cout << "# Testing setup:" << std::endl;
419  std::cout << " eps: " << epsilon << std::endl;
420  std::cout << " numeric: float" << std::endl;
421  retval = test<NumericT>(epsilon);
422  if ( retval == EXIT_SUCCESS )
423  std::cout << "# Test passed" << std::endl;
424  else
425  return retval;
426  }
427  std::cout << std::endl;
428  std::cout << "----------------------------------------------" << std::endl;
429  std::cout << std::endl;
430 
431 #ifdef VIENNACL_WITH_OPENCL
433 #endif
434  {
435  {
436  typedef double NumericT;
437  NumericT epsilon = 1.0E-13;
438  std::cout << "# Testing setup:" << std::endl;
439  std::cout << " eps: " << epsilon << std::endl;
440  std::cout << " numeric: double" << std::endl;
441  retval = test<NumericT>(epsilon);
442  if ( retval == EXIT_SUCCESS )
443  std::cout << "# Test passed" << std::endl;
444  else
445  return retval;
446  }
447  std::cout << std::endl;
448  std::cout << "----------------------------------------------" << std::endl;
449  std::cout << std::endl;
450  }
451 #ifdef VIENNACL_WITH_OPENCL
452  else
453  std::cout << "No double precision support, skipping test..." << std::endl;
454 #endif
455 
456 
457  std::cout << std::endl;
458  std::cout << "------- Test completed --------" << std::endl;
459  std::cout << std::endl;
460 
461  return retval;
462 }
Sparse matrix class using a hybrid format composed of the ELL and CSR format for storing the nonzeros...
Definition: forwards.h:406
A reader and writer for the matrix market format is implemented here.
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
int test(Epsilon const &epsilon)
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
Some helper routines for reading/writing/printing scheduler expressions.
A tag class representing assignment.
Definition: forwards.h:81
void finish()
Synchronizes the execution. finish() will only return after all compute kernels (CUDA, OpenCL) have completed.
Definition: memory.hpp:54
void execute(statement const &s)
Definition: execute.hpp:279
viennacl::scalar< int > s2
viennacl::scalar< float > s1
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ocl::device const & current_device()
Convenience function for returning the active device in the current context.
Definition: backend.hpp:351
Implementation of the coordinate_matrix class.
float NumericT
Definition: bisect.cpp:40
viennacl::vector< float > v1
Implementation of the hyb_matrix class.
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
Definition: vector.hpp:841
Implementations of incomplete factorization preconditioners. Convenience header file.
Implementation of the compressed_matrix class.
bool double_support() const
ViennaCL convenience function: Returns true if the device supports double precision.
Definition: device.hpp:956
Implementation of the ell_matrix class.
Proxy classes for vectors.
viennacl::vector< int > v2
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
int main()
T norm_inf(std::vector< T, A > const &v1)
Definition: norm_inf.hpp:60
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) ...
A small collection of sequential random number generators.
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
float ScalarType
Definition: fft_1d.cpp:42
Provides the datastructures for dealing with a single statement such as 'x = y + z;'.
The main class for representing a statement such as x = inner_prod(y,z); at runtime.
Definition: forwards.h:502
ScalarType diff(ScalarType &s1, viennacl::scalar< ScalarType > &s2)
long read_matrix_market_file(MatrixT &mat, const char *file, long index_base=1)
Reads a sparse matrix from a file (MatrixMarket format)
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
Definition: vector.hpp:848
Common routines used within ILU-type preconditioners.
Implementation of the ViennaCL scalar class.
A sparse square matrix, where entries are stored as triplets (i,j, val), where i and j are the row an...