ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
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-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 
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
36 {
42 template<class NumericT, unsigned int AlignmentV>
43 class hankel_matrix
44 {
45 public:
48 
53  explicit hankel_matrix() {}
54 
61  explicit hankel_matrix(vcl_size_t rows, vcl_size_t cols) : elements_(rows, cols)
62  {
63  assert(rows == cols && bool("Hankel matrix must be square!"));
64  (void)cols; // avoid 'unused parameter' warning in optimized builds
65  }
66 
73  void resize(vcl_size_t sz, bool preserve = true)
74  {
75  elements_.resize(sz, preserve);
76  }
77 
82  handle_type const & handle() const { return elements_.handle(); }
83 
89  toeplitz_matrix<NumericT, AlignmentV> const & elements() const { return elements_; }
90 
94  vcl_size_t size1() const { return elements_.size1(); }
95 
99  vcl_size_t size2() const { return elements_.size2(); }
100 
106  vcl_size_t internal_size() const { return elements_.internal_size(); }
107 
115  entry_proxy<NumericT> operator()(unsigned int row_index, unsigned int col_index)
116  {
117  assert(row_index < size1() && col_index < size2() && bool("Invalid access"));
118 
119  return elements_(size1() - row_index - 1, col_index);
120  }
121 
129  {
130  elements_ += that.elements();
131  return *this;
132  }
133 
134 private:
135  hankel_matrix(hankel_matrix const &) {}
136  hankel_matrix & operator=(hankel_matrix const & t);
137 
138  toeplitz_matrix<NumericT, AlignmentV> elements_;
139 };
140 
147 template<typename NumericT, unsigned int AlignmentV>
148 void copy(std::vector<NumericT> const & cpu_vec, hankel_matrix<NumericT, AlignmentV> & gpu_mat)
149 {
150  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
151 
152  copy(cpu_vec, gpu_mat.elements());
153 }
154 
161 template<typename NumericT, unsigned int AlignmentV>
162 void copy(hankel_matrix<NumericT, AlignmentV> const & gpu_mat, std::vector<NumericT> & cpu_vec)
163 {
164  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && bool("Size mismatch"));
165 
166  copy(gpu_mat.elements(), cpu_vec);
167 }
168 
175 template<typename NumericT, unsigned int AlignmentV, typename MatrixT>
176 void copy(hankel_matrix<NumericT, AlignmentV> const & han_src, MatrixT& com_dst)
177 {
178  assert( (viennacl::traits::size1(com_dst) == han_src.size1()) && bool("Size mismatch") );
179  assert( (viennacl::traits::size2(com_dst) == han_src.size2()) && bool("Size mismatch") );
180 
181  vcl_size_t size = han_src.size1();
182  std::vector<NumericT> tmp(size * 2 - 1);
183  copy(han_src, tmp);
184 
185  for (vcl_size_t i = 0; i < size; i++)
186  for (vcl_size_t j = 0; j < size; j++)
187  com_dst(i, j) = tmp[i + j];
188 }
189 
196 template<typename NumericT, unsigned int AlignmentV, typename MatrixT>
197 void copy(MatrixT const & com_src, hankel_matrix<NumericT, AlignmentV>& han_dst)
198 {
199  assert( (han_dst.size1() == 0 || viennacl::traits::size1(com_src) == han_dst.size1()) && bool("Size mismatch") );
200  assert( (han_dst.size2() == 0 || viennacl::traits::size2(com_src) == han_dst.size2()) && bool("Size mismatch") );
201  assert( viennacl::traits::size2(com_src) == viennacl::traits::size1(com_src) && bool("Logic error: non-square Hankel matrix!") );
202 
204 
205  std::vector<NumericT> tmp(2*size - 1);
206 
207  for (vcl_size_t i = 0; i < size; i++)
208  tmp[i] = com_src(0, i);
209 
210  for (vcl_size_t i = 1; i < size; i++)
211  tmp[size + i - 1] = com_src(size - 1, i);
212 
213  viennacl::copy(tmp, han_dst);
214 }
215 
216 /*template<typename NumericT, unsigned int AlignmentV, unsigned int VECTOR_AlignmentV>
217  void prod_impl(hankel_matrix<NumericT, AlignmentV>& mat,
218  vector<NumericT, VECTOR_AlignmentV>& vec,
219  vector<NumericT, VECTOR_AlignmentV>& result)
220  {
221  prod_impl(mat.elements(), vec, result);
222  fft::reverse(result);
223  }*/
224 
225 template<class NumericT, unsigned int AlignmentV>
226 std::ostream & operator<<(std::ostream & s, hankel_matrix<NumericT, AlignmentV>& gpu_matrix)
227 {
228  vcl_size_t size = gpu_matrix.size1();
229  std::vector<NumericT> tmp(2*size - 1);
230  copy(gpu_matrix, tmp);
231  s << "[" << size << "," << size << "](";
232 
233  for (vcl_size_t i = 0; i < size; i++)
234  {
235  s << "(";
236  for (vcl_size_t j = 0; j < size; j++)
237  {
238  s << tmp[i + j];
239  //s << (int)i - (int)j;
240  if (j < (size - 1)) s << ",";
241  }
242  s << ")";
243  }
244  s << ")";
245  return s;
246 }
247 
248 //
249 // Specify available operations:
250 //
251 
254 namespace linalg
255 {
256 namespace detail
257 {
258  // x = A * y
259  template<typename T, unsigned int A>
260  struct op_executor<vector_base<T>, op_assign, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
261  {
262  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
263  {
264  // check for the special case x = A * x
265  if (viennacl::traits::handle(lhs) == viennacl::traits::handle(rhs.rhs()))
266  {
267  viennacl::vector<T> temp(lhs);
268  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
269  lhs = temp;
270  }
271  else
272  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), lhs);
273  }
274  };
275 
276  template<typename T, unsigned int A>
277  struct op_executor<vector_base<T>, op_inplace_add, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
278  {
279  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
280  {
281  viennacl::vector<T> temp(lhs);
282  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
283  lhs += temp;
284  }
285  };
286 
287  template<typename T, unsigned int A>
288  struct op_executor<vector_base<T>, op_inplace_sub, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> >
289  {
290  static void apply(vector_base<T> & lhs, vector_expression<const hankel_matrix<T, A>, const vector_base<T>, op_prod> const & rhs)
291  {
292  viennacl::vector<T> temp(lhs);
293  viennacl::linalg::prod_impl(rhs.lhs(), rhs.rhs(), temp);
294  lhs -= temp;
295  }
296  };
297 
298 
299  // x = A * vec_op
300  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
301  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> >
302  {
303  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)
304  {
305  viennacl::vector<T> temp(rhs.rhs());
306  viennacl::linalg::prod_impl(rhs.lhs(), temp, lhs);
307  }
308  };
309 
310  // x = A * vec_op
311  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
312  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> >
313  {
314  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)
315  {
316  viennacl::vector<T> temp(rhs.rhs());
317  viennacl::vector<T> temp_result(lhs);
318  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
319  lhs += temp_result;
320  }
321  };
322 
323  // x = A * vec_op
324  template<typename T, unsigned int A, typename LHS, typename RHS, typename OP>
325  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> >
326  {
327  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)
328  {
329  viennacl::vector<T> temp(rhs.rhs());
330  viennacl::vector<T> temp_result(lhs);
331  viennacl::linalg::prod_impl(rhs.lhs(), temp, temp_result);
332  lhs -= temp_result;
333  }
334  };
335 
336 
337 
338 } // namespace detail
339 } // namespace linalg
340 
342 }
343 #endif // VIENNACL_HANKEL_MATRIX_HPP
Implementations of operations using hankel_matrix. Experimental.
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
scalar< typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT< NumericT >::ResultType > value_type
hankel_matrix(vcl_size_t rows, vcl_size_t cols)
Creates the matrix with the given size.
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
vcl_size_t size2() const
Returns the number of columns of the matrix.
This file provides the forward declarations for the main types used within ViennaCL.
A Hankel matrix class.
Definition: forwards.h:412
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
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
toeplitz_matrix< NumericT, AlignmentV > & elements()
Returns an internal viennacl::toeplitz_matrix, which represents a Hankel matrix elements.
entry_proxy< NumericT > operator()(unsigned int row_index, unsigned int col_index)
Read-write access to a element of the matrix.
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:235
handle_type const & handle() const
Returns the OpenCL handle.
Definition: blas3.hpp:36
hankel_matrix()
The default constructor. Does not allocate any memory.
viennacl::backend::mem_handle handle_type
vcl_size_t size1() const
Returns the number of rows of the matrix.
std::size_t vcl_size_t
Definition: forwards.h:75
toeplitz_matrix< NumericT, AlignmentV > const & elements() const
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.
hankel_matrix< NumericT, AlignmentV > & operator+=(hankel_matrix< NumericT, AlignmentV > &that)
+= operation for Hankel matrices
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) ...
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
vcl_size_t internal_size() const
Returns the internal size of matrix representtion. Usually required for launching OpenCL kernels only...
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
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
Implementation of the toeplitz_matrix class for efficient manipulation of Toeplitz matrices...