ViennaCL - The Vienna Computing Library  1.5.2
compressed_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_COMPRESSED_MATRIX_HPP_
2 #define VIENNACL_COMPRESSED_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 <vector>
26 #include <list>
27 #include <map>
28 #include "viennacl/forwards.h"
29 #include "viennacl/vector.hpp"
30 
32 
33 #include "viennacl/tools/tools.hpp"
35 
36 namespace viennacl
37 {
38  namespace detail
39  {
40  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
41  void copy_impl(const CPU_MATRIX & cpu_matrix,
43  vcl_size_t nonzeros)
44  {
45  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
46  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
47 
48  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1);
49  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), nonzeros);
50  std::vector<SCALARTYPE> elements(nonzeros);
51 
52  vcl_size_t row_index = 0;
53  vcl_size_t data_index = 0;
54 
55  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
56  row_it != cpu_matrix.end1();
57  ++row_it)
58  {
59  row_buffer.set(row_index, data_index);
60  ++row_index;
61 
62  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
63  col_it != row_it.end();
64  ++col_it)
65  {
66  col_buffer.set(data_index, col_it.index2());
67  elements[data_index] = *col_it;
68  ++data_index;
69  }
70  data_index = viennacl::tools::align_to_multiple<vcl_size_t>(data_index, ALIGNMENT); //take care of alignment
71  }
72  row_buffer.set(row_index, data_index);
73 
74  gpu_matrix.set(row_buffer.get(),
75  col_buffer.get(),
76  &elements[0],
77  cpu_matrix.size1(),
78  cpu_matrix.size2(),
79  nonzeros);
80  }
81  }
82 
83  //provide copy-operation:
98  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
99  void copy(const CPU_MATRIX & cpu_matrix,
101  {
102  if ( cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0 )
103  {
104  //determine nonzeros:
105  vcl_size_t num_entries = 0;
106  for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
107  row_it != cpu_matrix.end1();
108  ++row_it)
109  {
110  vcl_size_t entries_per_row = 0;
111  for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
112  col_it != row_it.end();
113  ++col_it)
114  {
115  ++entries_per_row;
116  }
117  num_entries += viennacl::tools::align_to_multiple<vcl_size_t>(entries_per_row, ALIGNMENT);
118  }
119 
120  if (num_entries == 0) //we copy an empty matrix
121  num_entries = 1;
122 
123  //set up matrix entries:
124  viennacl::detail::copy_impl(cpu_matrix, gpu_matrix, num_entries);
125  }
126  }
127 
128 
129  //adapted for std::vector< std::map < > > argument:
135  template <typename SizeType, typename SCALARTYPE, unsigned int ALIGNMENT>
136  void copy(const std::vector< std::map<SizeType, SCALARTYPE> > & cpu_matrix,
138  {
139  vcl_size_t nonzeros = 0;
140  vcl_size_t max_col = 0;
141  for (vcl_size_t i=0; i<cpu_matrix.size(); ++i)
142  {
143  if (cpu_matrix[i].size() > 0)
144  nonzeros += ((cpu_matrix[i].size() - 1) / ALIGNMENT + 1) * ALIGNMENT;
145  if (cpu_matrix[i].size() > 0)
146  max_col = std::max<vcl_size_t>(max_col, (cpu_matrix[i].rbegin())->first);
147  }
148 
150  gpu_matrix,
151  nonzeros);
152  }
153 
154 #ifdef VIENNACL_WITH_UBLAS
155  template <typename ScalarType, typename F, vcl_size_t IB, typename IA, typename TA>
156  void copy(const boost::numeric::ublas::compressed_matrix<ScalarType, F, IB, IA, TA> & ublas_matrix,
158  {
159  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(ublas_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
160  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(ublas_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
161 
162  //we just need to copy the CSR arrays:
163  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), ublas_matrix.size1() + 1);
164  for (vcl_size_t i=0; i<=ublas_matrix.size1(); ++i)
165  row_buffer.set(i, ublas_matrix.index1_data()[i]);
166 
167  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), ublas_matrix.nnz());
168  for (vcl_size_t i=0; i<ublas_matrix.nnz(); ++i)
169  col_buffer.set(i, ublas_matrix.index2_data()[i]);
170 
171  gpu_matrix.set(row_buffer.get(),
172  col_buffer.get(),
173  &(ublas_matrix.value_data()[0]),
174  ublas_matrix.size1(),
175  ublas_matrix.size2(),
176  ublas_matrix.nnz());
177 
178  }
179 #endif
180 
181  #ifdef VIENNACL_WITH_EIGEN
182  template <typename SCALARTYPE, int flags, unsigned int ALIGNMENT>
183  void copy(const Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix,
184  compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix)
185  {
186  assert( (gpu_matrix.size1() == 0 || static_cast<vcl_size_t>(eigen_matrix.rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
187  assert( (gpu_matrix.size2() == 0 || static_cast<vcl_size_t>(eigen_matrix.cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
188 
189  std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(eigen_matrix.rows());
190 
191  for (int k=0; k < eigen_matrix.outerSize(); ++k)
192  for (typename Eigen::SparseMatrix<SCALARTYPE, flags>::InnerIterator it(eigen_matrix, k); it; ++it)
193  stl_matrix[it.row()][it.col()] = it.value();
194 
195  copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix, eigen_matrix.rows(), eigen_matrix.cols()), gpu_matrix);
196  }
197 #endif
198 
199 
200 #ifdef VIENNACL_WITH_MTL4
201  template <typename SCALARTYPE, unsigned int ALIGNMENT>
202  void copy(const mtl::compressed2D<SCALARTYPE> & cpu_matrix,
203  compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix)
204  {
205  assert( (gpu_matrix.size1() == 0 || static_cast<vcl_size_t>(cpu_matrix.num_rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
206  assert( (gpu_matrix.size2() == 0 || static_cast<vcl_size_t>(cpu_matrix.num_cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
207 
208  typedef mtl::compressed2D<SCALARTYPE> MatrixType;
209 
210  std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(cpu_matrix.num_rows());
211 
212  using mtl::traits::range_generator;
213  using mtl::traits::range::min;
214 
215  // Choose between row and column traversal
216  typedef typename min<range_generator<mtl::tag::row, MatrixType>,
217  range_generator<mtl::tag::col, MatrixType> >::type range_type;
218  range_type my_range;
219 
220  // Type of outer cursor
221  typedef typename range_type::type c_type;
222  // Type of inner cursor
223  typedef typename mtl::traits::range_generator<mtl::tag::nz, c_type>::type ic_type;
224 
225  // Define the property maps
226  typename mtl::traits::row<MatrixType>::type row(cpu_matrix);
227  typename mtl::traits::col<MatrixType>::type col(cpu_matrix);
228  typename mtl::traits::const_value<MatrixType>::type value(cpu_matrix);
229 
230  // Now iterate over the matrix
231  for (c_type cursor(my_range.begin(cpu_matrix)), cend(my_range.end(cpu_matrix)); cursor != cend; ++cursor)
232  for (ic_type icursor(mtl::begin<mtl::tag::nz>(cursor)), icend(mtl::end<mtl::tag::nz>(cursor)); icursor != icend; ++icursor)
233  stl_matrix[row(*icursor)][col(*icursor)] = value(*icursor);
234 
235  copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix, cpu_matrix.num_rows(), cpu_matrix.num_cols()), gpu_matrix);
236  }
237 #endif
238 
239 
240 
241 
242 
243 
244 
245  //
246  // gpu to cpu:
247  //
257  template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT>
259  CPU_MATRIX & cpu_matrix )
260  {
261  assert( (viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
262  assert( (viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
263 
264  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
265  {
266  //get raw data from memory:
267  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1);
268  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
269  std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
270 
271  //std::cout << "GPU->CPU, nonzeros: " << gpu_matrix.nnz() << std::endl;
272 
273  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
274  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
275  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0]));
276 
277  //fill the cpu_matrix:
278  vcl_size_t data_index = 0;
279  for (vcl_size_t row = 1; row <= gpu_matrix.size1(); ++row)
280  {
281  while (data_index < row_buffer[row])
282  {
283  if (col_buffer[data_index] >= gpu_matrix.size2())
284  {
285  std::cerr << "ViennaCL encountered invalid data at colbuffer[" << data_index << "]: " << col_buffer[data_index] << std::endl;
286  return;
287  }
288 
289  if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
290  cpu_matrix(row-1, static_cast<vcl_size_t>(col_buffer[data_index])) = elements[data_index];
291  ++data_index;
292  }
293  }
294  }
295  }
296 
297 
303  template <typename SCALARTYPE, unsigned int ALIGNMENT>
305  std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix)
306  {
307  tools::sparse_matrix_adapter<SCALARTYPE> temp(cpu_matrix, cpu_matrix.size(), cpu_matrix.size());
308  copy(gpu_matrix, temp);
309  }
310 
311 #ifdef VIENNACL_WITH_UBLAS
312  template <typename ScalarType, unsigned int ALIGNMENT, typename F, vcl_size_t IB, typename IA, typename TA>
314  boost::numeric::ublas::compressed_matrix<ScalarType> & ublas_matrix)
315  {
316  assert( (viennacl::traits::size1(ublas_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
317  assert( (viennacl::traits::size2(ublas_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
318 
319  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
320  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
321 
322  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
323  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
324 
325  ublas_matrix.clear();
326  ublas_matrix.reserve(gpu_matrix.nnz());
327 
328  ublas_matrix.set_filled(gpu_matrix.size1() + 1, gpu_matrix.nnz());
329 
330  for (vcl_size_t i=0; i<ublas_matrix.size1() + 1; ++i)
331  ublas_matrix.index1_data()[i] = row_buffer[i];
332 
333  for (vcl_size_t i=0; i<ublas_matrix.nnz(); ++i)
334  ublas_matrix.index2_data()[i] = col_buffer[i];
335 
336  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(ScalarType) * gpu_matrix.nnz(), &(ublas_matrix.value_data()[0]));
337 
338  }
339 #endif
340 
341 #ifdef VIENNACL_WITH_EIGEN
342  template <typename SCALARTYPE, int flags, unsigned int ALIGNMENT>
343  void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix,
344  Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix)
345  {
346  assert( (static_cast<vcl_size_t>(eigen_matrix.rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
347  assert( (static_cast<vcl_size_t>(eigen_matrix.cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
348 
349  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
350  {
351  //get raw data from memory:
352  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
353  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
354  std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
355 
356  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
357  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
358  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0]));
359 
360  eigen_matrix.setZero();
361  vcl_size_t data_index = 0;
362  for (vcl_size_t row = 1; row <= gpu_matrix.size1(); ++row)
363  {
364  while (data_index < row_buffer[row])
365  {
366  assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer"));
367  if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
368  eigen_matrix.insert(row-1, col_buffer[data_index]) = elements[data_index];
369  ++data_index;
370  }
371  }
372  }
373  }
374 #endif
375 
376 
377 
378 #ifdef VIENNACL_WITH_MTL4
379  template <typename SCALARTYPE, unsigned int ALIGNMENT>
380  void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix,
381  mtl::compressed2D<SCALARTYPE> & mtl4_matrix)
382  {
383  assert( (static_cast<vcl_size_t>(mtl4_matrix.num_rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
384  assert( (static_cast<vcl_size_t>(mtl4_matrix.num_cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
385 
386  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
387  {
388 
389  //get raw data from memory:
390  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
391  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
392  std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
393 
394  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
395  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
396  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0]));
397 
398  //set_to_zero(mtl4_matrix);
399  //mtl4_matrix.change_dim(gpu_matrix.size1(), gpu_matrix.size2());
400 
401  mtl::matrix::inserter< mtl::compressed2D<SCALARTYPE> > ins(mtl4_matrix);
402  vcl_size_t data_index = 0;
403  for (vcl_size_t row = 1; row <= gpu_matrix.size1(); ++row)
404  {
405  while (data_index < row_buffer[row])
406  {
407  assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer"));
408  if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
409  ins(row-1, col_buffer[data_index]) << typename mtl::Collection< mtl::compressed2D<SCALARTYPE> >::value_type(elements[data_index]);
410  ++data_index;
411  }
412  }
413  }
414  }
415 #endif
416 
417 
418 
419 
420 
422 
427  template<class SCALARTYPE, unsigned int ALIGNMENT /* see VCLForwards.h */>
429  {
430  public:
434 
436  compressed_matrix() : rows_(0), cols_(0), nonzeros_(0) {}
437 
446  : rows_(rows), cols_(cols), nonzeros_(nonzeros)
447  {
448  row_buffer_.switch_active_handle_id(ctx.memory_type());
449  col_buffer_.switch_active_handle_id(ctx.memory_type());
450  elements_.switch_active_handle_id(ctx.memory_type());
451 
452 #ifdef VIENNACL_WITH_OPENCL
453  if (ctx.memory_type() == OPENCL_MEMORY)
454  {
455  row_buffer_.opencl_handle().context(ctx.opencl_context());
456  col_buffer_.opencl_handle().context(ctx.opencl_context());
457  elements_.opencl_handle().context(ctx.opencl_context());
458  }
459 #endif
460  if (rows > 0)
461  {
463  }
464  if (nonzeros > 0)
465  {
467  viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * nonzeros, ctx);
468  }
469  }
470 
478  : rows_(rows), cols_(cols), nonzeros_(0)
479  {
480  row_buffer_.switch_active_handle_id(ctx.memory_type());
481  col_buffer_.switch_active_handle_id(ctx.memory_type());
482  elements_.switch_active_handle_id(ctx.memory_type());
483 
484 #ifdef VIENNACL_WITH_OPENCL
485  if (ctx.memory_type() == OPENCL_MEMORY)
486  {
487  row_buffer_.opencl_handle().context(ctx.opencl_context());
488  col_buffer_.opencl_handle().context(ctx.opencl_context());
489  elements_.opencl_handle().context(ctx.opencl_context());
490  }
491 #endif
492  if (rows > 0)
493  {
495  }
496  }
497 
498  explicit compressed_matrix(viennacl::context ctx) : rows_(0), cols_(0), nonzeros_(0)
499  {
500  row_buffer_.switch_active_handle_id(ctx.memory_type());
501  col_buffer_.switch_active_handle_id(ctx.memory_type());
502  elements_.switch_active_handle_id(ctx.memory_type());
503 
504 #ifdef VIENNACL_WITH_OPENCL
505  if (ctx.memory_type() == OPENCL_MEMORY)
506  {
507  row_buffer_.opencl_handle().context(ctx.opencl_context());
508  col_buffer_.opencl_handle().context(ctx.opencl_context());
509  elements_.opencl_handle().context(ctx.opencl_context());
510  }
511 #endif
512  }
513 
514 
515 #ifdef VIENNACL_WITH_OPENCL
516  explicit compressed_matrix(cl_mem mem_row_buffer, cl_mem mem_col_buffer, cl_mem mem_elements,
517  vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros) :
518  rows_(rows), cols_(cols), nonzeros_(nonzeros)
519  {
521  row_buffer_.opencl_handle() = mem_row_buffer;
522  row_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
523  row_buffer_.raw_size(sizeof(cl_uint) * (rows + 1));
524 
526  col_buffer_.opencl_handle() = mem_col_buffer;
527  col_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
528  col_buffer_.raw_size(sizeof(cl_uint) * nonzeros);
529 
531  elements_.opencl_handle() = mem_elements;
532  elements_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
533  elements_.raw_size(sizeof(SCALARTYPE) * nonzeros);
534  }
535 #endif
536 
537 
540  {
541  assert( (rows_ == 0 || rows_ == other.size1()) && bool("Size mismatch") );
542  assert( (cols_ == 0 || cols_ == other.size2()) && bool("Size mismatch") );
543 
544  rows_ = other.size1();
545  cols_ = other.size2();
546  nonzeros_ = other.nnz();
547 
548  viennacl::backend::typesafe_memory_copy<unsigned int>(other.row_buffer_, row_buffer_);
549  viennacl::backend::typesafe_memory_copy<unsigned int>(other.col_buffer_, col_buffer_);
550  viennacl::backend::typesafe_memory_copy<SCALARTYPE>(other.elements_, elements_);
551 
552  return *this;
553  }
554 
555 
565  void set(const void * row_jumper,
566  const void * col_buffer,
567  const SCALARTYPE * elements,
568  vcl_size_t rows,
569  vcl_size_t cols,
570  vcl_size_t nonzeros)
571  {
572  assert( (rows > 0) && bool("Error in compressed_matrix::set(): Number of rows must be larger than zero!"));
573  assert( (cols > 0) && bool("Error in compressed_matrix::set(): Number of columns must be larger than zero!"));
574  assert( (nonzeros > 0) && bool("Error in compressed_matrix::set(): Number of nonzeros must be larger than zero!"));
575  //std::cout << "Setting memory: " << cols + 1 << ", " << nonzeros << std::endl;
576 
577  //row_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
579 
580  //col_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
582 
583  //elements_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
584  viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * nonzeros, viennacl::traits::context(elements_), elements);
585 
586  nonzeros_ = nonzeros;
587  rows_ = rows;
588  cols_ = cols;
589  }
590 
592  void reserve(vcl_size_t new_nonzeros)
593  {
594  if (new_nonzeros > nonzeros_)
595  {
596  handle_type col_buffer_old;
597  handle_type elements_old;
598  viennacl::backend::memory_shallow_copy(col_buffer_, col_buffer_old);
599  viennacl::backend::memory_shallow_copy(elements_, elements_old);
600 
602  viennacl::backend::memory_create(col_buffer_, size_deducer.element_size() * new_nonzeros, viennacl::traits::context(col_buffer_));
603  viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * new_nonzeros, viennacl::traits::context(elements_));
604 
605  viennacl::backend::memory_copy(col_buffer_old, col_buffer_, 0, 0, size_deducer.element_size() * nonzeros_);
606  viennacl::backend::memory_copy(elements_old, elements_, 0, 0, sizeof(SCALARTYPE)* nonzeros_);
607 
608  nonzeros_ = new_nonzeros;
609  }
610  }
611 
618  void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve = true)
619  {
620  assert(new_size1 > 0 && new_size2 > 0 && bool("Cannot resize to zero size!"));
621 
622  if (new_size1 != rows_ || new_size2 != cols_)
623  {
624  std::vector<std::map<unsigned int, SCALARTYPE> > stl_sparse_matrix;
625  if (rows_ > 0)
626  {
627  if (preserve)
628  {
629  stl_sparse_matrix.resize(rows_);
630  viennacl::copy(*this, stl_sparse_matrix);
631  } else
632  stl_sparse_matrix[0][0] = 0;
633  } else {
634  stl_sparse_matrix.resize(new_size1);
635  stl_sparse_matrix[0][0] = 0; //enforces nonzero array sizes if matrix was initially empty
636  }
637 
638  stl_sparse_matrix.resize(new_size1);
639 
640  //discard entries with column index larger than new_size2
641  if (new_size2 < cols_ && rows_ > 0)
642  {
643  for (vcl_size_t i=0; i<stl_sparse_matrix.size(); ++i)
644  {
645  std::list<unsigned int> to_delete;
646  for (typename std::map<unsigned int, SCALARTYPE>::iterator it = stl_sparse_matrix[i].begin();
647  it != stl_sparse_matrix[i].end();
648  ++it)
649  {
650  if (it->first >= new_size2)
651  to_delete.push_back(it->first);
652  }
653 
654  for (std::list<unsigned int>::iterator it = to_delete.begin(); it != to_delete.end(); ++it)
655  stl_sparse_matrix[i].erase(*it);
656  }
657  }
658 
659  viennacl::copy(stl_sparse_matrix, *this);
660 
661  rows_ = new_size1;
662  cols_ = new_size2;
663  }
664  }
665 
668  {
669  assert( (i < rows_) && (j < cols_) && bool("compressed_matrix access out of bounds!"));
670 
671  vcl_size_t index = element_index(i, j);
672 
673  // check for element in sparsity pattern
674  if (index < nonzeros_)
675  return entry_proxy<SCALARTYPE>(index, elements_);
676 
677  // Element not found. Copying required. Very slow, but direct entry manipulation is painful anyway...
678  std::vector< std::map<unsigned int, SCALARTYPE> > cpu_backup(rows_);
679  tools::sparse_matrix_adapter<SCALARTYPE> adapted_cpu_backup(cpu_backup, rows_, cols_);
680  viennacl::copy(*this, adapted_cpu_backup);
681  cpu_backup[i][static_cast<unsigned int>(j)] = 0.0;
682  viennacl::copy(adapted_cpu_backup, *this);
683 
684  index = element_index(i, j);
685 
686  assert(index < nonzeros_);
687 
688  return entry_proxy<SCALARTYPE>(index, elements_);
689  }
690 
692  const vcl_size_t & size1() const { return rows_; }
694  const vcl_size_t & size2() const { return cols_; }
696  const vcl_size_t & nnz() const { return nonzeros_; }
697 
699  const handle_type & handle1() const { return row_buffer_; }
701  const handle_type & handle2() const { return col_buffer_; }
703  const handle_type & handle() const { return elements_; }
704 
706  handle_type & handle1() { return row_buffer_; }
708  handle_type & handle2() { return col_buffer_; }
710  handle_type & handle() { return elements_; }
711 
713  {
714  viennacl::backend::switch_memory_context<unsigned int>(row_buffer_, new_ctx);
715  viennacl::backend::switch_memory_context<unsigned int>(col_buffer_, new_ctx);
716  viennacl::backend::switch_memory_context<SCALARTYPE>(elements_, new_ctx);
717  }
718 
720  {
721  return row_buffer_.get_active_handle_id();
722  }
723 
724  private:
725 
726  vcl_size_t element_index(vcl_size_t i, vcl_size_t j)
727  {
728  //read row indices
729  viennacl::backend::typesafe_host_array<unsigned int> row_indices(row_buffer_, 2);
730  viennacl::backend::memory_read(row_buffer_, row_indices.element_size()*i, row_indices.element_size()*2, row_indices.get());
731 
732  //get column indices for row i:
733  viennacl::backend::typesafe_host_array<unsigned int> col_indices(col_buffer_, row_indices[1] - row_indices[0]);
734  viennacl::backend::memory_read(col_buffer_, col_indices.element_size()*row_indices[0], row_indices.element_size()*col_indices.size(), col_indices.get());
735 
736  //get entries for row i:
737  viennacl::backend::typesafe_host_array<SCALARTYPE> row_entries(elements_, row_indices[1] - row_indices[0]);
738  viennacl::backend::memory_read(elements_, sizeof(SCALARTYPE)*row_indices[0], sizeof(SCALARTYPE)*row_entries.size(), row_entries.get());
739 
740  for (vcl_size_t k=0; k<col_indices.size(); ++k)
741  {
742  if (col_indices[k] == j)
743  return row_indices[0] + k;
744  }
745 
746  // if not found, return index past the end of the matrix (cf. matrix.end() in the spirit of the STL)
747  return nonzeros_;
748  }
749 
750  // /** @brief Copy constructor is by now not available. */
751  //compressed_matrix(compressed_matrix const &);
752 
753 
754  vcl_size_t rows_;
755  vcl_size_t cols_;
756  vcl_size_t nonzeros_;
757  handle_type row_buffer_;
758  handle_type col_buffer_;
759  handle_type elements_;
760  };
761 
762 
763 
764  //
765  // Specify available operations:
766  //
767 
770  namespace linalg
771  {
772  namespace detail
773  {
774  // x = A * y
775  template <typename T, unsigned int A>
776  struct op_executor<vector_base<T>, op_assign, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
777  {
778  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
779  {
780  // check for the special case x = A * x
781  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
782  {
783  viennacl::vector<T> temp(lhs);
784  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
785  lhs = temp;
786  }
787  else
788  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
789  }
790  };
791 
792  template <typename T, unsigned int A>
793  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
794  {
795  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
796  {
797  viennacl::vector<T> temp(lhs);
798  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
799  lhs += temp;
800  }
801  };
802 
803  template <typename T, unsigned int A>
804  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
805  {
806  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
807  {
808  viennacl::vector<T> temp(lhs);
809  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
810  lhs -= temp;
811  }
812  };
813 
814 
815  // x = A * vec_op
816  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
817  struct op_executor<vector_base<T>, op_assign, vector_expression<const compressed_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
818  {
819  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
820  {
821  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
822  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
823  }
824  };
825 
826  // x = A * vec_op
827  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
828  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const compressed_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> >
829  {
830  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
831  {
832  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
833  viennacl::vector<T> temp_result(lhs);
834  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
835  lhs += temp_result;
836  }
837  };
838 
839  // x = A * vec_op
840  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
841  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const compressed_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
842  {
843  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
844  {
845  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
846  viennacl::vector<T> temp_result(lhs);
847  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
848  lhs -= temp_result;
849  }
850  };
851 
852  } // namespace detail
853  } // namespace linalg
854 
856 }
857 
858 #endif
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: compressed_matrix.hpp:432
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
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:172
void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve=true)
Resize the matrix.
Definition: compressed_matrix.hpp:618
handle_type & handle()
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:710
void switch_memory_context(viennacl::context new_ctx)
Definition: compressed_matrix.hpp:712
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
compressed_matrix(vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros=0, viennacl::context ctx=viennacl::context())
Construction of a compressed matrix with the supplied number of rows and columns. If the number of no...
Definition: compressed_matrix.hpp:445
A proxy class for entries in a vector.
vcl_size_t size_type
Definition: compressed_matrix.hpp:433
entry_proxy< SCALARTYPE > operator()(vcl_size_t i, vcl_size_t j)
Returns a reference to the (i,j)-th entry of the sparse matrix. If (i,j) does not exist (zero)...
Definition: compressed_matrix.hpp:667
void copy_impl(const CPU_MATRIX &cpu_matrix, compressed_compressed_matrix< SCALARTYPE > &gpu_matrix, vcl_size_t nonzero_rows, vcl_size_t nonzeros)
Definition: compressed_compressed_matrix.hpp:41
This file provides the forward declarations for the main types used within ViennaCL.
handle_type & handle2()
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:708
compressed_matrix()
Default construction of a compressed matrix. No memory is allocated.
Definition: compressed_matrix.hpp:436
handle_type & handle1()
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:706
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
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:701
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
memory_types
Definition: forwards.h:476
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:694
compressed_matrix & operator=(compressed_matrix const &other)
Assignment a compressed matrix from possibly another memory domain.
Definition: compressed_matrix.hpp:539
viennacl::backend::mem_handle handle_type
Definition: compressed_matrix.hpp:431
void set(vcl_size_t index, U value)
Definition: util.hpp:145
const vcl_size_t & nnz() const
Returns the number of nonzero entries.
Definition: compressed_matrix.hpp:696
vcl_size_t element_size(memory_types)
Definition: memory.hpp:299
viennacl::memory_types memory_context() const
Definition: compressed_matrix.hpp:719
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
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:91
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
Definition: forwards.h:480
Implementations of operations using sparse matrices.
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
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:692
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
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 ...
void reserve(vcl_size_t new_nonzeros)
Allocate memory for the supplied number of nonzeros in the matrix. Old values are preserved...
Definition: compressed_matrix.hpp:592
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
compressed_matrix(viennacl::context ctx)
Definition: compressed_matrix.hpp:498
compressed_matrix(vcl_size_t rows, vcl_size_t cols, viennacl::context ctx)
Construction of a compressed matrix with the supplied number of rows and columns. If the number of no...
Definition: compressed_matrix.hpp:477
Adapts a non-const sparse matrix type made up from std::vector<std::map<SizeType, SCALARTYPE> > to ba...
Definition: adapter.hpp:345
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:428
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 set(const void *row_jumper, const void *col_buffer, const SCALARTYPE *elements, vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros)
Sets the row, column and value arrays of the compressed matrix.
Definition: compressed_matrix.hpp:565
A proxy class for a single element of a vector or matrix. This proxy should not be noticed by end-use...
Definition: forwards.h:178
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
vcl_size_t element_size() const
Definition: util.hpp:159
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:703
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:699