ViennaCL - The Vienna Computing Library  1.5.1
ilut.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_DETAIL_ILUT_HPP_
2 #define VIENNACL_LINALG_DETAIL_ILUT_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 <vector>
26 #include <cmath>
27 #include <iostream>
28 #include "viennacl/forwards.h"
29 #include "viennacl/tools/tools.hpp"
30 
33 
35 
36 #include <map>
37 
38 namespace viennacl
39 {
40  namespace linalg
41  {
42 
45  class ilut_tag
46  {
47  public:
54  ilut_tag(unsigned int entries_per_row = 20,
55  double drop_tolerance = 1e-4,
56  bool with_level_scheduling = false) : entries_per_row_(entries_per_row), drop_tolerance_(drop_tolerance), use_level_scheduling_(with_level_scheduling) {}
57 
58  void set_drop_tolerance(double tol)
59  {
60  if (tol > 0)
61  drop_tolerance_ = tol;
62  }
63  double get_drop_tolerance() const { return drop_tolerance_; }
64 
65  void set_entries_per_row(unsigned int e)
66  {
67  if (e > 0)
68  entries_per_row_ = e;
69  }
70 
71  unsigned int get_entries_per_row() const { return entries_per_row_; }
72 
73  bool use_level_scheduling() const { return use_level_scheduling_; }
74  void use_level_scheduling(bool b) { use_level_scheduling_ = b; }
75 
76  private:
77  unsigned int entries_per_row_;
78  double drop_tolerance_;
79  bool use_level_scheduling_;
80  };
81 
82 
84  template <typename ScalarType, typename SizeType, typename SparseVector>
86  SizeType row,
87  SparseVector & w)
88  {
89  assert( (A.handle1().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
90  assert( (A.handle2().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
91  assert( (A.handle().get_active_handle_id() == viennacl::MAIN_MEMORY) && bool("System matrix must reside in main memory for ILU0") );
92 
93  ScalarType const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(A.handle());
94  unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.handle1());
95  unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.handle2());
96 
97  SizeType row_i_begin = static_cast<SizeType>(row_buffer[row]);
98  SizeType row_i_end = static_cast<SizeType>(row_buffer[row+1]);
99  ScalarType row_norm = 0;
100  for (SizeType buf_index_i = row_i_begin; buf_index_i < row_i_end; ++buf_index_i) //Note: We do not assume that the column indices within a row are sorted
101  {
102  ScalarType entry = elements[buf_index_i];
103  w[col_buffer[buf_index_i]] = entry;
104  row_norm += entry * entry;
105  }
106  return std::sqrt(row_norm);
107  }
108 
110  template <typename ScalarType, typename SizeType, typename SparseVector>
111  ScalarType setup_w(std::vector< std::map<SizeType, ScalarType> > const & A,
112  SizeType row,
113  SparseVector & w)
114  {
115  ScalarType row_norm = 0;
116  w = A[row];
117  for (typename std::map<SizeType, ScalarType>::const_iterator iter_w = w.begin(); iter_w != w.end(); ++iter_w)
118  row_norm += iter_w->second * iter_w->second;
119 
120  return std::sqrt(row_norm);
121  }
122 
123 
132  template<typename SparseMatrixType, typename ScalarType, typename SizeType>
133  void precondition(SparseMatrixType const & A,
134  std::vector< std::map<SizeType, ScalarType> > & output,
135  ilut_tag const & tag)
136  {
137  typedef std::map<SizeType, ScalarType> SparseVector;
138  typedef typename SparseVector::iterator SparseVectorIterator;
139  typedef typename std::map<SizeType, ScalarType>::const_iterator OutputRowConstIterator;
140  typedef std::multimap<ScalarType, std::pair<SizeType, ScalarType> > TemporarySortMap;
141 
142  assert(viennacl::traits::size1(A) == output.size() && bool("Output matrix size mismatch") );
143 
144  SparseVector w;
145  TemporarySortMap temp_map;
146 
147  for (SizeType i=0; i<viennacl::traits::size1(A); ++i) // Line 1
148  {
149  /* if (i%10 == 0)
150  std::cout << i << std::endl;*/
151 
152  //line 2: set up w
153  ScalarType row_norm = setup_w(A, i, w);
154  ScalarType tau_i = static_cast<ScalarType>(tag.get_drop_tolerance()) * row_norm;
155 
156  //line 3:
157  for (SparseVectorIterator w_k = w.begin(); w_k != w.end(); ++w_k)
158  {
159  SizeType k = w_k->first;
160  if (k >= i)
161  break;
162 
163  //line 4:
164  ScalarType a_kk = output[k][k];
165  if (a_kk == 0)
166  {
167  std::cerr << "ViennaCL: FATAL ERROR in ILUT(): Diagonal entry is zero in row " << k
168  << " while processing line " << i << "!" << std::endl;
169  throw "ILUT zero diagonal!";
170  }
171 
172  ScalarType w_k_entry = w_k->second / a_kk;
173  w_k->second = w_k_entry;
174 
175  //line 5: (dropping rule to w_k)
176  if ( std::fabs(w_k_entry) > tau_i)
177  {
178  //line 7:
179  for (OutputRowConstIterator u_k = output[k].begin(); u_k != output[k].end(); ++u_k)
180  {
181  if (u_k->first > k)
182  w[u_k->first] -= w_k_entry * u_k->second;
183  }
184  }
185  //else
186  // w.erase(k);
187 
188  } //for w_k
189 
190  //Line 10: Apply a dropping rule to w
191  //Sort entries which are kept
192  temp_map.clear();
193  for (SparseVectorIterator w_k = w.begin(); w_k != w.end(); ++w_k)
194  {
195  SizeType k = w_k->first;
196  ScalarType w_k_entry = w_k->second;
197 
198  ScalarType abs_w_k = std::fabs(w_k_entry);
199  if ( (abs_w_k > tau_i) || (k == i) )//do not drop diagonal element!
200  {
201 
202  if (abs_w_k == 0) // this can only happen for diagonal entry
203  throw "Triangular factor in ILUT singular!";
204 
205  temp_map.insert(std::make_pair(abs_w_k, std::make_pair(k, w_k_entry)));
206  }
207  }
208 
209  //Lines 10-12: write the largest p values to L and U
210  SizeType written_L = 0;
211  SizeType written_U = 0;
212  for (typename TemporarySortMap::reverse_iterator iter = temp_map.rbegin(); iter != temp_map.rend(); ++iter)
213  {
214  std::map<SizeType, ScalarType> & row_i = output[i];
215  SizeType j = (iter->second).first;
216  ScalarType w_j_entry = (iter->second).second;
217 
218  if (j < i) // Line 11: entry for L
219  {
220  if (written_L < tag.get_entries_per_row())
221  {
222  row_i[j] = w_j_entry;
223  ++written_L;
224  }
225  }
226  else if (j == i) // Diagonal entry is always kept
227  {
228  row_i[j] = w_j_entry;
229  }
230  else //Line 12: entry for U
231  {
232  if (written_U < tag.get_entries_per_row())
233  {
234  row_i[j] = w_j_entry;
235  ++written_U;
236  }
237  }
238  }
239 
240  w.clear(); //Line 13
241 
242  } //for i
243  }
244 
245 
248  template <typename MatrixType>
250  {
251  typedef typename MatrixType::value_type ScalarType;
252 
253  public:
254  ilut_precond(MatrixType const & mat, ilut_tag const & tag) : tag_(tag), LU(mat.size1(), mat.size2())
255  {
256  //initialize preconditioner:
257  //std::cout << "Start CPU precond" << std::endl;
258  init(mat);
259  //std::cout << "End CPU precond" << std::endl;
260  }
261 
262  template <typename VectorType>
263  void apply(VectorType & vec) const
264  {
265  //Note: Since vec can be a rather arbitrary vector type, we call the more generic version in the backend manually:
266  unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle1());
267  unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU.handle2());
268  ScalarType const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(LU.handle());
269 
270  viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec, LU.size2(), unit_lower_tag());
271  viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec, LU.size2(), upper_tag());
272  }
273 
274  private:
275  void init(MatrixType const & mat)
276  {
279  viennacl::switch_memory_context(temp, host_context);
280 
281  viennacl::copy(mat, temp);
282 
283  std::vector< std::map<unsigned int, ScalarType> > LU_temp(mat.size1());
284 
285  viennacl::linalg::precondition(temp, LU_temp, tag_);
286 
287  viennacl::switch_memory_context(LU, host_context);
288  viennacl::copy(LU_temp, LU);
289  }
290 
291  ilut_tag const & tag_;
293  };
294 
295 
300  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
301  class ilut_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
302  {
304 
305  public:
306  ilut_precond(MatrixType const & mat, ilut_tag const & tag) : tag_(tag), LU(mat.size1(), mat.size2())
307  {
308  //initialize preconditioner:
309  //std::cout << "Start GPU precond" << std::endl;
310  init(mat);
311  //std::cout << "End GPU precond" << std::endl;
312  }
313 
314  void apply(vector<ScalarType> & vec) const
315  {
317  {
318  if (tag_.use_level_scheduling())
319  {
320  //std::cout << "Using multifrontal on GPU..." << std::endl;
322  multifrontal_L_row_index_arrays_,
323  multifrontal_L_row_buffers_,
324  multifrontal_L_col_buffers_,
325  multifrontal_L_element_buffers_,
326  multifrontal_L_row_elimination_num_list_);
327 
328  vec = viennacl::linalg::element_div(vec, multifrontal_U_diagonal_);
329 
331  multifrontal_U_row_index_arrays_,
332  multifrontal_U_row_buffers_,
333  multifrontal_U_col_buffers_,
334  multifrontal_U_element_buffers_,
335  multifrontal_U_row_elimination_num_list_);
336  }
337  else
338  {
340  viennacl::context old_context = viennacl::traits::context(vec);
341  viennacl::switch_memory_context(vec, host_context);
344  viennacl::switch_memory_context(vec, old_context);
345  }
346  }
347  else //apply ILUT directly:
348  {
351  }
352  }
353 
354  private:
355  void init(MatrixType const & mat)
356  {
358  viennacl::switch_memory_context(LU, host_context);
359 
360  std::vector< std::map<unsigned int, ScalarType> > LU_temp(mat.size1());
361 
363  {
364  viennacl::linalg::precondition(mat, LU_temp, tag_);
365  }
366  else //we need to copy to CPU
367  {
368  viennacl::compressed_matrix<ScalarType> cpu_mat(mat.size1(), mat.size2());
369  viennacl::switch_memory_context(cpu_mat, host_context);
370 
371  cpu_mat = mat;
372 
373  viennacl::linalg::precondition(cpu_mat, LU_temp, tag_);
374  }
375 
376  viennacl::copy(LU_temp, LU);
377 
378  if (!tag_.use_level_scheduling())
379  return;
380 
381  //
382  // multifrontal part:
383  //
384 
385  viennacl::switch_memory_context(multifrontal_U_diagonal_, host_context);
386  multifrontal_U_diagonal_.resize(LU.size1(), false);
388 
390  multifrontal_U_diagonal_, //dummy
391  multifrontal_L_row_index_arrays_,
392  multifrontal_L_row_buffers_,
393  multifrontal_L_col_buffers_,
394  multifrontal_L_element_buffers_,
395  multifrontal_L_row_elimination_num_list_);
396 
397 
399  multifrontal_U_diagonal_,
400  multifrontal_U_row_index_arrays_,
401  multifrontal_U_row_buffers_,
402  multifrontal_U_col_buffers_,
403  multifrontal_U_element_buffers_,
404  multifrontal_U_row_elimination_num_list_);
405 
406  //
407  // Bring to device if necessary:
408  //
409 
410  // L:
411 
412  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_row_index_arrays_.begin();
413  it != multifrontal_L_row_index_arrays_.end();
414  ++it)
415  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
416 
417  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_row_buffers_.begin();
418  it != multifrontal_L_row_buffers_.end();
419  ++it)
420  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
421 
422  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_col_buffers_.begin();
423  it != multifrontal_L_col_buffers_.end();
424  ++it)
425  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
426 
427  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_element_buffers_.begin();
428  it != multifrontal_L_element_buffers_.end();
429  ++it)
430  viennacl::backend::switch_memory_context<ScalarType>(*it, viennacl::traits::context(mat));
431 
432 
433  // U:
434 
435  viennacl::switch_memory_context(multifrontal_U_diagonal_, viennacl::traits::context(mat));
436 
437  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_row_index_arrays_.begin();
438  it != multifrontal_U_row_index_arrays_.end();
439  ++it)
440  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
441 
442  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_row_buffers_.begin();
443  it != multifrontal_U_row_buffers_.end();
444  ++it)
445  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
446 
447  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_col_buffers_.begin();
448  it != multifrontal_U_col_buffers_.end();
449  ++it)
450  viennacl::backend::switch_memory_context<unsigned int>(*it, viennacl::traits::context(mat));
451 
452  for (typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_element_buffers_.begin();
453  it != multifrontal_U_element_buffers_.end();
454  ++it)
455  viennacl::backend::switch_memory_context<ScalarType>(*it, viennacl::traits::context(mat));
456 
457 
458  }
459 
460  ilut_tag const & tag_;
462 
463  std::list< viennacl::backend::mem_handle > multifrontal_L_row_index_arrays_;
464  std::list< viennacl::backend::mem_handle > multifrontal_L_row_buffers_;
465  std::list< viennacl::backend::mem_handle > multifrontal_L_col_buffers_;
466  std::list< viennacl::backend::mem_handle > multifrontal_L_element_buffers_;
467  std::list< vcl_size_t > multifrontal_L_row_elimination_num_list_;
468 
469  viennacl::vector<ScalarType> multifrontal_U_diagonal_;
470  std::list< viennacl::backend::mem_handle > multifrontal_U_row_index_arrays_;
471  std::list< viennacl::backend::mem_handle > multifrontal_U_row_buffers_;
472  std::list< viennacl::backend::mem_handle > multifrontal_U_col_buffers_;
473  std::list< viennacl::backend::mem_handle > multifrontal_U_element_buffers_;
474  std::list< vcl_size_t > multifrontal_U_row_elimination_num_list_;
475  };
476 
477  }
478 }
479 
480 
481 
482 
483 #endif
484 
485 
486 
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)
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
void apply(vector< ScalarType > &vec) const
Definition: ilut.hpp:314
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
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
void apply(VectorType &vec) const
Definition: ilut.hpp:263
ilut_precond(MatrixType const &mat, ilut_tag const &tag)
Definition: ilut.hpp:254
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
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
ScalarType setup_w(viennacl::compressed_matrix< ScalarType > const &A, SizeType row, SparseVector &w)
Dispatcher overload for extracting the row of nonzeros of a compressed matrix.
Definition: ilut.hpp:85
void set_drop_tolerance(double tol)
Definition: ilut.hpp:58
A tag class representing an upper triangular matrix.
Definition: forwards.h:708
double get_drop_tolerance() const
Definition: ilut.hpp:63
A tag for incomplete LU factorization with threshold (ILUT)
Definition: ilut.hpp:45
ilut_precond(MatrixType const &mat, ilut_tag const &tag)
Definition: ilut.hpp:306
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
ILUT preconditioner class, can be supplied to solve()-routines.
Definition: ilut.hpp:249
viennacl::memory_types memory_type() const
Definition: context.hpp:76
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.
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:41
bool use_level_scheduling() const
Definition: ilut.hpp:73
Common routines for single-threaded or OpenMP-enabled execution on CPU.
void set_entries_per_row(unsigned int e)
Definition: ilut.hpp:65
Definition: forwards.h:479
ilut_tag(unsigned int entries_per_row=20, double drop_tolerance=1e-4, bool with_level_scheduling=false)
The constructor.
Definition: ilut.hpp:54
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:713
unsigned int get_entries_per_row() const
Definition: ilut.hpp:71
const handle_type & handle() const
Returns the memory handle.
Definition: vector.hpp:878
void use_level_scheduling(bool b)
Definition: ilut.hpp:74
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 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