precondition_block_base.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_precondition_block_base_h
#define dealii_precondition_block_base_h


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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/householder.h>
#include <deal.II/lac/lapack_full_matrix.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

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

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

  enum Inversion
  {
    gauss_jordan,
    householder,
    svd
  };

  PreconditionBlockBase(bool      store_diagonals = false,
                        Inversion method          = gauss_jordan);

  ~PreconditionBlockBase() = default;

  void
  clear();

  void
  reinit(unsigned int nblocks,
         size_type    blocksize,
         bool         compress,
         Inversion    method = gauss_jordan);

  void
  inverses_computed(bool are_they);

  bool
  same_diagonal() const;

  bool
  store_diagonals() const;

  bool
  inverses_ready() const;

  unsigned int
  size() const;

  template <typename number2>
  void
  inverse_vmult(size_type              i,
                Vector<number2> &      dst,
                const Vector<number2> &src) const;

  template <typename number2>
  void
  inverse_Tvmult(size_type              i,
                 Vector<number2> &      dst,
                 const Vector<number2> &src) const;

  FullMatrix<number> &
  inverse(size_type i);

  Householder<number> &
  inverse_householder(size_type i);

  LAPACKFullMatrix<number> &
  inverse_svd(size_type i);

  const FullMatrix<number> &
  inverse(size_type i) const;

  const Householder<number> &
  inverse_householder(size_type i) const;

  const LAPACKFullMatrix<number> &
  inverse_svd(size_type i) const;

  FullMatrix<number> &
  diagonal(size_type i);

  const FullMatrix<number> &
  diagonal(size_type i) const;

  void
  log_statistics() const;

  std::size_t
  memory_consumption() const;

  DeclException0(ExcDiagonalsNotStored);

  DeclException0(ExcInverseNotAvailable);

protected:
  Inversion inversion;

private:
  unsigned int n_diagonal_blocks;

  std::vector<FullMatrix<number>> var_inverse_full;

  std::vector<Householder<number>> var_inverse_householder;

  std::vector<LAPACKFullMatrix<number>> var_inverse_svd;

  std::vector<FullMatrix<number>> var_diagonal;


  bool var_store_diagonals;

  bool var_same_diagonal;

  bool var_inverses_ready;
};

//----------------------------------------------------------------------//

template <typename number>
inline PreconditionBlockBase<number>::PreconditionBlockBase(bool      store,
                                                            Inversion method)
  : inversion(method)
  , n_diagonal_blocks(0)
  , var_store_diagonals(store)
  , var_same_diagonal(false)
  , var_inverses_ready(false)
{}


template <typename number>
inline void
PreconditionBlockBase<number>::clear()
{
  if (var_inverse_full.size() != 0)
    var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
  if (var_inverse_householder.size() != 0)
    var_inverse_householder.erase(var_inverse_householder.begin(),
                                  var_inverse_householder.end());
  if (var_inverse_svd.size() != 0)
    var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
  if (var_diagonal.size() != 0)
    var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
  var_same_diagonal  = false;
  var_inverses_ready = false;
  n_diagonal_blocks  = 0;
}

template <typename number>
inline void
PreconditionBlockBase<number>::reinit(unsigned int n,
                                      size_type    b,
                                      bool         compress,
                                      Inversion    method)
{
  inversion          = method;
  var_same_diagonal  = compress;
  var_inverses_ready = false;
  n_diagonal_blocks  = n;

  if (compress)
    {
      switch (inversion)
        {
          case gauss_jordan:
            var_inverse_full.resize(1);
            var_inverse_full[0].reinit(b, b);
            break;
          case householder:
            var_inverse_householder.resize(1);
            break;
          case svd:
            var_inverse_svd.resize(1);
            var_inverse_svd[0].reinit(b, b);
            break;
          default:
            Assert(false, ExcNotImplemented());
        }

      if (store_diagonals())
        {
          var_diagonal.resize(1);
          var_diagonal[0].reinit(b, b);
        }
    }
  else
    {
      // set the arrays to the right
      // size. we could do it like this:
      // var_inverse = vector<>(nblocks,FullMatrix<>())
      // but this would involve copying many
      // FullMatrix objects.
      //
      // the following is a neat trick which
      // avoids copying
      if (store_diagonals())
        {
          std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
          var_diagonal.swap(tmp);
        }

      switch (inversion)
        {
          case gauss_jordan:
            {
              std::vector<FullMatrix<number>> tmp(n, FullMatrix<number>(b));
              var_inverse_full.swap(tmp);
              break;
            }
          case householder:
            var_inverse_householder.resize(n);
            break;
          case svd:
            {
              std::vector<LAPACKFullMatrix<number>> tmp(
                n, LAPACKFullMatrix<number>(b));
              var_inverse_svd.swap(tmp);
              break;
            }
          default:
            Assert(false, ExcNotImplemented());
        }
    }
}


template <typename number>
inline unsigned int
PreconditionBlockBase<number>::size() const
{
  return n_diagonal_blocks;
}

template <typename number>
template <typename number2>
inline void
PreconditionBlockBase<number>::inverse_vmult(size_type              i,
                                             Vector<number2> &      dst,
                                             const Vector<number2> &src) const
{
  const size_type ii = same_diagonal() ? 0U : i;

  switch (inversion)
    {
      case gauss_jordan:
        AssertIndexRange(ii, var_inverse_full.size());
        var_inverse_full[ii].vmult(dst, src);
        break;
      case householder:
        AssertIndexRange(ii, var_inverse_householder.size());
        var_inverse_householder[ii].vmult(dst, src);
        break;
      case svd:
        AssertIndexRange(ii, var_inverse_svd.size());
        var_inverse_svd[ii].vmult(dst, src);
        break;
      default:
        Assert(false, ExcNotImplemented());
    }
}


template <typename number>
template <typename number2>
inline void
PreconditionBlockBase<number>::inverse_Tvmult(size_type              i,
                                              Vector<number2> &      dst,
                                              const Vector<number2> &src) const
{
  const size_type ii = same_diagonal() ? 0U : i;

  switch (inversion)
    {
      case gauss_jordan:
        AssertIndexRange(ii, var_inverse_full.size());
        var_inverse_full[ii].Tvmult(dst, src);
        break;
      case householder:
        AssertIndexRange(ii, var_inverse_householder.size());
        var_inverse_householder[ii].Tvmult(dst, src);
        break;
      case svd:
        AssertIndexRange(ii, var_inverse_svd.size());
        var_inverse_svd[ii].Tvmult(dst, src);
        break;
      default:
        Assert(false, ExcNotImplemented());
    }
}


template <typename number>
inline const FullMatrix<number> &
PreconditionBlockBase<number>::inverse(size_type i) const
{
  if (same_diagonal())
    return var_inverse_full[0];

  Assert(i < var_inverse_full.size(),
         ExcIndexRange(i, 0, var_inverse_full.size()));
  return var_inverse_full[i];
}


template <typename number>
inline const Householder<number> &
PreconditionBlockBase<number>::inverse_householder(size_type i) const
{
  if (same_diagonal())
    return var_inverse_householder[0];

  AssertIndexRange(i, var_inverse_householder.size());
  return var_inverse_householder[i];
}


template <typename number>
inline const LAPACKFullMatrix<number> &
PreconditionBlockBase<number>::inverse_svd(size_type i) const
{
  if (same_diagonal())
    return var_inverse_svd[0];

  AssertIndexRange(i, var_inverse_svd.size());
  return var_inverse_svd[i];
}


template <typename number>
inline const FullMatrix<number> &
PreconditionBlockBase<number>::diagonal(size_type i) const
{
  Assert(store_diagonals(), ExcDiagonalsNotStored());

  if (same_diagonal())
    return var_diagonal[0];

  Assert(i < var_diagonal.size(), ExcIndexRange(i, 0, var_diagonal.size()));
  return var_diagonal[i];
}


template <typename number>
inline FullMatrix<number> &
PreconditionBlockBase<number>::inverse(size_type i)
{
  Assert(var_inverse_full.size() != 0, ExcInverseNotAvailable());

  if (same_diagonal())
    return var_inverse_full[0];

  Assert(i < var_inverse_full.size(),
         ExcIndexRange(i, 0, var_inverse_full.size()));
  return var_inverse_full[i];
}


template <typename number>
inline Householder<number> &
PreconditionBlockBase<number>::inverse_householder(size_type i)
{
  Assert(var_inverse_householder.size() != 0, ExcInverseNotAvailable());

  if (same_diagonal())
    return var_inverse_householder[0];

  AssertIndexRange(i, var_inverse_householder.size());
  return var_inverse_householder[i];
}


template <typename number>
inline LAPACKFullMatrix<number> &
PreconditionBlockBase<number>::inverse_svd(size_type i)
{
  Assert(var_inverse_svd.size() != 0, ExcInverseNotAvailable());

  if (same_diagonal())
    return var_inverse_svd[0];

  AssertIndexRange(i, var_inverse_svd.size());
  return var_inverse_svd[i];
}


template <typename number>
inline FullMatrix<number> &
PreconditionBlockBase<number>::diagonal(size_type i)
{
  Assert(store_diagonals(), ExcDiagonalsNotStored());

  if (same_diagonal())
    return var_diagonal[0];

  Assert(i < var_diagonal.size(), ExcIndexRange(i, 0, var_diagonal.size()));
  return var_diagonal[i];
}


template <typename number>
inline bool
PreconditionBlockBase<number>::same_diagonal() const
{
  return var_same_diagonal;
}


template <typename number>
inline bool
PreconditionBlockBase<number>::store_diagonals() const
{
  return var_store_diagonals;
}


template <typename number>
inline void
PreconditionBlockBase<number>::inverses_computed(bool x)
{
  var_inverses_ready = x;
}


template <typename number>
inline bool
PreconditionBlockBase<number>::inverses_ready() const
{
  return var_inverses_ready;
}


template <typename number>
inline void
PreconditionBlockBase<number>::log_statistics() const
{
  deallog << "PreconditionBlockBase: " << size() << " blocks; ";

  if (inversion == svd)
    {
      unsigned int kermin = 100000000, kermax = 0;
      double       sigmin = 1.e300, sigmax = -1.e300;
      double       kappamin = 1.e300, kappamax = -1.e300;

      for (size_type b = 0; b < size(); ++b)
        {
          const LAPACKFullMatrix<number> &matrix = inverse_svd(b);
          size_type                       k      = 1;
          while (k <= matrix.n_cols() &&
                 matrix.singular_value(matrix.n_cols() - k) == 0)
            ++k;
          const double s0 = matrix.singular_value(0);
          const double sm = matrix.singular_value(matrix.n_cols() - k);
          const double co = sm / s0;

          if (kermin > k)
            kermin = k - 1;
          if (kermax < k)
            kermax = k - 1;
          if (s0 < sigmin)
            sigmin = s0;
          if (sm > sigmax)
            sigmax = sm;
          if (co < kappamin)
            kappamin = co;
          if (co > kappamax)
            kappamax = co;
        }
      deallog << "dim ker [" << kermin << ':' << kermax << "] sigma [" << sigmin
              << ':' << sigmax << "] kappa [" << kappamin << ':' << kappamax
              << ']' << std::endl;
    }
  else if (inversion == householder)
    {}
  else if (inversion == gauss_jordan)
    {}
  else
    {
      Assert(false, ExcNotImplemented());
    }
}


template <typename number>
inline std::size_t
PreconditionBlockBase<number>::memory_consumption() const
{
  std::size_t mem = sizeof(*this);
  for (size_type i = 0; i < var_inverse_full.size(); ++i)
    mem += MemoryConsumption::memory_consumption(var_inverse_full[i]);
  for (size_type i = 0; i < var_diagonal.size(); ++i)
    mem += MemoryConsumption::memory_consumption(var_diagonal[i]);
  return mem;
}


DEAL_II_NAMESPACE_CLOSE

#endif