All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sparse_matrix_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_SPARSE_MATRIX_OPERATIONS_HPP_
#define VIENNACL_LINALG_SPARSE_MATRIX_OPERATIONS_HPP_
/* =========================================================================
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 manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include "viennacl/forwards.h"
#include "viennacl/scalar.hpp"
#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/host_based/sparse_matrix_operations.hpp"
#ifdef VIENNACL_WITH_OPENCL
#include "viennacl/linalg/opencl/sparse_matrix_operations.hpp"
#endif
#ifdef VIENNACL_WITH_CUDA
#include "viennacl/linalg/cuda/sparse_matrix_operations.hpp"
#endif
namespace viennacl
{
namespace linalg
{
namespace detail
{
template<typename SparseMatrixType, typename SCALARTYPE, unsigned int VEC_ALIGNMENT>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value >::type
row_info(SparseMatrixType const & mat,
vector<SCALARTYPE, VEC_ALIGNMENT> & vec,
row_info_types info_selector)
{
switch (viennacl::traits::handle(mat).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::detail::row_info(mat, vec, info_selector);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::detail::row_info(mat, vec, info_selector);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::detail::row_info(mat, vec, info_selector);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
}
// A * x
template<typename SparseMatrixType, class ScalarType>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
prod_impl(const SparseMatrixType & mat,
const viennacl::vector_base<ScalarType> & vec,
ScalarType alpha,
viennacl::vector_base<ScalarType> & result,
ScalarType beta)
{
assert( (mat.size1() == result.size()) && bool("Size check failed for compressed matrix-vector product: size1(mat) != size(result)"));
assert( (mat.size2() == vec.size()) && bool("Size check failed for compressed matrix-vector product: size2(mat) != size(x)"));
switch (viennacl::traits::handle(mat).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::prod_impl(mat, vec, alpha, result, beta);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::prod_impl(mat, vec, alpha, result, beta);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::prod_impl(mat, vec, alpha, result, beta);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
// A * B
template<typename SparseMatrixType, class ScalarType>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
prod_impl(const SparseMatrixType & sp_mat,
const viennacl::matrix_base<ScalarType> & d_mat,
viennacl::matrix_base<ScalarType> & result)
{
assert( (sp_mat.size1() == result.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size1(sp_mat) != size1(result)"));
assert( (sp_mat.size2() == d_mat.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size2(sp_mat) != size1(d_mat)"));
switch (viennacl::traits::handle(sp_mat).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::prod_impl(sp_mat, d_mat, result);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::prod_impl(sp_mat, d_mat, result);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::prod_impl(sp_mat, d_mat, result);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
// A * transpose(B)
template<typename SparseMatrixType, class ScalarType>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
prod_impl(const SparseMatrixType & sp_mat,
const viennacl::matrix_expression<const viennacl::matrix_base<ScalarType>,
const viennacl::matrix_base<ScalarType>,
viennacl::op_trans>& d_mat,
viennacl::matrix_base<ScalarType> & result)
{
assert( (sp_mat.size1() == result.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size1(sp_mat) != size1(result)"));
assert( (sp_mat.size2() == d_mat.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size2(sp_mat) != size1(d_mat)"));
switch (viennacl::traits::handle(sp_mat).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::prod_impl(sp_mat, d_mat, result);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::prod_impl(sp_mat, d_mat, result);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::prod_impl(sp_mat, d_mat, result);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
// A * B with both A and B sparse
template<typename NumericT>
void
prod_impl(const viennacl::compressed_matrix<NumericT> & A,
const viennacl::compressed_matrix<NumericT> & B,
viennacl::compressed_matrix<NumericT> & C)
{
assert( (A.size2() == B.size1()) && bool("Size check failed for sparse matrix-matrix product: size2(A) != size1(B)"));
assert( (C.size1() == 0 || C.size1() == A.size1()) && bool("Size check failed for sparse matrix-matrix product: size1(A) != size1(C)"));
assert( (C.size2() == 0 || C.size2() == B.size2()) && bool("Size check failed for sparse matrix-matrix product: size2(B) != size2(B)"));
switch (viennacl::traits::handle(A).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::prod_impl(A, B, C);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::prod_impl(A, B, C);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::prod_impl(A, B, C);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
inplace_solve(const SparseMatrixType & mat,
viennacl::vector_base<ScalarType> & vec,
SOLVERTAG tag)
{
assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on compressed matrix: size1(mat) != size2(mat)"));
assert( (mat.size2() == vec.size()) && bool("Size check failed for compressed matrix-vector product: size2(mat) != size(x)"));
switch (viennacl::traits::handle(mat).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::inplace_solve(mat, vec, tag);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::inplace_solve(mat, vec, tag);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::inplace_solve(mat, vec, tag);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
inplace_solve(const matrix_expression<const SparseMatrixType, const SparseMatrixType, op_trans> & mat,
viennacl::vector_base<ScalarType> & vec,
SOLVERTAG tag)
{
assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on transposed compressed matrix: size1(mat) != size2(mat)"));
assert( (mat.size1() == vec.size()) && bool("Size check failed for transposed compressed matrix triangular solve: size1(mat) != size(x)"));
switch (viennacl::traits::handle(mat.lhs()).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::inplace_solve(mat, vec, tag);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::inplace_solve(mat, vec, tag);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::inplace_solve(mat, vec, tag);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
namespace detail
{
template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value>::type
block_inplace_solve(const matrix_expression<const SparseMatrixType, const SparseMatrixType, op_trans> & mat,
viennacl::backend::mem_handle const & block_index_array, vcl_size_t num_blocks,
viennacl::vector_base<ScalarType> const & mat_diagonal,
viennacl::vector_base<ScalarType> & vec,
SOLVERTAG tag)
{
assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on transposed compressed matrix: size1(mat) != size2(mat)"));
assert( (mat.size1() == vec.size()) && bool("Size check failed for transposed compressed matrix triangular solve: size1(mat) != size(x)"));
switch (viennacl::traits::handle(mat.lhs()).get_active_handle_id())
{
case viennacl::MAIN_MEMORY:
viennacl::linalg::host_based::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
break;
#ifdef VIENNACL_WITH_OPENCL
case viennacl::OPENCL_MEMORY:
viennacl::linalg::opencl::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
break;
#endif
#ifdef VIENNACL_WITH_CUDA
case viennacl::CUDA_MEMORY:
viennacl::linalg::cuda::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
break;
#endif
case viennacl::MEMORY_NOT_INITIALIZED:
throw memory_exception("not initialised!");
default:
throw memory_exception("not implemented");
}
}
}
} //namespace linalg
template<typename M1>
typename viennacl::enable_if<viennacl::is_any_sparse_matrix<M1>::value,
matrix_expression< const M1, const M1, op_trans>
>::type
trans(const M1 & mat)
{
return matrix_expression< const M1, const M1, op_trans>(mat, mat);
}
//free functions:
template<typename SCALARTYPE, typename SparseMatrixType>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
viennacl::vector<SCALARTYPE> >::type
operator+(viennacl::vector_base<SCALARTYPE> & result,
const viennacl::vector_expression< const SparseMatrixType, const viennacl::vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
{
assert(proxy.lhs().size1() == result.size() && bool("Dimensions for addition of sparse matrix-vector product to vector don't match!"));
vector<SCALARTYPE> temp(proxy.lhs().size1());
viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), temp);
result += temp;
return result;
}
template<typename SCALARTYPE, typename SparseMatrixType>
typename viennacl::enable_if< viennacl::is_any_sparse_matrix<SparseMatrixType>::value,
viennacl::vector<SCALARTYPE> >::type
operator-(viennacl::vector_base<SCALARTYPE> & result,
const viennacl::vector_expression< const SparseMatrixType, const viennacl::vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
{
assert(proxy.lhs().size1() == result.size() && bool("Dimensions for addition of sparse matrix-vector product to vector don't match!"));
vector<SCALARTYPE> temp(proxy.lhs().size1());
viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), temp);
result += temp;
return result;
}
} //namespace viennacl
#endif
viennacl::compressed_matrix::size2
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:929
viennacl::enable_if
Simple enable-if variant that uses the SFINAE pattern.
Definition: enable_if.hpp:30
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
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:376
viennacl::linalg::detail::row_info_types
row_info_types
Definition: forwards.h:837
viennacl::linalg::opencl::prod_impl
void prod_impl(const matrix_base< NumericT > &mat, bool trans_A, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:620
viennacl::memory_exception
Exception class in case of memory errors.
Definition: forwards.h:572
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
matrix.hpp
Implementation of the dense matrix class.
tools.hpp
Various little tools used here and there in ViennaCL.
sparse_matrix_operations.hpp
Implementations of operations using sparse matrices on the CPU using a single thread or OpenMP...
viennacl::linalg::host_based::inplace_solve
void inplace_solve(matrix_base< NumericT > const &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:130
viennacl::OPENCL_MEMORY
Definition: forwards.h:349
viennacl::matrix_base
Definition: matrix_def.hpp:103
viennacl::matrix_expression
Expression template class for representing a tree of expressions which ultimately result in a matrix...
Definition: forwards.h:341
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::linalg::cuda::detail::row_info
void row_info(compressed_matrix< NumericT, AligmentV > const &mat, vector_base< NumericT > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:107
viennacl::linalg::cuda::prod_impl
void prod_impl(const matrix_base< NumericT > &mat, bool mat_transpose, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:1464
viennacl::vector_expression
An expression template class that represents a binary operation that yields a vector.
Definition: forwards.h:239
viennacl::matrix_expression::size1
vcl_size_t size1() const
Returns the size of the result vector.
Definition: matrix.hpp:72
viennacl::operator-
viennacl::vector< NumericT > operator-(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 - A * v2', where A is a matrix.
Definition: matrix_operations.hpp:1200
viennacl::MEMORY_NOT_INITIALIZED
Definition: forwards.h:347
viennacl::CUDA_MEMORY
Definition: forwards.h:350
viennacl::linalg::host_based::prod_impl
void prod_impl(const matrix_base< NumericT > &mat, bool trans, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:993
sparse_matrix_operations.hpp
Implementations of operations using sparse matrices using CUDA.
viennacl::matrix_expression::size2
vcl_size_t size2() const
Definition: matrix.hpp:73
viennacl::vector_base
Common base class for dense vectors, vector ranges, and vector slices.
Definition: vector_def.hpp:104
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::linalg::host_based::detail::block_inplace_solve
void block_inplace_solve(const matrix_expression< const compressed_matrix< NumericT, AlignmentV >, const compressed_matrix< NumericT, AlignmentV >, op_trans > &L, viennacl::backend::mem_handle const &, vcl_size_t, vector_base< NumericT > const &, vector_base< NumericT > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:824
viennacl::vector
Definition: forwards.h:266
viennacl::linalg::cuda::detail::block_inplace_solve
void block_inplace_solve(const matrix_expression< const compressed_matrix< NumericT, AlignmentV >, const compressed_matrix< NumericT, AlignmentV >, op_trans > &L, viennacl::backend::mem_handle const &block_indices, vcl_size_t num_blocks, vector_base< NumericT > const &, vector_base< NumericT > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:970
viennacl::matrix_base::size1
size_type size1() const
Returns the number of rows.
Definition: matrix_def.hpp:224
sparse_matrix_operations.hpp
Implementations of operations using sparse matrices and OpenCL.
viennacl::MAIN_MEMORY
Definition: forwards.h:348
viennacl::linalg::opencl::inplace_solve
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for dense triangular systems. Matlab notation: A \ B.
Definition: direct_solve.hpp:77
viennacl::operator+
viennacl::vector< NumericT > operator+(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 + A * v2', where A is a matrix.
Definition: matrix_operations.hpp:1182
viennacl::op_prod
A tag class representing matrix-vector products and element-wise multiplications. ...
Definition: forwards.h:94
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::op_trans
A tag class representing transposed matrices.
Definition: forwards.h:220
viennacl::linalg::host_based::detail::row_info
void row_info(compressed_matrix< NumericT, AlignmentV > const &mat, vector_base< NumericT > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:53
viennacl::compressed_matrix< NumericT >
viennacl::linalg::prod_impl
void prod_impl(const matrix_base< NumericT > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:438
viennacl::traits::handle
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
viennacl::linalg::opencl::detail::block_inplace_solve
void block_inplace_solve(const matrix_expression< const compressed_matrix< NumericT, AlignmentV >, const compressed_matrix< NumericT, AlignmentV >, op_trans > &L, viennacl::backend::mem_handle const &block_indices, vcl_size_t num_blocks, vector_base< NumericT > const &, vector_base< NumericT > &x, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:483
viennacl::matrix_expression::lhs
LHS & lhs() const
Get left hand side operand.
Definition: matrix.hpp:66
viennacl::linalg::detail::row_info
viennacl::enable_if< viennacl::is_any_sparse_matrix< SparseMatrixType >::value >::type row_info(SparseMatrixType const &mat, vector< SCALARTYPE, VEC_ALIGNMENT > &vec, row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:50
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::linalg::opencl::detail::row_info
void row_info(compressed_matrix< NumericT, AlignmentV > const &A, vector_base< NumericT > &x, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:56
viennacl::backend::mem_handle::get_active_handle_id
memory_types get_active_handle_id() const
Returns an ID for the currently active memory buffer. Other memory buffers might contain old or no da...
Definition: mem_handle.hpp:118
viennacl::linalg::detail::block_inplace_solve
viennacl::enable_if< viennacl::is_any_sparse_matrix< SparseMatrixType >::value >::type block_inplace_solve(const matrix_expression< const SparseMatrixType, const SparseMatrixType, op_trans > &mat, viennacl::backend::mem_handle const &block_index_array, vcl_size_t num_blocks, viennacl::vector_base< ScalarType > const &mat_diagonal, viennacl::vector_base< ScalarType > &vec, SOLVERTAG tag)
Definition: sparse_matrix_operations.hpp:332