ViennaCL - The Vienna Computing Library  1.7.0
Free open-source GPU-accelerated linear algebra and solver library.
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-2015, 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 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 MatrixT1, typename MatrixT2>
43  void upper_inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
44  {
45  typedef typename MatrixT2::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 MatrixT1, typename MatrixT2>
68  void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & 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 MatrixT1, typename MatrixT2>
74  void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & 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 MatrixT1, typename MatrixT2>
83  void lower_inplace_solve_matrix(MatrixT1 & A, MatrixT2 & B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
84  {
85  typedef typename MatrixT2::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 MatrixT1, typename MatrixT2>
106  void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & 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 MatrixT1, typename MatrixT2>
112  void inplace_solve_matrix(MatrixT1 & A, MatrixT2 & 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 SolverTagT>
132  SolverTagT)
133 {
134  typedef NumericT value_type;
135 
136  value_type const * data_A = detail::extract_raw_pointer<value_type>(A);
137  value_type * data_B = detail::extract_raw_pointer<value_type>(B);
138 
139  vcl_size_t A_start1 = viennacl::traits::start1(A);
140  vcl_size_t A_start2 = viennacl::traits::start2(A);
143  //vcl_size_t A_size1 = viennacl::traits::size1(A);
144  vcl_size_t A_size2 = viennacl::traits::size2(A);
145  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(A);
146  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(A);
147 
148  vcl_size_t B_start1 = viennacl::traits::start1(B);
149  vcl_size_t B_start2 = viennacl::traits::start2(B);
152  //vcl_size_t B_size1 = viennacl::traits::size1(B);
153  vcl_size_t B_size2 = viennacl::traits::size2(B);
154  vcl_size_t B_internal_size1 = viennacl::traits::internal_size1(B);
155  vcl_size_t B_internal_size2 = viennacl::traits::internal_size2(B);
156 
157 
158  if (A.row_major() && B.row_major())
159  {
160  detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
161  detail::matrix_array_wrapper<value_type, row_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
162 
163  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
164  }
165  else if (A.row_major() && !B.row_major())
166  {
167  detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
168  detail::matrix_array_wrapper<value_type, column_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
169 
170  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
171  }
172  else if (!A.row_major() && B.row_major())
173  {
174  detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
175  detail::matrix_array_wrapper<value_type, row_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
176 
177  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
178  }
179  else
180  {
181  detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
182  detail::matrix_array_wrapper<value_type, column_major, false> wrapper_B(data_B, B_start1, B_start2, B_inc1, B_inc2, B_internal_size1, B_internal_size2);
183 
184  detail::inplace_solve_matrix(wrapper_A, wrapper_B, A_size2, B_size2, SolverTagT());
185  }
186 }
187 
188 
189 //
190 // Solve on vector
191 //
192 
193 namespace detail
194 {
195  //
196  // Upper solve:
197  //
198  template<typename MatrixT, typename VectorT>
199  void upper_inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, bool unit_diagonal)
200  {
201  typedef typename VectorT::value_type value_type;
202 
203  for (vcl_size_t i = 0; i < A_size; ++i)
204  {
205  vcl_size_t current_row = A_size - i - 1;
206 
207  for (vcl_size_t j = current_row + 1; j < A_size; ++j)
208  {
209  value_type A_element = A(current_row, j);
210  b(current_row) -= A_element * b(j);
211  }
212 
213  if (!unit_diagonal)
214  b(current_row) /= A(current_row, current_row);
215  }
216  }
217 
218  template<typename MatrixT, typename VectorT>
219  void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
220  {
221  upper_inplace_solve_vector(A, b, A_size, true);
222  }
223 
224  template<typename MatrixT, typename VectorT>
225  void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::upper_tag)
226  {
227  upper_inplace_solve_vector(A, b, A_size, false);
228  }
229 
230  //
231  // Lower solve:
232  //
233  template<typename MatrixT, typename VectorT>
234  void lower_inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, bool unit_diagonal)
235  {
236  typedef typename VectorT::value_type value_type;
237 
238  for (vcl_size_t i = 0; i < A_size; ++i)
239  {
240  for (vcl_size_t j = 0; j < i; ++j)
241  {
242  value_type A_element = A(i, j);
243  b(i) -= A_element * b(j);
244  }
245 
246  if (!unit_diagonal)
247  b(i) /= A(i, i);
248  }
249  }
250 
251  template<typename MatrixT, typename VectorT>
252  void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::unit_lower_tag)
253  {
254  lower_inplace_solve_vector(A, b, A_size, true);
255  }
256 
257  template<typename MatrixT, typename VectorT>
258  void inplace_solve_vector(MatrixT & A, VectorT & b, vcl_size_t A_size, viennacl::linalg::lower_tag)
259  {
260  lower_inplace_solve_vector(A, b, A_size, false);
261  }
262 
263 }
264 
265 template<typename NumericT, typename SolverTagT>
267  vector_base<NumericT> & vec,
268  SolverTagT)
269 {
270  typedef NumericT value_type;
271 
272  value_type const * data_A = detail::extract_raw_pointer<value_type>(mat);
273  value_type * data_v = detail::extract_raw_pointer<value_type>(vec);
274 
275  vcl_size_t A_start1 = viennacl::traits::start1(mat);
276  vcl_size_t A_start2 = viennacl::traits::start2(mat);
279  vcl_size_t A_size2 = viennacl::traits::size2(mat);
280  vcl_size_t A_internal_size1 = viennacl::traits::internal_size1(mat);
281  vcl_size_t A_internal_size2 = viennacl::traits::internal_size2(mat);
282 
285 
286  if (mat.row_major())
287  {
288  detail::matrix_array_wrapper<value_type const, row_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
289  detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
290 
291  detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SolverTagT());
292  }
293  else
294  {
295  detail::matrix_array_wrapper<value_type const, column_major, false> wrapper_A(data_A, A_start1, A_start2, A_inc1, A_inc2, A_internal_size1, A_internal_size2);
296  detail::vector_array_wrapper<value_type> wrapper_v(data_v, start1, inc1);
297 
298  detail::inplace_solve_vector(wrapper_A, wrapper_v, A_size2, SolverTagT());
299  }
300 }
301 
302 
303 } // namespace host_based
304 } // namespace linalg
305 } // namespace viennacl
306 
307 #endif
Helper class for accessing a strided subvector of a larger vector.
Definition: common.hpp:50
result_of::size_type< matrix_base< NumericT > >::type stride1(matrix_base< NumericT > const &s)
Definition: stride.hpp:55
Implementation of the dense matrix class.
vcl_size_t internal_size1(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per row of a ViennaCL matrix...
Definition: size.hpp:382
void lower_inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, bool unit_diagonal)
A tag class representing a lower triangular matrix.
Definition: forwards.h:849
void inplace_solve(matrix_base< NumericT > const &A, matrix_base< NumericT > &B, SolverTagT)
Direct inplace solver for triangular systems with multiple right hand sides, i.e. A \ B (MATLAB notat...
vcl_size_t internal_size2(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per column of a ViennaCL matrix...
Definition: size.hpp:390
result_of::size_type< viennacl::vector_base< T > >::type stride(viennacl::vector_base< T > const &s)
Definition: stride.hpp:45
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:65
void upper_inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, bool unit_diagonal)
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:201
float NumericT
Definition: bisect.cpp:40
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:34
void lower_inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
result_of::size_type< T >::type start2(T const &obj)
Definition: start.hpp:84
Helper array for accessing a strided submatrix embedded in a larger matrix.
Definition: common.hpp:73
Definition: blas3.hpp:36
A tag class representing an upper triangular matrix.
Definition: forwards.h:854
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:44
std::size_t vcl_size_t
Definition: forwards.h:75
Common routines for single-threaded or OpenMP-enabled execution on CPU.
result_of::size_type< matrix_base< NumericT > >::type stride2(matrix_base< NumericT > const &s)
Definition: stride.hpp:65
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
bool row_major() const
Definition: matrix_def.hpp:248
void upper_inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, bool unit_diagonal)
void inplace_solve_vector(MatrixT &A, VectorT &b, vcl_size_t A_size, viennacl::linalg::unit_upper_tag)
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:859
void inplace_solve_matrix(MatrixT1 &A, MatrixT2 &B, vcl_size_t A_size, vcl_size_t B_size, viennacl::linalg::unit_upper_tag)
A tag class representing an upper triangular matrix with unit diagonal.
Definition: forwards.h:864