ViennaCL - The Vienna Computing Library  1.2.0
compressed_matrix_operations.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_COMPRESSED_MATRIX_OPERATIONS_HPP_
2 #define VIENNACL_COMPRESSED_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"
32 
33 namespace viennacl
34 {
35  namespace linalg
36  {
37  // A * x
45  template<class SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
46  vector_expression<const compressed_matrix<SCALARTYPE, ALIGNMENT>,
47  const vector<SCALARTYPE, VECTOR_ALIGNMENT>,
50  {
53  op_prod >(mat, vec);
54  }
55 
65  template<class TYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
69  size_t NUM_THREADS = 0)
70  {
71  assert(mat.size1() == result.size());
72  assert(mat.size2() == vec.size());
73 
74  viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::compressed_matrix<TYPE, ALIGNMENT>::program_name(), "vec_mul");
75 
76  viennacl::ocl::enqueue(k(mat.handle1(), mat.handle2(), mat.handle(), vec, result, static_cast<cl_uint>(mat.size1())));
77  }
78 
84  template<typename SCALARTYPE, unsigned int MAT_ALIGNMENT, unsigned int VEC_ALIGNMENT>
86  {
87  viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::compressed_matrix<SCALARTYPE, MAT_ALIGNMENT>::program_name(), "lu_forward");
88  unsigned int threads = k.local_work_size();
89 
92  viennacl::ocl::local_mem(sizeof(int) * (threads+1)),
93  viennacl::ocl::local_mem(sizeof(SCALARTYPE) * threads),
94  vec, L.size1()));
95  }
96 
103  template<typename SCALARTYPE, unsigned int MAT_ALIGNMENT, unsigned int VEC_ALIGNMENT, typename TAG>
107  {
108  // do an inplace solve on the result vector:
110  result = vec;
111 
112  inplace_solve(L, result, tag);
113 
114  return result;
115  }
116 
117 
123  template<typename SCALARTYPE, unsigned int MAT_ALIGNMENT, unsigned int VEC_ALIGNMENT>
125  {
126  viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::compressed_matrix<SCALARTYPE, MAT_ALIGNMENT>::program_name(), "lu_backward");
127  unsigned int threads = k.local_work_size();
128 
131  viennacl::ocl::local_mem(sizeof(int) * (threads+2)),
132  viennacl::ocl::local_mem(sizeof(SCALARTYPE) * (threads+2)),
133  vec, U.size1()));
134  }
135 
142  template<typename SCALARTYPE, unsigned int MAT_ALIGNMENT, unsigned int VEC_ALIGNMENT, typename TAG>
145  viennacl::linalg::upper_tag const & tag)
146  {
147  // do an inplace solve on the result vector:
149  result = vec;
150 
151  inplace_solve(L, result, tag);
152 
153  return result;
154  }
155 
156 
157  } //namespace linalg
158 
159 
160 
161  //v = A * x
166  template <typename SCALARTYPE, unsigned int ALIGNMENT>
167  template <unsigned int MAT_ALIGNMENT>
171  viennacl::op_prod> & proxy)
172  {
173  // check for the special case x = A * x
174  if (proxy.rhs().handle() == this->handle())
175  {
176  viennacl::vector<SCALARTYPE, ALIGNMENT> result(proxy.rhs().size());
177  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
178  *this = result;
179  return *this;
180  }
181  else
182  {
183  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), *this);
184  return *this;
185  }
186  return *this;
187  }
188 
189  //v += A * x
194  template <typename SCALARTYPE, unsigned int ALIGNMENT>
195  template <unsigned int MAT_ALIGNMENT>
199  op_prod> & proxy)
200  {
201  vector<SCALARTYPE, ALIGNMENT> result(proxy.lhs().size1());
202  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
203  *this += result;
204  return *this;
205  }
206 
211  template <typename SCALARTYPE, unsigned int ALIGNMENT>
212  template <unsigned int MAT_ALIGNMENT>
216  op_prod> & proxy)
217  {
218  vector<SCALARTYPE, ALIGNMENT> result(proxy.lhs().size1());
219  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
220  *this -= result;
221  return *this;
222  }
223 
224 
225  //free functions:
230  template <typename SCALARTYPE, unsigned int ALIGNMENT>
231  template <unsigned int MAT_ALIGNMENT>
235  op_prod> & proxy)
236  {
237  assert(proxy.lhs().size1() == size());
239  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
240  result += *this;
241  return result;
242  }
243 
248  template <typename SCALARTYPE, unsigned int ALIGNMENT>
249  template <unsigned int MAT_ALIGNMENT>
253  op_prod> & proxy)
254  {
255  assert(proxy.lhs().size1() == size());
257  viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
258  result = *this - result;
259  return result;
260  }
261 
262 } //namespace viennacl
263 
264 
265 #endif