diagonal_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_diagonal_matrix_h
#define dealii_diagonal_matrix_h


#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>

DEAL_II_NAMESPACE_OPEN

template <typename VectorType = Vector<double>>
class DiagonalMatrix : public Subscriptor
{
public:
  using value_type = typename VectorType::value_type;
  using size_type  = typename VectorType::size_type;

  DiagonalMatrix() = default;

  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;

  void
  initialize_dof_vector(VectorType &dst) const;

  std::size_t
  memory_consumption() const;

private:
  VectorType diagonal;
};

/* ---------------------------------- Inline functions ------------------- */

#ifndef DOXYGEN

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;
}



template <typename VectorType>
void
DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
{
  dst = src;
  dst.scale(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);
}


#endif

DEAL_II_NAMESPACE_CLOSE

#endif