ViennaCL - The Vienna Computing Library  1.5.1
ilu0.hpp
Go to the documentation of this file.
1 
2 #ifndef VIENNACL_LINALG_DETAIL_ILU0_HPP_
3 #define VIENNACL_LINALG_DETAIL_ILU0_HPP_
4 
5 /* =========================================================================
6  Copyright (c) 2010-2014, Institute for Microelectronics,
7  Institute for Analysis and Scientific Computing,
8  TU Wien.
9  Portions of this software are copyright by UChicago Argonne, LLC.
10 
11  -----------------
12  ViennaCL - The Vienna Computing Library
13  -----------------
14 
15  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
16 
17  (A list of authors and contributors can be found in the PDF manual)
18 
19  License: MIT (X11), see file LICENSE in the base directory
20 ============================================================================= */
21 
39 #include <vector>
40 #include <cmath>
41 #include <iostream>
42 #include "viennacl/forwards.h"
43 #include "viennacl/tools/tools.hpp"
47 
49 
50 #include <map>
51 
52 namespace viennacl
53 {
54  namespace linalg
55  {
56 
59  class ilu0_tag
60  {
61  public:
62  ilu0_tag(bool with_level_scheduling = false) : use_level_scheduling_(with_level_scheduling) {}
63 
64  bool use_level_scheduling() const { return use_level_scheduling_; }
65  void use_level_scheduling(bool b) { use_level_scheduling_ = b; }
66 
67  private:
68  bool use_level_scheduling_;
69  };
70 
71 
78  template<typename ScalarType>
80  {
81  assert( (A.handle1().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
82  assert( (A.handle2().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
83  assert( (A.handle().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
84 
85  ScalarType * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(A.handle());
86  unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.handle1());
87  unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.handle2());
88 
89  // Note: Line numbers in the following refer to the algorithm in Saad's book
90 
91  for (vcl_size_t i=1; i<A.size1(); ++i) // Line 1
92  {
93  unsigned int row_i_begin = row_buffer[i];
94  unsigned int row_i_end = row_buffer[i+1];
95  for (unsigned int buf_index_k = row_i_begin; buf_index_k < row_i_end; ++buf_index_k) //Note: We do not assume that the column indices within a row are sorted
96  {
97  unsigned int k = col_buffer[buf_index_k];
98  if (k >= i)
99  continue; //Note: We do not assume that the column indices within a row are sorted
100 
101  unsigned int row_k_begin = row_buffer[k];
102  unsigned int row_k_end = row_buffer[k+1];
103 
104  // get a_kk:
105  ScalarType a_kk = 0;
106  for (unsigned int buf_index_akk = row_k_begin; buf_index_akk < row_k_end; ++buf_index_akk)
107  {
108  if (col_buffer[buf_index_akk] == k)
109  {
110  a_kk = elements[buf_index_akk];
111  break;
112  }
113  }
114 
115  ScalarType & a_ik = elements[buf_index_k];
116  a_ik /= a_kk; //Line 3
117 
118  for (unsigned int buf_index_j = row_i_begin; buf_index_j < row_i_end; ++buf_index_j) //Note: We do not assume that the column indices within a row are sorted
119  {
120  unsigned int j = col_buffer[buf_index_j];
121  if (j <= k)
122  continue;
123 
124  // determine a_kj:
125  ScalarType a_kj = 0;
126  for (unsigned int buf_index_akj = row_k_begin; buf_index_akj < row_k_end; ++buf_index_akj)
127  {
128  if (col_buffer[buf_index_akj] == j)
129  {
130  a_kk = elements[buf_index_akj];
131  break;
132  }
133  }
134 
135  //a_ij -= a_ik * a_kj
136  elements[buf_index_j] -= a_ik * a_kj; //Line 5
137  }
138  }
139  }
140 
141  }
142 
143 
146  template <typename MatrixType>
148  {
149  typedef typename MatrixType::value_type ScalarType;
150 
151  public:
152  ilu0_precond(MatrixType const & mat, ilu0_tag const & tag) : tag_(tag), LU()
153  {
154  //initialize preconditioner:
155  //std::cout << "Start CPU precond" << std::endl;
156  init(mat);
157  //std::cout << "End CPU precond" << std::endl;
158  }
159 
160  template <typename VectorType>
161  void apply(VectorType & vec) const
162  {
163  unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle1());
164  unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle2());
165  ScalarType const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(LU.handle());
166 
167  viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec, LU.size2(), unit_lower_tag());
168  viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec, LU.size2(), upper_tag());
169  }
170 
171  private:
172  void init(MatrixType const & mat)
173  {
175  viennacl::switch_memory_context(LU, host_context);
176 
177  viennacl::copy(mat, LU);
179  }
180 
181  ilu0_tag const & tag_;
182 
184  };
185 
186 
191  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
192  class ilu0_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
193  {
195 
196  public:
197  ilu0_precond(MatrixType const & mat, ilu0_tag const & tag) : tag_(tag), LU(mat.size1(), mat.size2())
198  {
199  //initialize preconditioner:
200  //std::cout << "Start GPU precond" << std::endl;
201  init(mat);
202  //std::cout << "End GPU precond" << std::endl;
203  }
204 
205  void apply(vector<ScalarType> & vec) const
206  {
209  {
210  if (tag_.use_level_scheduling())
211  {
212  //std::cout << "Using multifrontal on GPU..." << std::endl;
214  multifrontal_L_row_index_arrays_,
215  multifrontal_L_row_buffers_,
216  multifrontal_L_col_buffers_,
217  multifrontal_L_element_buffers_,
218  multifrontal_L_row_elimination_num_list_);
219 
220  vec = viennacl::linalg::element_div(vec, multifrontal_U_diagonal_);
221 
223  multifrontal_U_row_index_arrays_,
224  multifrontal_U_row_buffers_,
225  multifrontal_U_col_buffers_,
226  multifrontal_U_element_buffers_,
227  multifrontal_U_row_elimination_num_list_);
228  }
229  else
230  {
231  viennacl::context old_context = viennacl::traits::context(vec);
232  viennacl::switch_memory_context(vec, host_context);
235  viennacl::switch_memory_context(vec, old_context);
236  }
237  }
238  else //apply ILU0 directly on CPU
239  {
240  if (tag_.use_level_scheduling())
241  {
242  //std::cout << "Using multifrontal..." << std::endl;
244  multifrontal_L_row_index_arrays_,
245  multifrontal_L_row_buffers_,
246  multifrontal_L_col_buffers_,
247  multifrontal_L_element_buffers_,
248  multifrontal_L_row_elimination_num_list_);
249 
250  vec = viennacl::linalg::element_div(vec, multifrontal_U_diagonal_);
251 
253  multifrontal_U_row_index_arrays_,
254  multifrontal_U_row_buffers_,
255  multifrontal_U_col_buffers_,
256  multifrontal_U_element_buffers_,
257  multifrontal_U_row_elimination_num_list_);
258  }
259  else
260  {
263  }
264  }
265  }
266 
267  vcl_size_t levels() const { return multifrontal_L_row_index_arrays_.size(); }
268 
269  private:
270  void init(MatrixType const & mat)
271  {
273  viennacl::switch_memory_context(LU, host_context);
274  LU = mat;
276 
277  if (!tag_.use_level_scheduling())
278  return;
279 
280  // multifrontal part:
281  viennacl::switch_memory_context(multifrontal_U_diagonal_, host_context);
282  multifrontal_U_diagonal_.resize(LU.size1(), false);
284 
286  multifrontal_U_diagonal_, //dummy
287  multifrontal_L_row_index_arrays_,
288  multifrontal_L_row_buffers_,
289  multifrontal_L_col_buffers_,
290  multifrontal_L_element_buffers_,
291  multifrontal_L_row_elimination_num_list_);
292 
293 
295  multifrontal_U_diagonal_,
296  multifrontal_U_row_index_arrays_,
297  multifrontal_U_row_buffers_,
298  multifrontal_U_col_buffers_,
299  multifrontal_U_element_buffers_,
300  multifrontal_U_row_elimination_num_list_);
301 
302  //
303  // Bring to device if necessary:
304  //
305 
306  // L:
307  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_row_index_arrays_.begin();
308  it != multifrontal_L_row_index_arrays_.end();
309  ++it)
310  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
311 
312  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_row_buffers_.begin();
313  it != multifrontal_L_row_buffers_.end();
314  ++it)
315  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
316 
317  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_col_buffers_.begin();
318  it != multifrontal_L_col_buffers_.end();
319  ++it)
320  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
321 
322  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_element_buffers_.begin();
323  it != multifrontal_L_element_buffers_.end();
324  ++it)
325  viennacl::backend::switch_memory_context<ScalarType>(*it, viennacl::traits::context(mat));
326 
327 
328  // U:
329 
330  viennacl::switch_memory_context(multifrontal_U_diagonal_, viennacl::traits::context(mat));
331 
332  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_row_index_arrays_.begin();
333  it != multifrontal_U_row_index_arrays_.end();
334  ++it)
335  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
336 
337  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_row_buffers_.begin();
338  it != multifrontal_U_row_buffers_.end();
339  ++it)
340  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
341 
342  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_col_buffers_.begin();
343  it != multifrontal_U_col_buffers_.end();
344  ++it)
345  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
346 
347  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_element_buffers_.begin();
348  it != multifrontal_U_element_buffers_.end();
349  ++it)
350  viennacl::backend::switch_memory_context<ScalarType>(*it, viennacl::traits::context(mat));
351 
352  }
353 
354  ilu0_tag const & tag_;
356 
357  std::list< viennacl::backend::mem_handle > multifrontal_L_row_index_arrays_;
358  std::list< viennacl::backend::mem_handle > multifrontal_L_row_buffers_;
359  std::list< viennacl::backend::mem_handle > multifrontal_L_col_buffers_;
360  std::list< viennacl::backend::mem_handle > multifrontal_L_element_buffers_;
361  std::list< vcl_size_t > multifrontal_L_row_elimination_num_list_;
362 
363  viennacl::vector<ScalarType> multifrontal_U_diagonal_;
364  std::list< viennacl::backend::mem_handle > multifrontal_U_row_index_arrays_;
365  std::list< viennacl::backend::mem_handle > multifrontal_U_row_buffers_;
366  std::list< viennacl::backend::mem_handle > multifrontal_U_col_buffers_;
367  std::list< viennacl::backend::mem_handle > multifrontal_U_element_buffers_;
368  std::list< vcl_size_t > multifrontal_U_row_elimination_num_list_;
369 
370  };
371 
372  }
373 }
374 
375 
376 
377 
378 #endif
379 
380 
381 
void level_scheduling_setup_L(viennacl::compressed_matrix< ScalarType, ALIGNMENT > const &LU, vector< ScalarType > const &diagonal_LU, std::list< viennacl::backend::mem_handle > &row_index_arrays, std::list< viennacl::backend::mem_handle > &row_buffers, std::list< viennacl::backend::mem_handle > &col_buffers, std::list< viennacl::backend::mem_handle > &element_buffers, std::list< vcl_size_t > &row_elimination_num_list)
Definition: common.hpp:191
viennacl::vector_expression< const vector_base< T >, const vector_base< T >, op_element_binary< op_div > > element_div(vector_base< T > const &v1, vector_base< T > const &v2)
std::size_t vcl_size_t
Definition: forwards.h:58
void level_scheduling_setup_U(viennacl::compressed_matrix< ScalarType, ALIGNMENT > const &LU, vector< ScalarType > const &diagonal_LU, std::list< viennacl::backend::mem_handle > &row_index_arrays, std::list< viennacl::backend::mem_handle > &row_buffers, std::list< viennacl::backend::mem_handle > &col_buffers, std::list< viennacl::backend::mem_handle > &element_buffers, std::list< vcl_size_t > &row_elimination_num_list)
Definition: common.hpp:208
ILU0 preconditioner class, can be supplied to solve()-routines.
Definition: ilu0.hpp:147
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
Various little tools used here and there in ViennaCL.
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 for incomplete LU factorization with static pattern (ILU0)
Definition: ilu0.hpp:59
bool use_level_scheduling() const
Definition: ilu0.hpp:64
This file provides the forward declarations for the main types used within ViennaCL.
const handle_type & handle2() const
Returns the OpenCL handle to the column index array.
Definition: compressed_matrix.hpp:701
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
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:694
ilu0_tag(bool with_level_scheduling=false)
Definition: ilu0.hpp:62
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
void use_level_scheduling(bool b)
Definition: ilu0.hpp:65
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
A tag class representing an upper triangular matrix.
Definition: forwards.h:708
void apply(VectorType &vec) const
Definition: ilu0.hpp:161
Implementation of the compressed_matrix class.
void copy(std::vector< SCALARTYPE > &cpu_vec, circulant_matrix< SCALARTYPE, ALIGNMENT > &gpu_mat)
Copies a circulant matrix from the std::vector to the OpenCL device (either GPU or multi-core CPU) ...
Definition: circulant_matrix.hpp:150
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:692
ilu0_precond(MatrixType const &mat, ilu0_tag const &tag)
Definition: ilu0.hpp:152
void precondition(viennacl::compressed_matrix< ScalarType > &A, ilu0_tag const &)
Implementation of a ILU-preconditioner with static pattern. Optimized version for CSR matrices...
Definition: ilu0.hpp:79
Common routines used within ILU-type preconditioners.
ilu0_precond(MatrixType const &mat, ilu0_tag const &tag)
Definition: ilu0.hpp:197
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:41
Common routines for single-threaded or OpenMP-enabled execution on CPU.
Definition: forwards.h:479
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:713
const handle_type & handle() const
Returns the memory handle.
Definition: vector.hpp:878
void level_scheduling_substitute(vector< ScalarType > &vec, std::list< viennacl::backend::mem_handle > const &row_index_arrays, std::list< viennacl::backend::mem_handle > const &row_buffers, std::list< viennacl::backend::mem_handle > const &col_buffers, std::list< viennacl::backend::mem_handle > const &element_buffers, std::list< vcl_size_t > const &row_elimination_num_list)
Definition: common.hpp:224
void apply(vector< ScalarType > &vec) const
Definition: ilu0.hpp:205
Main interface routines for memory management.
void row_info(compressed_matrix< ScalarType, MAT_ALIGNMENT > const &mat, vector_base< ScalarType > &vec, viennacl::linalg::detail::row_info_types info_selector)
Definition: sparse_matrix_operations.hpp:47
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:703
const handle_type & handle1() const
Returns the OpenCL handle to the row index array.
Definition: compressed_matrix.hpp:699
void switch_memory_context(T &obj, viennacl::context new_ctx)
Generic convenience routine for migrating data of an object to a new memory domain.
Definition: memory.hpp:624