ViennaCL - The Vienna Computing Library  1.7.1
Free open-source GPU-accelerated linear algebra and solver library.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2016, 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 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 
31 #include "viennacl/linalg/sparse_matrix_operations.hpp"
32 
33 #include "viennacl/tools/tools.hpp"
34 #include "viennacl/tools/entry_proxy.hpp"
35 
36 #ifdef VIENNACL_WITH_UBLAS
37 #include <boost/numeric/ublas/matrix_sparse.hpp>
38 #endif
39 
40 namespace viennacl
41 {
42 namespace detail
43 {
44 
49  template<typename CPUMatrixT, typename NumericT, unsigned int AlignmentV>
50  void copy_impl(const CPUMatrixT & cpu_matrix,
51  compressed_matrix<NumericT, AlignmentV> & gpu_matrix,
52  vcl_size_t nonzeros)
53  {
54  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
55  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
56 
57  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1);
58  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), nonzeros);
59  std::vector<NumericT> elements(nonzeros);
60 
61  vcl_size_t row_index = 0;
62  vcl_size_t data_index = 0;
63 
64  for (typename CPUMatrixT::const_iterator1 row_it = cpu_matrix.begin1();
65  row_it != cpu_matrix.end1();
66  ++row_it)
67  {
68  row_buffer.set(row_index, data_index);
69  ++row_index;
70 
71  for (typename CPUMatrixT::const_iterator2 col_it = row_it.begin();
72  col_it != row_it.end();
73  ++col_it)
74  {
75  col_buffer.set(data_index, col_it.index2());
76  elements[data_index] = *col_it;
77  ++data_index;
78  }
79  data_index = viennacl::tools::align_to_multiple<vcl_size_t>(data_index, AlignmentV); //take care of alignment
80  }
81  row_buffer.set(row_index, data_index);
82 
83  gpu_matrix.set(row_buffer.get(),
84  col_buffer.get(),
85  &elements[0],
86  cpu_matrix.size1(),
87  cpu_matrix.size2(),
88  nonzeros);
89  }
90 }
91 
92 //
93 // host to device:
94 //
95 
96 //provide copy-operation:
111 template<typename CPUMatrixT, typename NumericT, unsigned int AlignmentV>
112 void copy(const CPUMatrixT & cpu_matrix,
113  compressed_matrix<NumericT, AlignmentV> & gpu_matrix )
114 {
115  if ( cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0 )
116  {
117  //determine nonzeros:
118  vcl_size_t num_entries = 0;
119  for (typename CPUMatrixT::const_iterator1 row_it = cpu_matrix.begin1();
120  row_it != cpu_matrix.end1();
121  ++row_it)
122  {
123  vcl_size_t entries_per_row = 0;
124  for (typename CPUMatrixT::const_iterator2 col_it = row_it.begin();
125  col_it != row_it.end();
126  ++col_it)
127  {
128  ++entries_per_row;
129  }
130  num_entries += viennacl::tools::align_to_multiple<vcl_size_t>(entries_per_row, AlignmentV);
131  }
132 
133  if (num_entries == 0) //we copy an empty matrix
134  num_entries = 1;
135 
136  //set up matrix entries:
137  viennacl::detail::copy_impl(cpu_matrix, gpu_matrix, num_entries);
138  }
139 }
140 
141 
142 //adapted for std::vector< std::map < > > argument:
148 template<typename SizeT, typename NumericT, unsigned int AlignmentV>
149 void copy(const std::vector< std::map<SizeT, NumericT> > & cpu_matrix,
150  compressed_matrix<NumericT, AlignmentV> & gpu_matrix )
151 {
152  vcl_size_t nonzeros = 0;
153  vcl_size_t max_col = 0;
154  for (vcl_size_t i=0; i<cpu_matrix.size(); ++i)
155  {
156  if (cpu_matrix[i].size() > 0)
157  nonzeros += ((cpu_matrix[i].size() - 1) / AlignmentV + 1) * AlignmentV;
158  if (cpu_matrix[i].size() > 0)
159  max_col = std::max<vcl_size_t>(max_col, (cpu_matrix[i].rbegin())->first);
160  }
161 
162  viennacl::detail::copy_impl(tools::const_sparse_matrix_adapter<NumericT, SizeT>(cpu_matrix, cpu_matrix.size(), max_col + 1),
163  gpu_matrix,
164  nonzeros);
165 }
166 
167 #ifdef VIENNACL_WITH_UBLAS
168 
172 template<typename ScalarType, typename F, vcl_size_t IB, typename IA, typename TA>
173 void copy(const boost::numeric::ublas::compressed_matrix<ScalarType, F, IB, IA, TA> & ublas_matrix,
174  viennacl::compressed_matrix<ScalarType, 1> & gpu_matrix)
175 {
176  assert( (gpu_matrix.size1() == 0 || viennacl::traits::size1(ublas_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
177  assert( (gpu_matrix.size2() == 0 || viennacl::traits::size2(ublas_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
178 
179  //we just need to copy the CSR arrays:
180  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), ublas_matrix.size1() + 1);
181  for (vcl_size_t i=0; i<=ublas_matrix.size1(); ++i)
182  row_buffer.set(i, ublas_matrix.index1_data()[i]);
183 
184  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), ublas_matrix.nnz());
185  for (vcl_size_t i=0; i<ublas_matrix.nnz(); ++i)
186  col_buffer.set(i, ublas_matrix.index2_data()[i]);
187 
188  gpu_matrix.set(row_buffer.get(),
189  col_buffer.get(),
190  &(ublas_matrix.value_data()[0]),
191  ublas_matrix.size1(),
192  ublas_matrix.size2(),
193  ublas_matrix.nnz());
194 
195 }
196 #endif
197 
198 #ifdef VIENNACL_WITH_ARMADILLO
199 
204 template<typename NumericT, unsigned int AlignmentV>
205 void copy(arma::SpMat<NumericT> const & arma_matrix,
206  viennacl::compressed_matrix<NumericT, AlignmentV> & vcl_matrix)
207 {
208  assert( (vcl_matrix.size1() == 0 || static_cast<vcl_size_t>(arma_matrix.n_rows) == vcl_matrix.size1()) && bool("Size mismatch") );
209  assert( (vcl_matrix.size2() == 0 || static_cast<vcl_size_t>(arma_matrix.n_cols) == vcl_matrix.size2()) && bool("Size mismatch") );
210 
211  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(vcl_matrix.handle1(), arma_matrix.n_rows + 1);
212  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(vcl_matrix.handle2(), arma_matrix.n_nonzero);
213  viennacl::backend::typesafe_host_array<NumericT > value_buffer(vcl_matrix.handle(), arma_matrix.n_nonzero);
214 
215  // Step 1: Count number of nonzeros in each row
216  for (vcl_size_t col=0; col < static_cast<vcl_size_t>(arma_matrix.n_cols); ++col)
217  {
218  vcl_size_t col_begin = static_cast<vcl_size_t>(arma_matrix.col_ptrs[col]);
219  vcl_size_t col_end = static_cast<vcl_size_t>(arma_matrix.col_ptrs[col+1]);
220  for (vcl_size_t i = col_begin; i < col_end; ++i)
221  {
222  unsigned int row = arma_matrix.row_indices[i];
223  row_buffer.set(row, row_buffer[row] + 1);
224  }
225  }
226 
227  // Step 2: Exclusive scan on row_buffer to obtain offsets
228  unsigned int offset = 0;
229  for (vcl_size_t i=0; i<row_buffer.size(); ++i)
230  {
231  unsigned int tmp = row_buffer[i];
232  row_buffer.set(i, offset);
233  offset += tmp;
234  }
235 
236  // Step 3: Fill data
237  std::vector<unsigned int> row_offsets(arma_matrix.n_rows);
238  for (vcl_size_t col=0; col < static_cast<vcl_size_t>(arma_matrix.n_cols); ++col)
239  {
240  vcl_size_t col_begin = static_cast<vcl_size_t>(arma_matrix.col_ptrs[col]);
241  vcl_size_t col_end = static_cast<vcl_size_t>(arma_matrix.col_ptrs[col+1]);
242  for (vcl_size_t i = col_begin; i < col_end; ++i)
243  {
244  unsigned int row = arma_matrix.row_indices[i];
245  col_buffer.set(row_buffer[row] + row_offsets[row], col);
246  value_buffer.set(row_buffer[row] + row_offsets[row], arma_matrix.values[i]);
247  row_offsets[row] += 1;
248  }
249  }
250 
251  vcl_matrix.set(row_buffer.get(), col_buffer.get(), reinterpret_cast<NumericT*>(value_buffer.get()),
252  arma_matrix.n_rows, arma_matrix.n_cols, arma_matrix.n_nonzero);
253 }
254 #endif
255 
256 #ifdef VIENNACL_WITH_EIGEN
257 
261 template<typename NumericT, int flags, unsigned int AlignmentV>
262 void copy(const Eigen::SparseMatrix<NumericT, flags> & eigen_matrix,
263  compressed_matrix<NumericT, AlignmentV> & gpu_matrix)
264 {
265  assert( (gpu_matrix.size1() == 0 || static_cast<vcl_size_t>(eigen_matrix.rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
266  assert( (gpu_matrix.size2() == 0 || static_cast<vcl_size_t>(eigen_matrix.cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
267 
268  std::vector< std::map<unsigned int, NumericT> > stl_matrix(eigen_matrix.rows());
269 
270  for (int k=0; k < eigen_matrix.outerSize(); ++k)
271  for (typename Eigen::SparseMatrix<NumericT, flags>::InnerIterator it(eigen_matrix, k); it; ++it)
272  stl_matrix[it.row()][it.col()] = it.value();
273 
274  copy(tools::const_sparse_matrix_adapter<NumericT>(stl_matrix, eigen_matrix.rows(), eigen_matrix.cols()), gpu_matrix);
275 }
276 #endif
277 
278 
279 #ifdef VIENNACL_WITH_MTL4
280 
284 template<typename NumericT, unsigned int AlignmentV>
285 void copy(const mtl::compressed2D<NumericT> & cpu_matrix,
286  compressed_matrix<NumericT, AlignmentV> & gpu_matrix)
287 {
288  assert( (gpu_matrix.size1() == 0 || static_cast<vcl_size_t>(cpu_matrix.num_rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
289  assert( (gpu_matrix.size2() == 0 || static_cast<vcl_size_t>(cpu_matrix.num_cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
290 
291  typedef mtl::compressed2D<NumericT> MatrixType;
292 
293  std::vector< std::map<unsigned int, NumericT> > stl_matrix(cpu_matrix.num_rows());
294 
295  using mtl::traits::range_generator;
296  using mtl::traits::range::min;
297 
298  // Choose between row and column traversal
299  typedef typename min<range_generator<mtl::tag::row, MatrixType>,
300  range_generator<mtl::tag::col, MatrixType> >::type range_type;
301  range_type my_range;
302 
303  // Type of outer cursor
304  typedef typename range_type::type c_type;
305  // Type of inner cursor
306  typedef typename mtl::traits::range_generator<mtl::tag::nz, c_type>::type ic_type;
307 
308  // Define the property maps
309  typename mtl::traits::row<MatrixType>::type row(cpu_matrix);
310  typename mtl::traits::col<MatrixType>::type col(cpu_matrix);
311  typename mtl::traits::const_value<MatrixType>::type value(cpu_matrix);
312 
313  // Now iterate over the matrix
314  for (c_type cursor(my_range.begin(cpu_matrix)), cend(my_range.end(cpu_matrix)); cursor != cend; ++cursor)
315  for (ic_type icursor(mtl::begin<mtl::tag::nz>(cursor)), icend(mtl::end<mtl::tag::nz>(cursor)); icursor != icend; ++icursor)
316  stl_matrix[row(*icursor)][col(*icursor)] = value(*icursor);
317 
318  copy(tools::const_sparse_matrix_adapter<NumericT>(stl_matrix, cpu_matrix.num_rows(), cpu_matrix.num_cols()), gpu_matrix);
319 }
320 #endif
321 
322 
323 
324 
325 
326 
327 
328 //
329 // device to host:
330 //
340 template<typename CPUMatrixT, typename NumericT, unsigned int AlignmentV>
341 void copy(const compressed_matrix<NumericT, AlignmentV> & gpu_matrix,
342  CPUMatrixT & cpu_matrix )
343 {
344  assert( (viennacl::traits::size1(cpu_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
345  assert( (viennacl::traits::size2(cpu_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
346 
347  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
348  {
349  //get raw data from memory:
350  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1);
351  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
352  std::vector<NumericT> elements(gpu_matrix.nnz());
353 
354  //std::cout << "GPU->CPU, nonzeros: " << gpu_matrix.nnz() << std::endl;
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(NumericT)* gpu_matrix.nnz(), &(elements[0]));
359 
360  //fill the cpu_matrix:
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  if (col_buffer[data_index] >= gpu_matrix.size2())
367  {
368  std::cerr << "ViennaCL encountered invalid data at colbuffer[" << data_index << "]: " << col_buffer[data_index] << std::endl;
369  return;
370  }
371 
372  if (std::fabs(elements[data_index]) > static_cast<NumericT>(0))
373  cpu_matrix(row-1, static_cast<vcl_size_t>(col_buffer[data_index])) = elements[data_index];
374  ++data_index;
375  }
376  }
377  }
378 }
379 
380 
386 template<typename NumericT, unsigned int AlignmentV>
387 void copy(const compressed_matrix<NumericT, AlignmentV> & gpu_matrix,
388  std::vector< std::map<unsigned int, NumericT> > & cpu_matrix)
389 {
390  assert( (cpu_matrix.size() == gpu_matrix.size1()) && bool("Size mismatch") );
391 
392  tools::sparse_matrix_adapter<NumericT> temp(cpu_matrix, gpu_matrix.size1(), gpu_matrix.size2());
393  copy(gpu_matrix, temp);
394 }
395 
396 #ifdef VIENNACL_WITH_UBLAS
397 
401 template<typename ScalarType, unsigned int AlignmentV, typename F, vcl_size_t IB, typename IA, typename TA>
402 void copy(viennacl::compressed_matrix<ScalarType, AlignmentV> const & gpu_matrix,
403  boost::numeric::ublas::compressed_matrix<ScalarType> & ublas_matrix)
404 {
405  assert( (viennacl::traits::size1(ublas_matrix) == gpu_matrix.size1()) && bool("Size mismatch") );
406  assert( (viennacl::traits::size2(ublas_matrix) == gpu_matrix.size2()) && bool("Size mismatch") );
407 
408  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
409  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
410 
411  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
412  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
413 
414  ublas_matrix.clear();
415  ublas_matrix.reserve(gpu_matrix.nnz());
416 
417  ublas_matrix.set_filled(gpu_matrix.size1() + 1, gpu_matrix.nnz());
418 
419  for (vcl_size_t i=0; i<ublas_matrix.size1() + 1; ++i)
420  ublas_matrix.index1_data()[i] = row_buffer[i];
421 
422  for (vcl_size_t i=0; i<ublas_matrix.nnz(); ++i)
423  ublas_matrix.index2_data()[i] = col_buffer[i];
424 
425  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(ScalarType) * gpu_matrix.nnz(), &(ublas_matrix.value_data()[0]));
426 
427 }
428 #endif
429 
430 #ifdef VIENNACL_WITH_ARMADILLO
431 
436 template<typename NumericT, unsigned int AlignmentV>
437 void copy(viennacl::compressed_matrix<NumericT, AlignmentV> & vcl_matrix,
438  arma::SpMat<NumericT> & arma_matrix)
439 {
440  assert( (static_cast<vcl_size_t>(arma_matrix.n_rows) == vcl_matrix.size1()) && bool("Size mismatch") );
441  assert( (static_cast<vcl_size_t>(arma_matrix.n_cols) == vcl_matrix.size2()) && bool("Size mismatch") );
442 
443  if ( vcl_matrix.size1() > 0 && vcl_matrix.size2() > 0 )
444  {
445  //get raw data from memory:
446  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(vcl_matrix.handle1(), vcl_matrix.size1() + 1);
447  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(vcl_matrix.handle2(), vcl_matrix.nnz());
448  viennacl::backend::typesafe_host_array<NumericT> elements (vcl_matrix.handle(), vcl_matrix.nnz());
449 
450  viennacl::backend::memory_read(vcl_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
451  viennacl::backend::memory_read(vcl_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
452  viennacl::backend::memory_read(vcl_matrix.handle(), 0, elements.raw_size(), elements.get());
453 
454  arma_matrix.zeros();
455  vcl_size_t data_index = 0;
456  for (vcl_size_t row = 1; row <= vcl_matrix.size1(); ++row)
457  {
458  while (data_index < row_buffer[row])
459  {
460  assert(col_buffer[data_index] < vcl_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer"));
461  if (elements[data_index] != static_cast<NumericT>(0.0))
462  arma_matrix(row-1, col_buffer[data_index]) = elements[data_index];
463  ++data_index;
464  }
465  }
466  }
467 }
468 #endif
469 
470 #ifdef VIENNACL_WITH_EIGEN
471 
472 template<typename NumericT, int flags, unsigned int AlignmentV>
473 void copy(compressed_matrix<NumericT, AlignmentV> & gpu_matrix,
474  Eigen::SparseMatrix<NumericT, flags> & eigen_matrix)
475 {
476  assert( (static_cast<vcl_size_t>(eigen_matrix.rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
477  assert( (static_cast<vcl_size_t>(eigen_matrix.cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
478 
479  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
480  {
481  //get raw data from memory:
482  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
483  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
484  std::vector<NumericT> elements(gpu_matrix.nnz());
485 
486  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
487  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
488  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(NumericT)* gpu_matrix.nnz(), &(elements[0]));
489 
490  eigen_matrix.setZero();
491  vcl_size_t data_index = 0;
492  for (vcl_size_t row = 1; row <= gpu_matrix.size1(); ++row)
493  {
494  while (data_index < row_buffer[row])
495  {
496  assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer"));
497  if (elements[data_index] != static_cast<NumericT>(0.0))
498  eigen_matrix.insert(row-1, col_buffer[data_index]) = elements[data_index];
499  ++data_index;
500  }
501  }
502  }
503 }
504 #endif
505 
506 
507 
508 #ifdef VIENNACL_WITH_MTL4
509 
510 template<typename NumericT, unsigned int AlignmentV>
511 void copy(compressed_matrix<NumericT, AlignmentV> & gpu_matrix,
512  mtl::compressed2D<NumericT> & mtl4_matrix)
513 {
514  assert( (static_cast<vcl_size_t>(mtl4_matrix.num_rows()) == gpu_matrix.size1()) && bool("Size mismatch") );
515  assert( (static_cast<vcl_size_t>(mtl4_matrix.num_cols()) == gpu_matrix.size2()) && bool("Size mismatch") );
516 
517  if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
518  {
519 
520  //get raw data from memory:
521  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1);
522  viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz());
523  std::vector<NumericT> elements(gpu_matrix.nnz());
524 
525  viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get());
526  viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get());
527  viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(NumericT)* gpu_matrix.nnz(), &(elements[0]));
528 
529  //set_to_zero(mtl4_matrix);
530  //mtl4_matrix.change_dim(gpu_matrix.size1(), gpu_matrix.size2());
531 
532  mtl::matrix::inserter< mtl::compressed2D<NumericT> > ins(mtl4_matrix);
533  vcl_size_t data_index = 0;
534  for (vcl_size_t row = 1; row <= gpu_matrix.size1(); ++row)
535  {
536  while (data_index < row_buffer[row])
537  {
538  assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer"));
539  if (elements[data_index] != static_cast<NumericT>(0.0))
540  ins(row-1, col_buffer[data_index]) << typename mtl::Collection< mtl::compressed2D<NumericT> >::value_type(elements[data_index]);
541  ++data_index;
542  }
543  }
544  }
545 }
546 #endif
547 
548 
549 
550 
551 
553 
558 template<class NumericT, unsigned int AlignmentV /* see VCLForwards.h */>
559 class compressed_matrix
560 {
561 public:
562  typedef viennacl::backend::mem_handle handle_type;
563  typedef scalar<typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT<NumericT>::ResultType> value_type;
564  typedef vcl_size_t size_type;
565 
567  compressed_matrix() : rows_(0), cols_(0), nonzeros_(0), row_block_num_(0) {}
568 
576  explicit compressed_matrix(vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros = 0, viennacl::context ctx = viennacl::context())
577  : rows_(rows), cols_(cols), nonzeros_(nonzeros), row_block_num_(0)
578  {
579  row_buffer_.switch_active_handle_id(ctx.memory_type());
580  col_buffer_.switch_active_handle_id(ctx.memory_type());
581  elements_.switch_active_handle_id(ctx.memory_type());
582  row_blocks_.switch_active_handle_id(ctx.memory_type());
583 
584 #ifdef VIENNACL_WITH_OPENCL
585  if (ctx.memory_type() == OPENCL_MEMORY)
586  {
587  row_buffer_.opencl_handle().context(ctx.opencl_context());
588  col_buffer_.opencl_handle().context(ctx.opencl_context());
589  elements_.opencl_handle().context(ctx.opencl_context());
590  row_blocks_.opencl_handle().context(ctx.opencl_context());
591  }
592 #endif
593  if (rows > 0)
594  {
595  viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * (rows + 1), ctx);
596  viennacl::vector_base<unsigned int> init_temporary(row_buffer_, size_type(rows+1), 0, 1);
597  init_temporary = viennacl::zero_vector<unsigned int>(size_type(rows+1), ctx);
598  }
599  if (nonzeros > 0)
600  {
601  viennacl::backend::memory_create(col_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * nonzeros, ctx);
602  viennacl::backend::memory_create(elements_, sizeof(NumericT) * nonzeros, ctx);
603  }
604  }
605 
612  explicit compressed_matrix(vcl_size_t rows, vcl_size_t cols, viennacl::context ctx)
613  : rows_(rows), cols_(cols), nonzeros_(0), row_block_num_(0)
614  {
615  row_buffer_.switch_active_handle_id(ctx.memory_type());
616  col_buffer_.switch_active_handle_id(ctx.memory_type());
617  elements_.switch_active_handle_id(ctx.memory_type());
618  row_blocks_.switch_active_handle_id(ctx.memory_type());
619 
620 #ifdef VIENNACL_WITH_OPENCL
621  if (ctx.memory_type() == OPENCL_MEMORY)
622  {
623  row_buffer_.opencl_handle().context(ctx.opencl_context());
624  col_buffer_.opencl_handle().context(ctx.opencl_context());
625  elements_.opencl_handle().context(ctx.opencl_context());
626  row_blocks_.opencl_handle().context(ctx.opencl_context());
627  }
628 #endif
629  if (rows > 0)
630  {
631  viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * (rows + 1), ctx);
632  viennacl::vector_base<unsigned int> init_temporary(row_buffer_, size_type(rows+1), 0, 1);
633  init_temporary = viennacl::zero_vector<unsigned int>(size_type(rows+1), ctx);
634  }
635  }
636 
641  explicit compressed_matrix(viennacl::context ctx) : rows_(0), cols_(0), nonzeros_(0), row_block_num_(0)
642  {
643  row_buffer_.switch_active_handle_id(ctx.memory_type());
644  col_buffer_.switch_active_handle_id(ctx.memory_type());
645  elements_.switch_active_handle_id(ctx.memory_type());
646  row_blocks_.switch_active_handle_id(ctx.memory_type());
647 
648 #ifdef VIENNACL_WITH_OPENCL
649  if (ctx.memory_type() == OPENCL_MEMORY)
650  {
651  row_buffer_.opencl_handle().context(ctx.opencl_context());
652  col_buffer_.opencl_handle().context(ctx.opencl_context());
653  elements_.opencl_handle().context(ctx.opencl_context());
654  row_blocks_.opencl_handle().context(ctx.opencl_context());
655  }
656 #endif
657  }
658 
659 
660 #ifdef VIENNACL_WITH_OPENCL
661 
670  explicit compressed_matrix(cl_mem mem_row_buffer, cl_mem mem_col_buffer, cl_mem mem_elements,
671  vcl_size_t rows, vcl_size_t cols, vcl_size_t nonzeros) :
672  rows_(rows), cols_(cols), nonzeros_(nonzeros), row_block_num_(0)
673  {
674  row_buffer_.switch_active_handle_id(viennacl::OPENCL_MEMORY);
675  row_buffer_.opencl_handle() = mem_row_buffer;
676  row_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
677  row_buffer_.raw_size(sizeof(cl_uint) * (rows + 1));
678 
679  col_buffer_.switch_active_handle_id(viennacl::OPENCL_MEMORY);
680  col_buffer_.opencl_handle() = mem_col_buffer;
681  col_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
682  col_buffer_.raw_size(sizeof(cl_uint) * nonzeros);
683 
684  elements_.switch_active_handle_id(viennacl::OPENCL_MEMORY);
685  elements_.opencl_handle() = mem_elements;
686  elements_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed.
687  elements_.raw_size(sizeof(NumericT) * nonzeros);
688 
689  //generate block information for CSR-adaptive:
690  generate_row_block_information();
691  }
692 #endif
693 
695  compressed_matrix(matrix_expression<const compressed_matrix, const compressed_matrix, op_prod> const & proxy)
696  : rows_(0), cols_(0), nonzeros_(0), row_block_num_(0)
697  {
698  viennacl::context ctx = viennacl::traits::context(proxy.lhs());
699 
700  row_buffer_.switch_active_handle_id(ctx.memory_type());
701  col_buffer_.switch_active_handle_id(ctx.memory_type());
702  elements_.switch_active_handle_id(ctx.memory_type());
703  row_blocks_.switch_active_handle_id(ctx.memory_type());
704 
705 #ifdef VIENNACL_WITH_OPENCL
706  if (ctx.memory_type() == OPENCL_MEMORY)
707  {
708  row_buffer_.opencl_handle().context(ctx.opencl_context());
709  col_buffer_.opencl_handle().context(ctx.opencl_context());
710  elements_.opencl_handle().context(ctx.opencl_context());
711  row_blocks_.opencl_handle().context(ctx.opencl_context());
712  }
713 #endif
714 
715  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
716  generate_row_block_information();
717  }
718 
720  compressed_matrix & operator=(compressed_matrix const & other)
721  {
722  assert( (rows_ == 0 || rows_ == other.size1()) && bool("Size mismatch") );
723  assert( (cols_ == 0 || cols_ == other.size2()) && bool("Size mismatch") );
724 
725  rows_ = other.size1();
726  cols_ = other.size2();
727  nonzeros_ = other.nnz();
728  row_block_num_ = other.row_block_num_;
729 
730  viennacl::backend::typesafe_memory_copy<unsigned int>(other.row_buffer_, row_buffer_);
731  viennacl::backend::typesafe_memory_copy<unsigned int>(other.col_buffer_, col_buffer_);
732  viennacl::backend::typesafe_memory_copy<unsigned int>(other.row_blocks_, row_blocks_);
733  viennacl::backend::typesafe_memory_copy<NumericT>(other.elements_, elements_);
734 
735  return *this;
736  }
737 
739  compressed_matrix & operator=(matrix_expression<const compressed_matrix, const compressed_matrix, op_prod> const & proxy)
740  {
741  assert( (rows_ == 0 || rows_ == proxy.lhs().size1()) && bool("Size mismatch") );
742  assert( (cols_ == 0 || cols_ == proxy.rhs().size2()) && bool("Size mismatch") );
743 
744  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
745  generate_row_block_information();
746 
747  return *this;
748  }
749 
750 
764  void set(const void * row_jumper,
765  const void * col_buffer,
766  const NumericT * elements,
767  vcl_size_t rows,
768  vcl_size_t cols,
769  vcl_size_t nonzeros)
770  {
771  assert( (rows > 0) && bool("Error in compressed_matrix::set(): Number of rows must be larger than zero!"));
772  assert( (cols > 0) && bool("Error in compressed_matrix::set(): Number of columns must be larger than zero!"));
773  assert( (nonzeros > 0) && bool("Error in compressed_matrix::set(): Number of nonzeros must be larger than zero!"));
774  //std::cout << "Setting memory: " << cols + 1 << ", " << nonzeros << std::endl;
775 
776  //row_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
777  viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>(row_buffer_).element_size() * (rows + 1), viennacl::traits::context(row_buffer_), row_jumper);
778 
779  //col_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
780  viennacl::backend::memory_create(col_buffer_, viennacl::backend::typesafe_host_array<unsigned int>(col_buffer_).element_size() * nonzeros, viennacl::traits::context(col_buffer_), col_buffer);
781 
782  //elements_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY);
783  viennacl::backend::memory_create(elements_, sizeof(NumericT) * nonzeros, viennacl::traits::context(elements_), elements);
784 
785  nonzeros_ = nonzeros;
786  rows_ = rows;
787  cols_ = cols;
788 
789  //generate block information for CSR-adaptive:
790  generate_row_block_information();
791  }
792 
794  void reserve(vcl_size_t new_nonzeros, bool preserve = true)
795  {
796  if (new_nonzeros > nonzeros_)
797  {
798  if (preserve)
799  {
800  handle_type col_buffer_old;
801  handle_type elements_old;
802  viennacl::backend::memory_shallow_copy(col_buffer_, col_buffer_old);
803  viennacl::backend::memory_shallow_copy(elements_, elements_old);
804 
805  viennacl::backend::typesafe_host_array<unsigned int> size_deducer(col_buffer_);
806  viennacl::backend::memory_create(col_buffer_, size_deducer.element_size() * new_nonzeros, viennacl::traits::context(col_buffer_));
807  viennacl::backend::memory_create(elements_, sizeof(NumericT) * new_nonzeros, viennacl::traits::context(elements_));
808 
809  viennacl::backend::memory_copy(col_buffer_old, col_buffer_, 0, 0, size_deducer.element_size() * nonzeros_);
810  viennacl::backend::memory_copy(elements_old, elements_, 0, 0, sizeof(NumericT)* nonzeros_);
811  }
812  else
813  {
814  viennacl::backend::typesafe_host_array<unsigned int> size_deducer(col_buffer_);
815  viennacl::backend::memory_create(col_buffer_, size_deducer.element_size() * new_nonzeros, viennacl::traits::context(col_buffer_));
816  viennacl::backend::memory_create(elements_, sizeof(NumericT) * new_nonzeros, viennacl::traits::context(elements_));
817  }
818 
819  nonzeros_ = new_nonzeros;
820  }
821  }
822 
829  void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve = true)
830  {
831  assert(new_size1 > 0 && new_size2 > 0 && bool("Cannot resize to zero size!"));
832 
833  if (new_size1 != rows_ || new_size2 != cols_)
834  {
835  if (!preserve)
836  {
837  viennacl::backend::typesafe_host_array<unsigned int> host_row_buffer(row_buffer_, new_size1 + 1);
838  viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * (new_size1 + 1), viennacl::traits::context(row_buffer_), host_row_buffer.get());
839  // faster version without initializing memory:
840  //viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * (new_size1 + 1), viennacl::traits::context(row_buffer_));
841  nonzeros_ = 0;
842  }
843  else
844  {
845  std::vector<std::map<unsigned int, NumericT> > stl_sparse_matrix;
846  if (rows_ > 0)
847  {
848  stl_sparse_matrix.resize(rows_);
849  viennacl::copy(*this, stl_sparse_matrix);
850  } else {
851  stl_sparse_matrix.resize(new_size1);
852  stl_sparse_matrix[0][0] = 0; //enforces nonzero array sizes if matrix was initially empty
853  }
854 
855  stl_sparse_matrix.resize(new_size1);
856 
857  //discard entries with column index larger than new_size2
858  if (new_size2 < cols_ && rows_ > 0)
859  {
860  for (vcl_size_t i=0; i<stl_sparse_matrix.size(); ++i)
861  {
862  std::list<unsigned int> to_delete;
863  for (typename std::map<unsigned int, NumericT>::iterator it = stl_sparse_matrix[i].begin();
864  it != stl_sparse_matrix[i].end();
865  ++it)
866  {
867  if (it->first >= new_size2)
868  to_delete.push_back(it->first);
869  }
870 
871  for (std::list<unsigned int>::iterator it = to_delete.begin(); it != to_delete.end(); ++it)
872  stl_sparse_matrix[i].erase(*it);
873  }
874  }
875 
876  viennacl::tools::sparse_matrix_adapter<NumericT> adapted_matrix(stl_sparse_matrix, new_size1, new_size2);
877  rows_ = new_size1;
878  cols_ = new_size2;
879  viennacl::copy(adapted_matrix, *this);
880  }
881 
882  rows_ = new_size1;
883  cols_ = new_size2;
884  }
885  }
886 
888  void clear()
889  {
890  viennacl::backend::typesafe_host_array<unsigned int> host_row_buffer(row_buffer_, rows_ + 1);
891  viennacl::backend::typesafe_host_array<unsigned int> host_col_buffer(col_buffer_, 1);
892  std::vector<NumericT> host_elements(1);
893 
894  viennacl::backend::memory_create(row_buffer_, host_row_buffer.element_size() * (rows_ + 1), viennacl::traits::context(row_buffer_), host_row_buffer.get());
895  viennacl::backend::memory_create(col_buffer_, host_col_buffer.element_size() * 1, viennacl::traits::context(col_buffer_), host_col_buffer.get());
896  viennacl::backend::memory_create(elements_, sizeof(NumericT) * 1, viennacl::traits::context(elements_), &(host_elements[0]));
897 
898  nonzeros_ = 0;
899  }
900 
902  entry_proxy<NumericT> operator()(vcl_size_t i, vcl_size_t j)
903  {
904  assert( (i < rows_) && (j < cols_) && bool("compressed_matrix access out of bounds!"));
905 
906  vcl_size_t index = element_index(i, j);
907 
908  // check for element in sparsity pattern
909  if (index < nonzeros_)
910  return entry_proxy<NumericT>(index, elements_);
911 
912  // Element not found. Copying required. Very slow, but direct entry manipulation is painful anyway...
913  std::vector< std::map<unsigned int, NumericT> > cpu_backup(rows_);
914  tools::sparse_matrix_adapter<NumericT> adapted_cpu_backup(cpu_backup, rows_, cols_);
915  viennacl::copy(*this, adapted_cpu_backup);
916  cpu_backup[i][static_cast<unsigned int>(j)] = 0.0;
917  viennacl::copy(adapted_cpu_backup, *this);
918 
919  index = element_index(i, j);
920 
921  assert(index < nonzeros_);
922 
923  return entry_proxy<NumericT>(index, elements_);
924  }
925 
927  const vcl_size_t & size1() const { return rows_; }
929  const vcl_size_t & size2() const { return cols_; }
931  const vcl_size_t & nnz() const { return nonzeros_; }
933  const vcl_size_t & blocks1() const { return row_block_num_; }
934 
936  const handle_type & handle1() const { return row_buffer_; }
938  const handle_type & handle2() const { return col_buffer_; }
940  const handle_type & handle3() const { return row_blocks_; }
942  const handle_type & handle() const { return elements_; }
943 
945  handle_type & handle1() { return row_buffer_; }
947  handle_type & handle2() { return col_buffer_; }
949  handle_type & handle3() { return row_blocks_; }
951  handle_type & handle() { return elements_; }
952 
957  void switch_memory_context(viennacl::context new_ctx)
958  {
959  viennacl::backend::switch_memory_context<unsigned int>(row_buffer_, new_ctx);
960  viennacl::backend::switch_memory_context<unsigned int>(col_buffer_, new_ctx);
961  viennacl::backend::switch_memory_context<unsigned int>(row_blocks_, new_ctx);
962  viennacl::backend::switch_memory_context<NumericT>(elements_, new_ctx);
963  }
964 
966  viennacl::memory_types memory_context() const
967  {
968  return row_buffer_.get_active_handle_id();
969  }
970 
971 private:
972 
974  vcl_size_t element_index(vcl_size_t i, vcl_size_t j)
975  {
976  //read row indices
977  viennacl::backend::typesafe_host_array<unsigned int> row_indices(row_buffer_, 2);
978  viennacl::backend::memory_read(row_buffer_, row_indices.element_size()*i, row_indices.element_size()*2, row_indices.get());
979 
980  //get column indices for row i:
981  viennacl::backend::typesafe_host_array<unsigned int> col_indices(col_buffer_, row_indices[1] - row_indices[0]);
982  viennacl::backend::memory_read(col_buffer_, col_indices.element_size()*row_indices[0], row_indices.element_size()*col_indices.size(), col_indices.get());
983 
984  for (vcl_size_t k=0; k<col_indices.size(); ++k)
985  {
986  if (col_indices[k] == j)
987  return row_indices[0] + k;
988  }
989 
990  // if not found, return index past the end of the matrix (cf. matrix.end() in the spirit of the STL)
991  return nonzeros_;
992  }
993 
994 public:
999  void generate_row_block_information()
1000  {
1001  viennacl::backend::typesafe_host_array<unsigned int> row_buffer(row_buffer_, rows_ + 1);
1002  viennacl::backend::memory_read(row_buffer_, 0, row_buffer.raw_size(), row_buffer.get());
1003 
1004  viennacl::backend::typesafe_host_array<unsigned int> row_blocks(row_buffer_, rows_ + 1);
1005 
1006  vcl_size_t num_entries_in_current_batch = 0;
1007 
1008  const vcl_size_t shared_mem_size = 1024; // number of column indices loaded to shared memory, number of floating point values loaded to shared memory
1009 
1010  row_block_num_ = 0;
1011  row_blocks.set(0, 0);
1012  for (vcl_size_t i=0; i<rows_; ++i)
1013  {
1014  vcl_size_t entries_in_row = vcl_size_t(row_buffer[i+1]) - vcl_size_t(row_buffer[i]);
1015  num_entries_in_current_batch += entries_in_row;
1016 
1017  if (num_entries_in_current_batch > shared_mem_size)
1018  {
1019  vcl_size_t rows_in_batch = i - row_blocks[row_block_num_];
1020  if (rows_in_batch > 0) // at least one full row is in the batch. Use current row in next batch.
1021  row_blocks.set(++row_block_num_, i--);
1022  else // row is larger than buffer in shared memory
1023  row_blocks.set(++row_block_num_, i+1);
1024  num_entries_in_current_batch = 0;
1025  }
1026  }
1027  if (num_entries_in_current_batch > 0)
1028  row_blocks.set(++row_block_num_, rows_);
1029 
1030  if (row_block_num_ > 0) //matrix might be empty...
1031  viennacl::backend::memory_create(row_blocks_,
1032  row_blocks.element_size() * (row_block_num_ + 1),
1033  viennacl::traits::context(row_buffer_), row_blocks.get());
1034 
1035  }
1036 
1037 private:
1038  // /** @brief Copy constructor is by now not available. */
1039  //compressed_matrix(compressed_matrix const &);
1040 
1041 private:
1042 
1043  vcl_size_t rows_;
1044  vcl_size_t cols_;
1045  vcl_size_t nonzeros_;
1046  vcl_size_t row_block_num_;
1047  handle_type row_buffer_;
1048  handle_type row_blocks_;
1049  handle_type col_buffer_;
1050  handle_type elements_;
1051 };
1052 
1058 template<typename NumericT, unsigned int AlignmentV>
1059 std::ostream & operator<<(std::ostream & os, compressed_matrix<NumericT, AlignmentV> const & A)
1060 {
1061  std::vector<std::map<unsigned int, NumericT> > tmp(A.size1());
1062  viennacl::copy(A, tmp);
1063  os << "compressed_matrix of size (" << A.size1() << ", " << A.size2() << ") with " << A.nnz() << " nonzeros:" << std::endl;
1064 
1065  for (vcl_size_t i=0; i<A.size1(); ++i)
1066  {
1067  for (typename std::map<unsigned int, NumericT>::const_iterator it = tmp[i].begin(); it != tmp[i].end(); ++it)
1068  os << " (" << i << ", " << it->first << ")\t" << it->second << std::endl;
1069  }
1070  return os;
1071 }
1072 
1073 //
1074 // Specify available operations:
1075 //
1076 
1079 namespace linalg
1080 {
1081 namespace detail
1082 {
1083  // x = A * y
1084  template<typename T, unsigned int A>
1085  struct op_executor<vector_base<T>, op_assign, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
1086  {
1087  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
1088  {
1089  // check for the special case x = A * x
1090  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
1091  {
1092  viennacl::vector<T> temp(lhs);
1093  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
1094  lhs = temp;
1095  }
1096  else
1097  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), lhs, T(0));
1098  }
1099  };
1100 
1101  template<typename T, unsigned int A>
1102  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
1103  {
1104  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
1105  {
1106  // check for the special case x += A * x
1107  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
1108  {
1109  viennacl::vector<T> temp(lhs);
1110  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
1111  lhs += temp;
1112  }
1113  else
1114  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), lhs, T(1));
1115  }
1116  };
1117 
1118  template<typename T, unsigned int A>
1119  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> >
1120  {
1121  static void apply(vector_base<T> & lhs, vector_expression<const compressed_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
1122  {
1123  // check for the special case x -= A * x
1124  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
1125  {
1126  viennacl::vector<T> temp(lhs);
1127  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(1), temp, T(0));
1128  lhs -= temp;
1129  }
1130  else
1131  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), T(-1), lhs, T(1));
1132  }
1133  };
1134 
1135 
1136  // x = A * vec_op
1137  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
1138  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> >
1139  {
1140  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)
1141  {
1142  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
1143  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
1144  }
1145  };
1146 
1147  // x = A * vec_op
1148  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
1149  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> >
1150  {
1151  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)
1152  {
1153  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
1154  viennacl::vector<T> temp_result(lhs);
1155  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
1156  lhs += temp_result;
1157  }
1158  };
1159 
1160  // x = A * vec_op
1161  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
1162  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> >
1163  {
1164  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)
1165  {
1166  viennacl::vector<T> temp(rhs.rhs(), viennacl::traits::context(rhs));
1167  viennacl::vector<T> temp_result(lhs);
1168  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
1169  lhs -= temp_result;
1170  }
1171  };
1172 
1173 } // namespace detail
1174 } // namespace linalg
1176 }
1177 
1178 #endif
viennacl::compressed_matrix::size2
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:929
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::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::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::compressed_matrix::switch_memory_context
void switch_memory_context(viennacl::context new_ctx)
Switches the memory context of the matrix.
Definition: compressed_matrix.hpp:957
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
entry_proxy.hpp
A proxy class for entries in a vector.
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::zero_vector
Definition: vector_def.hpp:93
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::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::compressed_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:942
viennacl::compressed_matrix::handle1
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:936
viennacl::compressed_matrix::nnz
const vcl_size_t & nnz() const
Returns the number of nonzero entries.
Definition: compressed_matrix.hpp:931
viennacl::backend::detail::element_size
vcl_size_t element_size(memory_types)
Definition: memory.hpp:299
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::compressed_matrix::generate_row_block_information
void generate_row_block_information()
Builds the row block information needed for fast sparse matrix-vector multiplications.
Definition: compressed_matrix.hpp:999
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::compressed_matrix::operator()
entry_proxy< NumericT > 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:902
viennacl::compressed_matrix::handle3
handle_type & handle3()
Returns the OpenCL handle to the row block array.
Definition: compressed_matrix.hpp:949
viennacl::compressed_matrix::handle
handle_type & handle()
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:951
sparse_matrix_operations.hpp
Implementations of operations using sparse matrices.
viennacl::compressed_matrix::handle2
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:938
viennacl::vector_base< unsigned int >
viennacl::compressed_matrix::clear
void clear()
Resets all entries in the matrix back to zero without changing the matrix size. Resets the sparsity p...
Definition: compressed_matrix.hpp:888
viennacl::tools::const_sparse_matrix_adapter
Adapts a constant sparse matrix type made up from std::vector > to basic ub...
Definition: adapter.hpp:183
viennacl::compressed_matrix::value_type
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< NumericT >::ResultType > value_type
Definition: compressed_matrix.hpp:563
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::compressed_matrix::size_type
vcl_size_t size_type
Definition: compressed_matrix.hpp:564
viennacl::compressed_matrix::compressed_matrix
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:612
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::detail::copy_impl
void copy_impl(const CPUMatrixT &cpu_matrix, compressed_compressed_matrix< NumericT > &gpu_matrix, vcl_size_t nonzero_rows, vcl_size_t nonzeros)
Definition: compressed_compressed_matrix.hpp:41
viennacl::compressed_matrix::memory_context
viennacl::memory_types memory_context() const
Returns the current memory context to determine whether the matrix is set up for OpenMP, OpenCL, or CUDA.
Definition: compressed_matrix.hpp:966
viennacl::matrix_expression::rhs
RHS & rhs() const
Get right hand side operand.
Definition: matrix.hpp:69
viennacl::compressed_matrix::reserve
void reserve(vcl_size_t new_nonzeros, bool preserve=true)
Allocate memory for the supplied number of nonzeros in the matrix. Old values are preserved...
Definition: compressed_matrix.hpp:794
viennacl::compressed_matrix::handle3
const handle_type & handle3() const
Returns the OpenCL handle to the row block array.
Definition: compressed_matrix.hpp:940
viennacl::compressed_matrix::handle1
handle_type & handle1()
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:945
viennacl::compressed_matrix::compressed_matrix
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:576
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::backend::memory_copy
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
viennacl::compressed_matrix::compressed_matrix
compressed_matrix(viennacl::context ctx)
Creates an empty compressed_matrix, but sets the respective context information.
Definition: compressed_matrix.hpp:641
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::compressed_matrix::set
void set(const void *row_jumper, const void *col_buffer, const NumericT *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:764
viennacl::backend::typesafe_host_array::set
void set(vcl_size_t index, U value)
Definition: util.hpp:115
viennacl::compressed_matrix::handle_type
viennacl::backend::mem_handle handle_type
Definition: compressed_matrix.hpp:562
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::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 > to basic u...
Definition: adapter.hpp:357
viennacl::compressed_matrix
A sparse square matrix in compressed sparse rows format.
Definition: compressed_matrix.hpp:559
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::linalg::detail::min
T min(const T &lhs, const T &rhs)
Minimum.
Definition: util.hpp:45
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::compressed_matrix::compressed_matrix
compressed_matrix(matrix_expression< const compressed_matrix, const compressed_matrix, op_prod > const &proxy)
Assignment a compressed matrix from the product of two compressed_matrix objects (C = A * B)...
Definition: compressed_matrix.hpp:695
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::memory_types
memory_types
Definition: forwards.h:345
viennacl::matrix_expression::lhs
LHS & lhs() const
Get left hand side operand.
Definition: matrix.hpp:66
viennacl::backend::typesafe_host_array::raw_size
vcl_size_t raw_size() const
Definition: util.hpp:111
viennacl::compressed_matrix::blocks1
const vcl_size_t & blocks1() const
Returns the internal number of row blocks for an adaptive SpMV.
Definition: compressed_matrix.hpp:933
viennacl::entry_proxy
A proxy class for a single element of a vector or matrix. This proxy should not be noticed by end-use...
Definition: forwards.h:233
viennacl::compressed_matrix::compressed_matrix
compressed_matrix()
Default construction of a compressed matrix. No memory is allocated.
Definition: compressed_matrix.hpp:567
viennacl::compressed_matrix::resize
void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve=true)
Resize the matrix.
Definition: compressed_matrix.hpp:829
op_assign
Definition: self_assign.cpp:132
viennacl::backend::memory_shallow_copy
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::compressed_matrix::operator=
compressed_matrix & operator=(compressed_matrix const &other)
Assignment a compressed matrix from possibly another memory domain.
Definition: compressed_matrix.hpp:720
viennacl::compressed_matrix::handle2
handle_type & handle2()
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:947
viennacl::compressed_matrix::operator=
compressed_matrix & operator=(matrix_expression< const compressed_matrix, const compressed_matrix, op_prod > const &proxy)
Assignment a compressed matrix from the product of two compressed_matrix objects (C = A * B)...
Definition: compressed_matrix.hpp:739
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