shifted_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__shifted_matrix_h
#define dealii__shifted_matrix_h


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

DEAL_II_NAMESPACE_OPEN

template<typename MatrixType>
class ShiftedMatrix
{
public:
  ShiftedMatrix (const MatrixType &A, const double sigma);

  void shift (const double sigma);

  double shift () const;

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

  template <typename VectorType>
  double residual (VectorType &dst, const VectorType &src, const VectorType &rhs) const;

private:
  SmartPointer<const MatrixType,ShiftedMatrix<MatrixType> > A;

  //    VectorType aux;
  double sigma;
};



template<typename MatrixType, class MatrixType2, class VectorType>
class ShiftedMatrixGeneralized
{
public:
  ShiftedMatrixGeneralized (const MatrixType  &A,
                            const MatrixType2 &M,
                            const double       sigma);

  void shift (const double sigma);

  double shift () const;

  void vmult (VectorType &dst, const VectorType &src) const;

  double residual (VectorType &dst, const VectorType &src, const VectorType &rhs) const;

private:
  SmartPointer<const MatrixType,ShiftedMatrixGeneralized<MatrixType,MatrixType2,VectorType> > A;
  SmartPointer<const MatrixType2,ShiftedMatrixGeneralized<MatrixType,MatrixType2,VectorType> > M;

  mutable VectorType aux;

  double sigma;
};


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

template <typename MatrixType>
inline
ShiftedMatrix<MatrixType>::ShiftedMatrix (const MatrixType &A,
                                          const double      sigma)
  :
  A(&A), sigma(sigma)
{}



template <typename MatrixType>
inline void
ShiftedMatrix<MatrixType>::shift (const double s)
{
  sigma = s;
}


template <typename MatrixType>
inline double
ShiftedMatrix<MatrixType>::shift () const
{
  return sigma;
}



template <typename MatrixType>
template <class VectorType>
inline void
ShiftedMatrix<MatrixType>::vmult (VectorType &dst, const VectorType &src) const
{
  A->vmult(dst, src);
  if (sigma != 0.)
    dst.add(sigma, src);
}


template <typename MatrixType>
template <class VectorType>
inline double
ShiftedMatrix<MatrixType>::residual (VectorType       &dst,
                                     const VectorType &src,
                                     const VectorType &rhs) const
{
  A->vmult(dst, src);
  if (sigma != 0.)
    dst.add(sigma, src);
  dst.sadd(-1.,1.,rhs);
  return dst.l2_norm ();
}


//---------------------------------------------------------------------------
template <typename MatrixType, class MatrixType2, class VectorType>
inline
ShiftedMatrixGeneralized<MatrixType, MatrixType2, VectorType>
::ShiftedMatrixGeneralized (const MatrixType  &A,
                            const MatrixType2 &M,
                            const double       sigma)
  :
  A(&A), M(&M), sigma(sigma)
{}


template <typename MatrixType, class MatrixType2, class VectorType>
inline void
ShiftedMatrixGeneralized<MatrixType, MatrixType2, VectorType>::shift (const double s)
{
  sigma = s;
}

template <typename MatrixType, class MatrixType2, class VectorType>
inline double
ShiftedMatrixGeneralized<MatrixType, MatrixType2, VectorType>::shift () const
{
  return sigma;
}


template <typename MatrixType, class MatrixType2, class VectorType>
inline void
ShiftedMatrixGeneralized<MatrixType, MatrixType2, VectorType>::vmult
(VectorType       &dst,
 const VectorType &src) const
{
  A->vmult(dst, src);
  if (sigma != 0.)
    {
      aux.reinit(dst);
      M->vmult(aux, src);
      dst.add(sigma, aux);
    }
}


template <typename MatrixType, class MatrixType2, class VectorType>
inline double
ShiftedMatrixGeneralized<MatrixType, MatrixType2, VectorType>::residual
(VectorType       &dst,
 const VectorType &src,
 const VectorType &rhs) const
{
  A->vmult(dst, src);
  if (sigma != 0.)
    {
      aux.reinit(dst);
      M->vmult(aux, src);
      dst.add(sigma, aux);
    }
  dst.sadd(-1.,1.,rhs);
  return dst.l2_norm ();
}


DEAL_II_NAMESPACE_CLOSE

#endif