sparse_mic.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_mic_h
#define dealii_sparse_mic_h

#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparse_decomposition.h>

DEAL_II_NAMESPACE_OPEN

template <typename number>
class SparseMIC : public SparseLUDecomposition<number>
{
public:
  typedef types::global_dof_index size_type;

  SparseMIC ();

  virtual ~SparseMIC();

  virtual void clear();

  typedef
  typename SparseLUDecomposition<number>::AdditionalData
  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;

  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_