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