operator.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_OPERATOR
#define MFEM_OPERATOR

#include "vector.hpp"

namespace mfem
{

class Operator
{
protected:
   int height, width;

public:
   explicit Operator(int s = 0) { height = width = s; }

   Operator(int h, int w) { height = h; width = w; }

   inline int Height() const { return height; }
   inline int NumRows() const { return height; }

   inline int Width() const { return width; }
   inline int NumCols() const { return width; }

   virtual void Mult(const Vector &x, Vector &y) const = 0;

   virtual void MultTranspose(const Vector &x, Vector &y) const
   { mfem_error ("Operator::MultTranspose() is not overloaded!"); }

   virtual Operator &GetGradient(const Vector &x) const
   {
      mfem_error("Operator::GetGradient() is not overloaded!");
      return const_cast<Operator &>(*this);
   }

   void PrintMatlab (std::ostream & out, int n = 0, int m = 0) const;

   virtual ~Operator() { }
};


class TimeDependentOperator : public Operator
{
protected:
   double t;

public:
   explicit TimeDependentOperator(int n = 0, double _t = 0.0)
      : Operator(n) { t = _t; }

   TimeDependentOperator(int h, int w, double _t = 0.0)
      : Operator(h, w) { t = _t; }

   virtual double GetTime() const { return t; }

   virtual void SetTime(const double _t) { t = _t; }

   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
   {
      mfem_error("TimeDependentOperator::ImplicitSolve() is not overloaded!");
   }

   virtual ~TimeDependentOperator() { }
};


class Solver : public Operator
{
public:
   bool iterative_mode;

   explicit Solver(int s = 0, bool iter_mode = false)
      : Operator(s) { iterative_mode = iter_mode; }

   Solver(int h, int w, bool iter_mode = false)
      : Operator(h, w) { iterative_mode = iter_mode; }

   virtual void SetOperator(const Operator &op) = 0;
};


class IdentityOperator : public Operator
{
public:
   explicit IdentityOperator(int n) : Operator(n) { }

   virtual void Mult(const Vector &x, Vector &y) const { y = x; }

   ~IdentityOperator() { }
};


class TransposeOperator : public Operator
{
private:
   const Operator &A;

public:
   TransposeOperator(const Operator *a)
      : Operator(a->Width(), a->Height()), A(*a) { }

   TransposeOperator(const Operator &a)
      : Operator(a.Width(), a.Height()), A(a) { }

   virtual void Mult(const Vector &x, Vector &y) const
   { A.MultTranspose(x, y); }

   virtual void MultTranspose(const Vector &x, Vector &y) const
   { A.Mult(x, y); }

   ~TransposeOperator() { }
};


class RAPOperator : public Operator
{
private:
   Operator & Rt;
   Operator & A;
   Operator & P;
   mutable Vector Px;
   mutable Vector APx;

public:
   RAPOperator(Operator &Rt_, Operator &A_, Operator &P_)
      : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_),
        Px(P.Height()), APx(A.Height()) { }

   virtual void Mult(const Vector & x, Vector & y) const
   { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }

   virtual void MultTranspose(const Vector & x, Vector & y) const
   { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }

   virtual ~RAPOperator() { }
};


class TripleProductOperator : public Operator
{
   Operator *A;
   Operator *B;
   Operator *C;
   bool ownA, ownB, ownC;
   mutable Vector t1, t2;

public:
   TripleProductOperator(Operator *A, Operator *B, Operator *C,
                         bool ownA, bool ownB, bool ownC)
      : Operator(A->Height(), C->Width())
      , A(A), B(B), C(C)
      , ownA(ownA), ownB(ownB), ownC(ownC)
      , t1(C->Height()), t2(B->Height())
   {}

   virtual void Mult(const Vector &x, Vector &y) const
   { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }

   virtual void MultTranspose(const Vector &x, Vector &y) const
   { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }

   virtual ~TripleProductOperator()
   {
      if (ownA) { delete A; }
      if (ownB) { delete B; }
      if (ownC) { delete C; }
   }
};

}

#endif