ViennaCL - The Vienna Computing Library  1.2.0
circulant_matrix.hpp
Go to the documentation of this file.
1 #ifndef _VIENNACL_CIRCULANT_MATRIX_HPP
2 #define _VIENNACL_CIRCULANT_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 
29 
30 #include "viennacl/fft.hpp"
31 
32 namespace viennacl
33 {
39  template<class SCALARTYPE, unsigned int ALIGNMENT>
41  {
42  public:
47  explicit circulant_matrix()
48  {
49  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
50  }
51 
58  explicit circulant_matrix(std::size_t rows, std::size_t cols) : elements_(rows)
59  {
60  assert(rows == cols && "Circulant matrix must be square!");
61  viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
62  }
63 
70  void resize(size_t sz, bool preserve = true)
71  {
72  elements_.resize(sz, preserve);
73  }
74 
79  viennacl::ocl::handle<cl_mem> handle() const { return elements_.handle(); }
80 
86  viennacl::vector<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }
87 
91  std::size_t size1() const { return elements_.size(); }
92 
96  std::size_t size2() const { return elements_.size(); }
97 
103  std::size_t internal_size() const { return elements_.internal_size(); }
104 
112  entry_proxy<SCALARTYPE> operator()(std::size_t row_index, std::size_t col_index)
113  {
114  int index = static_cast<int>(row_index) - static_cast<int>(col_index);
115 
116  assert(row_index < size1() && col_index < size2() && "Invalid access");
117 
118  while (index < 0)
119  index += size1();
120  return elements_[index];
121  }
122 
130  {
131  elements_ += that.elements();
132  return *this;
133  }
134 
135  private:
136  circulant_matrix(circulant_matrix const & t) {}
137  circulant_matrix & operator=(circulant_matrix const & t) {}
138 
140  };
141 
148  template <typename SCALARTYPE, unsigned int ALIGNMENT>
149  void copy(std::vector<SCALARTYPE>& cpu_vec, circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat)
150  {
151  assert(cpu_vec.size() == gpu_mat.size1() && "Size mismatch");
152  copy(cpu_vec, gpu_mat.elements());
153  }
154 
161  template <typename SCALARTYPE, unsigned int ALIGNMENT>
162  void copy(circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat, std::vector<SCALARTYPE>& cpu_vec)
163  {
164  assert(cpu_vec.size() == gpu_mat.size1() && "Size mismatch");
165  copy(gpu_mat.elements(), cpu_vec);
166  }
167 
174  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
175  void copy(circulant_matrix<SCALARTYPE, ALIGNMENT>& circ_src, MATRIXTYPE& com_dst) {
176  std::size_t size = circ_src.size1();
177  assert(size == com_dst.size1() && "Size mismatch");
178  assert(size == com_dst.size2() && "Size mismatch");
179  std::vector<SCALARTYPE> tmp(size);
180  copy(circ_src, tmp);
181 
182  for (std::size_t i = 0; i < size; i++) {
183  for (std::size_t j = 0; j < size; j++) {
184  int index = static_cast<int>(i) - static_cast<int>(j);
185  if (index < 0)
186  index = size + index;
187  com_dst(i, j) = tmp[index];
188  }
189  }
190  }
191 
198  template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
199  void copy(MATRIXTYPE& com_src, circulant_matrix<SCALARTYPE, ALIGNMENT>& circ_dst) {
200  std::size_t size = circ_dst.size1();
201  assert(size == com_src.size1() && "Size mismatch");
202  assert(size == com_src.size2() && "Size mismatch");
203 
204  std::vector<SCALARTYPE> tmp(size);
205 
206  for(std::size_t i = 0; i < size; i++) tmp[i] = com_src(i, 0);
207 
208  copy(tmp, circ_dst);
209  }
210 
211  /*namespace linalg
212  {
213  template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
214  void prod_impl(circulant_matrix<SCALARTYPE, ALIGNMENT> const & mat,
215  vector<SCALARTYPE, VECTOR_ALIGNMENT> const & vec,
216  vector<SCALARTYPE, VECTOR_ALIGNMENT>& result) {
217  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> circ(mat.elements().size() * 2);
218  fft::real_to_complex(mat.elements(), circ, mat.elements().size());
219 
220  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp(vec.size() * 2);
221  viennacl::vector<SCALARTYPE, VECTOR_ALIGNMENT> tmp2(vec.size() * 2);
222 
223  fft::real_to_complex(vec, tmp, vec.size());
224  fft::convolve(circ, tmp, tmp2);
225  fft::complex_to_real(tmp2, result, vec.size());
226  }
227  }*/
228 
234  template<class SCALARTYPE, unsigned int ALIGNMENT>
235  std::ostream & operator<<(std::ostream& s, circulant_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
236  {
237  std::size_t size = gpu_matrix.size1();
238  std::vector<SCALARTYPE> tmp(size);
239  copy(gpu_matrix, tmp);
240  s << "[" << size << "," << size << "](";
241 
242  for(std::size_t i = 0; i < size; i++) {
243  s << "(";
244  for(std::size_t j = 0; j < size; j++) {
245  int index = (int)i - (int)j;
246  if(index < 0) index = size + index;
247  s << tmp[index];
248  //s << index;
249  if(j < (size - 1)) s << ",";
250  }
251  s << ")";
252  }
253  s << ")";
254  return s;
255  }
256 }
257 
258 #endif // _VIENNACL_CIRCULANT_MATRIX_HPP