sparse_decomposition.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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:
  typedef typename SparseMatrix<number>::size_type size_type;

  virtual ~SparseLUDecomposition () = 0;

  virtual void clear();

  class AdditionalData
  {
  public:
    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