ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
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-2015, 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 
29 #include "viennacl/fft.hpp"
30 
32 
33 
34 namespace viennacl
35 {
36 
42 template<class NumericT, unsigned int AlignmentV>
43 class toeplitz_matrix
44 {
45 public:
48 
53  explicit toeplitz_matrix() {}
54 
60  explicit toeplitz_matrix(vcl_size_t rows, vcl_size_t cols) : elements_(rows * 2)
61  {
62  assert(rows == cols && bool("Toeplitz matrix must be square!"));
63  (void)cols; // avoid 'unused parameter' warning in optimized builds
64  }
65 
66 
73  void resize(vcl_size_t sz, bool preserve = true)
74  {
75  elements_.resize(sz * 2, preserve);
76  }
77 
82  handle_type const & handle() const { return elements_.handle(); }
83 
89  viennacl::vector<NumericT, AlignmentV> const & elements() const { return elements_; }
90 
91 
95  vcl_size_t size1() const { return elements_.size() / 2; }
96 
100  vcl_size_t size2() const { return elements_.size() / 2; }
101 
107  vcl_size_t internal_size() const { return elements_.internal_size(); }
108 
109 
118  {
119  assert(row_index < size1() && col_index < size2() && bool("Invalid access"));
120 
121  long index = static_cast<long>(col_index) - static_cast<long>(row_index);
122 
123  if (index < 0)
124  index = -index;
125  else if
126  (index > 0) index = 2 * static_cast<long>(size1()) - index;
127  return elements_[vcl_size_t(index)];
128  }
129 
130 
138  {
139  elements_ += that.elements();
140  return *this;
141  }
142 
143 private:
144  toeplitz_matrix(toeplitz_matrix const &) {}
145  toeplitz_matrix & operator=(toeplitz_matrix const & t);
146 
147 
149 };
150 
157 template<typename NumericT, unsigned int AlignmentV>
158 void copy(std::vector<NumericT> const & cpu_vec, toeplitz_matrix<NumericT, AlignmentV>& gpu_mat)
159 {
160  assert( (gpu_mat.size1() == 0 || (gpu_mat.size1() * 2 - 1) == cpu_vec.size()) && bool("Size mismatch"));
161 
162  vcl_size_t size = gpu_mat.size1();
163  std::vector<NumericT> rvrs(cpu_vec.size());
164  std::copy(cpu_vec.begin(), cpu_vec.end(), rvrs.begin());
165  std::reverse(rvrs.begin(), rvrs.end());
166 
167  std::vector<NumericT> tmp(size * 2);
168  typedef typename std::vector<NumericT>::difference_type difference_type;
169  std::copy(rvrs.begin() + difference_type(size) - 1, rvrs.end(), tmp.begin());
170  std::copy(rvrs.begin(), rvrs.begin() + difference_type(size) - 1, tmp.begin() + difference_type(size) + 1);
171  tmp[size] = 0.0;
172  copy(tmp, gpu_mat.elements());
173 }
174 
181 template<typename NumericT, unsigned int AlignmentV>
182 void copy(toeplitz_matrix<NumericT, AlignmentV> const & gpu_mat, std::vector<NumericT> & cpu_vec)
183 {
184  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
185 
186  vcl_size_t size = gpu_mat.size1();
187  std::vector<NumericT> tmp(size * 2);
188  copy(gpu_mat.elements(), tmp);
189  std::reverse(tmp.begin(), tmp.end());
190 
191  typedef typename std::vector<NumericT>::difference_type difference_type;
192  std::copy(tmp.begin(), tmp.begin() + difference_type(size) - 1, cpu_vec.begin() + difference_type(size));
193  std::copy(tmp.begin() + difference_type(size), tmp.end(), cpu_vec.begin());
194 
195 }
196 
203 template<typename NumericT, unsigned int AlignmentV, typename MatrixT>
204 void copy(toeplitz_matrix<NumericT, AlignmentV> const & tep_src, MatrixT & com_dst)
205 {
206  assert(tep_src.size1() == viennacl::traits::size1(com_dst) && bool("Size mismatch"));
207  assert(tep_src.size2() == viennacl::traits::size2(com_dst) && bool("Size mismatch"));
208 
209  vcl_size_t size = tep_src.size1();
210  std::vector<NumericT> tmp(tep_src.size1() * 2 - 1);
211  copy(tep_src, tmp);
212 
213  for (vcl_size_t i = 0; i < size; i++)
214  for (vcl_size_t j = 0; j < size; j++)
215  com_dst(i, j) = tmp[static_cast<vcl_size_t>(static_cast<int>(j) - static_cast<int>(i) + static_cast<int>(size) - 1)];
216 }
217 
224 template<typename NumericT, unsigned int AlignmentV, typename MatrixT>
225 void copy(MatrixT const & com_src, toeplitz_matrix<NumericT, AlignmentV>& tep_dst)
226 {
227  assert( (tep_dst.size1() == 0 || tep_dst.size1() == viennacl::traits::size1(com_src)) && bool("Size mismatch"));
228  assert( (tep_dst.size2() == 0 || tep_dst.size2() == viennacl::traits::size2(com_src)) && bool("Size mismatch"));
229 
230  vcl_size_t size = tep_dst.size1();
231 
232  std::vector<NumericT> tmp(2*size - 1);
233 
234  for (long i = static_cast<long>(size) - 1; i >= 0; i--)
235  tmp[size - static_cast<vcl_size_t>(i) - 1] = com_src(static_cast<vcl_size_t>(i), 0);
236 
237  for (vcl_size_t i = 1; i < size; i++)
238  tmp[size + i - 1] = com_src(0, i);
239 
240  copy(tmp, tep_dst);
241 }
242 
243 /*template<typename NumericT, unsigned int AlignmentV, unsigned int VECTOR_AlignmentV>
244  void prod_impl(toeplitz_matrix<NumericT, AlignmentV>& mat,
245  vector<NumericT, VECTOR_AlignmentV>& vec,
246  vector<NumericT, VECTOR_AlignmentV>& result) {
247  viennacl::vector<NumericT, VECTOR_AlignmentV> tep(mat.elements().size() * 2);
248  fft::real_to_complex(mat.elements(), tep, mat.elements().size());
249 
250  viennacl::vector<NumericT, VECTOR_AlignmentV> tmp(vec.size() * 4);
251  viennacl::vector<NumericT, VECTOR_AlignmentV> tmp2(vec.size() * 4);
252 
253  tmp.clear();
254  copy(vec, tmp);
255  fft::real_to_complex(tmp, tmp2, vec.size() * 2);
256  fft::convolve(tep, tmp2, tmp);
257  fft::complex_to_real(tmp, tmp2, vec.size() * 2);
258  copy(tmp2.begin(), tmp2.begin() + vec.size(), result.begin());
259  }*/
260 
266 template<class NumericT, unsigned int AlignmentV>
267 std::ostream & operator<<(std::ostream & s, toeplitz_matrix<NumericT, AlignmentV>& gpu_matrix)
268 {
269  vcl_size_t size = gpu_matrix.size1();
270  std::vector<NumericT> tmp(2*size - 1);
271  copy(gpu_matrix, tmp);
272  s << "[" << size << "," << size << "](";
273 
274  for (vcl_size_t i = 0; i < size; i++)
275  {
276  s << "(";
277  for (vcl_size_t j = 0; j < size; j++)
278  {
279  s << tmp[vcl_size_t(static_cast<int>(j) - static_cast<int>(i) + static_cast<int>(size - 1))];
280  //s << (int)i - (int)j;
281  if (j < (size - 1)) s << ",";
282  }
283  s << ")";
284  }
285  s << ")";
286  return s;
287 }
288 
289 //
290 // Specify available operations:
291 //
292 
295 namespace linalg
296 {
297 namespace detail
298 {
299  // x = A * y
300  template<typename T, unsigned int A>
301  struct op_executor<vector_base<T>, op_assign, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
302  {
303  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
304  {
305  // check for the special case x = A * x
306  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
307  {
308  viennacl::vector<T> temp(lhs);
309  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
310  lhs = temp;
311  }
312  else
313  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
314  }
315  };
316 
317  template<typename T, unsigned int A>
318  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
319  {
320  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
321  {
322  viennacl::vector<T> temp(lhs);
323  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
324  lhs += temp;
325  }
326  };
327 
328  template<typename T, unsigned int A>
329  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> >
330  {
331  static void apply(vector_base<T> & lhs, vector_expression<const toeplitz_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
332  {
333  viennacl::vector<T> temp(lhs);
334  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
335  lhs -= temp;
336  }
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_assign, vector_expression<const toeplitz_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 toeplitz_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::linalg::prod_impl(rhs.lhs(), temp, lhs);
348  }
349  };
350 
351  // x = A * vec_op
352  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
353  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> >
354  {
355  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)
356  {
357  viennacl::vector<T> temp(rhs.rhs());
358  viennacl::vector<T> temp_result(lhs);
359  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
360  lhs += temp_result;
361  }
362  };
363 
364  // x = A * vec_op
365  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
366  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> >
367  {
368  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)
369  {
370  viennacl::vector<T> temp(rhs.rhs());
371  viennacl::vector<T> temp_result(lhs);
372  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
373  lhs -= temp_result;
374  }
375  };
376 
377 } // namespace detail
378 } // namespace linalg
379 
382 }
383 
384 #endif // VIENNACL_TOEPLITZ_MATRIX_HPP
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
Implementations of operations using toeplitz_matrix. Experimental.
handle_type const & handle() const
Returns the OpenCL handle.
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
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:201
vcl_size_t internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
toeplitz_matrix< NumericT, AlignmentV > & operator+=(toeplitz_matrix< NumericT, AlignmentV > &that)
+= operation for Toeplitz matrices
toeplitz_matrix()
The default constructor. Does not allocate any memory.
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
void resize(vcl_size_t sz, bool preserve=true)
Resizes the matrix. Existing entries can be preserved.
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:235
Definition: blas3.hpp:36
vcl_size_t size1() const
Returns the number of rows of the matrix.
viennacl::backend::mem_handle handle_type
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< NumericT >::ResultType > value_type
vcl_size_t size2() const
Returns the number of columns of the matrix.
std::size_t vcl_size_t
Definition: forwards.h:75
Implementations of the OpenCL backend, where all contexts are stored in.
void copy(MatrixT const &com_src, toeplitz_matrix< NumericT, AlignmentV > &tep_dst)
Copies a the matrix-like object to the Toeplitz matrix from the OpenCL device (either GPU or multi-co...
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) ...
viennacl::vector< NumericT, AlignmentV > & elements()
Returns an internal viennacl::vector, which represents a Toeplitz matrix elements.
viennacl::vector< NumericT, AlignmentV > const & elements() const
entry_proxy< NumericT > operator()(vcl_size_t row_index, vcl_size_t col_index)
Read-write access to a single element of the matrix.
All routines related to the Fast Fourier Transform. Experimental.
A Toeplitz matrix class.
Definition: forwards.h:415
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
void reverse(viennacl::vector_base< NumericT > &in)
Reverse vector to oposite order and save it in input vector.
void prod_impl(const matrix_base< NumericT > &mat, const vector_base< NumericT > &vec, vector_base< NumericT > &result)
Carries out matrix-vector multiplication.
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
toeplitz_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
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