1 #ifndef VIENNACL_LINALG_DETAIL_ILUT_HPP_
2 #define VIENNACL_LINALG_DETAIL_ILUT_HPP_
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) {}
61 drop_tolerance_ = tol;
77 unsigned int entries_per_row_;
78 double drop_tolerance_;
79 bool use_level_scheduling_;
84 template <
typename ScalarType,
typename SizeType,
typename SparseVector>
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());
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)
102 ScalarType entry = elements[buf_index_i];
103 w[col_buffer[buf_index_i]] = entry;
104 row_norm += entry * entry;
106 return std::sqrt(row_norm);
110 template <
typename ScalarType,
typename SizeType,
typename SparseVector>
111 ScalarType
setup_w(std::vector< std::map<SizeType, ScalarType> >
const & A,
115 ScalarType row_norm = 0;
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;
120 return std::sqrt(row_norm);
132 template<
typename SparseMatrixType,
typename ScalarType,
typename SizeType>
134 std::vector< std::map<SizeType, ScalarType> > & output,
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;
145 TemporarySortMap temp_map;
153 ScalarType row_norm =
setup_w(A, i, w);
157 for (SparseVectorIterator w_k = w.begin(); w_k != w.end(); ++w_k)
159 SizeType k = w_k->first;
164 ScalarType a_kk = output[k][k];
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!";
172 ScalarType w_k_entry = w_k->second / a_kk;
173 w_k->second = w_k_entry;
176 if ( std::fabs(w_k_entry) > tau_i)
179 for (OutputRowConstIterator u_k = output[k].begin(); u_k != output[k].end(); ++u_k)
182 w[u_k->first] -= w_k_entry * u_k->second;
193 for (SparseVectorIterator w_k = w.begin(); w_k != w.end(); ++w_k)
195 SizeType k = w_k->first;
196 ScalarType w_k_entry = w_k->second;
198 ScalarType abs_w_k = std::fabs(w_k_entry);
199 if ( (abs_w_k > tau_i) || (k == i) )
203 throw "Triangular factor in ILUT singular!";
205 temp_map.insert(std::make_pair(abs_w_k, std::make_pair(k, w_k_entry)));
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)
214 std::map<SizeType, ScalarType> & row_i = output[i];
215 SizeType j = (iter->second).first;
216 ScalarType w_j_entry = (iter->second).second;
222 row_i[j] = w_j_entry;
228 row_i[j] = w_j_entry;
234 row_i[j] = w_j_entry;
248 template <
typename MatrixType>
251 typedef typename MatrixType::value_type ScalarType;
262 template <
typename VectorType>
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());
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());
275 void init(MatrixType
const & mat)
283 std::vector< std::map<unsigned int, ScalarType> > LU_temp(mat.size1());
291 ilut_tag
const & tag_;
300 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT>
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_);
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_);
355 void init(MatrixType
const & mat)
360 std::vector< std::map<unsigned int, ScalarType> > LU_temp(mat.size1());
386 multifrontal_U_diagonal_.resize(LU.size1(),
false);
390 multifrontal_U_diagonal_,
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_);
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_);
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();
417 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_row_buffers_.begin();
418 it != multifrontal_L_row_buffers_.end();
422 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_col_buffers_.begin();
423 it != multifrontal_L_col_buffers_.end();
427 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_L_element_buffers_.begin();
428 it != multifrontal_L_element_buffers_.end();
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();
442 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_row_buffers_.begin();
443 it != multifrontal_U_row_buffers_.end();
447 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_col_buffers_.begin();
448 it != multifrontal_U_col_buffers_.end();
452 for (
typename std::list< viennacl::backend::mem_handle >::iterator it = multifrontal_U_element_buffers_.begin();
453 it != multifrontal_U_element_buffers_.end();
460 ilut_tag
const & tag_;
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_;
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_;
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
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
Definition: forwards.h:696
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