ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
sparse_prod.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 //
25 // *** System
26 //
27 #include <iostream>
28 #include <vector>
29 #include <map>
30 
31 
32 //
33 // *** ViennaCL
34 //
35 #include "viennacl/scalar.hpp"
37 #include "viennacl/linalg/prod.hpp"
38 
40 
41 //
42 // -------------------------------------------------------------
43 //
44 
45 /* Routine for computing the relative difference of two matrices. 1 is returned if the sparsity patterns do not match. */
46 template<typename IndexT, typename NumericT, typename MatrixT>
47 NumericT diff(std::vector<std::map<IndexT, NumericT> > const & stl_A,
48  MatrixT & vcl_A)
49 {
51 
52  NumericT error = NumericT(-1.0);
53 
54  NumericT const * vcl_A_elements = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(vcl_A.handle());
55  unsigned int const * vcl_A_row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(vcl_A.handle1());
56  unsigned int const * vcl_A_col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(vcl_A.handle2());
57 
58 
59  /* Simultaneously compare the sparsity patterns of both matrices against each other. */
60 
61  unsigned int const * vcl_A_current_col_ptr = vcl_A_col_buffer;
62  NumericT const * vcl_A_current_val_ptr = vcl_A_elements;
63 
64  for (std::size_t row = 0; row < stl_A.size(); ++row)
65  {
66  if (vcl_A_current_col_ptr != vcl_A_col_buffer + vcl_A_row_buffer[row])
67  {
68  std::cerr << "Sparsity pattern mismatch detected: Start of row out of sync!" << std::endl;
69  std::cerr << " STL row: " << row << std::endl;
70  std::cerr << " ViennaCL col ptr is: " << vcl_A_current_col_ptr << std::endl;
71  std::cerr << " ViennaCL col ptr should: " << vcl_A_col_buffer + vcl_A_row_buffer[row] << std::endl;
72  std::cerr << " ViennaCL col ptr value: " << *vcl_A_current_col_ptr << std::endl;
73  return NumericT(1.0);
74  }
75 
76  //std::cout << "Row " << row_it.index1() << ": " << std::endl;
77  for (typename std::map<IndexT, NumericT>::const_iterator col_it = stl_A[row].begin();
78  col_it != stl_A[row].end();
79  ++col_it, ++vcl_A_current_col_ptr, ++vcl_A_current_val_ptr)
80  {
81  if (col_it->first != std::size_t(*vcl_A_current_col_ptr))
82  {
83  std::cerr << "Sparsity pattern mismatch detected!" << std::endl;
84  std::cerr << " STL row: " << row << std::endl;
85  std::cerr << " STL col: " << col_it->first << std::endl;
86  std::cerr << " ViennaCL row entries: " << vcl_A_row_buffer[row] << ", " << vcl_A_row_buffer[row + 1] << std::endl;
87  std::cerr << " ViennaCL entry in row: " << vcl_A_current_col_ptr - (vcl_A_col_buffer + vcl_A_row_buffer[row]) << std::endl;
88  std::cerr << " ViennaCL col: " << *vcl_A_current_col_ptr << std::endl;
89  return NumericT(1.0);
90  }
91 
92  // compute relative error (we know for sure that the uBLAS matrix only carries nonzero entries:
93  NumericT current_error = std::fabs(col_it->second - *vcl_A_current_val_ptr) / std::max(std::fabs(col_it->second), std::fabs(*vcl_A_current_val_ptr));
94 
95  if (current_error > 0.1)
96  {
97  std::cerr << "Value mismatch detected!" << std::endl;
98  std::cerr << " STL row: " << row << std::endl;
99  std::cerr << " STL col: " << col_it->first << std::endl;
100  std::cerr << " STL value: " << col_it->second << std::endl;
101  std::cerr << " ViennaCL value: " << *vcl_A_current_val_ptr << std::endl;
102  return NumericT(1.0);
103  }
104 
105  if (current_error > error)
106  error = current_error;
107  }
108  }
109 
110  return error;
111 }
112 
113 template<typename IndexT, typename NumericT>
114 void prod(std::vector<std::map<IndexT, NumericT> > const & stl_A,
115  std::vector<std::map<IndexT, NumericT> > const & stl_B,
116  std::vector<std::map<IndexT, NumericT> > & stl_C)
117 {
118  for (std::size_t i=0; i<stl_A.size(); ++i)
119  for (typename std::map<IndexT, NumericT>::const_iterator it_A = stl_A[i].begin(); it_A != stl_A[i].end(); ++it_A)
120  {
121  IndexT row_B = it_A->first;
122  for (typename std::map<IndexT, NumericT>::const_iterator it_B = stl_B[row_B].begin(); it_B != stl_B[row_B].end(); ++it_B)
123  stl_C[i][it_B->first] += it_A->second * it_B->second;
124  }
125 }
126 
127 
128 //
129 // -------------------------------------------------------------
130 //
131 template< typename NumericT, typename Epsilon >
132 int test(Epsilon const& epsilon)
133 {
134  int retval = EXIT_SUCCESS;
135 
137 
138  std::size_t N = 210;
139  std::size_t K = 300;
140  std::size_t M = 420;
141  std::size_t nnz_row = 40;
142  // --------------------------------------------------------------------------
143  std::vector<std::map<unsigned int, NumericT> > stl_A(N);
144  std::vector<std::map<unsigned int, NumericT> > stl_B(K);
145  std::vector<std::map<unsigned int, NumericT> > stl_C(N);
146 
147  for (std::size_t i=0; i<stl_A.size(); ++i)
148  for (std::size_t j=0; j<nnz_row; ++j)
149  stl_A[i][static_cast<unsigned int>(randomNumber() * NumericT(K))] = NumericT(1.0) + NumericT();
150 
151  for (std::size_t i=0; i<stl_B.size(); ++i)
152  for (std::size_t j=0; j<nnz_row; ++j)
153  stl_B[i][static_cast<unsigned int>(randomNumber() * NumericT(M))] = NumericT(1.0) + NumericT();
154 
155 
159 
160  viennacl::tools::sparse_matrix_adapter<NumericT> adapted_stl_A(stl_A, N, K);
161  viennacl::tools::sparse_matrix_adapter<NumericT> adapted_stl_B(stl_B, K, M);
162  viennacl::copy(adapted_stl_A, vcl_A);
163  viennacl::copy(adapted_stl_B, vcl_B);
164 
165  // --------------------------------------------------------------------------
166  std::cout << "Testing products: STL" << std::endl;
167  prod(stl_A, stl_B, stl_C);
168 
169  std::cout << "Testing products: compressed_matrix" << std::endl;
170  vcl_C = viennacl::linalg::prod(vcl_A, vcl_B);
171 
172  if ( std::fabs(diff(stl_C, vcl_C)) > epsilon )
173  {
174  std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_C)" << std::endl;
175  std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
176  retval = EXIT_FAILURE;
177  }
178 
180  if ( std::fabs(diff(stl_C, vcl_D)) > epsilon )
181  {
182  std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_D)" << std::endl;
183  std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
184  retval = EXIT_FAILURE;
185  }
186 
188  if ( std::fabs(diff(stl_C, vcl_E)) > epsilon )
189  {
190  std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_E)" << std::endl;
191  std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
192  retval = EXIT_FAILURE;
193  }
194 
195  // --------------------------------------------------------------------------
196  return retval;
197 }
198 //
199 // -------------------------------------------------------------
200 //
201 int main()
202 {
203  std::cout << std::endl;
204  std::cout << "----------------------------------------------" << std::endl;
205  std::cout << "----------------------------------------------" << std::endl;
206  std::cout << "## Test :: Sparse Matrix Product" << std::endl;
207  std::cout << "----------------------------------------------" << std::endl;
208  std::cout << "----------------------------------------------" << std::endl;
209  std::cout << std::endl;
210 
211  int retval = EXIT_SUCCESS;
212 
213  std::cout << std::endl;
214  std::cout << "----------------------------------------------" << std::endl;
215  std::cout << std::endl;
216  {
217  typedef float NumericT;
218  NumericT epsilon = static_cast<NumericT>(1E-4);
219  std::cout << "# Testing setup:" << std::endl;
220  std::cout << " eps: " << epsilon << std::endl;
221  std::cout << " numeric: float" << std::endl;
222  retval = test<NumericT>(epsilon);
223  if ( retval == EXIT_SUCCESS )
224  std::cout << "# Test passed" << std::endl;
225  else
226  return retval;
227  }
228  std::cout << std::endl;
229  std::cout << "----------------------------------------------" << std::endl;
230  std::cout << std::endl;
231 
232 #ifdef VIENNACL_WITH_OPENCL
234 #endif
235  {
236  {
237  typedef double NumericT;
238  NumericT epsilon = 1.0E-12;
239  std::cout << "# Testing setup:" << std::endl;
240  std::cout << " eps: " << epsilon << std::endl;
241  std::cout << " numeric: double" << std::endl;
242  retval = test<NumericT>(epsilon);
243  if ( retval == EXIT_SUCCESS )
244  std::cout << "# Test passed" << std::endl;
245  else
246  return retval;
247  }
248  std::cout << std::endl;
249  std::cout << "----------------------------------------------" << std::endl;
250  std::cout << std::endl;
251  }
252 #ifdef VIENNACL_WITH_OPENCL
253  else
254  std::cout << "No double precision support, skipping test..." << std::endl;
255 #endif
256 
257 
258  std::cout << std::endl;
259  std::cout << "------- Test completed --------" << std::endl;
260  std::cout << std::endl;
261 
262  return retval;
263 }
int main()
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
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
float NumericT
Definition: bisect.cpp:40
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
int test(Epsilon const &epsilon)
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
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
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:900
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)
NumericT diff(std::vector< std::map< IndexT, NumericT > > const &stl_A, MatrixT &vcl_A)
Definition: sparse_prod.cpp:47
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.
Adapts a non-const sparse matrix type made up from std::vector > to basic u...
Definition: adapter.hpp:357
Implementation of the ViennaCL scalar class.
void switch_memory_context(T &obj, viennacl::context new_ctx)
Generic convenience routine for migrating data of an object to a new memory domain.
Definition: memory.hpp:622