sparse_decomposition.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 2023 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_sparse_decomposition_h
#define dealii_sparse_decomposition_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/lac/sparse_matrix.h>
 
#include <cmath>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename number>
class SparseLUDecomposition : protected SparseMatrix<number>,
                              public virtual Subscriptor
{
protected:
  SparseLUDecomposition();
 
public:
  using size_type = typename SparseMatrix<number>::size_type;
 
  virtual ~SparseLUDecomposition() override = 0;
 
  virtual void
  clear() override;
 
  class AdditionalData
  {
  public:
    explicit AdditionalData(const double       strengthen_diagonal   = 0.,
                            const unsigned int extra_off_diagonals   = 0,
                            const bool         use_previous_sparsity = false,
                            const SparsityPattern *use_this_sparsity = nullptr);
 
    double strengthen_diagonal;
 
    unsigned int extra_off_diagonals;
 
    bool use_previous_sparsity;
 
    const SparsityPattern *use_this_sparsity;
  };
 
  template <typename somenumber>
  void
  initialize(const SparseMatrix<somenumber> &matrix,
             const AdditionalData            parameters);
 
  bool
  empty() const;
 
  size_type
  m() const;
 
  size_type
  n() const;
 
  template <class OutVector, class InVector>
  void
  vmult_add(OutVector &dst, const InVector &src) const;
 
  template <class OutVector, class InVector>
  void
  Tvmult_add(OutVector &dst, const InVector &src) const;
 
  virtual std::size_t
  memory_consumption() const;
 
  DeclException1(ExcInvalidStrengthening,
                 double,
                 << "The strengthening parameter " << arg1
                 << " is not greater or equal than zero!");
protected:
  template <typename somenumber>
  void
  copy_from(const SparseMatrix<somenumber> &matrix);
 
  virtual void
  strengthen_diagonal_impl();
 
  virtual number
  get_strengthen_diagonal(const number rowsum, const size_type row) const;
 
  double strengthen_diagonal;
 
  std::vector<const size_type *> prebuilt_lower_bound;
 
  void
  prebuild_lower_bound();
 
private:
  SparsityPattern *own_sparsity;
};
 
//---------------------------------------------------------------------------
 
#ifndef DOXYGEN
 
template <typename number>
inline number
SparseLUDecomposition<number>::get_strengthen_diagonal(
  const number /*rowsum*/,
  const size_type /*row*/) const
{
  return strengthen_diagonal;
}
 
 
 
template <typename number>
inline bool
SparseLUDecomposition<number>::empty() const
{
  return SparseMatrix<number>::empty();
}
 
 
template <typename number>
inline typename SparseLUDecomposition<number>::size_type
SparseLUDecomposition<number>::m() const
{
  return SparseMatrix<number>::m();
}
 
 
template <typename number>
inline typename SparseLUDecomposition<number>::size_type
SparseLUDecomposition<number>::n() const
{
  return SparseMatrix<number>::n();
}
 
// Note: This function is required for full compatibility with
// the LinearOperator class. ::MatrixInterfaceWithVmultAdd
// picks up the vmult_add function in the protected SparseMatrix
// base class.
template <typename number>
template <class OutVector, class InVector>
inline void
SparseLUDecomposition<number>::vmult_add(OutVector &     dst,
                                         const InVector &src) const
{
  OutVector tmp;
  tmp.reinit(dst);
  this->vmult(tmp, src);
  dst += tmp;
}
 
// Note: This function is required for full compatibility with
// the LinearOperator class. ::MatrixInterfaceWithVmultAdd
// picks up the vmult_add function in the protected SparseMatrix
// base class.
template <typename number>
template <class OutVector, class InVector>
inline void
SparseLUDecomposition<number>::Tvmult_add(OutVector &     dst,
                                          const InVector &src) const
{
  OutVector tmp;
  tmp.reinit(dst);
  this->Tvmult(tmp, src);
  dst += tmp;
}
 
//---------------------------------------------------------------------------
 
 
template <typename number>
SparseLUDecomposition<number>::AdditionalData::AdditionalData(
  const double           strengthen_diag,
  const unsigned int     extra_off_diag,
  const bool             use_prev_sparsity,
  const SparsityPattern *use_this_spars)
  : strengthen_diagonal(strengthen_diag)
  , extra_off_diagonals(extra_off_diag)
  , use_previous_sparsity(use_prev_sparsity)
  , use_this_sparsity(use_this_spars)
{}
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_sparse_decomposition_h