precondition_block_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2016 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__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/subscriptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/memory_consumption.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:
  typedef types::global_dof_index size_type;

  enum Inversion
  {
    gauss_jordan,
    householder,
    svd
  };

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

  ~PreconditionBlockBase();

  void clear();

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

  void inverses_computed(bool are_they);

  void set_same_diagonal ();

  bool same_diagonal () const;

  bool store_diagonals() const;

  bool inverses_ready () const;

  bool empty () const;

  unsigned int size() const;

  number el(size_type i, size_type j) 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