sparse_vanka.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_sparse_vanka_h
#define dealii_sparse_vanka_h



#include <deal.II/base/config.h>

#include <deal.II/base/multithread_info.h>
#include <deal.II/base/smartpointer.h>

#include <map>
#include <vector>

DEAL_II_NAMESPACE_OPEN

template <typename number>
class FullMatrix;
template <typename number>
class SparseMatrix;
template <typename number>
class Vector;

template <typename number>
class SparseVanka;
template <typename number>
class SparseBlockVanka;

template <typename number>
class SparseVanka
{
public:
  using size_type = types::global_dof_index;

  SparseVanka();

  SparseVanka(const SparseMatrix<number> &M,
              const std::vector<bool> &   selected,
              const bool                  conserve_memory = false,
              const unsigned int n_threads = MultithreadInfo::n_threads());

  ~SparseVanka();

  class AdditionalData
  {
  public:
    AdditionalData(const std::vector<bool> &selected,
                   const bool               conserve_memory = false,
                   const unsigned int n_threads = MultithreadInfo::n_threads());

    const std::vector<bool> &selected;

    const bool conserve_mem;

    const unsigned int n_threads;
  };


  void
  initialize(const SparseMatrix<number> &M,
             const AdditionalData &      additional_data);

  template <typename number2>
  void
  vmult(Vector<number2> &dst, const Vector<number2> &src) const;

  template <typename number2>
  void
  Tvmult(Vector<number2> &dst, const Vector<number2> &src) const;

  size_type
  m() const;

  size_type
  n() const;

protected:
  template <typename number2>
  void
  apply_preconditioner(Vector<number2> &              dst,
                       const Vector<number2> &        src,
                       const std::vector<bool> *const dof_mask = nullptr) const;

  std::size_t
  memory_consumption() const;

private:
  SmartPointer<const SparseMatrix<number>, SparseVanka<number>> matrix;

  bool conserve_mem;

  const std::vector<bool> *selected;

  unsigned int n_threads;

  mutable std::vector<SmartPointer<FullMatrix<float>, SparseVanka<number>>>
    inverses;

  size_type _m;

  size_type _n;

  void
  compute_inverses();

  void
  compute_inverses(const size_type begin, const size_type end);

  void
  compute_inverse(const size_type row, std::vector<size_type> &local_indices);

  template <typename T>
  friend class SparseBlockVanka;
};



template <typename number>
class SparseBlockVanka : public SparseVanka<number>
{
public:
  using size_type = types::global_dof_index;

  enum BlockingStrategy
  {
    index_intervals,
    adaptive
  };

  SparseBlockVanka(const SparseMatrix<number> &M,
                   const std::vector<bool> &   selected,
                   const unsigned int          n_blocks,
                   const BlockingStrategy      blocking_strategy,
                   const bool                  conserve_memory = false,
                   const unsigned int n_threads = MultithreadInfo::n_threads());

  template <typename number2>
  void
  vmult(Vector<number2> &dst, const Vector<number2> &src) const;

  std::size_t
  memory_consumption() const;

private:
  const unsigned int n_blocks;

  std::vector<std::vector<bool>> dof_masks;

  void
  compute_dof_masks(const SparseMatrix<number> &M,
                    const std::vector<bool> &   selected,
                    const BlockingStrategy      blocking_strategy);
};

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

#ifndef DOXYGEN

template <typename number>
inline typename SparseVanka<number>::size_type
SparseVanka<number>::m() const
{
  Assert(_m != 0, ExcNotInitialized());
  return _m;
}

template <typename number>
inline typename SparseVanka<number>::size_type
SparseVanka<number>::n() const
{
  Assert(_n != 0, ExcNotInitialized());
  return _n;
}

template <typename number>
template <typename number2>
inline void
SparseVanka<number>::Tvmult(Vector<number2> & /*dst*/,
                            const Vector<number2> & /*src*/) const
{
  AssertThrow(false, ExcNotImplemented());
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif