sparse_mic.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 2018 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_mic_h
#define dealii_sparse_mic_h
 
#include <deal.II/lac/sparse_decomposition.h>
#include <deal.II/lac/sparse_matrix.h>
 
DEAL_II_NAMESPACE_OPEN
 
template <typename number>
class SparseMIC : public SparseLUDecomposition<number>
{
public:
  using size_type = types::global_dof_index;
 
  SparseMIC();
 
  virtual ~SparseMIC() override;
 
  virtual void
  clear() override;
 
  using AdditionalData = typename SparseLUDecomposition<number>::AdditionalData;
 
  template <typename somenumber>
  void
  initialize(const SparseMatrix<somenumber> &matrix,
             const AdditionalData &          parameters = AdditionalData());
 
  template <typename somenumber>
  void
  vmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
 
  template <typename somenumber>
  void
  Tvmult(Vector<somenumber> &dst, const Vector<somenumber> &src) const;
 
  std::size_t
  memory_consumption() const override;
 
  DeclException0(ExcStrengthenDiagonalTooSmall);
  DeclException1(ExcInvalidStrengthening,
                 double,
                 << "The strengthening parameter " << arg1
                 << " is not greater or equal than zero!");
  DeclException2(ExcDecompositionNotStable,
                 int,
                 double,
                 << "The diagonal element (" << arg1 << "," << arg1 << ") is "
                 << arg2 << ", but must be positive");
 
private:
  std::vector<number> diag;
 
  std::vector<number> inv_diag;
 
  std::vector<number> inner_sums;
 
  number
  get_rowsum(const size_type row) const;
};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_