ViennaCL - The Vienna Computing Library  1.2.0
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-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 
24 #include "viennacl/forwards.h"
25 #include "viennacl/vector.hpp"
26 #include "viennacl/ocl/context.hpp"
27 
28 #include "viennacl/fft.hpp"
29 
31 
32 
33 namespace viennacl {
34 
40  template<class SCALARTYPE, unsigned int ALIGNMENT>
42  {
43  public:
44 
49  explicit toeplitz_matrix()
50  {
51  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
52  }
53 
59  explicit toeplitz_matrix(std::size_t rows, std::size_t cols) : elements_(rows * 2)
60  {
61  assert(rows == cols && "Toeplitz matrix must be square!");
62  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
63  }
64 
65 
72  void resize(size_t sz, bool preserve = true) {
73  elements_.resize(sz * 2, preserve);
74  }
75 
80  viennacl::ocl::handle<cl_mem> handle() const { return elements_.handle(); }
81 
87  viennacl::vector<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
88 
89 
93  std::size_t size1() const { return elements_.size() / 2; }
94 
98  std::size_t size2() const { return elements_.size() / 2; }
99 
105  std::size_t internal_size() const { return elements_.internal_size(); }
106 
107 
115  entry_proxy<SCALARTYPE> operator()(std::size_t row_index, std::size_t col_index)
116  {
117  assert(row_index < size1() && col_index < size2() && "Invalid access");
118 
119  int index = static_cast<int>(col_index) - static_cast<int>(row_index);
120 
121  if (index < 0)
122  index = -index;
123  else if
124  (index > 0) index = 2 * 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 & t) {}
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  std::size_t size = gpu_mat.size1();
158  assert((size * 2 - 1) == cpu_vec.size() && "Size mismatch");
159  std::vector<SCALARTYPE> rvrs(cpu_vec.size());
160  std::copy(cpu_vec.begin(), cpu_vec.end(), rvrs.begin());
161  std::reverse(rvrs.begin(), rvrs.end());
162 
163  std::vector<SCALARTYPE> tmp(size * 2);
164  std::copy(rvrs.begin() + size - 1, rvrs.end(), tmp.begin());
165  std::copy(rvrs.begin(), rvrs.begin() + size - 1, tmp.begin() + size + 1);
166  tmp[size] = 0.0;
167  copy(tmp, gpu_mat.elements());
168  }
169 
176  template <typename SCALARTYPE, unsigned int ALIGNMENT>
177  void copy(toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & gpu_mat, std::vector<SCALARTYPE> & cpu_vec)
178  {
179  std::size_t size = gpu_mat.size1();
180  assert((size * 2 - 1) == cpu_vec.size() && "Size mismatch");
181  std::vector<SCALARTYPE> tmp(size * 2);
182  copy(gpu_mat.elements(), tmp);
183  std::reverse(tmp.begin(), tmp.end());
184 
185  std::copy(tmp.begin(), tmp.begin() + size - 1, cpu_vec.begin() + size);
186  std::copy(tmp.begin() + size, tmp.end(), cpu_vec.begin());
187 
188  }
189 
196  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
197  void copy(toeplitz_matrix<SCALARTYPE, ALIGNMENT> const & tep_src, MATRIXTYPE & com_dst)
198  {
199  std::size_t size = tep_src.size1();
200  assert(size == com_dst.size1() && "Size mismatch");
201  assert(size == com_dst.size2() && "Size mismatch");
202  std::vector<SCALARTYPE> tmp(tep_src.size1() * 2 - 1);
203  copy(tep_src, tmp);
204 
205  for(std::size_t i = 0; i < size; i++)
206  for(std::size_t j = 0; j < size; j++)
207  com_dst(i, j) = tmp[static_cast<int>(j) - static_cast<int>(i) + static_cast<int>(size) - 1];
208  }
209 
216  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
217  void copy(MATRIXTYPE const & com_src, toeplitz_matrix<SCALARTYPE, ALIGNMENT>& tep_dst)
218  {
219  std::size_t size = tep_dst.size1();
220  assert(size == com_src.size1() && "Size mismatch");
221  assert(size == com_src.size2() && "Size mismatch");
222 
223  std::vector<SCALARTYPE> tmp(2*size - 1);
224 
225  for(int i = size - 1; i >= 0; i--)
226  tmp[size - i - 1] = com_src(i, 0);
227 
228  for(std::size_t i = 1; i < size; i++)
229  tmp[size + i - 1] = com_src(0, i);
230 
231  copy(tmp, tep_dst);
232  }
233 
234  /*template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
235  void prod_impl(toeplitz_matrix<SCALARTYPE, ALIGNMENT>& mat,
236  vector<SCALARTYPE, VECTOR_ALIGNMENT>& vec,
237  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result) {
238  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tep(mat.elements().size() * 2);
239  fft::real_to_complex(mat.elements(), tep, mat.elements().size());
240 
241  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp(vec.size() * 4);
242  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp2(vec.size() * 4);
243 
244  tmp.clear();
245  copy(vec, tmp);
246  fft::real_to_complex(tmp, tmp2, vec.size() * 2);
247  fft::convolve(tep, tmp2, tmp);
248  fft::complex_to_real(tmp, tmp2, vec.size() * 2);
249  copy(tmp2.begin(), tmp2.begin() + vec.size(), result.begin());
250  }*/
251 
257  template<class SCALARTYPE, unsigned int ALIGNMENT>
258  std::ostream & operator<<(std::ostream & s, toeplitz_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
259  {
260  std::size_t size = gpu_matrix.size1();
261  std::vector<SCALARTYPE> tmp(2*size - 1);
262  copy(gpu_matrix, tmp);
263  s << "[" << size << "," << size << "](";
264 
265  for(std::size_t i = 0; i < size; i++) {
266  s << "(";
267  for(std::size_t j = 0; j < size; j++) {
268  s << tmp[(int)j - (int)i + (int)size - 1];
269  //s << (int)i - (int)j;
270  if(j < (size - 1)) s << ",";
271  }
272  s << ")";
273  }
274  s << ")";
275  return s;
276  }
277 
278 }
279 
280 #endif // _VIENNACL_TOEPLITZ_MATRIX_HPP