Skip to content

Commit

Permalink
Merge pull request #14244 from peterrum/diagonal_ladv
Browse files Browse the repository at this point in the history 
DiagonalMatrix: specialize for LA:d:V
  • Loading branch information
@kronbichler
kronbichler committed Sep 7, 2022
2 parents ef74666 + 76fc8de commit 6dff4a6
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 42 deletions.
52 changes: 50 additions & 2 deletions include/deal.II/lac/diagonal_matrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,20 @@

DEAL_II_NAMESPACE_OPEN

// forward declarations
#ifndef DOXYGEN
template <typename number>
class Vector;
namespace LinearAlgebra
{
namespace distributed
{
template <typename, typename>
class Vector;
} // namespace distributed
} // namespace LinearAlgebra
#endif

/**
* This class represents a <i>n x n</i> diagonal matrix based on a vector of
* size <i>n</i>. The matrix-vector products are realized by @p
Expand Down Expand Up @@ -386,12 +400,46 @@ DiagonalMatrix<VectorType>::add(const size_type i,



namespace internal
{
namespace DiagonalMatrix
{
template <typename VectorType>
void
assign_and_scale(VectorType & dst,
const VectorType &src,
const VectorType &diagonal)
{
dst = src;
dst.scale(diagonal);
}

template <typename Number>
void
assign_and_scale(
LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> & dst,
const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &src,
const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
&diagonal)
{
const auto dst_ptr = dst.begin();
const auto src_ptr = src.begin();
const auto diagonal_ptr = diagonal.begin();

DEAL_II_OPENMP_SIMD_PRAGMA
for (unsigned int i = 0; i < src.locally_owned_size(); ++i)
dst_ptr[i] = src_ptr[i] * diagonal_ptr[i];
}
} // namespace DiagonalMatrix
} // namespace internal



template <typename VectorType>
void
DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
{
dst = src;
dst.scale(diagonal);
internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
}


Expand Down
83 changes: 43 additions & 40 deletions include/deal.II/lac/precondition.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -542,7 +542,8 @@ namespace internal
typename PreconditionerType>
constexpr bool has_vmult_with_std_functions =
is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value &&
std::is_same<PreconditionerType,
dealii::DiagonalMatrix<VectorType>>::value &&
(std::is_same<VectorType,
dealii::Vector<typename VectorType::value_type>>::value ||
std::is_same<
Expand Down Expand Up @@ -1175,19 +1176,19 @@ namespace internal
// 2) specialized implementation for inverse-diagonal preconditioner
template <typename MatrixType,
typename VectorType,
std::enable_if_t<
!IsBlockVector<VectorType>::value &&
!has_vmult_with_std_functions<MatrixType,
VectorType,
DiagonalMatrix<VectorType>>,
VectorType> * = nullptr>
std::enable_if_t<!IsBlockVector<VectorType>::value &&
!has_vmult_with_std_functions<
MatrixType,
VectorType,
dealii::DiagonalMatrix<VectorType>>,
VectorType> * = nullptr>
void
step_operations(const MatrixType & A,
const DiagonalMatrix<VectorType> &preconditioner,
VectorType & dst,
const VectorType & src,
const double relaxation,
VectorType & tmp,
step_operations(const MatrixType & A,
const dealii::DiagonalMatrix<VectorType> &preconditioner,
VectorType & dst,
const VectorType & src,
const double relaxation,
VectorType & tmp,
VectorType &,
const unsigned int i,
const bool transposed)
Expand Down Expand Up @@ -1247,19 +1248,19 @@ namespace internal
// matrix that accepts ranges
template <typename MatrixType,
typename VectorType,
std::enable_if_t<
!IsBlockVector<VectorType>::value &&
has_vmult_with_std_functions<MatrixType,
VectorType,
DiagonalMatrix<VectorType>>,
VectorType> * = nullptr>
std::enable_if_t<!IsBlockVector<VectorType>::value &&
has_vmult_with_std_functions<
MatrixType,
VectorType,
dealii::DiagonalMatrix<VectorType>>,
VectorType> * = nullptr>
void
step_operations(const MatrixType & A,
const DiagonalMatrix<VectorType> &preconditioner,
VectorType & dst,
const VectorType & src,
const double relaxation,
VectorType & tmp,
step_operations(const MatrixType & A,
const dealii::DiagonalMatrix<VectorType> &preconditioner,
VectorType & dst,
const VectorType & src,
const double relaxation,
VectorType & tmp,
VectorType &,
const unsigned int i,
const bool transposed)
Expand Down Expand Up @@ -2871,15 +2872,16 @@ namespace internal
// selection for diagonal matrix around deal.II vector
template <typename Number>
inline void
vector_updates(const ::dealii::Vector<Number> & rhs,
const DiagonalMatrix<::dealii::Vector<Number>> &jacobi,
const unsigned int iteration_index,
const double factor1,
const double factor2,
::dealii::Vector<Number> &solution_old,
::dealii::Vector<Number> &temp_vector1,
::dealii::Vector<Number> &,
::dealii::Vector<Number> &solution)
vector_updates(
const ::dealii::Vector<Number> & rhs,
const dealii::DiagonalMatrix<::dealii::Vector<Number>> &jacobi,
const unsigned int iteration_index,
const double factor1,
const double factor2,
::dealii::Vector<Number> & solution_old,
::dealii::Vector<Number> & temp_vector1,
::dealii::Vector<Number> &,
::dealii::Vector<Number> &solution)
{
VectorUpdater<Number> upd(rhs.begin(),
jacobi.get_vector().begin(),
Expand Down Expand Up @@ -2909,7 +2911,7 @@ namespace internal
inline void
vector_updates(
const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
const DiagonalMatrix<
const dealii::DiagonalMatrix<
LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
const unsigned int iteration_index,
const double factor1,
Expand Down Expand Up @@ -3053,13 +3055,14 @@ namespace internal
template <typename MatrixType, typename VectorType>
inline void
initialize_preconditioner(
const MatrixType & matrix,
std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
const MatrixType & matrix,
std::shared_ptr<dealii::DiagonalMatrix<VectorType>> &preconditioner)
{
if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
{
if (preconditioner.get() == nullptr)
preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
preconditioner =
std::make_shared<dealii::DiagonalMatrix<VectorType>>();

Assert(
preconditioner->m() == 0,
Expand Down Expand Up @@ -3438,8 +3441,8 @@ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
// We do not need the second temporary vector in case we have a
// DiagonalMatrix as preconditioner and use deal.II's own vectors
using NumberType = typename VectorType::value_type;
if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
false ||
if (std::is_same<PreconditionerType,
dealii::DiagonalMatrix<VectorType>>::value == false ||
(std::is_same<VectorType, dealii::Vector<NumberType>>::value == false &&
((std::is_same<VectorType,
LinearAlgebra::distributed::
Expand Down

0 comments on commit 6dff4a6

Please sign in to comment.