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
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-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 "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 NumericT, unsigned int AlignmentV>
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<NumericT, AlignmentV> 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_[static_cast<vcl_size_t>(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 NumericT, unsigned int AlignmentV>
150 void copy(std::vector<NumericT>& cpu_vec, circulant_matrix<NumericT, AlignmentV>& 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 NumericT, unsigned int AlignmentV>
163 void copy(circulant_matrix<NumericT, AlignmentV>& gpu_mat, std::vector<NumericT>& 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 NumericT, unsigned int AlignmentV, typename MatrixT>
176 void copy(circulant_matrix<NumericT, AlignmentV>& circ_src, MatrixT& com_dst)
177 {
178  vcl_size_t size = circ_src.size1();
179  assert(size == viennacl::traits::size1(com_dst) && bool("Size mismatch"));
180  assert(size == viennacl::traits::size2(com_dst) && bool("Size mismatch"));
181  std::vector<NumericT> tmp(size);
182  copy(circ_src, tmp);
183 
184  for (vcl_size_t i = 0; i < size; i++)
185  {
186  for (vcl_size_t j = 0; j < size; j++)
187  {
188  long index = static_cast<long>(i) - static_cast<long>(j);
189  if (index < 0)
190  index += static_cast<long>(size);
191  com_dst(i, j) = tmp[static_cast<vcl_size_t>(index)];
192  }
193  }
194 }
195 
202 template<typename NumericT, unsigned int AlignmentV, typename MatrixT>
203 void copy(MatrixT& com_src, circulant_matrix<NumericT, AlignmentV>& circ_dst)
204 {
205  assert( (circ_dst.size1() == 0 || circ_dst.size1() == viennacl::traits::size1(com_src)) && bool("Size mismatch"));
206  assert( (circ_dst.size2() == 0 || circ_dst.size2() == viennacl::traits::size2(com_src)) && bool("Size mismatch"));
207 
209 
210  std::vector<NumericT> tmp(size);
211 
212  for (vcl_size_t i = 0; i < size; i++) tmp[i] = com_src(i, 0);
213 
214  copy(tmp, circ_dst);
215 }
216 
217 /*namespace linalg
218  {
219  template<typename NumericT, unsigned int AlignmentV, unsigned int VECTOR_AlignmentV>
220  void prod_impl(circulant_matrix<NumericT, AlignmentV> const & mat,
221  vector<NumericT, VECTOR_AlignmentV> const & vec,
222  vector<NumericT, VECTOR_AlignmentV>& result) {
223  viennacl::vector<NumericT, VECTOR_AlignmentV> circ(mat.elements().size() * 2);
224  fft::real_to_complex(mat.elements(), circ, mat.elements().size());
225 
226  viennacl::vector<NumericT, VECTOR_AlignmentV> tmp(vec.size() * 2);
227  viennacl::vector<NumericT, VECTOR_AlignmentV> tmp2(vec.size() * 2);
228 
229  fft::real_to_complex(vec, tmp, vec.size());
230  fft::convolve(circ, tmp, tmp2);
231  fft::complex_to_real(tmp2, result, vec.size());
232  }
233  }*/
234 
240 template<class NumericT, unsigned int AlignmentV>
241 std::ostream & operator<<(std::ostream& s, circulant_matrix<NumericT, AlignmentV>& gpu_matrix)
242 {
243  vcl_size_t size = gpu_matrix.size1();
244  std::vector<NumericT> tmp(size);
245  copy(gpu_matrix, tmp);
246  s << "[" << size << "," << size << "](";
247 
248  for (vcl_size_t i = 0; i < size; i++)
249  {
250  s << "(";
251  for (vcl_size_t j = 0; j < size; j++)
252  {
253  long index = static_cast<long>(i) - static_cast<long>(j);
254  if (index < 0) index = static_cast<long>(size) + index;
255  s << tmp[vcl_size_t(index)];
256  //s << index;
257  if (j < (size - 1)) s << ",";
258  }
259  s << ")";
260  }
261  s << ")";
262  return s;
263 }
264 
265 //
266 // Specify available operations:
267 //
268 
271 namespace linalg
272 {
273 namespace detail
274 {
275  // x = A * y
276  template<typename T, unsigned int A>
277  struct op_executor<vector_base<T>, op_assign, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
278  {
279  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
280  {
281  // check for the special case x = A * x
282  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
283  {
284  viennacl::vector<T> temp(lhs);
285  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
286  lhs = temp;
287  }
288  else
289  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
290  }
291  };
292 
293  template<typename T, unsigned int A>
294  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
295  {
296  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
297  {
298  viennacl::vector<T> temp(lhs);
299  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
300  lhs += temp;
301  }
302  };
303 
304  template<typename T, unsigned int A>
305  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> >
306  {
307  static void apply(vector_base<T> & lhs, vector_expression<const circulant_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
308  {
309  viennacl::vector<T> temp(lhs);
310  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
311  lhs -= temp;
312  }
313  };
314 
315 
316  // x = A * vec_op
317  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
318  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> >
319  {
320  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)
321  {
322  viennacl::vector<T> temp(rhs.rhs());
323  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
324  }
325  };
326 
327  // x = A * vec_op
328  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
329  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> >
330  {
331  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)
332  {
333  viennacl::vector<T> temp(rhs.rhs());
334  viennacl::vector<T> temp_result(lhs);
335  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
336  lhs += temp_result;
337  }
338  };
339 
340  // x = A * vec_op
341  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
342  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> >
343  {
344  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)
345  {
346  viennacl::vector<T> temp(rhs.rhs());
347  viennacl::vector<T> temp_result(lhs);
348  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
349  lhs -= temp_result;
350  }
351  };
352 
353 } // namespace detail
354 } // namespace linalg
355 
357 }
358 
359 #endif // VIENNACL_CIRCULANT_MATRIX_HPP
viennacl::vector< NumericT, AlignmentV > const & elements() const
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
viennacl::vector< NumericT, AlignmentV > & elements()
Returns an internal viennacl::vector, which represents a circulant matrix elements.
Implementations of operations using circulant_matrix. Experimental.
viennacl::backend::mem_handle handle_type
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
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< NumericT >::ResultType > value_type
This file provides the forward declarations for the main types used within ViennaCL.
void resize(vcl_size_t sz, bool preserve=true)
Resizes the matrix. Existing entries can be preserved.
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
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
vcl_size_t internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
vcl_size_t size1() const
Returns the number of rows of the matrix.
void resize(size_type new_size, bool preserve=true)
Resizes the allocated memory for the vector. Pads the memory to be a multiple of 'AlignmentV'.
Definition: vector.hpp:1046
std::size_t vcl_size_t
Definition: forwards.h:75
entry_proxy< NumericT > operator()(vcl_size_t row_index, vcl_size_t col_index)
Read-write access to a single element of the matrix.
circulant_matrix()
The default constructor. Does not allocate any memory.
Implementations of the OpenCL backend, where all contexts are stored in.
handle_type const & handle() const
Returns the OpenCL handle.
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
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) ...
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
All routines related to the Fast Fourier Transform. Experimental.
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
A Circulant matrix class.
circulant_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
vcl_size_t size2() const
Returns the number of columns of the matrix.
void prod_impl(const matrix_base< NumericT > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
size_type internal_size() const
Returns the internal length of the vector, which is given by size() plus the extra memory due to padd...
Definition: vector_def.hpp:120
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
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
const handle_type & handle() const
Returns the memory handle.
Definition: vector_def.hpp:128
circulant_matrix< NumericT, AlignmentV > & operator+=(circulant_matrix< NumericT, AlignmentV > &that)
+= operation for circulant matrices