ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
spmdm.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 //
24 // include necessary system headers
25 //
26 #include <iostream>
27 #include <cmath>
28 
29 //
30 // ublas includes
31 //
32 #include <boost/numeric/ublas/io.hpp>
33 #include <boost/numeric/ublas/triangular.hpp>
34 #include <boost/numeric/ublas/matrix_sparse.hpp>
35 #include <boost/numeric/ublas/matrix.hpp>
36 #include <boost/numeric/ublas/matrix_proxy.hpp>
37 #include <boost/numeric/ublas/operation.hpp>
38 #include <boost/numeric/ublas/operation_sparse.hpp>
39 #include <boost/numeric/ublas/io.hpp>
40 #include <boost/numeric/ublas/lu.hpp>
41 
42 // Must be set if you want to use ViennaCL algorithms on ublas objects
43 #define VIENNACL_WITH_UBLAS 1
44 //#define VIENNACL_WITH_OPENCL 1
45 //#define VIENNACL_WITH_CUDA 1
46 //#define VIENNACL_DEBUG_KERNEL 1
47 //#define VIENNACL_BUILD_INFO 1
48 
49 //
50 // ViennaCL includes
51 //
52 #include "viennacl/scalar.hpp"
53 #include "viennacl/vector.hpp"
54 #include "viennacl/matrix.hpp"
55 #include "viennacl/linalg/direct_solve.hpp"
56 #include "viennacl/compressed_matrix.hpp"
57 #include "viennacl/coordinate_matrix.hpp"
58 #include "viennacl/ell_matrix.hpp"
59 #include "viennacl/hyb_matrix.hpp"
60 #include "viennacl/linalg/prod.hpp" //generic matrix-vector product
61 #include "viennacl/linalg/norm_2.hpp" //generic l2-norm for vectors
62 #include "viennacl/io/matrix_market.hpp"
63 #include "viennacl/tools/random.hpp"
64 
65 
66 using namespace boost::numeric;
67 
68 template< typename ScalarType >
69 int check_matrices(const ublas::matrix< ScalarType >& ref_mat, const ublas::matrix< ScalarType >& mat, ScalarType eps) {
70 
71  std::size_t size1, size2;
72  size1 = ref_mat.size1(); size2 = ref_mat.size2();
73  if ( (size1 != mat.size1()) || (size2 != mat.size2()) )
74  return EXIT_FAILURE;
75 
76  for (unsigned int i = 0; i < size1; i++)
77  for (unsigned int j = 0; j < size2; j++)
78  {
79  ScalarType rel_error = std::abs(ref_mat(i,j) - mat(i,j)) / std::max(std::abs(ref_mat(i,j)), std::abs(mat(i,j)));
80  if ( rel_error > eps ) {
81  std::cout << "ERROR: Verification failed at (" << i <<", "<< j << "): "
82  << " Expected: " << ref_mat(i,j) << ", got: " << mat(i,j) << " (relative error: " << rel_error << ")" << std::endl;
83  return EXIT_FAILURE;
84  }
85  }
86 
87  std::cout << "Everything went well!" << std::endl;
88  return EXIT_SUCCESS;
89 }
90 
91 template<typename NumericT, typename ResultLayoutT, typename FactorLayoutT>
92 int test(NumericT epsilon)
93 {
94  int retVal = EXIT_SUCCESS;
95 
96  viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
97 
98  ublas::compressed_matrix<NumericT> ublas_lhs;
99 
100  if (viennacl::io::read_matrix_market_file(ublas_lhs, "../examples/testdata/mat65k.mtx") == EXIT_FAILURE)
101  {
102  std::cout << "Error reading Matrix file" << std::endl;
103  return EXIT_FAILURE;
104  }
105 
106  // add some extra weight to diagonal in order to avoid issues with round-off errors
107  for (std::size_t i=0; i<ublas_lhs.size1(); ++i)
108  ublas_lhs(i,i) *= NumericT(1.5);
109 
110  std::size_t cols_rhs = 1;
111 
112  viennacl::compressed_matrix<NumericT> compressed_lhs;
113  viennacl::ell_matrix<NumericT> ell_lhs;
114  viennacl::coordinate_matrix<NumericT> coo_lhs;
115  viennacl::hyb_matrix<NumericT> hyb_lhs;
116 
117  ublas::matrix<NumericT> ublas_result;
118  viennacl::matrix<NumericT, ResultLayoutT> result;
119 
120  viennacl::copy( ublas_lhs, compressed_lhs);
121  viennacl::copy( ublas_lhs, ell_lhs);
122  viennacl::copy( ublas_lhs, coo_lhs);
123  viennacl::copy( ublas_lhs, hyb_lhs);
124 
125  ublas::matrix<NumericT> ublas_rhs1(ublas_lhs.size2(), cols_rhs);
126  viennacl::matrix<NumericT, FactorLayoutT> rhs1(ublas_lhs.size2(), cols_rhs);
127 
128  ublas::matrix<NumericT> ublas_rhs2;
129  viennacl::matrix<NumericT, FactorLayoutT> rhs2;
130 
131  ublas::matrix<NumericT> temp(ublas_rhs1.size1(), cols_rhs);
132 
133  for (unsigned int i = 0; i < ublas_rhs1.size1(); i++)
134  for (unsigned int j = 0; j < ublas_rhs1.size2(); j++)
135  ublas_rhs1(i,j) = NumericT(0.5) + NumericT(0.1) * randomNumber();
136  viennacl::copy( ublas_rhs1, rhs1);
137 
138  ublas_rhs2 = ublas::trans( ublas_rhs1);
139  viennacl::copy( ublas_rhs2, rhs2);
140 
141  /* gold result */
142  ublas_result = ublas::prod( ublas_lhs, ublas_rhs1);
143 
144  /******************************************************************/
145  std::cout << "Testing compressed(CSR) lhs * dense rhs" << std::endl;
146  result = viennacl::linalg::prod( compressed_lhs, rhs1);
147 
148  temp.clear();
149  viennacl::copy( result, temp);
150  retVal = check_matrices(ublas_result, temp, epsilon);
151 
152  /******************************************************************/
153  std::cout << "Testing compressed(ELL) lhs * dense rhs" << std::endl;
154  result.clear();
155  result = viennacl::linalg::prod( ell_lhs, rhs1);
156 
157  temp.clear();
158  viennacl::copy( result, temp);
159  check_matrices(ublas_result, temp, epsilon);
160 
161  /******************************************************************/
162 
163  std::cout << "Testing compressed(COO) lhs * dense rhs" << std::endl;
164  result.clear();
165  result = viennacl::linalg::prod( coo_lhs, rhs1);
166 
167  temp.clear();
168  viennacl::copy( result, temp);
169  check_matrices(ublas_result, temp, epsilon);
170 
171  /******************************************************************/
172 
173  std::cout << "Testing compressed(HYB) lhs * dense rhs" << std::endl;
174  result.clear();
175  result = viennacl::linalg::prod( hyb_lhs, rhs1);
176 
177  temp.clear();
178  viennacl::copy( result, temp);
179  check_matrices(ublas_result, temp, epsilon);
180 
181  /******************************************************************/
182 
183  /* gold result */
184  ublas_result = ublas::prod( ublas_lhs, ublas::trans(ublas_rhs2));
185 
186  /******************************************************************/
187  std::cout << std::endl << "Testing compressed(CSR) lhs * transposed dense rhs:" << std::endl;
188  result.clear();
189  result = viennacl::linalg::prod( compressed_lhs, viennacl::trans(rhs2));
190 
191  temp.clear();
192  viennacl::copy( result, temp);
193  retVal = check_matrices(ublas_result, temp, epsilon);
194 
195  /******************************************************************/
196  std::cout << "Testing compressed(ELL) lhs * transposed dense rhs" << std::endl;
197  result.clear();
198  result = viennacl::linalg::prod( ell_lhs, viennacl::trans(rhs2));
199 
200  temp.clear();
201  viennacl::copy( result, temp);
202  check_matrices(ublas_result, temp, epsilon);
203 
204  /******************************************************************/
205  std::cout << "Testing compressed(COO) lhs * transposed dense rhs" << std::endl;
206  result.clear();
207  result = viennacl::linalg::prod( coo_lhs, viennacl::trans(rhs2));
208 
209  temp.clear();
210  viennacl::copy( result, temp);
211  check_matrices(ublas_result, temp, epsilon);
212 
213  /******************************************************************/
214 
215  std::cout << "Testing compressed(HYB) lhs * transposed dense rhs" << std::endl;
216  result.clear();
217  result = viennacl::linalg::prod( hyb_lhs, viennacl::trans(rhs2));
218 
219  temp.clear();
220  viennacl::copy( result, temp);
221  check_matrices(ublas_result, temp, epsilon);
222 
223  /******************************************************************/
224  if (retVal == EXIT_SUCCESS) {
225  std::cout << "Tests passed successfully" << std::endl;
226  }
227 
228  return retVal;
229 }
230 
231 //
232 // -------------------------------------------------------------
233 //
234 int main()
235 {
236  std::cout << std::endl;
237  std::cout << "----------------------------------------------" << std::endl;
238  std::cout << "----------------------------------------------" << std::endl;
239  std::cout << "## Test :: Sparse-Dense Matrix Multiplication" << std::endl;
240  std::cout << "----------------------------------------------" << std::endl;
241  std::cout << "----------------------------------------------" << std::endl;
242  std::cout << std::endl;
243 
244  int retval = EXIT_SUCCESS;
245 
246  std::cout << std::endl;
247  std::cout << "----------------------------------------------" << std::endl;
248  std::cout << std::endl;
249  {
250  typedef float NumericT;
251  NumericT epsilon = static_cast<NumericT>(1E-4);
252  std::cout << "# Testing setup:" << std::endl;
253  std::cout << " eps: " << epsilon << std::endl;
254  std::cout << " numeric: float" << std::endl;
255  std::cout << " layout: row-major, row-major" << std::endl;
256  retval = test<NumericT, viennacl::row_major, viennacl::row_major>(epsilon);
257  if ( retval == EXIT_SUCCESS )
258  std::cout << "# Test passed" << std::endl;
259  else
260  return retval;
261 
262  std::cout << "# Testing setup:" << std::endl;
263  std::cout << " eps: " << epsilon << std::endl;
264  std::cout << " numeric: float" << std::endl;
265  std::cout << " layout: row-major, column-major" << std::endl;
266  retval = test<NumericT, viennacl::row_major, viennacl::column_major>(epsilon);
267  if ( retval == EXIT_SUCCESS )
268  std::cout << "# Test passed" << std::endl;
269  else
270  return retval;
271 
272  std::cout << "# Testing setup:" << std::endl;
273  std::cout << " eps: " << epsilon << std::endl;
274  std::cout << " numeric: float" << std::endl;
275  std::cout << " layout: column-major, row-major" << std::endl;
276  retval = test<NumericT, viennacl::column_major, viennacl::row_major>(epsilon);
277  if ( retval == EXIT_SUCCESS )
278  std::cout << "# Test passed" << std::endl;
279  else
280  return retval;
281 
282  std::cout << "# Testing setup:" << std::endl;
283  std::cout << " eps: " << epsilon << std::endl;
284  std::cout << " numeric: float" << std::endl;
285  std::cout << " layout: column-major, column-major" << std::endl;
286  retval = test<NumericT, viennacl::column_major, viennacl::column_major>(epsilon);
287  if ( retval == EXIT_SUCCESS )
288  std::cout << "# Test passed" << std::endl;
289  else
290  return retval;
291 
292  }
293  std::cout << std::endl;
294  std::cout << "----------------------------------------------" << std::endl;
295  std::cout << std::endl;
296 
297 #ifdef VIENNACL_WITH_OPENCL
298  if ( viennacl::ocl::current_device().double_support() )
299 #endif
300  {
301  {
302  typedef double NumericT;
303  NumericT epsilon = 1.0E-12;
304  std::cout << "# Testing setup:" << std::endl;
305  std::cout << " eps: " << epsilon << std::endl;
306  std::cout << " numeric: double" << std::endl;
307  std::cout << " layout: row-major, row-major" << std::endl;
308  retval = test<NumericT, viennacl::row_major, viennacl::row_major>(epsilon);
309  if ( retval == EXIT_SUCCESS )
310  std::cout << "# Test passed" << std::endl;
311  else
312  return retval;
313 
314  std::cout << "# Testing setup:" << std::endl;
315  std::cout << " eps: " << epsilon << std::endl;
316  std::cout << " numeric: double" << std::endl;
317  std::cout << " layout: row-major, column-major" << std::endl;
318  retval = test<NumericT, viennacl::row_major, viennacl::column_major>(epsilon);
319  if ( retval == EXIT_SUCCESS )
320  std::cout << "# Test passed" << std::endl;
321  else
322  return retval;
323 
324  std::cout << "# Testing setup:" << std::endl;
325  std::cout << " eps: " << epsilon << std::endl;
326  std::cout << " numeric: double" << std::endl;
327  std::cout << " layout: column-major, row-major" << std::endl;
328  retval = test<NumericT, viennacl::column_major, viennacl::row_major>(epsilon);
329  if ( retval == EXIT_SUCCESS )
330  std::cout << "# Test passed" << std::endl;
331  else
332  return retval;
333 
334  std::cout << "# Testing setup:" << std::endl;
335  std::cout << " eps: " << epsilon << std::endl;
336  std::cout << " numeric: double" << std::endl;
337  std::cout << " layout: column-major, column-major" << std::endl;
338  retval = test<NumericT, viennacl::column_major, viennacl::column_major>(epsilon);
339  if ( retval == EXIT_SUCCESS )
340  std::cout << "# Test passed" << std::endl;
341  else
342  return retval;
343  }
344  std::cout << std::endl;
345  std::cout << "----------------------------------------------" << std::endl;
346  std::cout << std::endl;
347  }
348 #ifdef VIENNACL_WITH_OPENCL
349  else
350  std::cout << "No double precision support, skipping test..." << std::endl;
351 #endif
352 
353 
354  std::cout << std::endl;
355  std::cout << "------- Test completed --------" << std::endl;
356  std::cout << std::endl;
357 
358  return retval;
359 }
360 
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
check_matrices
int check_matrices(const ublas::matrix< ScalarType > &ref_mat, const ublas::matrix< ScalarType > &mat, ScalarType eps)
Definition: spmdm.cpp:69
matrix_market.hpp
A reader and writer for the matrix market format is implemented here.
norm_2.hpp
Generic interface for the l^2-norm. See viennacl/linalg/vector_operations.hpp for implementations...
viennacl::trans
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:374
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...
matrix.hpp
Implementation of the dense matrix class.
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
main
int main()
Definition: spmdm.cpp:234
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
viennacl::traits::size2
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:201
coordinate_matrix.hpp
Implementation of the coordinate_matrix class.
NumericT
float NumericT
Definition: bisect.cpp:40
test
int test(NumericT epsilon)
Definition: spmdm.cpp:92
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::matrix_base< NumericT >::clear
void clear()
Resets all entries to zero.
compressed_matrix.hpp
Implementation of the compressed_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
ell_matrix.hpp
Implementation of the ell_matrix class.
direct_solve.hpp
Implementations of dense direct solvers are found here.
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::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.
ScalarType
float ScalarType
Definition: fft_1d.cpp:42
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
scalar.hpp
Implementation of the ViennaCL scalar class.
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