ViennaCL - The Vienna Computing Library  1.5.1
sparse_matrix_operations.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_SPARSE_MATRIX_OPERATIONS_HPP_
2 #define VIENNACL_LINALG_SPARSE_MATRIX_OPERATIONS_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 "viennacl/forwards.h"
26 #include "viennacl/scalar.hpp"
27 #include "viennacl/vector.hpp"
28 #include "viennacl/matrix.hpp"
29 #include "viennacl/tools/tools.hpp"
31 
32 #ifdef VIENNACL_WITH_OPENCL
34 #endif
35 
36 #ifdef VIENNACL_WITH_CUDA
38 #endif
39 
40 namespace viennacl
41 {
42  namespace linalg
43  {
44 
45  namespace detail
46  {
47 
48  template<typename SparseMatrixType, typename SCALARTYPE, unsigned int VEC_ALIGNMENT>
50  row_info(SparseMatrixType const & mat,
52  row_info_types info_selector)
53  {
54  switch (viennacl::traits::handle(mat).get_active_handle_id())
55  {
57  viennacl::linalg::host_based::detail::row_info(mat, vec, info_selector);
58  break;
59 #ifdef VIENNACL_WITH_OPENCL
61  viennacl::linalg::opencl::detail::row_info(mat, vec, info_selector);
62  break;
63 #endif
64 #ifdef VIENNACL_WITH_CUDA
66  viennacl::linalg::cuda::detail::row_info(mat, vec, info_selector);
67  break;
68 #endif
70  throw memory_exception("not initialised!");
71  default:
72  throw memory_exception("not implemented");
73  }
74  }
75 
76  }
77 
78 
79 
80  // A * x
81 
90  template<typename SparseMatrixType, class ScalarType>
92  prod_impl(const SparseMatrixType & mat,
95  {
96  assert( (mat.size1() == result.size()) && bool("Size check failed for compressed matrix-vector product: size1(mat) != size(result)"));
97  assert( (mat.size2() == vec.size()) && bool("Size check failed for compressed matrix-vector product: size2(mat) != size(x)"));
98 
100  {
102  viennacl::linalg::host_based::prod_impl(mat, vec, result);
103  break;
104 #ifdef VIENNACL_WITH_OPENCL
106  viennacl::linalg::opencl::prod_impl(mat, vec, result);
107  break;
108 #endif
109 #ifdef VIENNACL_WITH_CUDA
111  viennacl::linalg::cuda::prod_impl(mat, vec, result);
112  break;
113 #endif
115  throw memory_exception("not initialised!");
116  default:
117  throw memory_exception("not implemented");
118  }
119  }
120 
121 
122  // A * B
131  template<typename SparseMatrixType, class ScalarType, typename F1, typename F2>
133  prod_impl(const SparseMatrixType & sp_mat,
136  {
137  assert( (sp_mat.size1() == result.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size1(sp_mat) != size1(result)"));
138  assert( (sp_mat.size2() == d_mat.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size2(sp_mat) != size1(d_mat)"));
139 
141  {
143  viennacl::linalg::host_based::prod_impl(sp_mat, d_mat, result);
144  break;
145 #ifdef VIENNACL_WITH_OPENCL
147  viennacl::linalg::opencl::prod_impl(sp_mat, d_mat, result);
148  break;
149 #endif
150 #ifdef VIENNACL_WITH_CUDA
152  viennacl::linalg::cuda::prod_impl(sp_mat, d_mat, result);
153  break;
154 #endif
156  throw memory_exception("not initialised!");
157  default:
158  throw memory_exception("not implemented");
159  }
160  }
161 
162  // A * transpose(B)
171  template<typename SparseMatrixType, class ScalarType, typename F1, typename F2>
173  prod_impl(const SparseMatrixType & sp_mat,
176  viennacl::op_trans>& d_mat,
178  {
179  assert( (sp_mat.size1() == result.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size1(sp_mat) != size1(result)"));
180  assert( (sp_mat.size2() == d_mat.size1()) && bool("Size check failed for compressed matrix - dense matrix product: size2(sp_mat) != size1(d_mat)"));
181 
183  {
185  viennacl::linalg::host_based::prod_impl(sp_mat, d_mat, result);
186  break;
187 #ifdef VIENNACL_WITH_OPENCL
189  viennacl::linalg::opencl::prod_impl(sp_mat, d_mat, result);
190  break;
191 #endif
192 #ifdef VIENNACL_WITH_CUDA
194  viennacl::linalg::cuda::prod_impl(sp_mat, d_mat, result);
195  break;
196 #endif
198  throw memory_exception("not initialised!");
199  default:
200  throw memory_exception("not implemented");
201  }
202  }
203 
210  template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
212  inplace_solve(const SparseMatrixType & mat,
214  SOLVERTAG tag)
215  {
216  assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on compressed matrix: size1(mat) != size2(mat)"));
217  assert( (mat.size2() == vec.size()) && bool("Size check failed for compressed matrix-vector product: size2(mat) != size(x)"));
218 
220  {
223  break;
224 #ifdef VIENNACL_WITH_OPENCL
227  break;
228 #endif
229 #ifdef VIENNACL_WITH_CUDA
232  break;
233 #endif
235  throw memory_exception("not initialised!");
236  default:
237  throw memory_exception("not implemented");
238  }
239  }
240 
241 
248  template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
252  SOLVERTAG tag)
253  {
254  assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on transposed compressed matrix: size1(mat) != size2(mat)"));
255  assert( (mat.size1() == vec.size()) && bool("Size check failed for transposed compressed matrix triangular solve: size1(mat) != size(x)"));
256 
257  switch (viennacl::traits::handle(mat.lhs()).get_active_handle_id())
258  {
261  break;
262 #ifdef VIENNACL_WITH_OPENCL
265  break;
266 #endif
267 #ifdef VIENNACL_WITH_CUDA
270  break;
271 #endif
273  throw memory_exception("not initialised!");
274  default:
275  throw memory_exception("not implemented");
276  }
277  }
278 
279 
280 
281  namespace detail
282  {
283 
284  template<typename SparseMatrixType, class ScalarType, typename SOLVERTAG>
287  viennacl::backend::mem_handle const & block_index_array, vcl_size_t num_blocks,
288  viennacl::vector_base<ScalarType> const & mat_diagonal,
290  SOLVERTAG tag)
291  {
292  assert( (mat.size1() == mat.size2()) && bool("Size check failed for triangular solve on transposed compressed matrix: size1(mat) != size2(mat)"));
293  assert( (mat.size1() == vec.size()) && bool("Size check failed for transposed compressed matrix triangular solve: size1(mat) != size(x)"));
294 
295  switch (viennacl::traits::handle(mat.lhs()).get_active_handle_id())
296  {
298  viennacl::linalg::host_based::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
299  break;
300  #ifdef VIENNACL_WITH_OPENCL
302  viennacl::linalg::opencl::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
303  break;
304  #endif
305  #ifdef VIENNACL_WITH_CUDA
307  viennacl::linalg::cuda::detail::block_inplace_solve(mat, block_index_array, num_blocks, mat_diagonal, vec, tag);
308  break;
309  #endif
311  throw memory_exception("not initialised!");
312  default:
313  throw memory_exception("not implemented");
314  }
315  }
316 
317 
318  }
319 
320 
321 
322  } //namespace linalg
323 
324 
326  template<typename M1>
328  matrix_expression< const M1, const M1, op_trans>
329  >::type
330  trans(const M1 & mat)
331  {
333  }
334 
335  //free functions:
341  template <typename SCALARTYPE, typename SparseMatrixType>
345  const viennacl::vector_expression< const SparseMatrixType, const viennacl::vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
346  {
347  assert(proxy.lhs().size1() == result.size() && bool("Dimensions for addition of sparse matrix-vector product to vector don't match!"));
348  vector<SCALARTYPE> temp(proxy.lhs().size1());
349  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), temp);
350  result += temp;
351  return result;
352  }
353 
359  template <typename SCALARTYPE, typename SparseMatrixType>
363  const viennacl::vector_expression< const SparseMatrixType, const viennacl::vector_base<SCALARTYPE>, viennacl::op_prod> & proxy)
364  {
365  assert(proxy.lhs().size1() == result.size() && bool("Dimensions for addition of sparse matrix-vector product to vector don't match!"));
366  vector<SCALARTYPE> temp(proxy.lhs().size1());
367  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), temp);
368  result += temp;
369  return result;
370  }
371 
372 } //namespace viennacl
373 
374 
375 #endif
Simple enable-if variant that uses the SFINAE pattern.
Definition: enable_if.hpp:29
size_type size1() const
Returns the number of rows.
Definition: matrix.hpp:625
std::size_t vcl_size_t
Definition: forwards.h:58
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:859
Definition: forwards.h:478
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:330
Exception class in case of memory errors.
Definition: forwards.h:485
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for dense triangular systems. Matlab notation: A \ B.
Definition: direct_solve.hpp:54
Implementation of the dense matrix class.
Various little tools used here and there in ViennaCL.
Implementations of operations using sparse matrices on the CPU using a single thread or OpenMP...
A dense matrix class.
Definition: forwards.h:290
void row_info(compressed_matrix< ScalarType, MAT_ALIGNMENT > const &mat, vector_base< ScalarType > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:100
Expression template class for representing a tree of expressions which ultimately result in a matrix...
Definition: forwards.h:283
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:297
This file provides the forward declarations for the main types used within ViennaCL.
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:737
An expression template class that represents a binary operation that yields a vector.
Definition: forwards.h:181
Definition: forwards.h:481
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:130
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
Definition: forwards.h:480
Implementations of operations using sparse matrices using CUDA.
vcl_size_t size1() const
Returns the size of the result vector.
Definition: matrix.hpp:180
vcl_size_t size2() const
Definition: matrix.hpp:181
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
void block_inplace_solve(const matrix_expression< const compressed_matrix< ScalarType, MAT_ALIGNMENT >, const compressed_matrix< ScalarType, MAT_ALIGNMENT >, op_trans > &L, viennacl::backend::mem_handle const &block_indices, vcl_size_t num_blocks, vector_base< ScalarType > const &, vector_base< ScalarType > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:613
viennacl::vector< NumericT > operator-(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 - A * v2', where A is a matrix.
Definition: matrix_operations.hpp:858
void block_inplace_solve(const matrix_expression< const compressed_matrix< ScalarType, MAT_ALIGNMENT >, const compressed_matrix< ScalarType, MAT_ALIGNMENT >, op_trans > &L, viennacl::backend::mem_handle const &block_indices, vcl_size_t num_blocks, vector_base< ScalarType > const &, vector_base< ScalarType > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:289
Implementations of operations using sparse matrices and OpenCL.
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for dense triangular systems. Matlab notation: A \ B.
Definition: direct_solve.hpp:75
void block_inplace_solve(const matrix_expression< const compressed_matrix< ScalarType, MAT_ALIGNMENT >, const compressed_matrix< ScalarType, MAT_ALIGNMENT >, op_trans > &L, viennacl::backend::mem_handle const &, vcl_size_t, vector_base< ScalarType > const &, vector_base< ScalarType > &vec, viennacl::linalg::unit_lower_tag)
Definition: sparse_matrix_operations.hpp:620
A tag class representing matrix-vector products and element-wise multiplications. ...
Definition: forwards.h:76
row_info_types
Definition: forwards.h:691
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
Definition: forwards.h:479
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
A tag class representing transposed matrices.
Definition: forwards.h:165
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:547
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
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
LHS & lhs() const
Get left hand side operand.
Definition: matrix.hpp:174
viennacl::enable_if< viennacl::is_any_sparse_matrix< SparseMatrixType >::value >::type row_info(SparseMatrixType const &mat, vector< SCALARTYPE, VEC_ALIGNMENT > &vec, row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:50
void row_info(compressed_matrix< SCALARTYPE, MAT_ALIGNMENT > const &mat, vector_base< SCALARTYPE > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:53
Implementation of the ViennaCL scalar class.
viennacl::vector< NumericT > operator+(const vector_base< NumericT > &v1, const vector_expression< const matrix_base< NumericT, F >, const vector_base< NumericT >, op_prod > &proxy)
Implementation of the operation 'result = v1 + A * v2', where A is a matrix.
Definition: matrix_operations.hpp:840
void row_info(compressed_matrix< ScalarType, MAT_ALIGNMENT > const &mat, vector_base< ScalarType > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:47
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:1364
viennacl::enable_if< viennacl::is_any_sparse_matrix< SparseMatrixType >::value >::type block_inplace_solve(const matrix_expression< const SparseMatrixType, const SparseMatrixType, op_trans > &mat, viennacl::backend::mem_handle const &block_index_array, vcl_size_t num_blocks, viennacl::vector_base< ScalarType > const &mat_diagonal, viennacl::vector_base< ScalarType > &vec, SOLVERTAG tag)
Definition: sparse_matrix_operations.hpp:286