ViennaCL - The Vienna Computing Library  1.5.2
ell_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_ELL_MATRIX_HPP_
2 #define VIENNACL_ELL_MATRIX_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2014, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8  Portions of this software are copyright by UChicago Argonne, LLC.
9 
10  -----------------
11  ViennaCL - The Vienna Computing Library
12  -----------------
13 
14  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
15 
16  (A list of authors and contributors can be found in the PDF manual)
17 
18  License: MIT (X11), see file LICENSE in the base directory
19 ============================================================================= */
20 
28 #include "viennacl/forwards.h"
29 #include "viennacl/vector.hpp"
30 
31 #include "viennacl/tools/tools.hpp"
32 
34 
35 namespace viennacl
36 {
52  template<typename SCALARTYPE, unsigned int ALIGNMENT /* see forwards.h for default argument */>
53  class ell_matrix
54  {
55  public:
59 
60  ell_matrix() : rows_(0), cols_(0), maxnnz_(0) {}
61 
62  ell_matrix(viennacl::context ctx) : rows_(0), cols_(0), maxnnz_(0)
63  {
64  coords_.switch_active_handle_id(ctx.memory_type());
65  elements_.switch_active_handle_id(ctx.memory_type());
66 
67 #ifdef VIENNACL_WITH_OPENCL
68  if (ctx.memory_type() == OPENCL_MEMORY)
69  {
70  coords_.opencl_handle().context(ctx.opencl_context());
71  elements_.opencl_handle().context(ctx.opencl_context());
72  }
73 #endif
74  }
75 
76  public:
77  vcl_size_t internal_size1() const { return viennacl::tools::align_to_multiple<vcl_size_t>(rows_, ALIGNMENT); }
78  vcl_size_t internal_size2() const { return viennacl::tools::align_to_multiple<vcl_size_t>(cols_, ALIGNMENT); }
79 
80  vcl_size_t size1() const { return rows_; }
81  vcl_size_t size2() const { return cols_; }
82 
83  vcl_size_t internal_maxnnz() const {return viennacl::tools::align_to_multiple<vcl_size_t>(maxnnz_, ALIGNMENT); }
84  vcl_size_t maxnnz() const { return maxnnz_; }
85 
86  vcl_size_t nnz() const { return rows_ * maxnnz_; }
88 
89  handle_type & handle() { return elements_; }
90  const handle_type & handle() const { return elements_; }
91 
92  handle_type & handle2() { return coords_; }
93  const handle_type & handle2() const { return coords_; }
94 
95  #if defined(_MSC_VER) && _MSC_VER < 1500 //Visual Studio 2005 needs special treatment
96  template <typename CPU_MATRIX>
97  friend void copy(const CPU_MATRIX & cpu_matrix, ell_matrix & gpu_matrix );
98  #else
99  template <typename CPU_MATRIX, typename T, unsigned int ALIGN>
100  friend void copy(const CPU_MATRIX & cpu_matrix, ell_matrix<T, ALIGN> & gpu_matrix );
101  #endif
102 
103  private:
104  vcl_size_t rows_;
105  vcl_size_t cols_;
106  vcl_size_t maxnnz_;
107 
108  handle_type coords_;
109  handle_type elements_;
110  };
111 
112  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
113  void copy(const CPU_MATRIX& cpu_matrix, ell_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix )
114  {
115  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
116  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
117 
118  if(cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0)
119  {
120  //determine max capacity for row
121  vcl_size_t max_entries_per_row = 0;
122  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1(); row_it != cpu_matrix.end1(); ++row_it)
123  {
124  vcl_size_t num_entries = 0;
125  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
126  {
127  ++num_entries;
128  }
129 
130  max_entries_per_row = std::max(max_entries_per_row, num_entries);
131  }
132 
133  //setup GPU matrix
134  gpu_matrix.maxnnz_ = max_entries_per_row;
135  gpu_matrix.rows_ = cpu_matrix.size1();
136  gpu_matrix.cols_ = cpu_matrix.size2();
137 
138  vcl_size_t nnz = gpu_matrix.internal_nnz();
139 
141  std::vector<SCALARTYPE> elements(nnz, 0);
142 
143  // std::cout << "ELL_MATRIX copy " << gpu_matrix.maxnnz_ << " " << gpu_matrix.rows_ << " " << gpu_matrix.cols_ << " "
144  // << gpu_matrix.internal_maxnnz() << "\n";
145 
146  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1(); row_it != cpu_matrix.end1(); ++row_it)
147  {
148  vcl_size_t data_index = 0;
149 
150  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin(); col_it != row_it.end(); ++col_it)
151  {
152  coords.set(gpu_matrix.internal_size1() * data_index + col_it.index1(), col_it.index2());
153  elements[gpu_matrix.internal_size1() * data_index + col_it.index1()] = *col_it;
154  //std::cout << *col_it << "\n";
155  data_index++;
156  }
157  }
158 
159  viennacl::backend::memory_create(gpu_matrix.handle2(), coords.raw_size(), traits::context(gpu_matrix.handle2()), coords.get());
160  viennacl::backend::memory_create(gpu_matrix.handle(), sizeof(SCALARTYPE) * elements.size(), traits::context(gpu_matrix.handle()), &(elements[0]));
161  }
162  }
163 
164  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
165  void copy(const ell_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix, CPU_MATRIX& cpu_matrix)
166  {
167  assert( (viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
168  assert( (viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
169 
170  if(gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0)
171  {
172  std::vector<SCALARTYPE> elements(gpu_matrix.internal_nnz());
174 
175  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE) * elements.size(), &(elements[0]));
176  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, coords.raw_size(), coords.get());
177 
178  for(vcl_size_t row = 0; row < gpu_matrix.size1(); row++)
179  {
180  for(vcl_size_t ind = 0; ind < gpu_matrix.internal_maxnnz(); ind++)
181  {
182  vcl_size_t offset = gpu_matrix.internal_size1() * ind + row;
183 
184  if(elements[offset] == static_cast<SCALARTYPE>(0.0))
185  continue;
186 
187  if(coords[offset] >= gpu_matrix.size2())
188  {
189  std::cerr << "ViennaCL encountered invalid data " << offset << " " << ind << " " << row << " " << coords[offset] << " " << gpu_matrix.size2() << std::endl;
190  return;
191  }
192 
193  cpu_matrix(row, coords[offset]) = elements[offset];
194  }
195  }
196  }
197  }
198 
199 
200  //
201  // Specify available operations:
202  //
203 
206  namespace linalg
207  {
208  namespace detail
209  {
210  // x = A * y
211  template <typename T, unsigned int A>
212  struct op_executor<vector_base<T>, op_assign, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
213  {
214  static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
215  {
216  // check for the special case x = A * x
217  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
218  {
219  viennacl::vector<T> temp(lhs);
220  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
221  lhs = temp;
222  }
223  else
224  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
225  }
226  };
227 
228  template <typename T, unsigned int A>
229  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
230  {
231  static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
232  {
233  viennacl::vector<T> temp(lhs);
234  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
235  lhs += temp;
236  }
237  };
238 
239  template <typename T, unsigned int A>
240  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> >
241  {
242  static void apply(vector_base<T> & lhs, vector_expression<const ell_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
243  {
244  viennacl::vector<T> temp(lhs);
245  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
246  lhs -= temp;
247  }
248  };
249 
250 
251  // x = A * vec_op
252  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
253  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> >
254  {
255  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)
256  {
257  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
258  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
259  }
260  };
261 
262  // x = A * vec_op
263  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
264  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> >
265  {
266  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)
267  {
268  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
269  viennacl::vector<T> temp_result(lhs);
270  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
271  lhs += temp_result;
272  }
273  };
274 
275  // x = A * vec_op
276  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
277  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> >
278  {
279  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)
280  {
281  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
282  viennacl::vector<T> temp_result(lhs);
283  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
284  lhs -= temp_result;
285  }
286  };
287 
288  } // namespace detail
289  } // namespace linalg
290 
292 }
293 
294 #endif
295 
296 
Helper class implementing an array on the host. Default case: No conversion necessary.
Definition: util.hpp:95
vcl_size_t internal_nnz() const
Definition: ell_matrix.hpp:87
std::size_t vcl_size_t
Definition: forwards.h:58
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:172
vcl_size_t size() const
Definition: util.hpp:163
const handle_type & handle() const
Definition: ell_matrix.hpp:90
Various little tools used here and there in ViennaCL.
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
vcl_size_t size1() const
Definition: ell_matrix.hpp:80
vcl_size_t size2() const
Definition: ell_matrix.hpp:81
This file provides the forward declarations for the main types used within ViennaCL.
friend void copy(const CPU_MATRIX &cpu_matrix, ell_matrix< T, ALIGN > &gpu_matrix)
vcl_size_t internal_maxnnz() const
Definition: ell_matrix.hpp:83
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
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:245
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
vcl_size_t internal_size2() const
Definition: ell_matrix.hpp:78
Definition: forwards.h:480
ell_matrix(viennacl::context ctx)
Definition: ell_matrix.hpp:62
const handle_type & handle2() const
Definition: ell_matrix.hpp:93
Sparse matrix class using the ELLPACK format for storing the nonzeros.
Definition: ell_matrix.hpp:53
handle_type & handle()
Definition: ell_matrix.hpp:89
Implementations of operations using sparse matrices.
viennacl::backend::mem_handle handle_type
Definition: ell_matrix.hpp:56
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &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::memory_types memory_type() const
Definition: context.hpp:76
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
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
vcl_size_t nnz() const
Definition: ell_matrix.hpp:86
ell_matrix()
Definition: ell_matrix.hpp:60
vcl_size_t raw_size() const
Returns the number of bytes of the currently active buffer.
Definition: mem_handle.hpp:203
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:41
vcl_size_t maxnnz() const
Definition: ell_matrix.hpp:84
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
handle_type & handle2()
Definition: ell_matrix.hpp:92
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: ell_matrix.hpp:57
vcl_size_t size_type
Definition: ell_matrix.hpp:58
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
void prod_impl(const matrix_base< NumericT, F > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
Definition: matrix_operations.hpp:350
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
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:94
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
vcl_size_t internal_size1() const
Definition: ell_matrix.hpp:77