1 #ifndef VIENNACL_LINALG_DETAIL_BLOCK_ILU_HPP_
2 #define VIENNACL_LINALG_DETAIL_BLOCK_ILU_HPP_
42 template <
typename VectorType,
typename ValueType,
typename SizeType = vcl_
size_t>
52 ) : vec_(v), start_(start_index), size_(vec_size) {}
56 assert(index < size_ &&
bool(
"Index out of bounds!"));
57 return vec_[start_ + index];
62 assert(index < size_ &&
bool(
"Index out of bounds!"));
63 return vec_[start_ + index];
66 SizeType
size()
const {
return size_; }
81 template <
typename ScalarType>
93 ScalarType
const * A_elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(A.
handle());
94 unsigned int const * A_row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.
handle1());
95 unsigned int const * A_col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(A.
handle2());
97 ScalarType * output_elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(diagonal_block_A.
handle());
98 unsigned int * output_row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(diagonal_block_A.
handle1());
99 unsigned int * output_col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(diagonal_block_A.
handle2());
104 unsigned int buffer_col_start = A_row_buffer[
row];
105 unsigned int buffer_col_end = A_row_buffer[
row+1];
107 output_row_buffer[
row - start_index] =
static_cast<unsigned int>(output_counter);
109 for (
unsigned int buf_index = buffer_col_start; buf_index < buffer_col_end; ++buf_index)
111 unsigned int col = A_col_buffer[buf_index];
112 if (col < start_index)
115 if (col >= static_cast<unsigned int>(stop_index))
118 output_col_buffer[output_counter] =
static_cast<unsigned int>(col - start_index);
119 output_elements[output_counter] = A_elements[buf_index];
122 output_row_buffer[
row - start_index + 1] =
static_cast<unsigned int>(output_counter);
134 template <
typename MatrixType,
typename ILUTag>
137 typedef typename MatrixType::value_type ScalarType;
146 ) : tag_(tag), LU_blocks(num_blocks)
150 block_indices_.resize(num_blocks);
153 vcl_size_t start_index = ( i * mat.size1()) / num_blocks;
154 vcl_size_t stop_index = ((i+1) * mat.size1()) / num_blocks;
156 block_indices_[i] = std::pair<vcl_size_t, vcl_size_t>(start_index, stop_index);
167 index_vector_type
const & block_boundaries
168 ) : tag_(tag), block_indices_(block_boundaries), LU_blocks(block_boundaries.
size())
177 template <
typename VectorType>
180 for (
vcl_size_t i=0; i<block_indices_.size(); ++i)
184 unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU_blocks[i].handle1());
185 unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(LU_blocks[i].handle2());
186 ScalarType
const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(LU_blocks[i].handle());
188 viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec_range, LU_blocks[i].size2(),
unit_lower_tag());
189 viennacl::linalg::host_based::detail::csr_inplace_solve<ScalarType>(row_buffer, col_buffer, elements, vec_range, LU_blocks[i].size2(),
upper_tag());
195 void init(MatrixType
const & A)
202 unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(mat.handle1());
204 #ifdef VIENNACL_WITH_OPENMP
205 #pragma omp parallel for
207 for (
long i=0; i<static_cast<long>(block_indices_.size()); ++i)
210 vcl_size_t block_size = block_indices_[i].second - block_indices_[i].first;
211 vcl_size_t block_nnz = row_buffer[block_indices_[i].second] - row_buffer[block_indices_[i].first];
218 preconditioner_dispatch(mat_block, LU_blocks[i], tag_);
235 std::vector< std::map<unsigned int, ScalarType> > temp(mat_block.
size1());
243 index_vector_type block_indices_;
244 std::vector< viennacl::compressed_matrix<ScalarType> > LU_blocks;
255 template <
typename ScalarType,
unsigned int MAT_ALIGNMENT,
typename ILUTag>
269 block_indices_(num_blocks),
274 LU_blocks(num_blocks)
277 block_indices_.resize(num_blocks);
283 block_indices_[i] = std::pair<vcl_size_t, vcl_size_t>(start_index, stop_index);
294 index_vector_type
const & block_boundaries
296 block_indices_(block_boundaries),
301 LU_blocks(block_boundaries.
size())
326 void init(MatrixType
const & A)
333 unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(mat.handle1());
335 #ifdef VIENNACL_WITH_OPENMP
336 #pragma omp parallel for
338 for (
long i=0; i<static_cast<long>(block_indices_.size()); ++i)
341 vcl_size_t block_size = block_indices_[i].second - block_indices_[i].first;
342 vcl_size_t block_nnz = row_buffer[block_indices_[i].second] - row_buffer[block_indices_[i].first];
349 preconditioner_dispatch(mat_block, LU_blocks[i], tag_);
361 for (
vcl_size_t i=0; i<block_indices_.size(); ++i)
363 block_indices_uint.
set(2*i, block_indices_[i].first);
364 block_indices_uint.set(2*i + 1, block_indices_[i].second);
369 blocks_to_device(mat.size1());
376 std::vector< std::map<unsigned int, ScalarType> > L_transposed(matrix_size);
377 std::vector< std::map<unsigned int, ScalarType> > U_transposed(matrix_size);
378 std::vector<ScalarType> entries_D(matrix_size);
383 for (
vcl_size_t block_index = 0; block_index < LU_blocks.size(); ++block_index)
385 MatrixType
const & current_block = LU_blocks[block_index];
387 unsigned int const * row_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(current_block.handle1());
388 unsigned int const * col_buffer = viennacl::linalg::host_based::detail::extract_raw_pointer<unsigned int>(current_block.handle2());
389 ScalarType
const * elements = viennacl::linalg::host_based::detail::extract_raw_pointer<ScalarType>(current_block.handle());
391 vcl_size_t block_start = block_indices_[block_index].first;
396 unsigned int buffer_col_start = row_buffer[
row];
397 unsigned int buffer_col_end = row_buffer[
row+1];
399 for (
unsigned int buf_index = buffer_col_start; buf_index < buffer_col_end; ++buf_index)
401 unsigned int col = col_buffer[buf_index];
404 L_transposed[col + block_start][
static_cast<unsigned int>(
row + block_start)] = elements[buf_index];
406 entries_D[
row + block_start] = elements[buf_index];
408 U_transposed[col + block_start][
static_cast<unsigned int>(
row + block_start)] = elements[buf_index];
416 tools::const_sparse_matrix_adapter<ScalarType, unsigned int> adapted_L_transposed(L_transposed, matrix_size, matrix_size);
417 tools::const_sparse_matrix_adapter<ScalarType, unsigned int> adapted_U_transposed(U_transposed, matrix_size, matrix_size);
435 std::vector< std::map<unsigned int, ScalarType> > temp(mat_block.
size1());
444 index_vector_type block_indices_;
450 std::vector< MatrixType > LU_blocks;
Helper class implementing an array on the host. Default case: No conversion necessary.
Definition: util.hpp:95
std::size_t vcl_size_t
Definition: forwards.h:58
viennacl::enable_if< viennacl::is_any_sparse_matrix< M1 >::value, matrix_expression< const M1, const M1, op_trans > >::type trans(const M1 &mat)
Returns an expression template class representing a transposed matrix.
Definition: sparse_matrix_operations.hpp:330
ilu_vector_range(VectorType &v, SizeType start_index, SizeType vec_size)
Definition: block_ilu.hpp:49
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
This file provides the forward declarations for the main types used within ViennaCL.
std::vector< std::pair< vcl_size_t, vcl_size_t > > index_vector_type
Definition: block_ilu.hpp:140
Implementations of incomplete factorization preconditioners with static nonzero pattern.
A block ILU preconditioner class, can be supplied to solve()-routines.
Definition: block_ilu.hpp:135
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
ValueType & operator[](SizeType index)
Definition: block_ilu.hpp:60
void set(vcl_size_t index, U value)
Definition: util.hpp:145
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
block_ilu_precond(MatrixType const &mat, ILUTag const &tag, vcl_size_t num_blocks=8)
Definition: block_ilu.hpp:265
Main namespace in ViennaCL. Holds all the basic types such as vector, matrix, etc. and defines operations upon them.
Definition: cpu_ram.hpp:29
void apply(vector< ScalarType > &vec) const
Definition: block_ilu.hpp:310
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:144
A tag class representing an upper triangular matrix.
Definition: forwards.h:708
A tag for incomplete LU factorization with threshold (ILUT)
Definition: ilut.hpp:45
std::vector< std::pair< vcl_size_t, vcl_size_t > > index_vector_type
Definition: block_ilu.hpp:262
block_ilu_precond(MatrixType const &mat, ILUTag const &tag, vcl_size_t num_blocks=8)
Definition: block_ilu.hpp:143
void apply(VectorType &vec) const
Definition: block_ilu.hpp:178
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
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
block_ilu_precond(MatrixType const &mat, ILUTag const &tag, index_vector_type const &block_boundaries)
Definition: block_ilu.hpp:292
block_ilu_precond(MatrixType const &mat, ILUTag const &tag, index_vector_type const &block_boundaries)
Definition: block_ilu.hpp:165
ValueType & operator()(SizeType index)
Definition: block_ilu.hpp:54
Implementations of an incomplete factorization preconditioner with threshold (ILUT) ...
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:41
Definition: forwards.h:479
A tag class representing a lower triangular matrix with unit diagonal.
Definition: forwards.h:713
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:62
void extract_block_matrix(viennacl::compressed_matrix< ScalarType > const &A, viennacl::compressed_matrix< ScalarType > &diagonal_block_A, vcl_size_t start_index, vcl_size_t stop_index)
Extracts a diagonal block from a larger system matrix.
Definition: block_ilu.hpp:82
void memory_create(mem_handle &handle, vcl_size_t size_in_bytes, viennacl::context const &ctx, const void *host_ptr=NULL)
Creates an array of the specified size. If the second argument is provided, the buffer is initialized...
Definition: memory.hpp:87
Helper range class for representing a subvector of a larger buffer.
Definition: block_ilu.hpp:43
SizeType size() const
Definition: block_ilu.hpp:66
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
viennacl::enable_if< viennacl::is_any_sparse_matrix< SparseMatrixType >::value >::type block_inplace_solve(const matrix_expression< const SparseMatrixType, const SparseMatrixType, op_trans > &mat, viennacl::backend::mem_handle const &block_index_array, vcl_size_t num_blocks, viennacl::vector_base< ScalarType > const &mat_diagonal, viennacl::vector_base< ScalarType > &vec, SOLVERTAG tag)
Definition: sparse_matrix_operations.hpp:286
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