diagonal_matrix.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2016 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_diagonal_matrix_h
#define dealii_diagonal_matrix_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>
 
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
 
template <typename VectorType = Vector<double>>
class DiagonalMatrix : public EnableObserverPointer
{
public:
  using value_type = typename VectorType::value_type;
  using size_type  = typename VectorType::size_type;
 
  DiagonalMatrix() = default;
 
  explicit DiagonalMatrix(const VectorType &vec);
 
  void
  reinit(const VectorType &vec);
 
  void
  compress(VectorOperation::values operation);
 
  VectorType &
  get_vector();
 
  void
  clear();
 
  const VectorType &
  get_vector() const;
 
  size_type
  m() const;
 
  size_type
  n() const;
 
  value_type
  operator()(const size_type i, const size_type j) const;
 
  value_type &
  operator()(const size_type i, const size_type j);
 
  template <typename number2>
  void
  add(const size_type  row,
      const size_type  n_cols,
      const size_type *col_indices,
      const number2   *values,
      const bool       elide_zero_values      = true,
      const bool       col_indices_are_sorted = false);
 
  void
  add(const size_type i, const size_type j, const value_type value);
 
  void
  vmult(VectorType &dst, const VectorType &src) const;
 
  void
  Tvmult(VectorType &dst, const VectorType &src) const;
 
  void
  vmult_add(VectorType &dst, const VectorType &src) const;
 
  void
  Tvmult_add(VectorType &dst, const VectorType &src) const;
 
  value_type
  apply(const unsigned int index, const value_type src) const;
 
  void
  apply_to_subrange(const unsigned int begin_range,
                    const unsigned int end_range,
                    const value_type  *src_pointer_to_current_range,
                    value_type        *dst_pointer_to_current_range) const;
 
  void
  initialize_dof_vector(VectorType &dst) const;
 
  std::size_t
  memory_consumption() const;
 
private:
  VectorType diagonal;
};
 
/* ---------------------------------- Inline functions ------------------- */
 
#ifndef DOXYGEN
 
template <typename VectorType>
DiagonalMatrix<VectorType>::DiagonalMatrix(const VectorType &vec)
  : diagonal(vec)
{}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::clear()
{
  diagonal.reinit(0);
}
 
 
 
template <typename VectorType>
std::size_t
DiagonalMatrix<VectorType>::memory_consumption() const
{
  return diagonal.memory_consumption();
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::reinit(const VectorType &vec)
{
  diagonal = vec;
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::initialize_dof_vector(VectorType &dst) const
{
  dst.reinit(diagonal);
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::compress(VectorOperation::values operation)
{
  diagonal.compress(operation);
}
 
 
 
template <typename VectorType>
VectorType &
DiagonalMatrix<VectorType>::get_vector()
{
  return diagonal;
}
 
 
 
template <typename VectorType>
const VectorType &
DiagonalMatrix<VectorType>::get_vector() const
{
  return diagonal;
}
 
 
 
template <typename VectorType>
typename VectorType::size_type
DiagonalMatrix<VectorType>::m() const
{
  return diagonal.size();
}
 
 
 
template <typename VectorType>
typename VectorType::size_type
DiagonalMatrix<VectorType>::n() const
{
  return diagonal.size();
}
 
 
 
template <typename VectorType>
typename VectorType::value_type
DiagonalMatrix<VectorType>::operator()(const size_type i,
                                       const size_type j) const
{
  Assert(i == j, ExcIndexRange(j, i, i + 1));
  (void)j;
  return diagonal(i);
}
 
 
 
template <typename VectorType>
typename VectorType::value_type &
DiagonalMatrix<VectorType>::operator()(const size_type i, const size_type j)
{
  Assert(i == j, ExcIndexRange(j, i, i + 1));
  (void)j;
  return diagonal(i);
}
 
 
 
template <typename VectorType>
template <typename number2>
void
DiagonalMatrix<VectorType>::add(const size_type  row,
                                const size_type  n_cols,
                                const size_type *col_indices,
                                const number2   *values,
                                const bool,
                                const bool)
{
  for (size_type i = 0; i < n_cols; ++i)
    if (col_indices[i] == row)
      diagonal(row) += values[i];
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::add(const size_type  i,
                                const size_type  j,
                                const value_type value)
{
  if (i == j)
    diagonal(i) += value;
}
 
 
 
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)
    {
      auto *const       dst_ptr      = dst.begin();
      const auto *const src_ptr      = src.begin();
      const auto *const 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
{
  internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
{
  vmult(dst, src);
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::vmult_add(VectorType       &dst,
                                      const VectorType &src) const
{
  VectorType tmp(src);
  tmp.scale(diagonal);
  dst += tmp;
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::Tvmult_add(VectorType       &dst,
                                       const VectorType &src) const
{
  vmult_add(dst, src);
}
 
 
 
template <typename VectorType>
typename VectorType::value_type
DiagonalMatrix<VectorType>::apply(const unsigned int index,
                                  const value_type   src) const
{
  AssertIndexRange(index, diagonal.locally_owned_elements().n_elements());
  return diagonal.local_element(index) * src;
}
 
 
 
template <typename VectorType>
void
DiagonalMatrix<VectorType>::apply_to_subrange(
  const unsigned int begin_range,
  const unsigned int end_range,
  const value_type  *src_pointer_to_current_range,
  value_type        *dst_pointer_to_current_range) const
{
  AssertIndexRange(begin_range,
                   diagonal.locally_owned_elements().n_elements() + 1);
  AssertIndexRange(end_range,
                   diagonal.locally_owned_elements().n_elements() + 1);
 
  const value_type  *diagonal_entry = diagonal.begin() + begin_range;
  const unsigned int length         = end_range - begin_range;
 
  DEAL_II_OPENMP_SIMD_PRAGMA
  for (unsigned int i = 0; i < length; ++i)
    dst_pointer_to_current_range[i] =
      diagonal_entry[i] * src_pointer_to_current_range[i];
}
 
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif