compositions.hh

Go to the documentation of this file.

 
#ifndef compositions_hh
#define compositions_hh
 
#include "basics/outputOperator.hh"
#include "basics/exceptions.hh"
#include "space/space.hh"
#include "function/basis.hh"
#include "function/vector.hh"
#include "basics/debug.hh"
 
// debugging
#define ComposeConstr_D 0
#define ComposeAppl_D 0
#define SchurComplContr_D 0
 
namespace concepts {
 
  // forward declaration
  template<class F>
  class SparseMatrix;
 
  template<class F>
  class DenseMatrix;
 
  template<class F>
  class Matrix;
 
  template<class F>
  class HashedSparseMatrix;
 
  // ************************************************************** Operator **
 
  template<class F>
  class Operator : public virtual OutputOperator {
  public:
    typedef F type;
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;
 
    Operator(uint dimX, uint dimY) : dimX_(dimX), dimY_(dimY) {}
    virtual ~Operator() {}
 
    virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY, Function<c_type>& fncX);
 
    virtual void operator()();
 
    virtual inline const uint dimX() const { return dimX_; }
 
    virtual inline const uint dimY() const { return dimY_; }
 
    virtual void show_messages() { }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    uint dimX_, dimY_;
  };
 
  // *********************************************************** VecOperator **
 
  template<class F>
  class VecOperator : public Operator<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;
 
    VecOperator(uint dimX, uint dimY) : Operator<F>(dimY, dimX) {}
 
    virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX);
 
    virtual void operator()(const Vector<r_type>& fncY, Vector<F>& fncX);
    virtual void operator()(const Vector<c_type>& fncY, Vector<c_type>& fncX);
 
    void operator()(const Matrix<r_type>& mX, Matrix<F>& mY);
    void operator()(const Matrix<c_type>& mX, Matrix<c_type>& mY);
    void operator()();
  protected:
    virtual std::ostream& info(std::ostream& os) const;
    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX) = 0;
    virtual void apply_() = 0;
  };
 
  // *************************************************************** Compose **
 
  template<class F, class H = F>
  class Compose : public Operator<F> {
  public:
    inline Compose(Operator<F>& A, Operator<H>& B);
 
    inline void operator()(const Function<F>& fncY, Function<F>& fncX);
 
    bool collapse(Matrix<F>& dest) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Operator<F>& A_;
    Operator<H>& B_;
    Vector<F> f_;
 
    void collapse_(DenseMatrix<F>& A, DenseMatrix<H>& B, Matrix<F>& dest) const;
  };
 
  template<class F, class H>
  Compose<F,H>::Compose(Operator<F>& A, Operator<H>& B)
    : Operator<F>(A.dimX(), B.dimY()), A_(A), B_(B), f_(B.dimX()) {
    DEBUGL(ComposeConstr_D, 
           "A is (" << A.dimX() << "x" << A.dimY() << ") *" <<
           "B is (" << B.dimX() << "x" << B.dimY() << ")");
    conceptsAssert(A.dimY() == B.dimX(), Assertion());
  }
 
  template<class F, class H>
  void Compose<F,H>::operator()(const Function<F>& fncY, Function<F>& fncX)
  {
    DEBUGL(ComposeAppl_D, *this);
    conceptsAssert3(fncY.dim() == this->dimY_, Assertion(),
                    "fncY.dim = " << fncY.dim() << ", dimY = " << this->dimY());
    conceptsAssert3(fncX.dim() == this->dimX(), Assertion(),
                    "fncX.dim = " << fncX.dim() << ", dimX = " << this->dimX());
    DEBUGL(ComposeAppl_D, "y = " << fncY);
    B_(fncY, f_);
    DEBUGL(ComposeAppl_D, "f = " << f_);
    A_(f_, fncX);
    DEBUGL(ComposeAppl_D, "x = " << fncX);
  }
 
  // ************************************************************** Multiple **
 
  template<class F>
  class Multiple : public Operator<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;
 
    Multiple(Operator<F>& A, F a)
      : Operator<F>(A.dimX(), A.dimY()), A_(A), a_(a) {}
 
    inline void operator()(const Function<r_type>& fncY, Function<F>& fncX)
    { 
      apply_(fncY, fncX);
    }
    inline void operator()(const Function<c_type>& fncY,
                           Function<c_type>& fncX)  { 
      apply_(fncY, fncX);
    }
 
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Operator<F>& A_;
    F a_;
 
    template<class H, class I>
    inline void apply_(const Function<H>& fncY, Function<I>& fncX);
  };
 
  template<class F>
  template<class H, class I>
  void Multiple<F>:: apply_(const Function<H>& fncY, Function<I>& fncX) {
    conceptsAssert(fncY.dim() == this->dimY(), Assertion());
    conceptsAssert(fncX.dim() == this->dimX(), Assertion());
    A_(fncY, fncX);
    if (a_ != 1.0)
      fncX *= a_;
  }
 
  // ***************************************************************** LiCoI **
 
  template<class F>
  class LiCoI : public Operator<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;
 
    LiCoI(Operator<F>& A, F a = 1.0, F b = 1.0)
      : Operator<F>(A.dimX(), A.dimY()), A_(A), a_(a), b_(b) {}
 
    virtual void operator()(const Function<r_type>& fncY, Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX);
 
    bool collapse(Matrix<F>& dest, const F fact = 1.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Operator<F>& A_;
    F a_;
    F b_;
  };
 
  // ****************************************************************** LiCo **
 
  template<class F>
  class LiCo : public Operator<F> {
  public:
    typedef typename Realtype<F>::type r_type;
    typedef typename Cmplxtype<F>::type c_type;
 
    LiCo(Operator<F>& A, Operator<F>& B, F a = 1.0, F b = 1.0) :
      Operator<F>(A.dimX(), A.dimY()), A_(A), B_(B), a_(a), b_(b) {
      conceptsAssert(A.dimX() == B.dimX(), Assertion());
      conceptsAssert(A.dimY() == B.dimY(), Assertion());
    }
 
    inline void operator()(const Function<r_type>& fncY, Function<F>& fncX)
    {
      apply_(fncY, fncX);
    }
    inline void operator()(const Function<c_type>& fncY, Function<c_type>& fncX)
    {
      apply_(fncY, fncX);
    }
 
    bool collapse(Matrix<F>& dest, const F fact = 1.0) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Operator<F>& A_;
    Operator<F>& B_;
    F a_;
    F b_;
 
    template<class H, class I>
    inline void apply_(const Function<H>& fncY, Function<I>& fncX);
  };
 
  template<class F>
  template<class H, class I>
  void LiCo<F>::apply_(const Function<H>& fncY, Function<I>& fncX) {
    conceptsAssert(this->dimX() == fncX.dim(), Assertion());
    conceptsAssert(this->dimY() == fncY.dim(), Assertion());
    A_(fncY, fncX);
    fncX *= a_;
    Vector<I> f(this->dimX());   // Intermediate vector
    B_(fncY, f);
    fncX.add(f, b_);
  }
 
  // ************************************************************ SchurCompl **
 
  template<class F>
  class SchurCompl : public LiCo<F> {
  public:
    SchurCompl(Operator<F>& A_II_inv, Operator<F>& A_IB,
               Operator<F>& A_BI,     Operator<F>& A_BB);
    virtual ~SchurCompl() {
      delete D_; delete C_;
    }
  private:
    Compose<F> *C_, *D_;
  };
 
  template<class F>
  SchurCompl<F>::SchurCompl(Operator<F>& A_II_inv, Operator<F>& A_IB,
                            Operator<F>& A_BI,     Operator<F>& A_BB)
    : LiCo<F>(A_BB, 
              *(D_ = new Compose<F>(A_BI, 
                                    *(C_ = new Compose<F>(A_II_inv, A_IB)))
                ), 1.0, -1.0)
  {
    DEBUGL(SchurComplContr_D, "A_II_inv = " << 
           (DenseMatrix<F>(SparseMatrix<F>(A_II_inv, true))));
    DEBUGL(SchurComplContr_D, "A_IB = " << 
           (DenseMatrix<F>(SparseMatrix<F>(A_IB, true))));
    DEBUGL(SchurComplContr_D, "A_BI = " << 
           (DenseMatrix<F>(SparseMatrix<F>(A_BI, true))));
    DEBUGL(SchurComplContr_D, "A_BB = " << 
           (DenseMatrix<F>(SparseMatrix<F>(A_BB, true))));
    DEBUGL(SchurComplContr_D, "C_ = AII^-1 * AIB = " << 
           (DenseMatrix<F>(SparseMatrix<F>(*C_, true))));
    DEBUGL(SchurComplContr_D, "D_ = ABI * AII^-1 * AIB = " << 
           (DenseMatrix<F>(SparseMatrix<F>(*D_, true))));
    DEBUGL(SchurComplContr_D, "S_ = A_BB - A_BI * A_II^-1 * A_IB = " << 
           (DenseMatrix<F>(SparseMatrix<F>(*this, true))));
  }
 
} // namespace concepts
 
#endif // compositions_hh