vectorsMatrices.hh

Go to the documentation of this file.

 
#ifndef vectMatr_hh
#define vectMatr_hh
 
#include <iostream>
#include <cmath>
#include <cstring>
#include <type_traits>
 
#include "functionCompiler.hh"
#include "exceptions.hh"
#include "typedefs.hh"
#include "vectorsMatricesForward.hh"
#include "debug.hh"
#include <type_traits>
 
// debugging Mapping
#define MappingAppll_D 0
#define MappingAppl_D 0
 
#define MappingConstructorHack 1
 
namespace concepts {
 
  template<typename F, typename G>
  void memorycpy(F *dest, const G *src, size_t n) {
    for (size_t i = 0; i < n; ++i) *dest++ = *src++;
  }
 
  // ***************************************************************** Point **
 
  template<class F, uint dim>
  class Point {
  public:
    Point() { F* d = data_; for(uint i = dim; i--;) *d++ = 0; }
    template<class G>
    Point(const Point<G, dim>& p) {
      memorycpy(data_, (const G*)p, dim);
    }
    template<uint dim2>
    Point(const Point<F, dim2>& p) {
      if (dim2<dim) {
        F* d = data_+dim2; for(uint i = dim-dim2; i--;) *d++ = 0;
      }
      memorycpy(data_, (const F*)p, ((dim < dim2) ? dim : dim2));
    }
    Point(F x) { F* d = data_; for(uint i = dim; i--;) *d++ = x; }
    Point(const F e0, const F e1, const F e2 = 0) {
      uint i = 0;
      data_[i++] = e0;
      if (dim > 1)
        data_[i++] = e1;
      if (dim > 2)
        data_[i] = e2;
    }
    Point(const uchar *pgm, const Real x, const Real y, const Real z = 0.0);
 
    template<uint dim2>
    Point<F, dim>& operator=(const Point<F, dim2>& b) {
      if (this != &b)
        memorycpy(data_, b.data_, ((dim < dim2) ? dim : dim2));
      return *this;
    }
 
    const F& operator[](const uint i) const {
      conceptsAssert(i < dim, concepts::Assertion());
      return data_[i]; }
    F& operator[](const uint i) {
      conceptsAssert(i < dim, concepts::Assertion());
      return data_[i]; }
 
    operator F*() { return data_; }
    operator const F*() const { return data_; }
 
    Point<F, dim>& operator+=(const Point<F, dim>& p) {
      F* d = data_; const F* dp = p.data_;
      for(uint i = dim; i--;) *d++ += *dp++;
      return *this;
    }
    Point<F, dim> operator+ (const Point<F, dim>& p) const {
      Point<F, dim> res(*this);
      return res += p;
    }
    Point<F, dim>& operator-=(const Point<F, dim>& p) {
      F* d = data_; const F* dp = p.data_;
      for(uint i = dim; i--;) *d++ -= *dp++;
      return *this;
    }
    Point<F, dim> operator- (const Point<F, dim>& p) const {
      Point<F, dim> res(*this);
      return res -= p;
    }
 
    F operator*(const Point<F, dim>& b) const;
    F operator*(const UnitNd<dim>& b) const;
    F dot(const Point<F, dim>& b) const;
 
    template<class G>
    Point<F, dim>& operator*=(const G n) {
      F* d = data_; 
      for(uint i = dim; i--;) *d++ *= n;
      return *this;
    }
 
    template<class G>
    Point<F, dim>& operator/=(const G n) {
      F* d = data_; for(uint i = dim; i--;) *d++ /= n;
      return *this;
    }
 
    Point<Cmplx, dim> operator*(const Cmplx n) const {
      Point<Cmplx, dim> res(*this);
      return res *= n;
    }
    Point<F, dim> operator*(const Real n) const {
      Point<F, dim> res(*this);
      return res *= n;
    }
    
    Point<Cmplx, dim> operator/(const Cmplx n) const {
      Point<Cmplx, dim> res(*this);
      return res /= n;
    }
    Point<F, dim> operator/(const Real n) const {
      Point<F, dim> res(*this);
      return res /= n;
    }
 
    Point<F, dim>& scale(const Point<F, dim>& n) {
      F* d = data_; const F* e = n.data_;
      for(uint i = dim; i--;) *d++ *= *e++;
      return *this;
    }
 
    Point<F, 3> operator^(const Point<F, 3>& b) const;
    F operator^(const Point<F, 2>& b) const;
 
    Real l2() const { return std::sqrt(l2_2()); }
    Real l2_2() const;
    Real linfty() const;
 
    void lincomb(const Point<F, dim>& a, const Point<F, dim>& b,
                 const F ca = 1.0, const F cb = 1.0);
 
    void max(const Point<F, dim>& a, const Point<F, dim>& b);
 
    void min(const Point<F, dim>& a, const Point<F, dim>& b);
 
    Point<F, dim>& ortho(const Point<F, dim>& a);
    Point<F, dim>& ortho();
    Point<F, dim> ortho() const;
 
    std::ostream& info(std::ostream& os) const;
  private:
    F data_[dim];
  };
 
  template<class F, uint dim>
  std::ostream& operator<<(std::ostream& os, const Point<F, dim>& p) {
    return p.info(os);
  }
 
  // set the real type of a vector to the real type of the components
  template<typename F, uint dim>
  struct Realtype<Point<F,dim> > {
    typedef typename Realtype<F>::type type;
  };
 
  // set the data type of a vector to the data type of the components
  template<typename F, uint dim>
  struct Datatype<Point<F,dim> > {
    //typedef typename Datatype<F>::type type;
    typedef F type;
  };
 
  template <class F, uint dim>
  inline bool operator==(const Point<F,dim>& x, const Point<F,dim>& y)
  {
    const F* d = (const F*)x; const F* e = (const F*)y;
    for (uint i = dim; i--;)
      if (*d++ != *e++) return false;
    return true;
  }
 
  template <class F, uint dim>
  inline Point<typename Combtype<F,Real>::type,dim> operator*(const Real x, const Point<F,dim>& y)
  {
    return y * x;
  }
 
  template <class F, uint dim>
  inline Point<typename Combtype<F,Cmplx>::type,dim> operator*(const Cmplx x, const Point<F,dim>& y)
  {
    return y * x;
  }
 
  template<uint dim>
  Cmplx operator*(const Point<Cmplx, dim>& a, const Point<Real, dim>& b) {
    Cmplx res = 0.0;
    const Cmplx* d = (const Cmplx*)a; const Real* db = (const Real*)b;
    for(uint i = 2; i--;) 
      res += (*d++) * (*db++);
    return res;
  }
 
  template<uint dim>
  Cmplx operator*(const Point<Real, dim>& a, const Point<Cmplx, dim>& b) {
    Cmplx res = 0.0;
    const Real* d = (const Real*)a; const Cmplx* db = (const Cmplx*)b;
    for(uint i = dim; i--;) 
      res += (*d++) * (*db++);
    return res;
  }
 
  // **************************************************************** UnitNd **
 
  template <uint dim>
  class UnitNd : public Point<Real, dim> {
    Real len_;
 
  public:
    UnitNd() : Point<Real, dim>() {len_ = 0.0;}
    inline UnitNd(const Point<Real, dim>& u);
    inline UnitNd(Real x, Real y, Real z);
 
    Real length() const {return len_;}
  };
 
  template <uint dim>
  inline UnitNd<dim>::UnitNd(const Point<Real, dim>& u)
    : Point<Real, dim>(u) {
    len_ = this->l2();  *this *= (1.0 / len_);
  }
 
  template <uint dim>
  inline UnitNd<dim>::UnitNd(Real x, Real y, Real z)
    : Point<Real, dim>(x, y, z) {
    len_ = this->l2();  *this *= (1.0 / len_);
  }
 
  // ************************************************ Mapping<F, DimY, DimX> **  
 
  template<class F, uint DimY, uint DimX=DimY>
  class Mapping {
  public:
    Mapping() {}
    
    Mapping(const Mapping<F, DimY, DimX>& m) {
      memorycpy(data_, m.data_, DimX*DimY);
    }
 
    template <class H>
    Mapping(const Mapping<H, DimY, DimX>& m) {
      for (uint i=0; i < DimY; ++i) {
        for (uint j=0; j < DimX; ++j) {
          data_[lin_(i,j)] = m(i,j);
        }
      }
    }
 
#if MappingConstructorHack
    template <class G>
    Mapping(G x) { F* d = data_; for(uint i = DimX*DimY; i--;) *d++ = x; }
#else
    Mapping(F x) { F* d = data_; for(uint i = DimX*DimY; i--;) *d++ = x; }
#endif
    
#if MappingConstructorHack
    template<class... FieldTs>
    Mapping(FieldTs... data) : data_{data...} {};
#else
    template<class... FieldTs>   
    Mapping(std::enable_if<std::is_floating_point<FieldTs>::value, FieldTs>::type... data) : data_{data...} {
      static_assert(sizeof...(FieldTs)==DimX*DimY,"Wrong number of arguments.");    
    }
#endif
    
    //     /// Constructor for a 2D mapping
    //     Mapping(const F e0, const F e1,
    //             const F e2, const F e3);
    
    //     /// Constructor for a 3D mapping
    //     Mapping(const F e0, const F e1, const F e2,
    //             const F e3, const F e4, const F e5,
    //             const F e6, const F e7, const F e8);
    
    //     /// Constructor for a 6D mapping
    //     Mapping(const F e0, const F e1, const F e2,
    //             const F e3, const F e4, const F e5,
    //             const F e6, const F e7, const F e8, 
    //             const F e9, const F e10, const F e11, 
    //             const F e12, const F e13, const F e14, 
    //             const F e15, const F e16, const F e17, 
    //             const F e18, const F e19, const F e20, 
    //             const F e21, const F e22, const F e23, 
    //             const F e24, const F e25, const F e26, 
    //             const F e27, const F e28, const F e29, 
    //             const F e30, const F e31, const F e32, 
    //             const F e33, const F e34, const F e35);
    
    Mapping(const Point<F, DimY> first);
    
    Mapping(const Point<F, DimY> first, 
            const Point<F, DimY> second);
    
    Mapping(const Point<F, DimY> first,
            const Point<F, DimY> second,
            const Point<F, DimY> third);
 
    Mapping(const UnitNd<DimY>& n);
    
    
    F determinant() const;
    
    Mapping<F, DimY, DimX> adjugate() const;
 
    Mapping<F, DimY-1, DimX-1> subMatrix(uint i, uint j) const;
 
    F trace() const;
    
    Mapping<F, DimX, DimY> inverse() const;
    
    const F& operator()(uint i, uint j) const {
      conceptsAssert(i < DimY, Assertion());
      conceptsAssert(j < DimX, Assertion());
      DEBUGL(MappingAppll_D, '(' << i << ", " << j << ") is " << data_[i*DimY+j]
             << " at " << &(data_[i*DimX+j]));
      return data_[i*DimX+j];
    }
 
    F& operator()(uint i, uint j) {
      conceptsAssert(i < DimY, Assertion());
      conceptsAssert(j < DimX, Assertion());
      DEBUGL(MappingAppll_D, '(' << i << ", " << j << ") is " << data_[i*DimY+j]
             << " at " << &(data_[i*DimX+j]));
      return data_[i*DimX+j];
    }
 
    template<class H>
    Point<typename Combtype<F,H>::type, DimY> operator*(const Point<H, DimX>& p) const;
 
    template<class H>
    Point<H, DimY> mapPoint(const Point<H, DimX>& b) const {
      return this->operator*(b);
    }
 
    template<class G, uint DimZ>
    Mapping<typename Combtype<G,F>::type, DimY, DimZ> operator*(const Mapping<G, DimX, DimZ>& m) const;
    
    Mapping<F, DimY, DimX>& operator*=(const F n);
    
    Mapping<F, DimY, DimX>& operator*=(const Mapping<F, DimY, DimX>& m);
    
    template<class G>
    Mapping<typename Combtype<G,F>::type, DimY, DimX> operator*(
                const G n) const;
 
    Mapping<F, DimY, DimX>& operator+=(const Mapping<F, DimY, DimX>& M);
    
    Mapping<F, DimY, DimX>& scaleRows(const Point<F, DimY>& s);
    
    Mapping<F, DimY, DimX>& scaleCols(const Point<F, DimX>& s);
    
    Mapping<F, DimX, DimY> transpose() const;
    
    Mapping<F, DimY, DimY> prodTranspose() const;
    
    Point<F, DimY> column(const uint i) const;
 
    void zeros() { F* d = data_;
      for (uint i = 0; i < DimX*DimY; ++i) *d++ = 0; }
    
    template <class H>
    void mapTranspose(const Point<H, DimY>& y, Point<H, DimX>& x) const;
    
    void rank1Product(const Point<F, DimX> x,
                      const Point<F, DimY> y);
    
    template<class H, class J, uint dimy, uint dimx>
    void operator()(const Point<H, dimy>& y, Point<J, dimx>& x) const;
    
    std::ostream& info(std::ostream& os) const;
 
  private: 
    F data_[DimX*DimY];
    
    const uint lin_(uint i, uint j) const { return i*DimX+j; }
  };
 
  // Move to cc file?
  
  template<class F, uint DimY, uint DimX>
  template<class H>
  Point<typename Combtype<F,H>::type, DimY> Mapping<F, DimY, DimX>::operator*(const Point<H, DimX>& p) const
  {
    Point<typename Combtype<F,H>::type, DimY> res;
    this->operator()(p, res);
    return res;
  }
  
  template<class F, uint DimY, uint DimX>
  template<class H>
  inline void Mapping<F, DimY, DimX>::mapTranspose(const Point<H, DimY>& y,
               Point<H, DimX>& x) const 
  {
    H* dr = (H*)x;
    for(uint i = 0; i < DimX; i++, dr++) {
      const H* dp = (const H*)y;
      const F* dm = data_ + i;
      *dr = 0;
      for(uint j = DimY; j--; ++dp) {
        *dr += (*dm) * (*dp);
        dm += DimX;
      }
    }
  }
    
  template<class F, uint DimY, uint DimX>
  template<class H, class J, uint dimy, uint dimx>
  void Mapping<F, DimY, DimX>::operator()(const Point<H, dimy>& y, Point<J, dimx>& x) const
  {
    DEBUGL(MappingAppl_D, *this << " * " << y << " -> " << x);
    if (dimx == dimy && (void*)&x == (void*)&y) {
      
      J*       dx = (J*)x;
      const F* d  = this->data_;
      const H* yp = (const H*)y;
      H yy[dimy];
      for (uint k = 0; k < dimy; k++)  yy[k] = yp[k];
      const H* dy = yy;
      for(uint i = 0; i < (DimY < dimx ? DimY : dimx); ++dx, ++i) {
        dy = yy;
        *dx = 0;
        for(uint j = 0; j < (DimX < dimy ? DimX : dimy); ++j) {
          *dx += d[j] * (*dy++);
        }
        d += DimX;
      }
      if (DimY < dimx) {
        for(uint j = DimY; j < (dimx < dimy ? dimx : dimy); ++j)
          *dx++ = *dy++;
        for(uint k = 0; k < dimx - (DimX < dimy ? dimy : DimX); ++k)
          *dx++ = 0;
      }
      
    } // if (dimx == dimy && &x == &y)
    else {
      DEBUGL(MappingAppl_D, "x != y");
      J*       dx = (J*)x;
      const F* d  = this->data_;
      const H* dy = (const H*)y;
      for(uint i = 0; i < (DimY < dimx ? DimY : dimx); ++dx, ++i) {
        dy = (const H*)y;
        *dx = 0;
        for(uint j = 0; j < (DimX < dimy ? DimX : dimy); ++j) {
          *dx += d[j] * (*dy++);
        }
        DEBUGL(MappingAppl_D, "x = " << x);
        d += DimX;
      }
      if (DimY < dimx) {
        for(uint j = DimY; j < (dimx < dimy ? dimx : dimy); ++j)
          *dx++ = *dy++;
        for(uint k = 0; k < dimx - (DimX < dimy ? dimy : DimX); ++k)
          *dx++ = 0;
      }
      
    } // else (dimx == dimy && &x == &y)
  }
 
  template<class F, uint DimY, uint DimX>
  std::ostream& operator<<(std::ostream& os, const Mapping<F, DimY, DimX>& m) {
    return m.info(os);
  }
 
  // set the real type of a matrix to the real type of the components
  template<typename F, uint DimX, uint DimY>
  struct Realtype<Mapping<F, DimY, DimX> > {
    typedef typename Realtype<F>::type type;
  };
 
  // set the data type of a matrix to the data type of the components
  template<typename F, uint DimX, uint DimY>
  struct Datatype<Mapping<F, DimY, DimX> > {
    typedef typename Datatype<F>::type type;
  };
  
  typedef Mapping<Real, 2> MapReal2d;
  typedef Mapping<Real, 3> MapReal3d;
  typedef Mapping<Cmplx, 2> MapCmplx2d;
  typedef Mapping<Cmplx, 3> MapCmplx3d;  
  
  
  // ******************************************************** GeneralMapping **
 
  template<class F, uint dim>
  struct GeneralMapping {
    typedef typename concepts::Mapping<F,dim> Type;
  };
  
  template<class F>
  struct GeneralMapping<F,1> {
    typedef F Type;
  };
  
 
  // **************************************************************** number **
 
  template<class F, uint dim>
  struct number<Point<F,dim> > {
    static inline std::string name() {
      return number<F>::name() + char('0'+dim) + 'd';
    }
  };
 
  template<class F, uint dim>
  struct number<Mapping<F,dim> > {
    static inline std::string name() {
      return "Map" + number<F>::name() + char('0'+dim) + 'd';
    }
  };
 
  // ********************************************************** GeneralPoint **
 
  template<class F, uint CoordDim>
  struct GeneralPoint : std::conditional<CoordDim==1,
                                         F,
                                         concepts::Point<F,CoordDim>> {};
 
  // ************************************************************ Coordinate **
  template<uint CoordDim>
  struct Coordinate : GeneralPoint<Real, CoordDim> {};
 
  // ******************************************************* CoordinateParam **
  template<int CoordDim>
  struct CoordinateParam : std::conditional<CoordDim==1,
                                            const typename Coordinate<1>::type,
                                            const typename Coordinate<CoordDim>::type &> {};
 
  
} // namespace concepts
 
namespace std {
 
  template<class F, uint dim>
  inline const concepts::Point<F,dim> conj(const concepts::Point<F,dim>& v) { 
    concepts::Point<F,dim> res;
    F* d = (F*)res; const F* e = (const F*)v;
    for(uint i = dim; i--;) *d++ += conj(*e++);
    return res;
  }
 
  template<class F, uint dim>
  inline const concepts::Real norm(const concepts::Point<F,dim>& v) { 
    return v.l2_2();
  }
 
  // ******************************************************************** arg **
 
  concepts::Real arg(const concepts::Point<concepts::Real,2>& p);
 
} // namespace std
 
#endif // vectMatr_hh