All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sparse_prod.cpp
Go to the documentation of this file.
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the PDF manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
//
// *** System
//
#include <iostream>
#include <vector>
#include <map>
//
// *** ViennaCL
//
#include "viennacl/scalar.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/tools/random.hpp"
//
// -------------------------------------------------------------
//
/* Routine for computing the relative difference of two matrices. 1 is returned if the sparsity patterns do not match. */
template<typename IndexT, typename NumericT, typename MatrixT>
NumericT diff(std::vector<std::map<IndexT, NumericT> > const & stl_A,
MatrixT & vcl_A)
{
viennacl::switch_memory_context(vcl_A, viennacl::context(viennacl::MAIN_MEMORY));
NumericT error = NumericT(-1.0);
NumericT const * vcl_A_elements = viennacl::linalg::host_based::detail::extract_raw_pointer<NumericT>(vcl_A.handle());
unsigned int const * vcl_A_row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(vcl_A.handle1());
unsigned int const * vcl_A_col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(vcl_A.handle2());
/* Simultaneously compare the sparsity patterns of both matrices against each other. */
unsigned int const * vcl_A_current_col_ptr = vcl_A_col_buffer;
NumericT const * vcl_A_current_val_ptr = vcl_A_elements;
for (std::size_t row = 0; row < stl_A.size(); ++row)
{
if (vcl_A_current_col_ptr != vcl_A_col_buffer + vcl_A_row_buffer[row])
{
std::cerr << "Sparsity pattern mismatch detected: Start of row out of sync!" << std::endl;
std::cerr << " STL row: " << row << std::endl;
std::cerr << " ViennaCL col ptr is: " << vcl_A_current_col_ptr << std::endl;
std::cerr << " ViennaCL col ptr should: " << vcl_A_col_buffer + vcl_A_row_buffer[row] << std::endl;
std::cerr << " ViennaCL col ptr value: " << *vcl_A_current_col_ptr << std::endl;
return NumericT(1.0);
}
//std::cout << "Row " << row_it.index1() << ": " << std::endl;
for (typename std::map<IndexT, NumericT>::const_iterator col_it = stl_A[row].begin();
col_it != stl_A[row].end();
++col_it, ++vcl_A_current_col_ptr, ++vcl_A_current_val_ptr)
{
if (col_it->first != std::size_t(*vcl_A_current_col_ptr))
{
std::cerr << "Sparsity pattern mismatch detected!" << std::endl;
std::cerr << " STL row: " << row << std::endl;
std::cerr << " STL col: " << col_it->first << std::endl;
std::cerr << " ViennaCL row entries: " << vcl_A_row_buffer[row] << ", " << vcl_A_row_buffer[row + 1] << std::endl;
std::cerr << " ViennaCL entry in row: " << vcl_A_current_col_ptr - (vcl_A_col_buffer + vcl_A_row_buffer[row]) << std::endl;
std::cerr << " ViennaCL col: " << *vcl_A_current_col_ptr << std::endl;
return NumericT(1.0);
}
// compute relative error (we know for sure that the uBLAS matrix only carries nonzero entries:
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));
if (current_error > 0.1)
{
std::cerr << "Value mismatch detected!" << std::endl;
std::cerr << " STL row: " << row << std::endl;
std::cerr << " STL col: " << col_it->first << std::endl;
std::cerr << " STL value: " << col_it->second << std::endl;
std::cerr << " ViennaCL value: " << *vcl_A_current_val_ptr << std::endl;
return NumericT(1.0);
}
if (current_error > error)
error = current_error;
}
}
return error;
}
template<typename IndexT, typename NumericT>
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)
{
for (std::size_t i=0; i<stl_A.size(); ++i)
for (typename std::map<IndexT, NumericT>::const_iterator it_A = stl_A[i].begin(); it_A != stl_A[i].end(); ++it_A)
{
IndexT row_B = it_A->first;
for (typename std::map<IndexT, NumericT>::const_iterator it_B = stl_B[row_B].begin(); it_B != stl_B[row_B].end(); ++it_B)
stl_C[i][it_B->first] += it_A->second * it_B->second;
}
}
//
// -------------------------------------------------------------
//
template< typename NumericT, typename Epsilon >
int test(Epsilon const& epsilon)
{
int retval = EXIT_SUCCESS;
viennacl::tools::uniform_random_numbers<NumericT> randomNumber;
std::size_t N = 210;
std::size_t K = 300;
std::size_t M = 420;
std::size_t nnz_row = 40;
// --------------------------------------------------------------------------
std::vector<std::map<unsigned int, NumericT> > stl_A(N);
std::vector<std::map<unsigned int, NumericT> > stl_B(K);
std::vector<std::map<unsigned int, NumericT> > stl_C(N);
for (std::size_t i=0; i<stl_A.size(); ++i)
for (std::size_t j=0; j<nnz_row; ++j)
stl_A[i][static_cast<unsigned int>(randomNumber() * NumericT(K))] = NumericT(1.0) + NumericT();
for (std::size_t i=0; i<stl_B.size(); ++i)
for (std::size_t j=0; j<nnz_row; ++j)
stl_B[i][static_cast<unsigned int>(randomNumber() * NumericT(M))] = NumericT(1.0) + NumericT();
viennacl::compressed_matrix<NumericT> vcl_A(N, K);
viennacl::compressed_matrix<NumericT> vcl_B(K, M);
viennacl::compressed_matrix<NumericT> vcl_C;
viennacl::tools::sparse_matrix_adapter<NumericT> adapted_stl_A(stl_A, N, K);
viennacl::tools::sparse_matrix_adapter<NumericT> adapted_stl_B(stl_B, K, M);
viennacl::copy(adapted_stl_A, vcl_A);
viennacl::copy(adapted_stl_B, vcl_B);
// --------------------------------------------------------------------------
std::cout << "Testing products: STL" << std::endl;
prod(stl_A, stl_B, stl_C);
std::cout << "Testing products: compressed_matrix" << std::endl;
vcl_C = viennacl::linalg::prod(vcl_A, vcl_B);
if ( std::fabs(diff(stl_C, vcl_C)) > epsilon )
{
std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_C)" << std::endl;
std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
retval = EXIT_FAILURE;
}
viennacl::compressed_matrix<NumericT> vcl_D = viennacl::linalg::prod(vcl_A, vcl_B);
if ( std::fabs(diff(stl_C, vcl_D)) > epsilon )
{
std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_D)" << std::endl;
std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
retval = EXIT_FAILURE;
}
viennacl::compressed_matrix<NumericT> vcl_E(viennacl::linalg::prod(vcl_A, vcl_B));
if ( std::fabs(diff(stl_C, vcl_E)) > epsilon )
{
std::cout << "# Error at operation: matrix-matrix product with compressed_matrix (vcl_E)" << std::endl;
std::cout << " diff: " << std::fabs(diff(stl_C, vcl_C)) << std::endl;
retval = EXIT_FAILURE;
}
// --------------------------------------------------------------------------
return retval;
}
//
// -------------------------------------------------------------
//
int main()
{
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "## Test :: Sparse Matrix Product" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
int retval = EXIT_SUCCESS;
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
{
typedef float NumericT;
NumericT epsilon = static_cast<NumericT>(1E-4);
std::cout << "# Testing setup:" << std::endl;
std::cout << " eps: " << epsilon << std::endl;
std::cout << " numeric: float" << std::endl;
retval = test<NumericT>(epsilon);
if ( retval == EXIT_SUCCESS )
std::cout << "# Test passed" << std::endl;
else
return retval;
}
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
#ifdef VIENNACL_WITH_OPENCL
if ( viennacl::ocl::current_device().double_support() )
#endif
{
{
typedef double NumericT;
NumericT epsilon = 1.0E-12;
std::cout << "# Testing setup:" << std::endl;
std::cout << " eps: " << epsilon << std::endl;
std::cout << " numeric: double" << std::endl;
retval = test<NumericT>(epsilon);
if ( retval == EXIT_SUCCESS )
std::cout << "# Test passed" << std::endl;
else
return retval;
}
std::cout << std::endl;
std::cout << "----------------------------------------------" << std::endl;
std::cout << std::endl;
}
#ifdef VIENNACL_WITH_OPENCL
else
std::cout << "No double precision support, skipping test..." << std::endl;
#endif
std::cout << std::endl;
std::cout << "------- Test completed --------" << std::endl;
std::cout << std::endl;
return retval;
}
main
int main()
Definition: sparse_prod.cpp:201
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
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
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::context
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
test
int test(Epsilon const &epsilon)
Definition: sparse_prod.cpp:132
viennacl::tools::uniform_random_numbers
Random number generator for returning uniformly distributed values in the closed interval [0...
Definition: random.hpp:44
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
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::MAIN_MEMORY
Definition: forwards.h:348
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
diff
NumericT diff(std::vector< std::map< IndexT, NumericT > > const &stl_A, MatrixT &vcl_A)
Definition: sparse_prod.cpp:47
viennacl::copy
void copy(std::vector< NumericT > &cpu_vec, circulant_matrix< NumericT, AlignmentV > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
random.hpp
A small collection of sequential random number generators.
viennacl::tools::sparse_matrix_adapter
Adapts a non-const sparse matrix type made up from std::vector<std::map<SizeT, NumericT> > to basic u...
Definition: adapter.hpp:357
viennacl::compressed_matrix< NumericT >
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::switch_memory_context
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