All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
chow_patel_ilu.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_DETAIL_CHOW_PATEL_ILU_HPP_
#define VIENNACL_LINALG_DETAIL_CHOW_PATEL_ILU_HPP_
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the PDF manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <vector>
#include <cmath>
#include <iostream>
#include "viennacl/forwards.h"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/detail/ilu/common.hpp"
#include "viennacl/linalg/ilu_operations.hpp"
#include "viennacl/linalg/prod.hpp"
#include "viennacl/backend/memory.hpp"
namespace viennacl
{
namespace linalg
{
class chow_patel_tag
{
public:
chow_patel_tag(vcl_size_t num_sweeps = 3, vcl_size_t num_jacobi_iters = 2) : sweeps_(num_sweeps), jacobi_iters_(num_jacobi_iters) {}
vcl_size_t sweeps() const { return sweeps_; }
void sweeps(vcl_size_t num) { sweeps_ = num; }
vcl_size_t jacobi_iters() const { return jacobi_iters_; }
void jacobi_iters(vcl_size_t num) { jacobi_iters_ = num; }
private:
vcl_size_t sweeps_;
vcl_size_t jacobi_iters_;
};
namespace detail
{
template<typename NumericT>
void precondition(viennacl::compressed_matrix<NumericT> const & A,
viennacl::compressed_matrix<NumericT> & L,
viennacl::vector<NumericT> & diag_L,
viennacl::compressed_matrix<NumericT> & L_trans,
chow_patel_tag const & tag)
{
// make sure L and U have correct dimensions:
L.resize(A.size1(), A.size2(), false);
// initialize L and U from values in A:
viennacl::linalg::extract_L(A, L);
// diagonally scale values from A in L:
viennacl::linalg::icc_scale(A, L);
viennacl::vector<NumericT> aij_L(L.nnz(), viennacl::traits::context(A));
viennacl::backend::memory_copy(L.handle(), aij_L.handle(), 0, 0, sizeof(NumericT) * L.nnz());
// run sweeps:
for (vcl_size_t i=0; i<tag.sweeps(); ++i)
viennacl::linalg::icc_chow_patel_sweep(L, aij_L);
// transpose L to obtain L_trans:
viennacl::linalg::ilu_transpose(L, L_trans);
// form (I - D_L^{-1}L) and (I - D_U^{-1} U), with U := L_trans
viennacl::linalg::ilu_form_neumann_matrix(L, diag_L);
viennacl::linalg::ilu_form_neumann_matrix(L_trans, diag_L);
}
template<typename NumericT>
void precondition(viennacl::compressed_matrix<NumericT> const & A,
viennacl::compressed_matrix<NumericT> & L,
viennacl::vector<NumericT> & diag_L,
viennacl::compressed_matrix<NumericT> & U,
viennacl::vector<NumericT> & diag_U,
chow_patel_tag const & tag)
{
// make sure L and U have correct dimensions:
L.resize(A.size1(), A.size2(), false);
U.resize(A.size1(), A.size2(), false);
// initialize L and U from values in A:
viennacl::linalg::extract_LU(A, L, U);
// diagonally scale values from A in L and U:
viennacl::linalg::ilu_scale(A, L, U);
// transpose storage layout of U from CSR to CSC via transposition
viennacl::compressed_matrix<NumericT> U_trans;
viennacl::linalg::ilu_transpose(U, U_trans);
// keep entries of a_ij for the sweeps
viennacl::vector<NumericT> aij_L (L.nnz(), viennacl::traits::context(A));
viennacl::vector<NumericT> aij_U_trans(U_trans.nnz(), viennacl::traits::context(A));
viennacl::backend::memory_copy( L.handle(), aij_L.handle(), 0, 0, sizeof(NumericT) * L.nnz());
viennacl::backend::memory_copy(U_trans.handle(), aij_U_trans.handle(), 0, 0, sizeof(NumericT) * U_trans.nnz());
// run sweeps:
for (vcl_size_t i=0; i<tag.sweeps(); ++i)
viennacl::linalg::ilu_chow_patel_sweep(L, aij_L, U_trans, aij_U_trans);
// transpose U_trans back:
viennacl::linalg::ilu_transpose(U_trans, U);
// form (I - D_L^{-1}L) and (I - D_U^{-1} U)
viennacl::linalg::ilu_form_neumann_matrix(L, diag_L);
viennacl::linalg::ilu_form_neumann_matrix(U, diag_U);
}
}
template<typename MatrixT>
class chow_patel_icc_precond
{
// only works with compressed_matrix!
typedef typename MatrixT::CHOW_PATEL_ICC_ONLY_WORKS_WITH_COMPRESSED_MATRIX error_type;
};
template<typename NumericT, unsigned int AlignmentV>
class chow_patel_icc_precond< viennacl::compressed_matrix<NumericT, AlignmentV> >
{
public:
chow_patel_icc_precond(viennacl::compressed_matrix<NumericT, AlignmentV> const & A, chow_patel_tag const & tag)
: tag_(tag),
L_(0, 0, 0, viennacl::traits::context(A)),
diag_L_(A.size1(), viennacl::traits::context(A)),
L_trans_(0, 0, 0, viennacl::traits::context(A)),
x_k_(A.size1(), viennacl::traits::context(A)),
b_(A.size1(), viennacl::traits::context(A))
{
viennacl::linalg::detail::precondition(A, L_, diag_L_, L_trans_, tag_);
}
template<typename VectorT>
void apply(VectorT & vec) const
{
//
// y = L^{-1} b through Jacobi iteration y_{k+1} = (I - D^{-1}L)y_k + D^{-1}x
//
b_ = viennacl::linalg::element_div(vec, diag_L_);
x_k_ = b_;
for (unsigned int i=0; i<tag_.jacobi_iters(); ++i)
{
vec = viennacl::linalg::prod(L_, x_k_);
x_k_ = vec + b_;
}
//
// x = U^{-1} y through Jacobi iteration x_{k+1} = (I - D^{-1}L^T)x_k + D^{-1}b
//
b_ = viennacl::linalg::element_div(x_k_, diag_L_);
x_k_ = b_; // x_1 if x_0 \equiv 0
for (unsigned int i=0; i<tag_.jacobi_iters(); ++i)
{
vec = viennacl::linalg::prod(L_trans_, x_k_);
x_k_ = vec + b_;
}
// return result:
vec = x_k_;
}
private:
chow_patel_tag tag_;
viennacl::compressed_matrix<NumericT> L_;
viennacl::vector<NumericT> diag_L_;
viennacl::compressed_matrix<NumericT> L_trans_;
mutable viennacl::vector<NumericT> x_k_;
mutable viennacl::vector<NumericT> b_;
};
template<typename MatrixT>
class chow_patel_ilu_precond
{
// only works with compressed_matrix!
typedef typename MatrixT::CHOW_PATEL_ILU_ONLY_WORKS_WITH_COMPRESSED_MATRIX error_type;
};
template<typename NumericT, unsigned int AlignmentV>
class chow_patel_ilu_precond< viennacl::compressed_matrix<NumericT, AlignmentV> >
{
public:
chow_patel_ilu_precond(viennacl::compressed_matrix<NumericT, AlignmentV> const & A, chow_patel_tag const & tag)
: tag_(tag),
L_(0, 0, 0, viennacl::traits::context(A)),
diag_L_(A.size1(), viennacl::traits::context(A)),
U_(0, 0, 0, viennacl::traits::context(A)),
diag_U_(A.size1(), viennacl::traits::context(A)),
x_k_(A.size1(), viennacl::traits::context(A)),
b_(A.size1(), viennacl::traits::context(A))
{
viennacl::linalg::detail::precondition(A, L_, diag_L_, U_, diag_U_, tag_);
}
template<typename VectorT>
void apply(VectorT & vec) const
{
//
// y = L^{-1} b through Jacobi iteration y_{k+1} = (I - D^{-1}L)y_k + D^{-1}x
//
b_ = viennacl::linalg::element_div(vec, diag_L_);
x_k_ = b_;
for (unsigned int i=0; i<tag_.jacobi_iters(); ++i)
{
vec = viennacl::linalg::prod(L_, x_k_);
x_k_ = vec + b_;
}
//
// x = U^{-1} y through Jacobi iteration x_{k+1} = (I - D^{-1}U)x_k + D^{-1}b
//
b_ = viennacl::linalg::element_div(x_k_, diag_U_);
x_k_ = b_; // x_1 if x_0 \equiv 0
for (unsigned int i=0; i<tag_.jacobi_iters(); ++i)
{
vec = viennacl::linalg::prod(U_, x_k_);
x_k_ = vec + b_;
}
// return result:
vec = x_k_;
}
private:
chow_patel_tag tag_;
viennacl::compressed_matrix<NumericT> L_;
viennacl::vector<NumericT> diag_L_;
viennacl::compressed_matrix<NumericT> U_;
viennacl::vector<NumericT> diag_U_;
mutable viennacl::vector<NumericT> x_k_;
mutable viennacl::vector<NumericT> b_;
};
} // namespace linalg
} // namespace viennacl
#endif
viennacl::compressed_matrix::size2
const vcl_size_t & size2() const
Returns the number of columns.
Definition: compressed_matrix.hpp:929
viennacl::linalg::chow_patel_icc_precond< viennacl::compressed_matrix< NumericT, AlignmentV > >::apply
void apply(VectorT &vec) const
Preconditioner application: LL^Tx = b, computed via Ly = b, L^Tx = y using Jacobi iterations...
Definition: chow_patel_ilu.hpp:190
viennacl::linalg::element_div
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)
viennacl::linalg::chow_patel_tag
A tag for incomplete LU and incomplete Cholesky factorization with static pattern (Parallel-ILU0...
Definition: chow_patel_ilu.hpp:44
viennacl::linalg::chow_patel_tag::jacobi_iters
void jacobi_iters(vcl_size_t num)
Sets the number of Jacobi iterations for each triangular 'solve' when applying the preconditioner to ...
Definition: chow_patel_ilu.hpp:62
viennacl::compressed_matrix::size1
const vcl_size_t & size1() const
Returns the number of rows.
Definition: compressed_matrix.hpp:927
prod.hpp
Generic interface for matrix-vector and matrix-matrix products. See viennacl/linalg/vector_operations...
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::linalg::extract_L
void extract_L(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L)
Extracts the lower triangular part L from A.
Definition: ilu_operations.hpp:56
viennacl::linalg::ilu_chow_patel_sweep
void ilu_chow_patel_sweep(compressed_matrix< NumericT > &L, vector< NumericT > const &aij_L, compressed_matrix< NumericT > &U_trans, vector< NumericT > const &aij_U_trans)
Performs one nonlinear relaxation step in the Chow-Patel-ILU (cf. Algorithm 2 in paper) ...
Definition: ilu_operations.hpp:272
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::linalg::chow_patel_icc_precond
Parallel Chow-Patel ILU preconditioner class, can be supplied to solve()-routines.
Definition: chow_patel_ilu.hpp:158
viennacl::linalg::chow_patel_tag::sweeps
void sweeps(vcl_size_t num)
Sets the number of sweeps (i.e. number of nonlinear iterations) in the solver setup stage...
Definition: chow_patel_ilu.hpp:57
viennacl::compressed_matrix::handle
const handle_type & handle() const
Returns the OpenCL handle to the matrix entry array.
Definition: compressed_matrix.hpp:942
viennacl::compressed_matrix::nnz
const vcl_size_t & nnz() const
Returns the number of nonzero entries.
Definition: compressed_matrix.hpp:931
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::context
Represents a generic 'context' similar to an OpenCL context, but is backend-agnostic and thus also su...
Definition: context.hpp:39
ilu_operations.hpp
Implementations of specialized routines for the Chow-Patel parallel ILU preconditioner.
viennacl::linalg::chow_patel_tag::sweeps
vcl_size_t sweeps() const
Returns the number of sweeps (i.e. number of nonlinear iterations) in the solver setup stage...
Definition: chow_patel_ilu.hpp:55
viennacl::linalg::prod
VectorT prod(std::vector< std::vector< T, A1 >, A2 > const &matrix, VectorT const &vector)
Definition: prod.hpp:102
viennacl::linalg::chow_patel_tag::chow_patel_tag
chow_patel_tag(vcl_size_t num_sweeps=3, vcl_size_t num_jacobi_iters=2)
Constructor allowing to set the number of sweeps and Jacobi iterations.
Definition: chow_patel_ilu.hpp:52
viennacl::linalg::icc_scale
void icc_scale(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L)
Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L a...
Definition: ilu_operations.hpp:87
viennacl::linalg::ilu_scale
void ilu_scale(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L, compressed_matrix< NumericT > &U)
Scales the values extracted from A such that A' = DAD has unit diagonal. Updates values from A in L a...
Definition: ilu_operations.hpp:188
viennacl::linalg::extract_LU
void extract_LU(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &L, compressed_matrix< NumericT > &U)
Extracts the lower triangular part L and the upper triangular part U from A.
Definition: ilu_operations.hpp:156
viennacl::linalg::chow_patel_icc_precond< viennacl::compressed_matrix< NumericT, AlignmentV > >::chow_patel_icc_precond
chow_patel_icc_precond(viennacl::compressed_matrix< NumericT, AlignmentV > const &A, chow_patel_tag const &tag)
Definition: chow_patel_ilu.hpp:174
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::vector< NumericT >
viennacl::linalg::ilu_form_neumann_matrix
void ilu_form_neumann_matrix(compressed_matrix< NumericT > &R, vector< NumericT > &diag_R)
Extracts the lower triangular part L and the upper triangular part U from A.
Definition: ilu_operations.hpp:305
viennacl::linalg::ilu_transpose
void ilu_transpose(compressed_matrix< NumericT > const &A, compressed_matrix< NumericT > &B)
Transposition B <- A^T, where the aij-vector is permuted in the same way as the value array in A when...
Definition: ilu_operations.hpp:220
viennacl::backend::memory_copy
void memory_copy(mem_handle const &src_buffer, mem_handle &dst_buffer, vcl_size_t src_offset, vcl_size_t dst_offset, vcl_size_t bytes_to_copy)
Copies 'bytes_to_copy' bytes from address 'src_buffer + src_offset' to memory starting at address 'ds...
Definition: memory.hpp:140
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::linalg::icc_chow_patel_sweep
void icc_chow_patel_sweep(compressed_matrix< NumericT > &L, vector< NumericT > &aij_L)
Performs one nonlinear relaxation step in the Chow-Patel-ICC (cf. Algorithm 3 in paper, but for L rather than U)
Definition: ilu_operations.hpp:121
viennacl::compressed_matrix< NumericT >
common.hpp
Common routines used within ILU-type preconditioners.
viennacl::linalg::detail::precondition
void precondition(viennacl::compressed_matrix< NumericT > const &A, viennacl::compressed_matrix< NumericT > &L, viennacl::vector< NumericT > &diag_L, viennacl::compressed_matrix< NumericT > &L_trans, chow_patel_tag const &tag)
Implementation of the parallel ICC0 factorization, Algorithm 3 in Chow-Patel paper.
Definition: chow_patel_ilu.hpp:77
viennacl::compressed_matrix::resize
void resize(vcl_size_t new_size1, vcl_size_t new_size2, bool preserve=true)
Resize the matrix.
Definition: compressed_matrix.hpp:829
viennacl::linalg::chow_patel_ilu_precond
Parallel Chow-Patel ILU preconditioner class, can be supplied to solve()-routines.
Definition: chow_patel_ilu.hpp:236
memory.hpp
Main interface routines for memory management.
viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix< NumericT, AlignmentV > >::apply
void apply(VectorT &vec) const
Preconditioner application: LUx = b, computed via Ly = b, Ux = y using Jacobi iterations.
Definition: chow_patel_ilu.hpp:269
viennacl::linalg::chow_patel_tag::jacobi_iters
vcl_size_t jacobi_iters() const
Returns the number of Jacobi iterations (i.e. applications of x_{k+1} = (I - D^{-1}R)x_k + D^{-1} b) ...
Definition: chow_patel_ilu.hpp:60
viennacl::linalg::chow_patel_ilu_precond< viennacl::compressed_matrix< NumericT, AlignmentV > >::chow_patel_ilu_precond
chow_patel_ilu_precond(viennacl::compressed_matrix< NumericT, AlignmentV > const &A, chow_patel_tag const &tag)
Definition: chow_patel_ilu.hpp:252