ViennaCL - The Vienna Computing Library  1.5.2
toeplitz_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_TOEPLITZ_MATRIX_HPP
2 #define VIENNACL_TOEPLITZ_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 
29 #include "viennacl/fft.hpp"
30 
32 
33 
34 namespace viennacl {
35 
41  template<class SCALARTYPE, unsigned int ALIGNMENT>
42  class toeplitz_matrix
43  {
44  public:
47 
52  explicit toeplitz_matrix() {}
53 
59  explicit toeplitz_matrix(vcl_size_t rows, vcl_size_t cols) : elements_(rows * 2)
60  {
61  assert(rows == cols && bool("Toeplitz matrix must be square!"));
62  (void)cols; // avoid 'unused parameter' warning in optimized builds
63  }
64 
65 
72  void resize(vcl_size_t sz, bool preserve = true) {
73  elements_.resize(sz * 2, preserve);
74  }
75 
80  handle_type const & handle() const { return elements_.handle(); }
81 
87  viennacl::vector<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
88 
89 
93  vcl_size_t size1() const { return elements_.size() / 2; }
94 
98  vcl_size_t size2() const { return elements_.size() / 2; }
99 
105  vcl_size_t internal_size() const { return elements_.internal_size(); }
106 
107 
116  {
117  assert(row_index < size1() && col_index < size2() && bool("Invalid access"));
118 
119  long index = static_cast<long>(col_index) - static_cast<long>(row_index);
120 
121  if (index < 0)
122  index = -index;
123  else if
124  (index > 0) index = 2 * static_cast<long>(size1()) - index;
125  return elements_[index];
126  }
127 
128 
136  elements_ += that.elements();
137  return *this;
138  }
139 
140  private:
141  toeplitz_matrix(toeplitz_matrix const &) {}
142  toeplitz_matrix & operator=(toeplitz_matrix const & t);
143 
144 
146  };
147 
154  template <typename SCALARTYPE, unsigned int ALIGNMENT>
155  void copy(std::vector<SCALARTYPE> const & cpu_vec, toeplitz_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat)
156  {
157  assert( (gpu_mat.size1() == 0 || (gpu_mat.size1() * 2 - 1) == cpu_vec.size()) && bool("Size mismatch"));
158 
159  vcl_size_t size = gpu_mat.size1();
160  std::vector<SCALARTYPE> rvrs(cpu_vec.size());
161  std::copy(cpu_vec.begin(), cpu_vec.end(), rvrs.begin());
162  std::reverse(rvrs.begin(), rvrs.end());
163 
164  std::vector<SCALARTYPE> tmp(size * 2);
165  std::copy(rvrs.begin() + size - 1, rvrs.end(), tmp.begin());
166  std::copy(rvrs.begin(), rvrs.begin() + size - 1, tmp.begin() + size + 1);
167  tmp[size] = 0.0;
168  copy(tmp, gpu_mat.elements());
169  }
170 
177  template <typename SCALARTYPE, unsigned int ALIGNMENT>
178  void copy(toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & gpu_mat, std::vector<SCALARTYPE> & cpu_vec)
179  {
180  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
181 
182  vcl_size_t size = gpu_mat.size1();
183  std::vector<SCALARTYPE> tmp(size * 2);
184  copy(gpu_mat.elements(), tmp);
185  std::reverse(tmp.begin(), tmp.end());
186 
187  std::copy(tmp.begin(), tmp.begin() + size - 1, cpu_vec.begin() + size);
188  std::copy(tmp.begin() + size, tmp.end(), cpu_vec.begin());
189 
190  }
191 
198  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
199  void copy(toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & tep_src, MATRIXTYPE & com_dst)
200  {
201  assert(tep_src.size1() == viennacl::traits::size1(com_dst) && bool("Size mismatch"));
202  assert(tep_src.size2() == viennacl::traits::size2(com_dst) && bool("Size mismatch"));
203 
204  vcl_size_t size = tep_src.size1();
205  std::vector<SCALARTYPE> tmp(tep_src.size1() * 2 - 1);
206  copy(tep_src, tmp);
207 
208  for(vcl_size_t i = 0; i < size; i++)
209  for(vcl_size_t j = 0; j < size; j++)
210  com_dst(i, j) = tmp[static_cast<int>(j) - static_cast<int>(i) + static_cast<int>(size) - 1];
211  }
212 
219  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
220  void copy(MATRIXTYPE const & com_src, toeplitz_matrix<SCALARTYPE, ALIGNMENT>& tep_dst)
221  {
222  assert( (tep_dst.size1() == 0 || tep_dst.size1() == viennacl::traits::size1(com_src)) && bool("Size mismatch"));
223  assert( (tep_dst.size2() == 0 || tep_dst.size2() == viennacl::traits::size2(com_src)) && bool("Size mismatch"));
224 
225  vcl_size_t size = tep_dst.size1();
226 
227  std::vector<SCALARTYPE> tmp(2*size - 1);
228 
229  for(long i = static_cast<long>(size) - 1; i >= 0; i--)
230  tmp[size - i - 1] = com_src(i, 0);
231 
232  for(vcl_size_t i = 1; i < size; i++)
233  tmp[size + i - 1] = com_src(0, i);
234 
235  copy(tmp, tep_dst);
236  }
237 
238  /*template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
239  void prod_impl(toeplitz_matrix<SCALARTYPE, ALIGNMENT>& mat,
240  vector<SCALARTYPE, VECTOR_ALIGNMENT>& vec,
241  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result) {
242  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tep(mat.elements().size() * 2);
243  fft::real_to_complex(mat.elements(), tep, mat.elements().size());
244 
245  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp(vec.size() * 4);
246  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp2(vec.size() * 4);
247 
248  tmp.clear();
249  copy(vec, tmp);
250  fft::real_to_complex(tmp, tmp2, vec.size() * 2);
251  fft::convolve(tep, tmp2, tmp);
252  fft::complex_to_real(tmp, tmp2, vec.size() * 2);
253  copy(tmp2.begin(), tmp2.begin() + vec.size(), result.begin());
254  }*/
255 
261  template<class SCALARTYPE, unsigned int ALIGNMENT>
262  std::ostream & operator<<(std::ostream & s, toeplitz_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
263  {
264  vcl_size_t size = gpu_matrix.size1();
265  std::vector<SCALARTYPE> tmp(2*size - 1);
266  copy(gpu_matrix, tmp);
267  s << "[" << size << "," << size << "](";
268 
269  for(vcl_size_t i = 0; i < size; i++) {
270  s << "(";
271  for(vcl_size_t j = 0; j < size; j++) {
272  s << tmp[static_cast<int>(j) - static_cast<int>(i) + static_cast<int>(size - 1)];
273  //s << (int)i - (int)j;
274  if(j < (size - 1)) s << ",";
275  }
276  s << ")";
277  }
278  s << ")";
279  return s;
280  }
281 
282  //
283  // Specify available operations:
284  //
285 
288  namespace linalg
289  {
290  namespace detail
291  {
292  // x = A * y
293  template <typename T, unsigned int A>
294  struct op_executor<vector_base<T>, op_assign, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
295  {
296  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
297  {
298  // check for the special case x = A * x
299  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
300  {
301  viennacl::vector<T> temp(lhs);
302  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
303  lhs = temp;
304  }
305  else
306  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
307  }
308  };
309 
310  template <typename T, unsigned int A>
311  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
312  {
313  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
314  {
315  viennacl::vector<T> temp(lhs);
316  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
317  lhs += temp;
318  }
319  };
320 
321  template <typename T, unsigned int A>
322  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
323  {
324  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
325  {
326  viennacl::vector<T> temp(lhs);
327  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
328  lhs -= temp;
329  }
330  };
331 
332 
333  // x = A * vec_op
334  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
335  struct op_executor<vector_base<T>, op_assign, vector_expression<const toeplitz_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
336  {
337  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
338  {
339  viennacl::vector<T> temp(rhs.rhs());
340  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
341  }
342  };
343 
344  // x = A * vec_op
345  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
346  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const toeplitz_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> >
347  {
348  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
349  {
350  viennacl::vector<T> temp(rhs.rhs());
351  viennacl::vector<T> temp_result(lhs);
352  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
353  lhs += temp_result;
354  }
355  };
356 
357  // x = A * vec_op
358  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
359  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const toeplitz_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
360  {
361  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
362  {
363  viennacl::vector<T> temp(rhs.rhs());
364  viennacl::vector<T> temp_result(lhs);
365  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
366  lhs -= temp_result;
367  }
368  };
369 
370  } // namespace detail
371  } // namespace linalg
372 
375 }
376 
377 #endif // VIENNACL_TOEPLITZ_MATRIX_HPP
std::size_t vcl_size_t
Definition: forwards.h:58
Implementations of operations using toeplitz_matrix. Experimental.
viennacl::backend::mem_handle handle_type
Definition: toeplitz_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
viennacl::vector< SCALARTYPE, ALIGNMENT > & elements()
Returns an internal viennacl::vector, which represents a Toeplitz matrix elements.
Definition: toeplitz_matrix.hpp:86
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
This file provides the forward declarations for the main types used within ViennaCL.
vcl_size_t size2() const
Returns the number of columns of the matrix.
Definition: toeplitz_matrix.hpp:98
handle_type const & handle() const
Returns the OpenCL handle.
Definition: toeplitz_matrix.hpp:80
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
toeplitz_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
Definition: toeplitz_matrix.hpp:59
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
toeplitz_matrix< SCALARTYPE, ALIGNMENT > & operator+=(toeplitz_matrix< SCALARTYPE, ALIGNMENT > &that)
+= operation for Toeplitz matrices
Definition: toeplitz_matrix.hpp:135
A Toeplitz matrix class.
Definition: forwards.h:330
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
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: toeplitz_matrix.hpp:115
toeplitz_matrix()
The default constructor. Does not allocate any memory.
Definition: toeplitz_matrix.hpp:52
Implementations of the OpenCL backend, where all contexts are stored in.
void resize(vcl_size_t sz, bool preserve=true)
Resizes the matrix. Existing entries can be preserved.
Definition: toeplitz_matrix.hpp:72
viennacl::vector< SCALARTYPE, ALIGNMENT > const & elements() const
Definition: toeplitz_matrix.hpp:87
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
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 internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
Definition: toeplitz_matrix.hpp:105
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
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 copy(MATRIXTYPE const &com_src, toeplitz_matrix< SCALARTYPE, ALIGNMENT > &tep_dst)
Copies a the matrix-like object to the Toeplitz matrix from the OpenCL device (either GPU or multi-co...
Definition: toeplitz_matrix.hpp:220
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: toeplitz_matrix.hpp:46
vcl_size_t size1() const
Returns the number of rows of the matrix.
Definition: toeplitz_matrix.hpp:93