ViennaCL - The Vienna Computing Library  1.5.2
hankel_matrix.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_HANKEL_MATRIX_HPP
2 #define VIENNACL_HANKEL_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 
21 
26 #include "viennacl/forwards.h"
27 #include "viennacl/vector.hpp"
28 #include "viennacl/ocl/backend.hpp"
29 
31 #include "viennacl/fft.hpp"
32 
34 
35 namespace viennacl {
41  template<class SCALARTYPE, unsigned int ALIGNMENT>
42  class hankel_matrix
43  {
44  public:
47 
52  explicit hankel_matrix() {}
53 
60  explicit hankel_matrix(vcl_size_t rows, vcl_size_t cols) : elements_(rows, cols)
61  {
62  assert(rows == cols && bool("Hankel matrix must be square!"));
63  (void)cols; // avoid 'unused parameter' warning in optimized builds
64  }
65 
72  void resize(vcl_size_t sz, bool preserve = true)
73  {
74  elements_.resize(sz, preserve);
75  }
76 
81  handle_type const & handle() const { return elements_.handle(); }
82 
88  toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
89 
93  vcl_size_t size1() const { return elements_.size1(); }
94 
98  vcl_size_t size2() const { return elements_.size2(); }
99 
105  vcl_size_t internal_size() const { return elements_.internal_size(); }
106 
114  entry_proxy<SCALARTYPE> operator()(unsigned int row_index, unsigned int col_index)
115  {
116  assert(row_index < size1() && col_index < size2() && bool("Invalid access"));
117 
118  return elements_(size1() - row_index - 1, col_index);
119  }
120 
128  {
129  elements_ += that.elements();
130  return *this;
131  }
132 
133  private:
134  hankel_matrix(hankel_matrix const &) {}
135  hankel_matrix & operator=(hankel_matrix const & t);
136 
137  toeplitz_matrix<SCALARTYPE, ALIGNMENT> elements_;
138  };
139 
146  template <typename SCALARTYPE, unsigned int ALIGNMENT>
147  void copy(std::vector<SCALARTYPE> const & cpu_vec, hankel_matrix<SCALARTYPE, ALIGNMENT> & gpu_mat)
148  {
149  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
150 
151  copy(cpu_vec, gpu_mat.elements());
152  }
153 
160  template <typename SCALARTYPE, unsigned int ALIGNMENT>
161  void copy(hankel_matrix<SCALARTYPE, ALIGNMENT> const & gpu_mat, std::vector<SCALARTYPE> & cpu_vec)
162  {
163  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
164 
165  copy(gpu_mat.elements(), cpu_vec);
166  }
167 
174  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
175  void copy(hankel_matrix<SCALARTYPE, ALIGNMENT> const & han_src, MATRIXTYPE& com_dst)
176  {
177  assert( (viennacl::traits::size1(com_dst) == han_src.size1()) && bool("Size mismatch") );
178  assert( (viennacl::traits::size2(com_dst) == han_src.size2()) && bool("Size mismatch") );
179 
180  vcl_size_t size = han_src.size1();
181  std::vector<SCALARTYPE> tmp(size * 2 - 1);
182  copy(han_src, tmp);
183 
184  for (vcl_size_t i = 0; i < size; i++)
185  for (vcl_size_t j = 0; j < size; j++)
186  com_dst(i, j) = tmp[i + j];
187  }
188 
195  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
196  void copy(MATRIXTYPE const & com_src, hankel_matrix<SCALARTYPE, ALIGNMENT>& han_dst)
197  {
198  assert( (han_dst.size1() == 0 || viennacl::traits::size1(com_src) == han_dst.size1()) && bool("Size mismatch") );
199  assert( (han_dst.size2() == 0 || viennacl::traits::size2(com_src) == han_dst.size2()) && bool("Size mismatch") );
200  assert( viennacl::traits::size2(com_src) == viennacl::traits::size1(com_src) && bool("Logic error: non-square Hankel matrix!") );
201 
203 
204  std::vector<SCALARTYPE> tmp(2*size - 1);
205 
206  for (vcl_size_t i = 0; i < size; i++)
207  tmp[i] = com_src(0, i);
208 
209  for (vcl_size_t i = 1; i < size; i++)
210  tmp[size + i - 1] = com_src(size - 1, i);
211 
212  viennacl::copy(tmp, han_dst);
213  }
214 
215  /*template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
216  void prod_impl(hankel_matrix<SCALARTYPE, ALIGNMENT>& mat,
217  vector<SCALARTYPE, VECTOR_ALIGNMENT>& vec,
218  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result)
219  {
220  prod_impl(mat.elements(), vec, result);
221  fft::reverse(result);
222  }*/
223 
224  template<class SCALARTYPE, unsigned int ALIGNMENT>
225  std::ostream & operator<<(std::ostream & s, hankel_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
226  {
227  vcl_size_t size = gpu_matrix.size1();
228  std::vector<SCALARTYPE> tmp(2*size - 1);
229  copy(gpu_matrix, tmp);
230  s << "[" << size << "," << size << "](";
231 
232  for(vcl_size_t i = 0; i < size; i++) {
233  s << "(";
234  for(vcl_size_t j = 0; j < size; j++) {
235  s << tmp[i + j];
236  //s << (int)i - (int)j;
237  if(j < (size - 1)) s << ",";
238  }
239  s << ")";
240  }
241  s << ")";
242  return s;
243  }
244 
245  //
246  // Specify available operations:
247  //
248 
251  namespace linalg
252  {
253  namespace detail
254  {
255  // x = A * y
256  template <typename T, unsigned int A>
257  struct op_executor<vector_base<T>, op_assign, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
258  {
259  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
260  {
261  // check for the special case x = A * x
262  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
263  {
264  viennacl::vector<T> temp(lhs);
265  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
266  lhs = temp;
267  }
268  else
269  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
270  }
271  };
272 
273  template <typename T, unsigned int A>
274  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
275  {
276  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
277  {
278  viennacl::vector<T> temp(lhs);
279  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
280  lhs += temp;
281  }
282  };
283 
284  template <typename T, unsigned int A>
285  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
286  {
287  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
288  {
289  viennacl::vector<T> temp(lhs);
290  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
291  lhs -= temp;
292  }
293  };
294 
295 
296  // x = A * vec_op
297  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
298  struct op_executor<vector_base<T>, op_assign, vector_expression<const hankel_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
299  {
300  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
301  {
302  viennacl::vector<T> temp(rhs.rhs());
303  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
304  }
305  };
306 
307  // x = A * vec_op
308  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
309  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const hankel_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> >
310  {
311  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
312  {
313  viennacl::vector<T> temp(rhs.rhs());
314  viennacl::vector<T> temp_result(lhs);
315  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
316  lhs += temp_result;
317  }
318  };
319 
320  // x = A * vec_op
321  template <typename T, unsigned int A, typename LHS, typename RHS, typename OP>
322  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const hankel_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> >
323  {
324  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_expression<const LHS, const RHS, OP>, op_prod> const & rhs)
325  {
326  viennacl::vector<T> temp(rhs.rhs());
327  viennacl::vector<T> temp_result(lhs);
328  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
329  lhs -= temp_result;
330  }
331  };
332 
333 
334 
335  } // namespace detail
336  } // namespace linalg
337 
339 }
340 #endif // VIENNACL_HANKEL_MATRIX_HPP
vcl_size_t internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
Definition: hankel_matrix.hpp:105
Implementations of operations using hankel_matrix. Experimental.
hankel_matrix< SCALARTYPE, ALIGNMENT > & operator+=(hankel_matrix< SCALARTYPE, ALIGNMENT > &that)
+= operation for Hankel matrices
Definition: hankel_matrix.hpp:127
std::size_t vcl_size_t
Definition: forwards.h:58
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:172
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
entry_proxy< SCALARTYPE > operator()(unsigned int row_index, unsigned int col_index)
Read-write access to a element of the matrix.
Definition: hankel_matrix.hpp:114
This file provides the forward declarations for the main types used within ViennaCL.
A Hankel matrix class.
Definition: forwards.h:327
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
hankel_matrix()
The default constructor. Does not allocate any memory.
Definition: hankel_matrix.hpp:52
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
viennacl::backend::mem_handle handle_type
Definition: hankel_matrix.hpp:45
handle_type const & handle() const
Returns the OpenCL handle.
Definition: hankel_matrix.hpp:81
vcl_size_t size1() const
Returns the number of rows of the matrix.
Definition: hankel_matrix.hpp:93
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
toeplitz_matrix< SCALARTYPE, ALIGNMENT > & elements()
Returns an internal viennacl::toeplitz_matrix, which represents a Hankel matrix elements.
Definition: hankel_matrix.hpp:87
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
void resize(vcl_size_t sz, bool preserve=true)
Resizes the matrix. Existing entries can be preserved.
Definition: hankel_matrix.hpp:72
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< SCALARTYPE >::ResultType > value_type
Definition: hankel_matrix.hpp:46
Implementations of the OpenCL backend, where all contexts are stored in.
toeplitz_matrix< SCALARTYPE, ALIGNMENT > const & elements() const
Definition: hankel_matrix.hpp:88
hankel_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
Definition: hankel_matrix.hpp:60
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
All routines related to the Fast Fourier Transform. Experimental.
A Toeplitz matrix class.
Definition: forwards.h:330
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
vcl_size_t size2() const
Returns the number of columns of the matrix.
Definition: hankel_matrix.hpp:98
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:178
Implementation of the toeplitz_matrix class for efficient manipulation of Toeplitz matrices...