ViennaCL - The Vienna Computing Library  1.2.0
ilu.hpp
Go to the documentation of this file.
1 #ifndef VIENNACL_LINALG_ILU_HPP_
2 #define VIENNACL_LINALG_ILU_HPP_
3 
4 /* =========================================================================
5  Copyright (c) 2010-2011, Institute for Microelectronics,
6  Institute for Analysis and Scientific Computing,
7  TU Wien.
8 
9  -----------------
10  ViennaCL - The Vienna Computing Library
11  -----------------
12 
13  Project Head: Karl Rupp rupp@iue.tuwien.ac.at
14 
15  (A list of authors and contributors can be found in the PDF manual)
16 
17  License: MIT (X11), see file LICENSE in the base directory
18 ============================================================================= */
19 
24 #include <vector>
25 #include <cmath>
26 #include "viennacl/forwards.h"
27 #include "viennacl/tools/tools.hpp"
28 
29 #include <map>
30 
31 namespace viennacl
32 {
33  namespace linalg
34  {
35 
38  class ilut_tag
39  {
40  public:
46  ilut_tag(unsigned int entries_per_row = 20,
47  double drop_tolerance = 1e-4) : _entries_per_row(entries_per_row), _drop_tolerance(drop_tolerance) {};
48 
49  void set_drop_tolerance(double tol)
50  {
51  if (tol > 0)
52  _drop_tolerance = tol;
53  }
54  double get_drop_tolerance() const { return _drop_tolerance; }
55 
56  void set_entries_per_row(unsigned int e)
57  {
58  if (e > 0)
59  _entries_per_row = e;
60  }
61 
62  unsigned int get_entries_per_row() const { return _entries_per_row; }
63 
64  private:
65  unsigned int _entries_per_row;
66  double _drop_tolerance;
67  };
68 
69 
77  template <typename T>
78  void ilut_inc_row_iterator_to_row_index(T & row_iter, unsigned int k)
79  {
80  while (row_iter.index1() < k)
81  ++row_iter;
82  }
83 
91  template <typename ScalarType>
93  {
94  row_iter += k - row_iter.index1();
95  }
96 
104  template <typename ScalarType>
106  {
107  row_iter += k - row_iter.index1();
108  }
109 
110 
119  template<typename MatrixType, typename LUType>
120  void precondition(MatrixType const & input, LUType & output, ilut_tag const & tag)
121  {
122  typedef std::map<unsigned int, double> SparseVector;
123  typedef typename SparseVector::iterator SparseVectorIterator;
124  typedef typename MatrixType::const_iterator1 InputRowIterator; //iterate along increasing row index
125  typedef typename MatrixType::const_iterator2 InputColIterator; //iterate along increasing column index
126  typedef typename LUType::iterator1 OutputRowIterator; //iterate along increasing row index
127  typedef typename LUType::iterator2 OutputColIterator; //iterate along increasing column index
128 
129  output.clear();
130  assert(input.size1() == output.size1());
131  assert(input.size2() == output.size2());
132  output.resize(static_cast<unsigned int>(input.size1()), static_cast<unsigned int>(input.size2()), false);
133  SparseVector w;
134 
135  std::map<double, unsigned int> temp_map;
136 
137  for (InputRowIterator row_iter = input.begin1(); row_iter != input.end1(); ++row_iter)
138  {
139  /* if (i%10 == 0)
140  std::cout << i << std::endl;*/
141 
142  //line 2:
143  w.clear();
144  for (InputColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
145  w[static_cast<unsigned int>(col_iter.index2())] = *col_iter;
146 
147  //line 3:
148  OutputRowIterator row_iter_out = output.begin1();
149  for (SparseVectorIterator k = w.begin(); k != w.end();)
150  {
151  unsigned int index_k = k->first;
152  if (index_k >= static_cast<unsigned int>(row_iter.index1()))
153  break;
154 
155 
156  //while (row_iter_out.index1() < index_k)
157  // ++row_iter_out;
158  //if (row_iter_out.index1() < index_k)
159  // row_iter_out += index_k - row_iter_out.index1();
160  ilut_inc_row_iterator_to_row_index(row_iter_out, index_k);
161 
162  //line 4:
163  double temp = k->second / output(index_k, index_k);
164  if (output(index_k, index_k) == 0.0)
165  {
166  std::cerr << "ViennaCL: FATAL ERROR in ILUT(): Diagonal entry is zero in row " << index_k << "!" << std::endl;
167  }
168 
169  //line 5: (dropping rule to w_k)
170  if ( fabs(temp) > tag.get_drop_tolerance())
171  {
172  //line 7:
173  for (OutputColIterator j = row_iter_out.begin(); j != row_iter_out.end(); ++j)
174  {
175  if (j.index2() > index_k) //attention: manipulation of w(k->first) would invalidate iterator!
176  {
177  w[j.index2()] -= temp * *j;
178  }
179  }
180  ++k; //attention: manipulation of w(k->first) would invalidate iterator!
181  w[index_k] = temp;// - temp * A(index_k, index_k);
182  }
183  else
184  ++k;
185  } //for k
186 
187  //Line 10: Apply a dropping rule to w
188  //Step 1: Sort all entries:
189  temp_map.clear();
190  for (SparseVectorIterator k = w.begin(); k != w.end(); )
191  {
192  if (fabs(k->second) < tag.get_drop_tolerance())
193  {
194  long index = k->first;
195  ++k;
196  w.erase(index);
197  }
198  else
199  {
200  double temp = fabs(k->second);
201  while (temp_map.find(temp) != temp_map.end())
202  temp *= 1.00000001; //make entry slightly larger to maintain uniqueness of the entry
203  temp_map[temp] = k->first;
204  ++k;
205  }
206  }
207 
208  //Lines 10-12: write the largest p values to L and U
209  unsigned int written_L = 0;
210  unsigned int written_U = 0;
211  for (typename std::map<double, unsigned int>::reverse_iterator iter = temp_map.rbegin(); iter != temp_map.rend(); ++iter)
212  {
213  if (iter->second > static_cast<unsigned int>(row_iter.index1())) //entry for U
214  {
215  if (written_U < tag.get_entries_per_row())
216  {
217  output(static_cast<unsigned int>(row_iter.index1()), iter->second) = static_cast<typename LUType::value_type>(w[iter->second]);
218  ++written_U;
219  }
220  }
221  else if (iter->second == static_cast<unsigned int>(row_iter.index1()))
222  {
223  output(iter->second, iter->second) = static_cast<typename LUType::value_type>(w[static_cast<unsigned int>(row_iter.index1())]);
224  }
225  else //entry for L
226  {
227  if (written_L < tag.get_entries_per_row())
228  {
229  output(static_cast<unsigned int>(row_iter.index1()), iter->second) = static_cast<typename LUType::value_type>(w[iter->second]);
230  ++written_L;
231  }
232  }
233  }
234  } //for i
235  }
236 
237 
243  template<typename MatrixType, typename VectorType>
244  void ilu_inplace_solve(MatrixType const & mat, VectorType & vec, viennacl::linalg::unit_lower_tag)
245  {
246  typedef typename MatrixType::const_iterator1 InputRowIterator; //iterate along increasing row index
247  typedef typename MatrixType::const_iterator2 InputColIterator; //iterate along increasing column index
248 
249  for (InputRowIterator row_iter = mat.begin1(); row_iter != mat.end1(); ++row_iter)
250  {
251  for (InputColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
252  {
253  if (col_iter.index2() < col_iter.index1())
254  vec[col_iter.index1()] -= *col_iter * vec[col_iter.index2()];
255  }
256  }
257  }
258 
264  template<typename MatrixType, typename VectorType>
265  void ilu_inplace_solve(MatrixType const & mat, VectorType & vec, viennacl::linalg::upper_tag)
266  {
267  typedef typename MatrixType::const_reverse_iterator1 InputRowIterator; //iterate along increasing row index
268  typedef typename MatrixType::const_iterator2 InputColIterator; //iterate along increasing column index
269  typedef typename VectorType::value_type ScalarType;
270 
271  ScalarType diagonal_entry = 1.0;
272 
273  for (InputRowIterator row_iter = mat.rbegin1(); row_iter != mat.rend1(); ++row_iter)
274  {
275  for (InputColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
276  {
277  if (col_iter.index2() > col_iter.index1())
278  vec[col_iter.index1()] -= *col_iter * vec[col_iter.index2()];
279  if (col_iter.index2() == col_iter.index1())
280  diagonal_entry = *col_iter;
281  }
282  vec[row_iter.index1()] /= diagonal_entry;
283  }
284  }
285 
291  template<typename MatrixType, typename VectorType>
292  void ilu_lu_substitute(MatrixType const & mat, VectorType & vec)
293  {
294  ilu_inplace_solve(mat, vec, unit_lower_tag());
295  ilu_inplace_solve(mat, vec, upper_tag());
296  }
297 
298 
301  template <typename MatrixType>
303  {
304  typedef typename MatrixType::value_type ScalarType;
305 
306  public:
307  ilut_precond(MatrixType const & mat, ilut_tag const & tag) : _tag(tag), LU(mat.size1())
308  {
309  //initialize preconditioner:
310  //std::cout << "Start CPU precond" << std::endl;
311  init(mat);
312  //std::cout << "End CPU precond" << std::endl;
313  }
314 
315  template <typename VectorType>
316  void apply(VectorType & vec) const
317  {
319  viennacl::linalg::ilu_lu_substitute(LU_const_adapter, vec);
320  }
321 
322  private:
323  void init(MatrixType const & mat)
324  {
326  viennacl::linalg::precondition(mat, LU_adapter, _tag);
327  }
328 
329  ilut_tag const & _tag;
330  std::vector< std::map<unsigned int, ScalarType> > LU;
331  };
332 
333 
338  template <typename ScalarType, unsigned int MAT_ALIGNMENT>
339  class ilut_precond< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
340  {
342 
343  public:
344  ilut_precond(MatrixType const & mat, ilut_tag const & tag) : _tag(tag), LU(mat.size1())
345  {
346  //initialize preconditioner:
347  //std::cout << "Start GPU precond" << std::endl;
348  init(mat);
349  //std::cout << "End GPU precond" << std::endl;
350  }
351 
352  void apply(vector<ScalarType> & vec) const
353  {
354  copy(vec, temp_vec);
355  //lu_substitute(LU, vec);
357  viennacl::linalg::ilu_lu_substitute(LU_const_adapter, temp_vec);
358 
359  copy(temp_vec, vec);
360  }
361 
362  private:
363  void init(MatrixType const & mat)
364  {
365  std::vector< std::map<unsigned int, ScalarType> > temp(mat.size1());
366  //std::vector< std::map<unsigned int, ScalarType> > LU_cpu(mat.size1());
367 
368  //copy to cpu:
369  copy(mat, temp);
370 
373  viennacl::linalg::precondition(temp_adapter, LU_adapter, _tag);
374 
375  temp_vec.resize(mat.size1());
376 
377  //copy resulting preconditioner back to gpu:
378  //copy(LU_cpu, LU);
379  }
380 
381  ilut_tag const & _tag;
382  //MatrixType LU;
383  std::vector< std::map<unsigned int, ScalarType> > LU;
384  mutable std::vector<ScalarType> temp_vec;
385  };
386 
387  }
388 }
389 
390 
391 
392 
393 #endif
394 
395 
396