ViennaCL - The Vienna Computing Library  1.5.2
direct_solve.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DIRECT_SOLVE_HPP_
2 #define VIENNACL_LINALG_DIRECT_SOLVE_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 
25 #include "viennacl/forwards.h"
27 #include "viennacl/vector.hpp"
28 #include "viennacl/matrix.hpp"
30 
31 #ifdef VIENNACL_WITH_OPENCL
33 #endif
34 
35 #ifdef VIENNACL_WITH_CUDA
37 #endif
38 
39 namespace viennacl
40 {
41  namespace linalg
42  {
43 
44  //
45  // A \ B:
46  //
47 
53  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
55  {
56  assert( (viennacl::traits::size1(A) == viennacl::traits::size2(A)) && bool("Size check failed in inplace_solve(): size1(A) != size2(A)"));
57  assert( (viennacl::traits::size1(A) == viennacl::traits::size1(B)) && bool("Size check failed in inplace_solve(): size1(A) != size1(B)"));
58 
59  switch (viennacl::traits::handle(A).get_active_handle_id())
60  {
63  break;
64 #ifdef VIENNACL_WITH_OPENCL
66  viennacl::linalg::opencl::inplace_solve(A, B, SOLVERTAG());
67  break;
68 #endif
69 #ifdef VIENNACL_WITH_CUDA
71  viennacl::linalg::cuda::inplace_solve(A, B, SOLVERTAG());
72  break;
73 #endif
75  throw memory_exception("not initialised!");
76  default:
77  throw memory_exception("not implemented");
78  }
79  }
80 
86  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
89  SOLVERTAG)
90  {
91  assert( (viennacl::traits::size1(A) == viennacl::traits::size2(A)) && bool("Size check failed in inplace_solve(): size1(A) != size2(A)"));
92  assert( (viennacl::traits::size1(A) == viennacl::traits::size1(proxy_B)) && bool("Size check failed in inplace_solve(): size1(A) != size1(B^T)"));
93 
94  switch (viennacl::traits::handle(A).get_active_handle_id())
95  {
97  viennacl::linalg::host_based::inplace_solve(A, proxy_B, SOLVERTAG());
98  break;
99 #ifdef VIENNACL_WITH_OPENCL
101  viennacl::linalg::opencl::inplace_solve(A, proxy_B, SOLVERTAG());
102  break;
103 #endif
104 #ifdef VIENNACL_WITH_CUDA
106  viennacl::linalg::cuda::inplace_solve(A, proxy_B, SOLVERTAG());
107  break;
108 #endif
110  throw memory_exception("not initialised!");
111  default:
112  throw memory_exception("not implemented");
113  }
114  }
115 
116  //upper triangular solver for transposed lower triangular matrices
122  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
125  SOLVERTAG)
126  {
127  assert( (viennacl::traits::size1(proxy_A) == viennacl::traits::size2(proxy_A)) && bool("Size check failed in inplace_solve(): size1(A) != size2(A)"));
128  assert( (viennacl::traits::size1(proxy_A) == viennacl::traits::size1(B)) && bool("Size check failed in inplace_solve(): size1(A^T) != size1(B)"));
129 
130  switch (viennacl::traits::handle(proxy_A.lhs()).get_active_handle_id())
131  {
133  viennacl::linalg::host_based::inplace_solve(proxy_A, B, SOLVERTAG());
134  break;
135 #ifdef VIENNACL_WITH_OPENCL
137  viennacl::linalg::opencl::inplace_solve(proxy_A, B, SOLVERTAG());
138  break;
139 #endif
140 #ifdef VIENNACL_WITH_CUDA
142  viennacl::linalg::cuda::inplace_solve(proxy_A, B, SOLVERTAG());
143  break;
144 #endif
146  throw memory_exception("not initialised!");
147  default:
148  throw memory_exception("not implemented");
149  }
150  }
151 
157  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
160  SOLVERTAG)
161  {
162  assert( (viennacl::traits::size1(proxy_A) == viennacl::traits::size2(proxy_A)) && bool("Size check failed in inplace_solve(): size1(A) != size2(A)"));
163  assert( (viennacl::traits::size1(proxy_A) == viennacl::traits::size1(proxy_B)) && bool("Size check failed in inplace_solve(): size1(A^T) != size1(B^T)"));
164 
165  switch (viennacl::traits::handle(proxy_A.lhs()).get_active_handle_id())
166  {
168  viennacl::linalg::host_based::inplace_solve(proxy_A, proxy_B, SOLVERTAG());
169  break;
170 #ifdef VIENNACL_WITH_OPENCL
172  viennacl::linalg::opencl::inplace_solve(proxy_A, proxy_B, SOLVERTAG());
173  break;
174 #endif
175 #ifdef VIENNACL_WITH_CUDA
177  viennacl::linalg::cuda::inplace_solve(proxy_A, proxy_B, SOLVERTAG());
178  break;
179 #endif
181  throw memory_exception("not initialised!");
182  default:
183  throw memory_exception("not implemented");
184  }
185  }
186 
187  //
188  // A \ b
189  //
190 
191  template <typename NumericT, typename F, typename SOLVERTAG>
193  vector_base<NumericT> & vec,
194  SOLVERTAG)
195  {
196  assert( (mat.size1() == vec.size()) && bool("Size check failed in inplace_solve(): size1(A) != size(b)"));
197  assert( (mat.size2() == vec.size()) && bool("Size check failed in inplace_solve(): size2(A) != size(b)"));
198 
200  {
202  viennacl::linalg::host_based::inplace_solve(mat, vec, SOLVERTAG());
203  break;
204 #ifdef VIENNACL_WITH_OPENCL
206  viennacl::linalg::opencl::inplace_solve(mat, vec, SOLVERTAG());
207  break;
208 #endif
209 #ifdef VIENNACL_WITH_CUDA
211  viennacl::linalg::cuda::inplace_solve(mat, vec, SOLVERTAG());
212  break;
213 #endif
215  throw memory_exception("not initialised!");
216  default:
217  throw memory_exception("not implemented");
218  }
219  }
220 
226  template <typename NumericT, typename F, typename SOLVERTAG>
228  vector_base<NumericT> & vec,
229  SOLVERTAG)
230  {
231  assert( (proxy.lhs().size1() == vec.size()) && bool("Size check failed in inplace_solve(): size1(A) != size(b)"));
232  assert( (proxy.lhs().size2() == vec.size()) && bool("Size check failed in inplace_solve(): size2(A) != size(b)"));
233 
234  switch (viennacl::traits::handle(proxy.lhs()).get_active_handle_id())
235  {
237  viennacl::linalg::host_based::inplace_solve(proxy, vec, SOLVERTAG());
238  break;
239 #ifdef VIENNACL_WITH_OPENCL
241  viennacl::linalg::opencl::inplace_solve(proxy, vec, SOLVERTAG());
242  break;
243 #endif
244 #ifdef VIENNACL_WITH_CUDA
246  viennacl::linalg::cuda::inplace_solve(proxy, vec, SOLVERTAG());
247  break;
248 #endif
250  throw memory_exception("not initialised!");
251  default:
252  throw memory_exception("not implemented");
253  }
254  }
255 
257 
258 
265  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
267  const matrix_base<NumericT, F2> & B,
268  SOLVERTAG tag)
269  {
270  // do an inplace solve on the result vector:
271  matrix<NumericT, F2> result(B);
272 
273  inplace_solve(A, result, tag);
274 
275  return result;
276  }
277 
278 
280 
287  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
290  SOLVERTAG tag)
291  {
292  // do an inplace solve on the result vector:
293  matrix<NumericT, F2> result(proxy);
294 
295  inplace_solve(A, result, tag);
296 
297  return result;
298  }
299 
306  template <typename NumericT, typename F1, typename SOLVERTAG>
308  const vector_base<NumericT> & vec,
309  SOLVERTAG const & tag)
310  {
311  // do an inplace solve on the result vector:
312  vector<NumericT> result(vec);
313 
314  inplace_solve(mat, result, tag);
315 
316  return result;
317  }
318 
319 
321 
327  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
329  const matrix_base<NumericT, F2> & B,
330  SOLVERTAG tag)
331  {
332  // do an inplace solve on the result vector:
333  matrix<NumericT, F2> result(B);
334 
335  inplace_solve(proxy, result, tag);
336 
337  return result;
338  }
339 
340 
347  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
350  SOLVERTAG tag)
351  {
352  // do an inplace solve on the result vector:
353  matrix<NumericT, F2> result(proxy_B);
354 
355  inplace_solve(proxy_A, result, tag);
356 
357  return result;
358  }
359 
366  template <typename NumericT, typename F1, typename SOLVERTAG>
368  const vector_base<NumericT> & vec,
369  SOLVERTAG const & tag)
370  {
371  // do an inplace solve on the result vector:
372  vector<NumericT> result(vec);
373 
374  inplace_solve(proxy, result, tag);
375 
376  return result;
377  }
378 
379 
380  }
381 }
382 
383 #endif
Implementations of dense direct triangular solvers are found here.
size_type size1() const
Returns the number of rows.
Definition: matrix.hpp:625
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector.hpp:837
Definition: forwards.h:478
Exception class in case of memory errors.
Definition: forwards.h:485
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for dense triangular systems. Matlab notation: A \ B.
Definition: direct_solve.hpp:54
Implementation of the dense matrix class.
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
A dense matrix class.
Definition: forwards.h:290
Expression template class for representing a tree of expressions which ultimately result in a matrix...
Definition: forwards.h:283
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG tag)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:297
This file provides the forward declarations for the main types used within ViennaCL.
A dense matrix class.
Definition: forwards.h:293
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
size_type size2() const
Returns the number of columns.
Definition: matrix.hpp:627
Definition: forwards.h:481
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
Definition: direct_solve.hpp:130
memory_types get_active_handle_id() const
Returns an ID for the currently active memory buffer. Other memory buffers might contain old or no da...
Definition: mem_handle.hpp:91
Definition: forwards.h:480
Common base class for dense vectors, vector ranges, and vector slices.
Definition: forwards.h:205
A vector class representing a linear memory sequence on the GPU. Inspired by boost::numeric::ublas::v...
Definition: forwards.h:208
void inplace_solve(const matrix_base< NumericT, F1 > &A, matrix_base< NumericT, F2 > &B, SOLVERTAG)
Direct inplace solver for dense triangular systems. Matlab notation: A \ B.
Definition: direct_solve.hpp:75
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
Definition: forwards.h:479
Implementations of dense direct solvers using CUDA are found here.
A tag class representing transposed matrices.
Definition: forwards.h:165
viennacl::backend::mem_handle & handle(T &obj)
Returns the generic memory handle of an object. Non-const version.
Definition: handle.hpp:41
VectorType solve(const MatrixType &matrix, VectorType const &rhs, bicgstab_tag const &tag)
Implementation of the stabilized Bi-conjugate gradient solver.
Definition: bicgstab.hpp:92
Implementations of dense direct solvers are found here.
Simple enable-if variant that uses the SFINAE pattern.