1 #ifndef VIENNACL_LINALG_ILU_HPP_
2 #define VIENNACL_LINALG_ILU_HPP_
47 double drop_tolerance = 1e-4) : _entries_per_row(entries_per_row), _drop_tolerance(drop_tolerance) {};
52 _drop_tolerance = tol;
65 unsigned int _entries_per_row;
66 double _drop_tolerance;
80 while (row_iter.index1() < k)
91 template <
typename ScalarType>
94 row_iter += k - row_iter.index1();
104 template <
typename ScalarType>
107 row_iter += k - row_iter.index1();
119 template<
typename MatrixType,
typename LUType>
122 typedef std::map<unsigned int, double> SparseVector;
123 typedef typename SparseVector::iterator SparseVectorIterator;
124 typedef typename MatrixType::const_iterator1 InputRowIterator;
125 typedef typename MatrixType::const_iterator2 InputColIterator;
126 typedef typename LUType::iterator1 OutputRowIterator;
127 typedef typename LUType::iterator2 OutputColIterator;
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);
135 std::map<double, unsigned int> temp_map;
137 for (InputRowIterator row_iter = input.begin1(); row_iter != input.end1(); ++row_iter)
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;
148 OutputRowIterator row_iter_out = output.begin1();
149 for (SparseVectorIterator k = w.begin(); k != w.end();)
151 unsigned int index_k = k->first;
152 if (index_k >= static_cast<unsigned int>(row_iter.index1()))
163 double temp = k->second / output(index_k, index_k);
164 if (output(index_k, index_k) == 0.0)
166 std::cerr <<
"ViennaCL: FATAL ERROR in ILUT(): Diagonal entry is zero in row " << index_k <<
"!" << std::endl;
173 for (OutputColIterator j = row_iter_out.begin(); j != row_iter_out.end(); ++j)
175 if (j.index2() > index_k)
177 w[j.index2()] -= temp * *j;
190 for (SparseVectorIterator k = w.begin(); k != w.end(); )
194 long index = k->first;
200 double temp = fabs(k->second);
201 while (temp_map.find(temp) != temp_map.end())
203 temp_map[temp] = k->first;
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)
213 if (iter->second > static_cast<unsigned int>(row_iter.index1()))
217 output(static_cast<unsigned int>(row_iter.index1()), iter->second) =
static_cast<typename LUType::value_type
>(w[iter->second]);
221 else if (iter->second == static_cast<unsigned int>(row_iter.index1()))
223 output(iter->second, iter->second) =
static_cast<typename LUType::value_type
>(w[
static_cast<unsigned int>(row_iter.index1())]);
229 output(static_cast<unsigned int>(row_iter.index1()), iter->second) =
static_cast<typename LUType::value_type
>(w[iter->second]);
243 template<
typename MatrixType,
typename VectorType>
246 typedef typename MatrixType::const_iterator1 InputRowIterator;
247 typedef typename MatrixType::const_iterator2 InputColIterator;
249 for (InputRowIterator row_iter = mat.begin1(); row_iter != mat.end1(); ++row_iter)
251 for (InputColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
253 if (col_iter.index2() < col_iter.index1())
254 vec[col_iter.index1()] -= *col_iter * vec[col_iter.index2()];
264 template<
typename MatrixType,
typename VectorType>
267 typedef typename MatrixType::const_reverse_iterator1 InputRowIterator;
268 typedef typename MatrixType::const_iterator2 InputColIterator;
269 typedef typename VectorType::value_type ScalarType;
271 ScalarType diagonal_entry = 1.0;
273 for (InputRowIterator row_iter = mat.rbegin1(); row_iter != mat.rend1(); ++row_iter)
275 for (InputColIterator col_iter = row_iter.begin(); col_iter != row_iter.end(); ++col_iter)
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;
282 vec[row_iter.index1()] /= diagonal_entry;
291 template<
typename MatrixType,
typename VectorType>
301 template <
typename MatrixType>
304 typedef typename MatrixType::value_type ScalarType;
315 template <
typename VectorType>
323 void init(MatrixType
const & mat)
329 ilut_tag
const & _tag;
330 std::vector< std::map<unsigned int, ScalarType> > LU;
338 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
363 void init(MatrixType
const & mat)
365 std::vector< std::map<unsigned int, ScalarType> > temp(mat.size1());
375 temp_vec.resize(mat.size1());
381 ilut_tag
const & _tag;
383 std::vector< std::map<unsigned int, ScalarType> > LU;
384 mutable std::vector<ScalarType> temp_vec;