ViennaCL - The Vienna Computing Library  1.5.2
coordinate_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_COORDINATE_MATRIX_HPP_
2 #define VIENNACL_COORDINATE_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 
25 #include <map>
26 #include <vector>
27 #include <list>
28 
29 #include "viennacl/forwards.h"
30 #include "viennacl/vector.hpp"
31 
33 
34 namespace viennacl
35 {
36 
37 
38  //provide copy-operation:
46  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
47  void copy(const CPU_MATRIX & cpu_matrix,
49  {
50  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
51  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
52 
53  vcl_size_t group_num = 64;
54 
55  // Step 1: Determine nonzeros:
56  if ( cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0 )
57  {
58  vcl_size_t num_entries = 0;
59  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
60  row_it != cpu_matrix.end1();
61  ++row_it)
62  {
63  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
64  col_it != row_it.end();
65  ++col_it)
66  {
67  ++num_entries;
68  }
69  }
70 
71  // Step 2: Set up matrix data:
72  gpu_matrix.nonzeros_ = num_entries;
73  gpu_matrix.rows_ = cpu_matrix.size1();
74  gpu_matrix.cols_ = cpu_matrix.size2();
75 
76  viennacl::backend::typesafe_host_array<unsigned int> group_boundaries(gpu_matrix.handle3(), group_num + 1);
77  viennacl::backend::typesafe_host_array<unsigned int> coord_buffer(gpu_matrix.handle12(), 2*gpu_matrix.internal_nnz());
78  std::vector<SCALARTYPE> elements(gpu_matrix.internal_nnz());
79 
80  vcl_size_t data_index = 0;
81  vcl_size_t current_fraction = 0;
82 
83  group_boundaries.set(0, 0);
84  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
85  row_it != cpu_matrix.end1();
86  ++row_it)
87  {
88  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
89  col_it != row_it.end();
90  ++col_it)
91  {
92  coord_buffer.set(2*data_index, col_it.index1());
93  coord_buffer.set(2*data_index + 1, col_it.index2());
94  elements[data_index] = *col_it;
95  ++data_index;
96  }
97 
98  while (data_index > (current_fraction + 1) / static_cast<double>(group_num) * num_entries) //split data equally over 64 groups
99  group_boundaries.set(++current_fraction, data_index);
100  }
101 
102  //write end of last group:
103  group_boundaries.set(group_num, data_index);
104  //group_boundaries[1] = data_index; //for one compute unit
105 
106  //std::cout << "Group boundaries: " << std::endl;
107  //for (vcl_size_t i=0; i<group_boundaries.size(); ++i)
108  // std::cout << group_boundaries[i] << std::endl;
109 
110  viennacl::backend::memory_create(gpu_matrix.group_boundaries_, group_boundaries.raw_size(), traits::context(gpu_matrix.group_boundaries_), group_boundaries.get());
111  viennacl::backend::memory_create(gpu_matrix.coord_buffer_, coord_buffer.raw_size(), traits::context(gpu_matrix.coord_buffer_), coord_buffer.get());
112  viennacl::backend::memory_create(gpu_matrix.elements_, sizeof(SCALARTYPE)*elements.size(), traits::context(gpu_matrix.elements_), &(elements[0]));
113  }
114  }
115 
121  template <typename SCALARTYPE, unsigned int ALIGNMENT>
122  void copy(const std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix,
124  {
125  copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(cpu_matrix, cpu_matrix.size(), cpu_matrix.size()), gpu_matrix);
126  }
127 
128  //gpu to cpu:
138  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
140  CPU_MATRIX & cpu_matrix )
141  {
142  assert( (viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
143  assert( (viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
144 
145  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
146  {
147  //get raw data from memory:
148  viennacl::backend::typesafe_host_array<unsigned int> coord_buffer(gpu_matrix.handle12(), 2*gpu_matrix.nnz());
149  std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
150 
151  //std::cout << "GPU nonzeros: " << gpu_matrix.nnz() << std::endl;
152 
153  viennacl::backend::memory_read(gpu_matrix.handle12(), 0, coord_buffer.raw_size(), coord_buffer.get());
154  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE) * elements.size(), &(elements[0]));
155 
156  //fill the cpu_matrix:
157  for (vcl_size_t index = 0; index < gpu_matrix.nnz(); ++index)
158  cpu_matrix(coord_buffer[2*index], coord_buffer[2*index+1]) = elements[index];
159 
160  }
161  }
162 
168  template <typename SCALARTYPE, unsigned int ALIGNMENT>
170  std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix)
171  {
172  tools::sparse_matrix_adapter<SCALARTYPE> temp(cpu_matrix, gpu_matrix.size1(), gpu_matrix.size2());
173  copy(gpu_matrix, temp);
174  }
175 
176 
178 
185  template<class SCALARTYPE, unsigned int ALIGNMENT /* see forwards.h */ >
187  {
188  public:
192 
194  coordinate_matrix() : rows_(0), cols_(0), nonzeros_(0), group_num_(64) {}
195 
196  explicit coordinate_matrix(viennacl::context ctx) : rows_(0), cols_(0), nonzeros_(0), group_num_(64)
197  {
198  group_boundaries_.switch_active_handle_id(ctx.memory_type());
199  coord_buffer_.switch_active_handle_id(ctx.memory_type());
200  elements_.switch_active_handle_id(ctx.memory_type());
201 
202 #ifdef VIENNACL_WITH_OPENCL
203  if (ctx.memory_type() == OPENCL_MEMORY)
204  {
205  group_boundaries_.opencl_handle().context(ctx.opencl_context());
206  coord_buffer_.opencl_handle().context(ctx.opencl_context());
207  elements_.opencl_handle().context(ctx.opencl_context());
208  }
209 #endif
210  }
211 
220  rows_(rows), cols_(cols), nonzeros_(nonzeros)
221  {
222  if (nonzeros > 0)
223  {
226  viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * internal_nnz(), ctx);
227  }
228  else
229  {
230  group_boundaries_.switch_active_handle_id(ctx.memory_type());
231  coord_buffer_.switch_active_handle_id(ctx.memory_type());
232  elements_.switch_active_handle_id(ctx.memory_type());
233 
234  #ifdef VIENNACL_WITH_OPENCL
235  if (ctx.memory_type() == OPENCL_MEMORY)
236  {
237  group_boundaries_.opencl_handle().context(ctx.opencl_context());
238  coord_buffer_.opencl_handle().context(ctx.opencl_context());
239  elements_.opencl_handle().context(ctx.opencl_context());
240  }
241  #endif
242  }
243  }
244 
252  : rows_(rows), cols_(cols), nonzeros_(0)
253  {
254  group_boundaries_.switch_active_handle_id(ctx.memory_type());
255  coord_buffer_.switch_active_handle_id(ctx.memory_type());
256  elements_.switch_active_handle_id(ctx.memory_type());
257 
258 #ifdef VIENNACL_WITH_OPENCL
259  if (ctx.memory_type() == OPENCL_MEMORY)
260  {
261  group_boundaries_.opencl_handle().context(ctx.opencl_context());
262  coord_buffer_.opencl_handle().context(ctx.opencl_context());
263  elements_.opencl_handle().context(ctx.opencl_context());
264  }
265 #endif
266  }
267 
268 
270  void reserve(vcl_size_t new_nonzeros)
271  {
272  if (new_nonzeros > nonzeros_) //TODO: Do we need to initialize new memory with zero?
273  {
274  handle_type coord_buffer_old;
275  handle_type elements_old;
276  viennacl::backend::memory_shallow_copy(coord_buffer_, coord_buffer_old);
277  viennacl::backend::memory_shallow_copy(elements_, elements_old);
278 
279  vcl_size_t internal_new_nnz = viennacl::tools::align_to_multiple<vcl_size_t>(new_nonzeros, ALIGNMENT);
280  viennacl::backend::typesafe_host_array<unsigned int> size_deducer(coord_buffer_);
281  viennacl::backend::memory_create(coord_buffer_, size_deducer.element_size() * 2 * internal_new_nnz, viennacl::traits::context(coord_buffer_));
282  viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * internal_new_nnz, viennacl::traits::context(elements_));
283 
284  viennacl::backend::memory_copy(coord_buffer_old, coord_buffer_, 0, 0, size_deducer.element_size() * 2 * nonzeros_);
285  viennacl::backend::memory_copy(elements_old, elements_, 0, 0, sizeof(SCALARTYPE) * nonzeros_);
286 
287  nonzeros_ = new_nonzeros;
288  }
289  }
290 
297  void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve = true)
298  {
299  assert (new_size1 > 0 && new_size2 > 0);
300 
301  if (new_size1 < rows_ || new_size2 < cols_) //enlarge buffer
302  {
303  std::vector<std::map<unsigned int, SCALARTYPE> > stl_sparse_matrix;
304  if (rows_ > 0)
305  stl_sparse_matrix.resize(rows_);
306 
307  if (preserve && rows_ > 0)
308  viennacl::copy(*this, stl_sparse_matrix);
309 
310  stl_sparse_matrix.resize(new_size1);
311 
312  //std::cout << "Cropping STL matrix of size " << stl_sparse_matrix.size() << std::endl;
313  if (new_size2 < cols_ && rows_ > 0)
314  {
315  for (vcl_size_t i=0; i<stl_sparse_matrix.size(); ++i)
316  {
317  std::list<unsigned int> to_delete;
318  for (typename std::map<unsigned int, SCALARTYPE>::iterator it = stl_sparse_matrix[i].begin();
319  it != stl_sparse_matrix[i].end();
320  ++it)
321  {
322  if (it->first >= new_size2)
323  to_delete.push_back(it->first);
324  }
325 
326  for (std::list<unsigned int>::iterator it = to_delete.begin(); it != to_delete.end(); ++it)
327  stl_sparse_matrix[i].erase(*it);
328  }
329  //std::cout << "Cropping done..." << std::endl;
330  }
331 
332  rows_ = new_size1;
333  cols_ = new_size2;
334  viennacl::copy(stl_sparse_matrix, *this);
335  }
336 
337  rows_ = new_size1;
338  cols_ = new_size2;
339  }
340 
341 
343  vcl_size_t size1() const { return rows_; }
345  vcl_size_t size2() const { return cols_; }
347  vcl_size_t nnz() const { return nonzeros_; }
349  vcl_size_t internal_nnz() const { return viennacl::tools::align_to_multiple<vcl_size_t>(nonzeros_, ALIGNMENT); }
350 
352  const handle_type & handle12() const { return coord_buffer_; }
354  const handle_type & handle() const { return elements_; }
356  const handle_type & handle3() const { return group_boundaries_; }
357 
358  vcl_size_t groups() const { return group_num_; }
359 
360  #if defined(_MSC_VER) && _MSC_VER < 1500 //Visual Studio 2005 needs special treatment
361  template <typename CPU_MATRIX>
362  friend void copy(const CPU_MATRIX & cpu_matrix, coordinate_matrix & gpu_matrix );
363  #else
364  template <typename CPU_MATRIX, typename SCALARTYPE2, unsigned int ALIGNMENT2>
365  friend void copy(const CPU_MATRIX & cpu_matrix, coordinate_matrix<SCALARTYPE2, ALIGNMENT2> & gpu_matrix );
366  #endif
367 
368  private:
371 
373  coordinate_matrix & operator=(coordinate_matrix const &);
374 
375 
376  vcl_size_t rows_;
377  vcl_size_t cols_;
378  vcl_size_t nonzeros_;
379  vcl_size_t group_num_;
380  handle_type coord_buffer_;
381  handle_type elements_;
382  handle_type group_boundaries_;
383  };
384 
385 
386  //
387  // Specify available operations:
388  //
389 
392  namespace linalg
393  {
394  namespace detail
395  {
396  // x = A * y
397  template <typename T, unsigned int A>
398  struct op_executor<vector_base<T>, op_assign, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> >
399  {
400  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
401  {
402  // check for the special case x = A * x
403  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
404  {
405  viennacl::vector<T> temp(lhs);
406  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
407  lhs = temp;
408  }
409  else
410  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
411  }
412  };
413 
414  template <typename T, unsigned int A>
415  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> >
416  {
417  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
418  {
419  viennacl::vector<T> temp(lhs);
420  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
421  lhs += temp;
422  }
423  };
424 
425  template <typename T, unsigned int A>
426  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> >
427  {
428  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
429  {
430  viennacl::vector<T> temp(lhs);
431  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
432  lhs -= temp;
433  }
434  };
435 
436 
437  // x = A * vec_op
438  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
439  struct op_executor<vector_base<T>, op_assign, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
440  {
441  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
442  {
443  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
444  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
445  }
446  };
447 
448  // x += A * vec_op
449  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
450  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
451  {
452  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
453  {
454  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
455  viennacl::vector<T> temp_result(lhs);
456  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
457  lhs += temp_result;
458  }
459  };
460 
461  // x -= A * vec_op
462  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
463  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
464  {
465  static void apply(vector_base<T> & lhs, vector_expression<const coordinate_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
466  {
467  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
468  viennacl::vector<T> temp_result(lhs);
469  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
470  lhs -= temp_result;
471  }
472  };
473 
474  } // namespace detail
475  } // namespace linalg
476 
478 }
479 
480 #endif
Helper class implementing an array on the host. Default case: No conversion necessary.
Definition: util.hpp:95
std::size_t vcl_size_t
Definition: forwards.h:58
coordinate_matrix(vcl_size_t rows, vcl_size_t cols, viennacl::context ctx)
Construction of a coordinate matrix with the supplied number of rows and columns in the supplied cont...
Definition: coordinate_matrix.hpp:251
friend void copy(const CPU_MATRIX &cpu_matrix, coordinate_matrix< SCALARTYPE2, ALIGNMENT2 > &gpu_matrix)
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 size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:216
void reserve(vcl_size_t new_nonzeros)
Allocate memory for the supplied number of nonzeros in the matrix. Old values are preserved...
Definition: coordinate_matrix.hpp:270
vcl_size_t size2() const
Returns the number of columns.
Definition: coordinate_matrix.hpp:345
vcl_size_t size1() const
Returns the number of rows.
Definition: coordinate_matrix.hpp:343
This file provides the forward declarations for the main types used within ViennaCL.
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
void set(vcl_size_t index, U value)
Definition: util.hpp:145
vcl_size_t element_size(memory_types)
Definition: memory.hpp:299
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: coordinate_matrix.hpp:354
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
const handle_type & handle3() const
Returns the OpenCL handle to the group start index array.
Definition: coordinate_matrix.hpp:356
coordinate_matrix(viennacl::context ctx)
Definition: coordinate_matrix.hpp:196
Definition: forwards.h:480
vcl_size_t groups() const
Definition: coordinate_matrix.hpp:358
Implementations of operations using sparse matrices.
const handle_type & handle12() const
Returns the OpenCL handle to the (row, column) index array.
Definition: coordinate_matrix.hpp:352
coordinate_matrix(vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros=0, viennacl::context ctx=viennacl::context())
Construction of a coordinate matrix with the supplied number of rows and columns. If the number of no...
Definition: coordinate_matrix.hpp:219
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
Adapts a constant sparse matrix type made up from std::vector<std::map<SizeType, SCALARTYPE> > to bas...
Definition: adapter.hpp:176
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
coordinate_matrix()
Default construction of a coordinate matrix. No memory is allocated.
Definition: coordinate_matrix.hpp:194
void memory_copy(mem_handle const &src_buffer, mem_handle &dst_buffer, vcl_size_t src_offset, vcl_size_t dst_offset, vcl_size_t bytes_to_copy)
Copies 'bytes_to_copy' bytes from address 'src_buffer + src_offset' to memory starting at address 'ds...
Definition: memory.hpp:140
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
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
Adapts a non-const sparse matrix type made up from std::vector<std::map<SizeType, SCALARTYPE> > to ba...
Definition: adapter.hpp:345
vcl_size_t size_type
Definition: coordinate_matrix.hpp:191
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
void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve=true)
Resize the matrix.
Definition: coordinate_matrix.hpp:297
vcl_size_t internal_nnz() const
Returns the number of internal nonzero entries.
Definition: coordinate_matrix.hpp:349
void memory_shallow_copy(mem_handle const &src_buffer, mem_handle &dst_buffer)
A 'shallow' copy operation from an initialized buffer to an uninitialized buffer. The uninitialized b...
Definition: memory.hpp:177
viennacl::backend::mem_handle handle_type
Definition: coordinate_matrix.hpp:189
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:186
vcl_size_t nnz() const
Returns the number of nonzero entries.
Definition: coordinate_matrix.hpp:347
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: coordinate_matrix.hpp:190