All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ell_matrix.hpp
Go to the documentation of this file.
#ifndef VIENNACL_ELL_MATRIX_HPP_
#define VIENNACL_ELL_MATRIX_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/vector.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/sparse_matrix_operations.hpp"
namespace viennacl
{
template<typename NumericT, unsigned int AlignmentV /* see forwards.h for default argument */>
class ell_matrix
{
public:
typedef viennacl::backend::mem_handle handle_type;
typedef scalar<typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT<NumericT>::ResultType> value_type;
typedef vcl_size_t size_type;
ell_matrix() : rows_(0), cols_(0), maxnnz_(0) {}
ell_matrix(viennacl::context ctx) : rows_(0), cols_(0), maxnnz_(0)
{
coords_.switch_active_handle_id(ctx.memory_type());
elements_.switch_active_handle_id(ctx.memory_type());
#ifdef VIENNACL_WITH_OPENCL
if (ctx.memory_type() == OPENCL_MEMORY)
{
coords_.opencl_handle().context(ctx.opencl_context());
elements_.opencl_handle().context(ctx.opencl_context());
}
#endif
}
void clear()
{
maxnnz_ = 0;
viennacl::backend::typesafe_host_array<unsigned int> host_coords_buffer(coords_, internal_size1());
std::vector<NumericT> host_elements(internal_size1());
viennacl::backend::memory_create(coords_, host_coords_buffer.element_size() * internal_size1(), viennacl::traits::context(coords_), host_coords_buffer.get());
viennacl::backend::memory_create(elements_, sizeof(NumericT) * internal_size1(), viennacl::traits::context(elements_), &(host_elements[0]));
}
vcl_size_t internal_size1() const { return viennacl::tools::align_to_multiple<vcl_size_t>(rows_, AlignmentV); }
vcl_size_t internal_size2() const { return viennacl::tools::align_to_multiple<vcl_size_t>(cols_, AlignmentV); }
vcl_size_t size1() const { return rows_; }
vcl_size_t size2() const { return cols_; }
vcl_size_t internal_maxnnz() const {return viennacl::tools::align_to_multiple<vcl_size_t>(maxnnz_, AlignmentV); }
vcl_size_t maxnnz() const { return maxnnz_; }
vcl_size_t nnz() const { return rows_ * maxnnz_; }
vcl_size_t internal_nnz() const { return internal_size1() * internal_maxnnz(); }
handle_type & handle() { return elements_; }
const handle_type & handle() const { return elements_; }
handle_type & handle2() { return coords_; }
const handle_type & handle2() const { return coords_; }
#if defined(_MSC_VER) && _MSC_VER < 1500 //Visual Studio 2005 needs special treatment
template<typename CPUMatrixT>
friend void copy(const CPUMatrixT & cpu_matrix, ell_matrix & gpu_matrix );
#else
template<typename CPUMatrixT, typename T, unsigned int ALIGN>
friend void copy(const CPUMatrixT & cpu_matrix, ell_matrix<T, ALIGN> & gpu_matrix );
#endif
private:
vcl_size_t rows_;
vcl_size_t cols_;
vcl_size_t maxnnz_;
handle_type coords_;
handle_type elements_;
};
template<typename CPUMatrixT, typename NumericT, unsigned int AlignmentV>
void copy(const CPUMatrixT& cpu_matrix, ell_matrix<NumericT, AlignmentV>& gpu_matrix )
{
assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
if (cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0)
{
//determine max capacity for row
vcl_size_t max_entries_per_row = 0;
for (typename CPUMatrixT::const_iterator1 row_it = cpu_matrix.begin1(); row_it != cpu_matrix.end1(); ++row_it)
{
vcl_size_t num_entries = 0;
for (typename CPUMatrixT::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
++num_entries;
max_entries_per_row = std::max(max_entries_per_row, num_entries);
}
//setup GPU matrix
gpu_matrix.maxnnz_ = max_entries_per_row;
gpu_matrix.rows_ = cpu_matrix.size1();
gpu_matrix.cols_ = cpu_matrix.size2();
vcl_size_t nnz = gpu_matrix.internal_nnz();
viennacl::backend::typesafe_host_array<unsigned int> coords(gpu_matrix.handle2(), nnz);
std::vector<NumericT> elements(nnz, 0);
// std::cout << "ELL_MATRIX copy " << gpu_matrix.maxnnz_ << " " << gpu_matrix.rows_ << " " << gpu_matrix.cols_ << " "
// << gpu_matrix.internal_maxnnz() << "\n";
for (typename CPUMatrixT::const_iterator1 row_it = cpu_matrix.begin1(); row_it != cpu_matrix.end1(); ++row_it)
{
vcl_size_t data_index = 0;
for (typename CPUMatrixT::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
{
coords.set(gpu_matrix.internal_size1() * data_index + col_it.index1(), col_it.index2());
elements[gpu_matrix.internal_size1() * data_index + col_it.index1()] = *col_it;
//std::cout << *col_it << "\n";
data_index++;
}
}
viennacl::backend::memory_create(gpu_matrix.handle2(), coords.raw_size(), traits::context(gpu_matrix.handle2()), coords.get());
viennacl::backend::memory_create(gpu_matrix.handle(), sizeof(NumericT) * elements.size(), traits::context(gpu_matrix.handle()), &(elements[0]));
}
}
template<typename IndexT, typename NumericT, unsigned int AlignmentV>
void copy(std::vector< std::map<IndexT, NumericT> > const & cpu_matrix,
ell_matrix<NumericT, AlignmentV> & gpu_matrix)
{
vcl_size_t max_col = 0;
for (vcl_size_t i=0; i<cpu_matrix.size(); ++i)
{
if (cpu_matrix[i].size() > 0)
max_col = std::max<vcl_size_t>(max_col, (cpu_matrix[i].rbegin())->first);
}
viennacl::copy(tools::const_sparse_matrix_adapter<NumericT, IndexT>(cpu_matrix, cpu_matrix.size(), max_col + 1), gpu_matrix);
}
template<typename CPUMatrixT, typename NumericT, unsigned int AlignmentV>
void copy(const ell_matrix<NumericT, AlignmentV>& gpu_matrix, CPUMatrixT& cpu_matrix)
{
assert( (viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
assert( (viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
if (gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0)
{
std::vector<NumericT> elements(gpu_matrix.internal_nnz());
viennacl::backend::typesafe_host_array<unsigned int> coords(gpu_matrix.handle2(), gpu_matrix.internal_nnz());
viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(NumericT) * elements.size(), &(elements[0]));
viennacl::backend::memory_read(gpu_matrix.handle2(), 0, coords.raw_size(), coords.get());
for (vcl_size_t row = 0; row < gpu_matrix.size1(); row++)
{
for (vcl_size_t ind = 0; ind < gpu_matrix.internal_maxnnz(); ind++)
{
vcl_size_t offset = gpu_matrix.internal_size1() * ind + row;
NumericT val = elements[offset];
if (val <= 0 && val >= 0) // val == 0 without compiler warnings
continue;
if (coords[offset] >= gpu_matrix.size2())
{
std::cerr << "ViennaCL encountered invalid data " << offset << " " << ind << " " << row << " " << coords[offset] << " " << gpu_matrix.size2() << std::endl;
return;
}
cpu_matrix(row, coords[offset]) = val;
}
}
}
}
template<typename NumericT, unsigned int AlignmentV, typename IndexT>
void copy(const ell_matrix<NumericT, AlignmentV> & gpu_matrix,
std::vector< std::map<IndexT, NumericT> > & cpu_matrix)
{
if (cpu_matrix.size() == 0)
cpu_matrix.resize(gpu_matrix.size1());
assert(cpu_matrix.size() == gpu_matrix.size1() && bool("Matrix dimension mismatch!"));
tools::sparse_matrix_adapter<NumericT, IndexT> temp(cpu_matrix, gpu_matrix.size1(), gpu_matrix.size2());
viennacl::copy(gpu_matrix, temp);
}
//
// Specify available operations:
//
namespace linalg
{
namespace detail
{
// x = A * y
template<typename T, unsigned int A>
struct op_executor<vector_base<T>, op_assign, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
{
// check for the special case x = A * x
if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
{
viennacl::vector<T> temp(lhs);
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
lhs = temp;
}
else
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), lhs, T(0));
}
};
template<typename T, unsigned int A>
struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
{
// check for the special case x += A * x
if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
{
viennacl::vector<T> temp(lhs);
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
lhs += temp;
}
else
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), lhs, T(1));
}
};
template<typename T, unsigned int A>
struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
{
// check for the special case x -= A * x
if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
{
viennacl::vector<T> temp(lhs);
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
lhs -= temp;
}
else
viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(-1), lhs, T(1));
}
};
// x = A * vec_op
template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
struct op_executor<vector_base<T>, op_assign, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
{
viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
}
};
// x = A * vec_op
template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
{
viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
viennacl::vector<T> temp_result(lhs);
viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
lhs += temp_result;
}
};
// x = A * vec_op
template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
{
static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
{
viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
viennacl::vector<T> temp_result(lhs);
viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
lhs -= temp_result;
}
};
} // namespace detail
} // namespace linalg
}
#endif
viennacl::ell_matrix::value_type
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< NumericT >::ResultType > value_type
Definition: ell_matrix.hpp:57
viennacl::ell_matrix::size1
vcl_size_t size1() const
Definition: ell_matrix.hpp:91
viennacl::backend::typesafe_host_array
Helper class implementing an array on the host. Default case: No conversion necessary.
Definition: util.hpp:92
viennacl::backend::typesafe_host_array::element_size
vcl_size_t element_size() const
Definition: util.hpp:112
viennacl::ell_matrix::handle2
handle_type & handle2()
Definition: ell_matrix.hpp:103
viennacl::scalar
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
viennacl::ell_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: ell_matrix.hpp:77
viennacl::ell_matrix::handle2
const handle_type & handle2() const
Definition: ell_matrix.hpp:104
viennacl::backend::typesafe_host_array::size
vcl_size_t size() const
Definition: util.hpp:113
tools.hpp
Various little tools used here and there in ViennaCL.
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::OPENCL_MEMORY
Definition: forwards.h:349
viennacl::ell_matrix::nnz
vcl_size_t nnz() const
Definition: ell_matrix.hpp:97
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::ell_matrix::size2
vcl_size_t size2() const
Definition: ell_matrix.hpp:92
viennacl::ell_matrix::ell_matrix
ell_matrix()
Definition: ell_matrix.hpp:60
viennacl::backend::typesafe_host_array::get
void * get()
Definition: util.hpp:110
viennacl::backend::memory_read
void memory_read(mem_handle const &src_buffer, vcl_size_t src_offset, vcl_size_t bytes_to_read, void *ptr, bool async=false)
Reads data from a buffer back to main RAM.
Definition: memory.hpp:261
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
viennacl::ell_matrix::handle_type
viennacl::backend::mem_handle handle_type
Definition: ell_matrix.hpp:56
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
viennacl::ell_matrix::internal_size1
vcl_size_t internal_size1() const
Definition: ell_matrix.hpp:88
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::ell_matrix::handle
const handle_type & handle() const
Definition: ell_matrix.hpp:101
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::ell_matrix::internal_nnz
vcl_size_t internal_nnz() const
Definition: ell_matrix.hpp:98
viennacl::ell_matrix
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
viennacl::ell_matrix::copy
friend void copy(const CPUMatrixT &cpu_matrix, ell_matrix< T, ALIGN > &gpu_matrix)
sparse_matrix_operations.hpp
Implementations of operations using sparse matrices.
viennacl::tools::const_sparse_matrix_adapter
Adapts a constant sparse matrix type made up from std::vector<std::map<SizeT, NumericT> > to basic ub...
Definition: adapter.hpp:183
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::vector
Definition: forwards.h:266
viennacl::ell_matrix::maxnnz
vcl_size_t maxnnz() const
Definition: ell_matrix.hpp:95
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::context::memory_type
viennacl::memory_types memory_type() const
Definition: context.hpp:76
viennacl::ell_matrix::handle
handle_type & handle()
Definition: ell_matrix.hpp:100
viennacl::backend::mem_handle::switch_active_handle_id
void switch_active_handle_id(memory_types new_id)
Switches the currently active handle. If no support for that backend is provided, an exception is thr...
Definition: mem_handle.hpp:121
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
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
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::backend::mem_handle::raw_size
vcl_size_t raw_size() const
Returns the number of bytes of the currently active buffer.
Definition: mem_handle.hpp:230
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::backend::memory_create
void memory_create(mem_handle &handle, vcl_size_t size_in_bytes, viennacl::context const &ctx, const void *host_ptr=NULL)
Creates an array of the specified size. If the second argument is provided, the buffer is initialized...
Definition: memory.hpp:87
viennacl::ell_matrix::internal_size2
vcl_size_t internal_size2() const
Definition: ell_matrix.hpp:89
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::ell_matrix::size_type
vcl_size_t size_type
Definition: ell_matrix.hpp:58
viennacl::ell_matrix::internal_maxnnz
vcl_size_t internal_maxnnz() const
Definition: ell_matrix.hpp:94
op_assign
Definition: self_assign.cpp:132
viennacl::ell_matrix::ell_matrix
ell_matrix(viennacl::context ctx)
Definition: ell_matrix.hpp:62