ViennaCL - The Vienna Computing Library  1.5.1
direct_solve.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_HOST_BASED_DIRECT_SOLVE_HPP
2 #define VIENNACL_LINALG_HOST_BASED_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/vector.hpp"
26 #include "viennacl/matrix.hpp"
27 
29 
30 namespace viennacl
31 {
32  namespace linalg
33  {
34  namespace host_based
35  {
36 
37  namespace detail
38  {
39  //
40  // Upper solve:
41  //
42  template <typename MatrixType1, typename MatrixType2>
43  void upper_inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
44  {
45  typedef typename MatrixType2::value_type value_type;
46 
47  for (vcl_size_t i = 0; i < A_size; ++i)
48  {
49  vcl_size_t current_row = A_size - i - 1;
50 
51  for (vcl_size_t j = current_row + 1; j < A_size; ++j)
52  {
53  value_type A_element = A(current_row, j);
54  for (vcl_size_t k=0; k < B_size; ++k)
55  B(current_row, k) -= A_element * B(j, k);
56  }
57 
58  if (!unit_diagonal)
59  {
60  value_type A_diag = A(current_row, current_row);
61  for (vcl_size_t k=0; k < B_size; ++k)
62  B(current_row, k) /= A_diag;
63  }
64  }
65  }
66 
67  template <typename MatrixType1, typename MatrixType2>
68  void inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_upper_tag)
69  {
70  upper_inplace_solve_matrix(A, B, A_size, B_size, true);
71  }
72 
73  template <typename MatrixType1, typename MatrixType2>
74  void inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::upper_tag)
75  {
76  upper_inplace_solve_matrix(A, B, A_size, B_size, false);
77  }
78 
79  //
80  // Lower solve:
81  //
82  template <typename MatrixType1, typename MatrixType2>
83  void lower_inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
84  {
85  typedef typename MatrixType2::value_type value_type;
86 
87  for (vcl_size_t i = 0; i < A_size; ++i)
88  {
89  for (vcl_size_t j = 0; j < i; ++j)
90  {
91  value_type A_element = A(i, j);
92  for (vcl_size_t k=0; k < B_size; ++k)
93  B(i, k) -= A_element * B(j, k);
94  }
95 
96  if (!unit_diagonal)
97  {
98  value_type A_diag = A(i, i);
99  for (vcl_size_t k=0; k < B_size; ++k)
100  B(i, k) /= A_diag;
101  }
102  }
103  }
104 
105  template <typename MatrixType1, typename MatrixType2>
106  void inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_lower_tag)
107  {
108  lower_inplace_solve_matrix(A, B, A_size, B_size, true);
109  }
110 
111  template <typename MatrixType1, typename MatrixType2>
112  void inplace_solve_matrix(MatrixType1 & A, MatrixType2 & B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::lower_tag)
113  {
114  lower_inplace_solve_matrix(A, B, A_size, B_size, false);
115  }
116 
117  }
118 
119  //
120  // Note: By convention, all size checks are performed in the calling frontend. No need to double-check here.
121  //
122 
124 
129  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
131  {
132  typedef NumericT value_type;
133 
134  value_type const * data_A = detail::extract_raw_pointer<value_type>(A);
135  value_type * data_B = detail::extract_raw_pointer<value_type>(B);
136 
137  vcl_size_t A_start1 = viennacl::traits::start1(A);
138  vcl_size_t A_start2 = viennacl::traits::start2(A);
141  vcl_size_t A_size2 = viennacl::traits::size2(A);
142  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(A);
143  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(A);
144 
145  vcl_size_t B_start1 = viennacl::traits::start1(B);
146  vcl_size_t B_start2 = viennacl::traits::start2(B);
149  vcl_size_t B_size2 = viennacl::traits::size2(B);
150  vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(B);
151  vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(B);
152 
153 
154  detail::matrix_array_wrapper<value_type const, typename F1::orientation_category, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
155  detail::matrix_array_wrapper<value_type, typename F2::orientation_category, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
156 
157  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SOLVERTAG());
158  }
159 
165  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
168  SOLVERTAG)
169  {
170  typedef NumericT value_type;
171 
172  value_type const * data_A = detail::extract_raw_pointer<value_type>(A);
173  value_type * data_B = const_cast<value_type *>(detail::extract_raw_pointer<value_type>(proxy_B.lhs()));
174 
175  vcl_size_t A_start1 = viennacl::traits::start1(A);
176  vcl_size_t A_start2 = viennacl::traits::start2(A);
179  vcl_size_t A_size2 = viennacl::traits::size2(A);
180  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(A);
181  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(A);
182 
183  vcl_size_t B_start1 = viennacl::traits::start1(proxy_B.lhs());
184  vcl_size_t B_start2 = viennacl::traits::start2(proxy_B.lhs());
185  vcl_size_t B_inc1 = viennacl::traits::stride1(proxy_B.lhs());
186  vcl_size_t B_inc2 = viennacl::traits::stride2(proxy_B.lhs());
187  vcl_size_t B_size1 = viennacl::traits::size1(proxy_B.lhs());
188  vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(proxy_B.lhs());
189  vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(proxy_B.lhs());
190 
191 
192  detail::matrix_array_wrapper<value_type const, typename F1::orientation_category, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
193  detail::matrix_array_wrapper<value_type, typename F2::orientation_category, true> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
194 
195  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size1, SOLVERTAG());
196  }
197 
198  //upper triangular solver for transposed lower triangular matrices
204  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
207  SOLVERTAG)
208  {
209  typedef NumericT value_type;
210 
211  value_type const * data_A = detail::extract_raw_pointer<value_type>(proxy_A.lhs());
212  value_type * data_B = const_cast<value_type *>(detail::extract_raw_pointer<value_type>(B));
213 
214  vcl_size_t A_start1 = viennacl::traits::start1(proxy_A.lhs());
215  vcl_size_t A_start2 = viennacl::traits::start2(proxy_A.lhs());
216  vcl_size_t A_inc1 = viennacl::traits::stride1(proxy_A.lhs());
217  vcl_size_t A_inc2 = viennacl::traits::stride2(proxy_A.lhs());
218  vcl_size_t A_size2 = viennacl::traits::size2(proxy_A.lhs());
219  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(proxy_A.lhs());
220  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(proxy_A.lhs());
221 
222  vcl_size_t B_start1 = viennacl::traits::start1(B);
223  vcl_size_t B_start2 = viennacl::traits::start2(B);
226  vcl_size_t B_size2 = viennacl::traits::size2(B);
227  vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(B);
228  vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(B);
229 
230 
231  detail::matrix_array_wrapper<value_type const, typename F1::orientation_category, true> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
232  detail::matrix_array_wrapper<value_type, typename F2::orientation_category, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
233 
234  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SOLVERTAG());
235  }
236 
242  template <typename NumericT, typename F1, typename F2, typename SOLVERTAG>
245  SOLVERTAG)
246  {
247  typedef NumericT value_type;
248 
249  value_type const * data_A = detail::extract_raw_pointer<value_type>(proxy_A.lhs());
250  value_type * data_B = const_cast<value_type *>(detail::extract_raw_pointer<value_type>(proxy_B.lhs()));
251 
252  vcl_size_t A_start1 = viennacl::traits::start1(proxy_A.lhs());
253  vcl_size_t A_start2 = viennacl::traits::start2(proxy_A.lhs());
254  vcl_size_t A_inc1 = viennacl::traits::stride1(proxy_A.lhs());
255  vcl_size_t A_inc2 = viennacl::traits::stride2(proxy_A.lhs());
256  vcl_size_t A_size2 = viennacl::traits::size2(proxy_A.lhs());
257  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(proxy_A.lhs());
258  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(proxy_A.lhs());
259 
260  vcl_size_t B_start1 = viennacl::traits::start1(proxy_B.lhs());
261  vcl_size_t B_start2 = viennacl::traits::start2(proxy_B.lhs());
262  vcl_size_t B_inc1 = viennacl::traits::stride1(proxy_B.lhs());
263  vcl_size_t B_inc2 = viennacl::traits::stride2(proxy_B.lhs());
264  vcl_size_t B_size1 = viennacl::traits::size1(proxy_B.lhs());
265  vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(proxy_B.lhs());
266  vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(proxy_B.lhs());
267 
268 
269  detail::matrix_array_wrapper<value_type const, typename F1::orientation_category, true> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
270  detail::matrix_array_wrapper<value_type, typename F2::orientation_category, true> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
271 
272  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size1, SOLVERTAG());
273  }
274 
275  //
276  // Solve on vector
277  //
278 
279  namespace detail
280  {
281  //
282  // Upper solve:
283  //
284  template <typename MatrixType, typename VectorType>
285  void upper_inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, bool unit_diagonal)
286  {
287  typedef typename VectorType::value_type value_type;
288 
289  for (vcl_size_t i = 0; i < A_size; ++i)
290  {
291  vcl_size_t current_row = A_size - i - 1;
292 
293  for (vcl_size_t j = current_row + 1; j < A_size; ++j)
294  {
295  value_type A_element = A(current_row, j);
296  b(current_row) -= A_element * b(j);
297  }
298 
299  if (!unit_diagonal)
300  b(current_row) /= A(current_row, current_row);
301  }
302  }
303 
304  template <typename MatrixType, typename VectorType>
305  void inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
306  {
307  upper_inplace_solve_vector(A, b, A_size, true);
308  }
309 
310  template <typename MatrixType, typename VectorType>
311  void inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, viennacl::linalg::upper_tag)
312  {
313  upper_inplace_solve_vector(A, b, A_size, false);
314  }
315 
316  //
317  // Lower solve:
318  //
319  template <typename MatrixType, typename VectorType>
320  void lower_inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, bool unit_diagonal)
321  {
322  typedef typename VectorType::value_type value_type;
323 
324  for (vcl_size_t i = 0; i < A_size; ++i)
325  {
326  for (vcl_size_t j = 0; j < i; ++j)
327  {
328  value_type A_element = A(i, j);
329  b(i) -= A_element * b(j);
330  }
331 
332  if (!unit_diagonal)
333  b(i) /= A(i, i);
334  }
335  }
336 
337  template <typename MatrixType, typename VectorType>
338  void inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, viennacl::linalg::unit_lower_tag)
339  {
340  lower_inplace_solve_vector(A, b, A_size, true);
341  }
342 
343  template <typename MatrixType, typename VectorType>
344  void inplace_solve_vector(MatrixType & A, VectorType & b, vcl_size_t A_size, viennacl::linalg::lower_tag)
345  {
346  lower_inplace_solve_vector(A, b, A_size, false);
347  }
348 
349  }
350 
351  template <typename NumericT, typename F, typename SOLVERTAG>
353  vector_base<NumericT> & vec,
354  SOLVERTAG)
355  {
356  typedef NumericT value_type;
357 
358  value_type const * data_A = detail::extract_raw_pointer<value_type>(mat);
359  value_type * data_v = detail::extract_raw_pointer<value_type>(vec);
360 
361  vcl_size_t A_start1 = viennacl::traits::start1(mat);
362  vcl_size_t A_start2 = viennacl::traits::start2(mat);
365  vcl_size_t A_size2 = viennacl::traits::size2(mat);
366  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(mat);
367  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(mat);
368 
371 
372  detail::matrix_array_wrapper<value_type const, typename F::orientation_category, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
373  detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
374 
375  detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SOLVERTAG());
376  }
377 
378 
379 
385  template <typename NumericT, typename F, typename SOLVERTAG>
387  vector_base<NumericT> & vec,
388  SOLVERTAG)
389  {
390  typedef NumericT value_type;
391 
392  value_type const * data_A = detail::extract_raw_pointer<value_type>(proxy.lhs());
393  value_type * data_v = detail::extract_raw_pointer<value_type>(vec);
394 
395  vcl_size_t A_start1 = viennacl::traits::start1(proxy.lhs());
396  vcl_size_t A_start2 = viennacl::traits::start2(proxy.lhs());
397  vcl_size_t A_inc1 = viennacl::traits::stride1(proxy.lhs());
398  vcl_size_t A_inc2 = viennacl::traits::stride2(proxy.lhs());
399  vcl_size_t A_size2 = viennacl::traits::size2(proxy.lhs());
400  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(proxy.lhs());
401  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(proxy.lhs());
402 
405 
406  detail::matrix_array_wrapper<value_type const, typename F::orientation_category, true> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
407  detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
408 
409  detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SOLVERTAG());
410  }
411 
412 
413 
414  }
415  }
416 }
417 
418 #endif
Helper class for accessing a strided subvector of a larger vector.
Definition: common.hpp:49
std::size_t vcl_size_t
Definition: forwards.h:58
result_of::size_type< matrix_base< NumericT, F > >::type stride2(matrix_base< NumericT, F > const &s)
Definition: stride.hpp:68
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 tag class representing a lower triangular matrix.
Definition: forwards.h:703
void inplace_solve_vector(MatrixType &A, VectorType &b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
Definition: direct_solve.hpp:305
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
result_of::size_type< viennacl::vector_base< T > >::type stride(viennacl::vector_base< T > const &s)
Definition: stride.hpp:46
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:64
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
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
result_of::size_type< matrix_base< NumericT, F > >::type stride1(matrix_base< NumericT, F > const &s)
Definition: stride.hpp:57
result_of::size_type< T >::type start2(T const &obj)
Definition: start.hpp:83
Helper array for accessing a strided submatrix embedded in a larger matrix.
Definition: common.hpp:100
A tag class representing an upper triangular matrix.
Definition: forwards.h:708
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:43
Common base class for dense vectors, vector ranges, and vector slices.
Definition: forwards.h:205
vcl_size_t internal_size2(matrix_base< NumericT, F > const &mat)
Helper routine for obtaining the internal number of entries per column of a ViennaCL matrix...
Definition: size.hpp:287
void inplace_solve_matrix(MatrixType1 &A, MatrixType2 &B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_upper_tag)
Definition: direct_solve.hpp:68
Common routines for single-threaded or OpenMP-enabled execution on CPU.
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:713
void lower_inplace_solve_matrix(MatrixType1 &A, MatrixType2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
Definition: direct_solve.hpp:83
A tag class representing transposed matrices.
Definition: forwards.h:165
void lower_inplace_solve_vector(MatrixType &A, VectorType &b, vcl_size_t A_size, bool unit_diagonal)
Definition: direct_solve.hpp:320
vcl_size_t internal_size1(matrix_base< NumericT, F > const &mat)
Helper routine for obtaining the internal number of entries per row of a ViennaCL matrix...
Definition: size.hpp:279
void upper_inplace_solve_matrix(MatrixType1 &A, MatrixType2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
Definition: direct_solve.hpp:43
void upper_inplace_solve_vector(MatrixType &A, VectorType &b, vcl_size_t A_size, bool unit_diagonal)
Definition: direct_solve.hpp:285
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:718