schur_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2015 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__schur_matrix_h
#define dealii__schur_matrix_h

#include <deal.II/base/config.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/block_vector.h>
#include <vector>

DEAL_II_NAMESPACE_OPEN


template <class MA_inverse, class MB, class MDt, class MC>
class SchurMatrix : public Subscriptor
{
public:

  SchurMatrix(const MA_inverse &Ainv,
              const MB &B,
              const MDt &Dt,
              const MC &C,
              VectorMemory<BlockVector<double> > &mem,
              const std::vector<types::global_dof_index> &signature = std::vector<types::global_dof_index>(0));

  void prepare_rhs (BlockVector<double> &dst,
                    const BlockVector<double> &src) const;

  void vmult (BlockVector<double> &dst,
              const BlockVector<double> &src) const;

//  void Tmult(BlockVector<double>& dst, const BlockVector<double>& src) const;

  double residual (BlockVector<double> &dst,
                   const BlockVector<double> &src,
                   const BlockVector<double> &rhs) const;

  void postprocess (BlockVector<double> &dst,
                    const BlockVector<double> &src,
                    const BlockVector<double> &rhs) const;

  void debug_level(unsigned int l);
private:
  SchurMatrix (const SchurMatrix<MA_inverse, MB, MDt, MC> &);
  SchurMatrix &operator = (const SchurMatrix<MA_inverse, MB, MDt, MC> &);

  const SmartPointer<const MA_inverse,SchurMatrix<MA_inverse,MB,MDt,MC> > Ainv;
  const SmartPointer<const MB,SchurMatrix<MA_inverse,MB,MDt,MC> > B;
  const SmartPointer<const MDt,SchurMatrix<MA_inverse,MB,MDt,MC> > Dt;
  const SmartPointer<const MC,SchurMatrix<MA_inverse,MB,MDt,MC> > C;
  VectorMemory<BlockVector<double> > &mem;

  std::vector<types::global_dof_index> signature;

  unsigned int debug;
};

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

template <class MA_inverse, class MB, class MDt, class MC>
SchurMatrix<MA_inverse, MB, MDt, MC>
::SchurMatrix(const MA_inverse &Ainv,
              const MB &B,
              const MDt &Dt,
              const MC &C,
              VectorMemory<BlockVector<double> > &mem,
              const std::vector<types::global_dof_index> &signature)
  : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C),
    mem(mem),
    signature(signature),
    debug(0)
{
}


template <class MA_inverse, class MB, class MDt, class MC>
void
SchurMatrix<MA_inverse, MB, MDt, MC>
::debug_level(unsigned int l)
{
  debug = l;
}


template <class MA_inverse, class MB, class MDt, class MC>
void SchurMatrix<MA_inverse, MB, MDt, MC>
::vmult(BlockVector<double> &dst,
        const BlockVector<double> &src) const
{
  deallog.push("Schur");
  if (debug > 0)
    deallog << "src:" << src.l2_norm() << std::endl;

  C->vmult(dst, src);
  if (debug > 0)
    deallog << "C:" << dst.l2_norm() << std::endl;

  BlockVector<double> *h1 = mem.alloc();
  if (signature.size()>0)
    h1->reinit(signature);
  else
    h1->reinit(B->n_block_cols(), src.block(0).size());
  Dt->Tvmult(*h1,src);
  if (debug > 0)
    deallog << "Dt:" << h1->l2_norm() << std::endl;

  BlockVector<double> *h2 = mem.alloc();
  h2->reinit(*h1);
  Ainv->vmult(*h2, *h1);
  if (debug > 0)
    deallog << "Ainverse:" << h2->l2_norm() << std::endl;

  mem.free(h1);
  B->vmult_add(dst, *h2);
  if (debug > 0)
    deallog << "dst:" << dst.l2_norm() << std::endl;

  mem.free(h2);
  deallog.pop();
}


template <class MA_inverse, class MB, class MDt, class MC>
double SchurMatrix<MA_inverse, MB, MDt, MC>
::residual(BlockVector<double> &dst,
           const BlockVector<double> &src,
           const BlockVector<double> &rhs) const
{
  vmult(dst, src);
  dst *= -1.;
  dst += rhs;
  return dst.l2_norm();
}


template <class MA_inverse, class MB, class MDt, class MC>
void SchurMatrix<MA_inverse, MB, MDt, MC>
::prepare_rhs(BlockVector<double> &dst,
              const BlockVector<double> &src) const
{
  Assert (src.n_blocks() == B->n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), B->n_block_cols()));
  Assert (dst.n_blocks() == B->n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), B->n_block_rows()));

  deallog.push("Schur-prepare");
  if (debug > 0)
    deallog << "src:" << src.l2_norm() << std::endl;
  BlockVector<double> *h1 = mem.alloc();
  if (signature.size()>0)
    h1->reinit(signature);
  else
    h1->reinit(B->n_block_cols(), src.block(0).size());
  Ainv->vmult(*h1, src);
  if (debug > 0)
    deallog << "Ainverse:" << h1->l2_norm() << std::endl;
  B->vmult_add(dst, *h1);
  if (debug > 0)
    deallog << "dst:" << dst.l2_norm() << std::endl;
  mem.free(h1);
  deallog.pop();
}


template <class MA_inverse, class MB, class MDt, class MC>
void SchurMatrix<MA_inverse, MB, MDt, MC>
::postprocess(BlockVector<double> &dst,
              const BlockVector<double> &src,
              const BlockVector<double> &rhs) const
{
  Assert (dst.n_blocks() == B->n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), B->n_block_cols()));
  Assert (rhs.n_blocks() == B->n_block_cols(),
          ExcDimensionMismatch(rhs.n_blocks(), B->n_block_cols()));
  Assert (src.n_blocks() == B->n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), B->n_block_rows()));

  deallog.push("Schur-post");
  if (debug > 0)
    deallog << "src:" << src.l2_norm() << std::endl;
  BlockVector<double> *h1 = mem.alloc();
  if (signature.size()>0)
    h1->reinit(signature);
  else
    h1->reinit(B->n_block_cols(), src.block(0).size());
  Dt->Tvmult(*h1, src);
  if (debug > 0)
    deallog << "Dt:" << h1->l2_norm() << std::endl;
  h1->sadd(-1.,rhs);
  Ainv->vmult(dst,*h1);
  if (debug > 0)
    deallog << "dst:" << dst.l2_norm() << std::endl;
  mem.free(h1);
  deallog.pop();
}


DEAL_II_NAMESPACE_CLOSE

#endif