ViennaCL - The Vienna Computing Library  1.2.0
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-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
20 
25 #include "viennacl/forwards.h"
26 #include "viennacl/vector.hpp"
27 #include "viennacl/ocl/context.hpp"
28 
30 #include "viennacl/fft.hpp"
31 
33 
34 namespace viennacl {
40  template<class SCALARTYPE, unsigned int ALIGNMENT>
42  {
43  public:
48  explicit hankel_matrix()
49  {
50  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
51  }
52 
59  explicit hankel_matrix(std::size_t rows, std::size_t cols) : elements_(rows, cols)
60  {
61  assert(rows == cols && "Hankel matrix must be square!");
62  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
63  }
64 
71  void resize(size_t sz, bool preserve = true)
72  {
73  elements_.resize(sz, preserve);
74  }
75 
80  viennacl::ocl::handle<cl_mem> handle() const { return elements_.handle(); }
81 
87  toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
88 
92  std::size_t size1() const { return elements_.size1(); }
93 
97  std::size_t size2() const { return elements_.size2(); }
98 
104  std::size_t internal_size() const { return elements_.internal_size(); }
105 
113  entry_proxy<SCALARTYPE> operator()(unsigned int row_index, unsigned int col_index)
114  {
115  assert(row_index < size1() && col_index < size2() && "Invalid access");
116 
117  return elements_(size1() - row_index - 1, col_index);
118  }
119 
127  {
128  elements_ += that.elements();
129  return *this;
130  }
131 
132  private:
133  hankel_matrix(hankel_matrix const & t) {}
134  hankel_matrix & operator=(hankel_matrix const & t) {}
135 
136  toeplitz_matrix<SCALARTYPE, ALIGNMENT> elements_;
137  };
138 
145  template <typename SCALARTYPE, unsigned int ALIGNMENT>
146  void copy(std::vector<SCALARTYPE> const & cpu_vec, hankel_matrix<SCALARTYPE, ALIGNMENT> & gpu_mat)
147  {
148  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && "Size mismatch");
149 
150  copy(cpu_vec, gpu_mat.elements());
151  }
152 
159  template <typename SCALARTYPE, unsigned int ALIGNMENT>
160  void copy(hankel_matrix<SCALARTYPE, ALIGNMENT> const & gpu_mat, std::vector<SCALARTYPE> & cpu_vec)
161  {
162  assert((gpu_mat.size1() * 2 - 1) == cpu_vec.size() && "Size mismatch");
163 
164  copy(gpu_mat.elements(), cpu_vec);
165  }
166 
173  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
174  void copy(hankel_matrix<SCALARTYPE, ALIGNMENT> const & han_src, MATRIXTYPE& com_dst)
175  {
176  std::size_t size = han_src.size1();
177  assert(size == com_dst.size1() && "Size mismatch");
178  assert(size == com_dst.size2() && "Size mismatch");
179  std::vector<SCALARTYPE> tmp(size * 2 - 1);
180  copy(han_src, tmp);
181 
182  for (std::size_t i = 0; i < size; i++)
183  for (std::size_t j = 0; j < size; j++)
184  com_dst(i, j) = tmp[i + j];
185  }
186 
193  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
194  void copy(MATRIXTYPE const & com_src, hankel_matrix<SCALARTYPE, ALIGNMENT>& han_dst)
195  {
196  std::size_t size = han_dst.size1();
197  assert(size == com_src.size1() && "Size mismatch");
198  assert(size == com_src.size2() && "Size mismatch");
199 
200  std::vector<SCALARTYPE> tmp(2*size - 1);
201 
202  for (std::size_t i = 0; i < size; i++)
203  tmp[i] = com_src(0, i);
204 
205  for (std::size_t i = 1; i < size; i++)
206  tmp[size + i - 1] = com_src(size - 1, i);
207 
208  viennacl::copy(tmp, han_dst);
209  }
210 
211  /*template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
212  void prod_impl(hankel_matrix<SCALARTYPE, ALIGNMENT>& mat,
213  vector<SCALARTYPE, VECTOR_ALIGNMENT>& vec,
214  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result)
215  {
216  prod_impl(mat.elements(), vec, result);
217  fft::reverse(result);
218  }*/
219 
220  template<class SCALARTYPE, unsigned int ALIGNMENT>
221  std::ostream & operator<<(std::ostream & s, hankel_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
222  {
223  std::size_t size = gpu_matrix.size1();
224  std::vector<SCALARTYPE> tmp(2*size - 1);
225  copy(gpu_matrix, tmp);
226  s << "[" << size << "," << size << "](";
227 
228  for(std::size_t i = 0; i < size; i++) {
229  s << "(";
230  for(std::size_t j = 0; j < size; j++) {
231  s << tmp[i + j];
232  //s << (int)i - (int)j;
233  if(j < (size - 1)) s << ",";
234  }
235  s << ")";
236  }
237  s << ")";
238  return s;
239  }
240 }
241 #endif // _VIENNACL_HANKEL_MATRIX_HPP