permutation.hh

Go to the documentation of this file.

 
#ifndef permutation_hh
#define permutation_hh
 
#include "matrix.hh"
#include "basics/exceptions.hh"
 
// debugging
#include "basics/debug.hh"
 
#include <type_traits>
#include <typeinfo>
 
#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;
    typedef typename std::conditional<std::is_same<typename Realtype<F>::type, F>::value ,
      typename Realtype<F>::type, typename Cmplxtype<F>::type >::type d_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