ViennaCL - The Vienna Computing Library  1.5.2
circulant_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_CIRCULANT_MATRIX_HPP
2 #define VIENNACL_CIRCULANT_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 "viennacl/forwards.h"
26 #include "viennacl/vector.hpp"
27 #include "viennacl/ocl/backend.hpp"
28 
30 
31 #include "viennacl/fft.hpp"
32 
33 namespace viennacl
34 {
40  template<class SCALARTYPE, unsigned int ALIGNMENT>
42  {
43  public:
46 
51  explicit circulant_matrix() {}
52 
59  explicit circulant_matrix(vcl_size_t rows, vcl_size_t cols) : elements_(rows)
60  {
61  assert(rows == cols && bool("Circulant matrix must be square!"));
62  (void)cols; // avoid 'unused parameter' warning in optimized builds
63  }
64 
71  void resize(vcl_size_t sz, bool preserve = true)
72  {
73  elements_.resize(sz, preserve);
74  }
75 
80  handle_type const & handle() const { return elements_.handle(); }
81 
87  viennacl::vector<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
88 
92  vcl_size_t size1() const { return elements_.size(); }
93 
97  vcl_size_t size2() const { return elements_.size(); }
98 
104  vcl_size_t internal_size() const { return elements_.internal_size(); }
105 
114  {
115  long index = static_cast<long>(row_index) - static_cast<long>(col_index);
116 
117  assert(row_index < size1() && col_index < size2() && bool("Invalid access"));
118 
119  while (index < 0)
120  index += static_cast<long>(size1());
121  return elements_[index];
122  }
123 
131  {
132  elements_ += that.elements();
133  return *this;
134  }
135 
136  private:
138  circulant_matrix & operator=(circulant_matrix const & t);
139 
141  };
142 
149  template <typename SCALARTYPE, unsigned int ALIGNMENT>
150  void copy(std::vector<SCALARTYPE>& cpu_vec, circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat)
151  {
152  assert( (gpu_mat.size1() == 0 || cpu_vec.size() == gpu_mat.size1()) && bool("Size mismatch"));
153  copy(cpu_vec, gpu_mat.elements());
154  }
155 
162  template <typename SCALARTYPE, unsigned int ALIGNMENT>
163  void copy(circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat, std::vector<SCALARTYPE>& cpu_vec)
164  {
165  assert(cpu_vec.size() == gpu_mat.size1() && bool("Size mismatch"));
166  copy(gpu_mat.elements(), cpu_vec);
167  }
168 
175  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
176  void copy(circulant_matrix<SCALARTYPE, ALIGNMENT>& circ_src, MATRIXTYPE& com_dst) {
177  vcl_size_t size = circ_src.size1();
178  assert(size == viennacl::traits::size1(com_dst) && bool("Size mismatch"));
179  assert(size == viennacl::traits::size2(com_dst) && bool("Size mismatch"));
180  std::vector<SCALARTYPE> tmp(size);
181  copy(circ_src, tmp);
182 
183  for (vcl_size_t i = 0; i < size; i++) {
184  for (vcl_size_t j = 0; j < size; j++) {
185  long index = static_cast<long>(i) - static_cast<long>(j);
186  if (index < 0)
187  index = static_cast<long>(size + index);
188  com_dst(i, j) = tmp[index];
189  }
190  }
191  }
192 
199  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
200  void copy(MATRIXTYPE& com_src, circulant_matrix<SCALARTYPE, ALIGNMENT>& circ_dst)
201  {
202  assert( (circ_dst.size1() == 0 || circ_dst.size1() == viennacl::traits::size1(com_src)) && bool("Size mismatch"));
203  assert( (circ_dst.size2() == 0 || circ_dst.size2() == viennacl::traits::size2(com_src)) && bool("Size mismatch"));
204 
206 
207  std::vector<SCALARTYPE> tmp(size);
208 
209  for(vcl_size_t i = 0; i < size; i++) tmp[i] = com_src(i, 0);
210 
211  copy(tmp, circ_dst);
212  }
213 
214  /*namespace linalg
215  {
216  template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
217  void prod_impl(circulant_matrix<SCALARTYPE, ALIGNMENT> const & mat,
218  vector<SCALARTYPE, VECTOR_ALIGNMENT> const & vec,
219  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result) {
220  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> circ(mat.elements().size() * 2);
221  fft::real_to_complex(mat.elements(), circ, mat.elements().size());
222 
223  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp(vec.size() * 2);
224  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp2(vec.size() * 2);
225 
226  fft::real_to_complex(vec, tmp, vec.size());
227  fft::convolve(circ, tmp, tmp2);
228  fft::complex_to_real(tmp2, result, vec.size());
229  }
230  }*/
231 
237  template<class SCALARTYPE, unsigned int ALIGNMENT>
238  std::ostream & operator<<(std::ostream& s, circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
239  {
240  vcl_size_t size = gpu_matrix.size1();
241  std::vector<SCALARTYPE> tmp(size);
242  copy(gpu_matrix, tmp);
243  s << "[" << size << "," << size << "](";
244 
245  for(vcl_size_t i = 0; i < size; i++) {
246  s << "(";
247  for(vcl_size_t j = 0; j < size; j++) {
248  long index = static_cast<long>(i) - static_cast<long>(j);
249  if(index < 0) index = static_cast<long>(size) + index;
250  s << tmp[index];
251  //s << index;
252  if(j < (size - 1)) s << ",";
253  }
254  s << ")";
255  }
256  s << ")";
257  return s;
258  }
259 
260  //
261  // Specify available operations:
262  //
263 
266  namespace linalg
267  {
268  namespace detail
269  {
270  // x = A * y
271  template <typename T, unsigned int A>
272  struct op_executor<vector_base<T>, op_assign, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
273  {
274  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
275  {
276  // check for the special case x = A * x
277  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
278  {
279  viennacl::vector<T> temp(lhs);
280  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
281  lhs = temp;
282  }
283  else
284  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
285  }
286  };
287 
288  template <typename T, unsigned int A>
289  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
290  {
291  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
292  {
293  viennacl::vector<T> temp(lhs);
294  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
295  lhs += temp;
296  }
297  };
298 
299  template <typename T, unsigned int A>
300  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
301  {
302  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
303  {
304  viennacl::vector<T> temp(lhs);
305  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
306  lhs -= temp;
307  }
308  };
309 
310 
311  // x = A * vec_op
312  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
313  struct op_executor<vector_base<T>, op_assign, vector_expression<const circulant_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
314  {
315  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
316  {
317  viennacl::vector<T> temp(rhs.rhs());
318  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
319  }
320  };
321 
322  // x = A * vec_op
323  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
324  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const circulant_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> >
325  {
326  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
327  {
328  viennacl::vector<T> temp(rhs.rhs());
329  viennacl::vector<T> temp_result(lhs);
330  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
331  lhs += temp_result;
332  }
333  };
334 
335  // x = A * vec_op
336  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
337  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const circulant_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
338  {
339  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
340  {
341  viennacl::vector<T> temp(rhs.rhs());
342  viennacl::vector<T> temp_result(lhs);
343  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
344  lhs -= temp_result;
345  }
346  };
347 
348  } // namespace detail
349  } // namespace linalg
350 
352 }
353 
354 #endif // VIENNACL_CIRCULANT_MATRIX_HPP
std::size_t vcl_size_t
Definition: forwards.h:58
circulant_matrix< SCALARTYPE, ALIGNMENT > & operator+=(circulant_matrix< SCALARTYPE, ALIGNMENT > &that)
+= operation for circulant matrices
Definition: circulant_matrix.hpp:130
Implementations of operations using circulant_matrix. Experimental.
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: circulant_matrix.hpp:45
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:172
entry_proxy< SCALARTYPE > operator()(vcl_size_t row_index, vcl_size_t col_index)
Read-write access to a single element of the matrix.
Definition: circulant_matrix.hpp:113
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
handle_type const & handle() const
Returns the OpenCL handle.
Definition: circulant_matrix.hpp:80
viennacl::backend::mem_handle handle_type
Definition: circulant_matrix.hpp:44
This file provides the forward declarations for the main types used within ViennaCL.
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
circulant_matrix()
The default constructor. Does not allocate any memory.
Definition: circulant_matrix.hpp:51
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
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
viennacl::vector< SCALARTYPE, ALIGNMENT > & elements()
Returns an internal viennacl::vector, which represents a circulant matrix elements.
Definition: circulant_matrix.hpp:86
circulant_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
Definition: circulant_matrix.hpp:59
viennacl::vector< SCALARTYPE, ALIGNMENT > const & elements() const
Definition: circulant_matrix.hpp:87
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
Implementations of the OpenCL backend, where all contexts are stored in.
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
vcl_size_t size2() const
Returns the number of columns of the matrix.
Definition: circulant_matrix.hpp:97
vcl_size_t internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
Definition: circulant_matrix.hpp:104
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
All routines related to the Fast Fourier Transform. Experimental.
vcl_size_t size1() const
Returns the number of rows of the matrix.
Definition: circulant_matrix.hpp:92
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
A Circulant matrix class.
Definition: circulant_matrix.hpp:41
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
void resize(vcl_size_t sz, bool preserve=true)
Resizes the matrix. Existing entries can be preserved.
Definition: circulant_matrix.hpp:71