Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DiagonalMatrix: specialize for LA:d:V #14244

Merged
merged 1 commit into from
Sep 7, 2022
Merged

DiagonalMatrix: specialize for LA:d:V #14244

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
52 changes: 50 additions & 2 deletions include/deal.II/lac/diagonal_matrix.h
View file
Open in desktop
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
Open in desktop
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