operator/permutation.hh

Go to the documentation of this file.

// permutation and restriction/extension operator

#ifndef permutation_hh
#define permutation_hh

#include "matrix.hh"
#include "basics/exceptions.hh"

// debugging
#include "basics/debug.hh"

#define PermutationAppl_D    0

#define trivExtRestrConstr_D 0
#define trivExtRestrAppl_D   0

namespace concepts {

  // forward declaration
  template<typename F>
  class TrivExtendRestrict;

  // *********************************************************** Permutation **

  template<typename F>
  class Permutation : public Matrix<F> {
  public:
    typedef typename Realtype<F>::type  r_type;
    typedef typename Cmplxtype<F>::type c_type;

    template<class G>
    Permutation(const Space<G>& space,
                const Array<int>& perm, bool transpose = false);

    Permutation(uint dim, const Array<int>& perm, bool transpose = false);
    Permutation(const Permutation<F>& perm, bool transpose = false);

    virtual void operator()(const Function<r_type>& fncY,
                            Function<F>& fncX);
    virtual void operator()(const Function<c_type>& fncY,
                            Function<c_type>& fncX);

    int operator[](const uint i) const { return p_[i]; }

    virtual F operator()(const uint i, const uint j) const;
    virtual F& operator()(const uint i, const uint j);

    virtual void transpMult(const Vector<r_type>& fncY,
                            Vector<F>& fncX);
    virtual void transpMult(const Vector<c_type>& fncY,
                            Vector<c_type>& fncX);

    void composeRestr(const TrivExtendRestrict<F>& restr,
                      Matrix<F>& dest) const;

    void convertToMatrix(Matrix<F>& dest) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    Array<int> p_;
  };

  // **************************************************** TrivExtendRestrict **

  template<typename F>
  class TrivExtendRestrict : public concepts::Operator<F> {
  public:
    template<class G>
    TrivExtendRestrict(const concepts::Space<G>& spcX,
                       const concepts::Space<G>& spcY, bool extend = false)
      : Operator<F>(spcX.dim(), spcY.dim()), extend_(extend)
    {
      DEBUGL(trivExtRestrConstr_D, "extend = " << extend
             << ", dimX = " << this->dimX() << ", dimY = " << this->dimY());
#ifdef DEBUG
      if (extend_) {
        conceptsAssert(spcX.dim() >= spcY.dim(), concepts::Assertion());
      } else {
        conceptsAssert(spcX.dim() <= spcY.dim(), concepts::Assertion());
      }
#endif
    }

    TrivExtendRestrict(uint dimX, uint dimY, bool extend = false)
      : Operator<F>(dimX, dimY), extend_(extend)
    {
      DEBUGL(trivExtRestrConstr_D, "extend = " << extend
             << ", dimX = " << this->dimX() << ", dimY = " << this->dimY());
#ifdef DEBUG
      if (extend_) {
        conceptsAssert(dimX >= dimY, concepts::Assertion());
      } else {
        conceptsAssert(dimX <= dimY, concepts::Assertion());
      }
#endif
    }

    virtual void operator()(const concepts::Function<F>& fncY,
                            concepts::Function<F>& fncX);
    bool extend() const { return extend_; }

    void convertToMatrix(Matrix<F>& dest) const;
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    const bool extend_;
  };

  template<typename F>
  void TrivExtendRestrict<F>::operator()
    (const concepts::Function<F>& fncY, concepts::Function<F>& fncX) {
    conceptsAssert(this->dimX() == fncX.dim(), concepts::Assertion());
    conceptsAssert(this->dimY() == fncY.dim(), concepts::Assertion());
    DEBUGL(trivExtRestrAppl_D, "y = " << fncY);
    const uint diff = this->dimX() > this->dimY() ? 
      this->dimX() - this->dimY() : this->dimY() - this->dimX();
    DEBUGL(trivExtRestrAppl_D, "diff = " << diff << ", extend = " << extend_);
    if (extend_) {
      for (uint i = 0; i < diff; ++i) fncX(i) = 0;
      for (uint i = 0; i < this->dimY(); ++i) fncX(i+diff) = fncY(i);
    } else {
      for (uint i = 0; i < this->dimX(); ++i) fncX(i) = fncY(i+diff);
    }
    DEBUGL(trivExtRestrAppl_D, "x = " << fncX);
  }

} // namespace concepts

#endif // permutation_hh