ViennaCL - The Vienna Computing Library  1.2.0
circulant_matrix_operations.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_CIRCULANT_MATRIX_OPERATIONS_HPP_
2 #define VIENNACL_LINALG_CIRCULANT_MATRIX_OPERATIONS_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/ocl/device.hpp"
26 #include "viennacl/ocl/handle.hpp"
27 #include "viennacl/ocl/kernel.hpp"
28 #include "viennacl/scalar.hpp"
29 #include "viennacl/vector.hpp"
30 #include "viennacl/tools/tools.hpp"
31 #include "viennacl/fft.hpp"
32 //#include "viennacl/linalg/kernels/coordinate_matrix_kernels.h"
33 
34 namespace viennacl
35 {
36  namespace linalg
37  {
38 
39 
40  // A * x
48  template<class SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
49  vector_expression<const circulant_matrix<SCALARTYPE, ALIGNMENT>,
50  const vector<SCALARTYPE, VECTOR_ALIGNMENT>,
53  {
56  op_prod >(mat, vec);
57  }
58 
59  // A * x
68  template<class SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
73  size_t NUM_THREADS)
74  {
77  viennacl::op_prod >(mat, vec);
78  }
79 
88  template<class SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
92  {
93  assert(mat.size1() == result.size());
94  assert(mat.size2() == vec.size());
95  //result.clear();
96 
97  //std::cout << "prod(circulant_matrix" << ALIGNMENT << ", vector) called with internal_nnz=" << mat.internal_nnz() << std::endl;
98 
100  viennacl::detail::fft::real_to_complex(mat.elements(), circ, mat.elements().size());
101 
104 
105  viennacl::detail::fft::real_to_complex(vec, tmp, vec.size());
106  viennacl::linalg::convolve(circ, tmp, tmp2);
107  viennacl::detail::fft::complex_to_real(tmp2, result, vec.size());
108 
109  }
110 
111  } //namespace linalg
112 
113 
114 
119  template <typename SCALARTYPE, unsigned int ALIGNMENT>
120  template <unsigned int MAT_ALIGNMENT>
124  viennacl::op_prod> & proxy)
125  {
126  // check for the special case x = A * x
127  if (proxy.rhs().handle() == this->handle())
128  {
129  viennacl::vector<SCALARTYPE, ALIGNMENT> result(proxy.rhs().size());
130  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
131  *this = result;
132  return *this;
133  }
134  else
135  {
136  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
137  return *this;
138  }
139  return *this;
140  }
141 
142  //v += A * x
147  template <typename SCALARTYPE, unsigned int ALIGNMENT>
148  template <unsigned int MAT_ALIGNMENT>
152  op_prod> & proxy)
153  {
154  vector<SCALARTYPE, ALIGNMENT> result(proxy.lhs().size1());
155  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
156  *this += result;
157  return *this;
158  }
159 
164  template <typename SCALARTYPE, unsigned int ALIGNMENT>
165  template <unsigned int MAT_ALIGNMENT>
169  op_prod> & proxy)
170  {
171  vector<SCALARTYPE, ALIGNMENT> result(proxy.get_lhs().size1());
172  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
173  *this -= result;
174  return *this;
175  }
176 
177 
178  //free functions:
183  template <typename SCALARTYPE, unsigned int ALIGNMENT>
184  template <unsigned int MAT_ALIGNMENT>
188  op_prod> & proxy)
189  {
190  assert(proxy.get_lhs().size1() == size());
192  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
193  result += *this;
194  return result;
195  }
196 
201  template <typename SCALARTYPE, unsigned int ALIGNMENT>
202  template <unsigned int MAT_ALIGNMENT>
206  op_prod> & proxy)
207  {
208  assert(proxy.get_lhs().size1() == size());
210  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
211  result = *this - result;
212  return result;
213  }
214 
215 } //namespace viennacl
216 
217 
218 #endif