basics/vectorsMatrices.hh

Go to the documentation of this file.

/* vectors and matrices for 2D and 3D
*/

#ifndef vectMatr_hh
#define vectMatr_hh

#include <iostream>
#include <cmath>
#include <cstring>

#include "functionCompiler.hh"
#include "exceptions.hh"
#include "typedefs.hh"
#include "vectorsMatricesForward.hh"
#include "debug.hh"

// debugging Mapping
#define MappingAppll_D 0
#define MappingAppl_D 0

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, int 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<int 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) {
          data_[0] = e0;
          if (dim > 1)
            data_[1] = e1;
          if (dim > 2)
            data_[2] = e2;
        }
        Point(const uchar *pgm, const Real x, const Real y, const Real z = 0.0);

        template<int 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 int i) const {
          conceptsAssert(i < dim, concepts::Assertion());
          return data_[i]; }
          F& operator[](const int 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 {
              return *this * b;
            }

            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;
              }
            template<class G>
            Point<typename Combtype<F,G>::type, dim> operator*(const G n) const {
              Point<typename Combtype<F,G>::type, 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;
            }

            F cross(const Point<F, dim>& b) const {
              return data_[0]*b[1] - data_[1]*b[0];
            }
            Point<F, dim> operator^(const Point<F, dim>& 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, int 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, int 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, int dim>
  struct Datatype<Point<F,dim> > {
    //typedef typename Datatype<F>::type type;
    typedef F type;
  };

  template <class F, int 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, int dim>
  inline Point<typename Combtype<F,Real>::type,dim> operator*(const Real x, const Point<F,dim>& y)
  {
    return y * x;
  }

  template <class F, int dim>
  inline Point<typename Combtype<F,Cmplx>::type,dim> operator*(const Cmplx x, const Point<F,dim>& y)
  {
    return y * x;
  }

  template<int 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<int 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++) * conj(*db++);
    return res;
  }

  // ********************************************************** GeneralPoint **

  template<class F, uint dim>
  struct GeneralPoint {
    typedef typename concepts::Point<F,dim> Type;
  };
  
  template<class F>
  struct GeneralPoint<F,1> {
    typedef F Type;
  };


  // **************************************************************** UnitNd **

  template <int 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 <int dim>
    inline UnitNd<dim>::UnitNd(const Point<Real, dim>& u)
    : Point<Real, dim>(u) {
      len_ = this->l2();  *this *= (1.0 / len_);
    }

  template <int 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 **

  template<class F, int dim>
  class Mapping {
  public:
    Mapping() {}
    Mapping(const Mapping<F, dim>& m) {
      memorycpy(data_, m.data_, dim*dim);
    }

    template <class H>
    Mapping(const Mapping<H, dim>& m) {
      for(int i=0; i < dim*dim; ++i)
        data_[i] = m.data_[i];
    }

    // FIXME: this is not really nice
    friend class Mapping<Cmplx, dim>;
    friend class Mapping<Real, dim>;

    Mapping(const F e0, const F e1,
            const F e2, const F e3);
    Mapping(const Point<F, dim> first, const Point<F, dim> second);
    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);
    Mapping(const Point<F, dim> first,
            const Point<F, dim> second,
            const Point<F, dim> third);
    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 UnitNd<dim>& n);

    const F& operator()(uint i, uint j) const {
      conceptsAssert(i < dim, Assertion());
      conceptsAssert(j < dim, Assertion());
      DEBUGL(MappingAppll_D, '(' << i << ", " << j << ") is " << data_[i*dim+j]
             << " at " << &(data_[i*dim+j]));
      return data_[i*dim+j];
    }

    F& operator()(uint i, uint j) {
      conceptsAssert(i < dim, Assertion());
      conceptsAssert(j < dim, Assertion());
      DEBUGL(MappingAppll_D, '(' << i << ", " << j << ") is " << data_[i*dim+j]
             << " at " << &(data_[i*dim+j]));
      return data_[i*dim+j];
    }


    Mapping(F x) { F* d = data_; for(uint i = dim*dim; i--;) *d++ = x; }

    template<class H>
    Point<typename Combtype<F,H>::type, dim> operator*(const Point<H, dim>& p) const;

    template<class H>
    Point<H, dim> mapPoint(const Point<H, dim>& b) const {
      return this->operator*(b);
    }

    template<class G>
    Mapping<typename Combtype<F,G>::type, dim> 
    operator*(const Mapping<G, dim>& m) const;
    Mapping<F, dim>& operator*=(const F n);
    Mapping<F, dim>& operator*=(const Mapping<F, dim>& m);
    template<class G>
    Mapping<typename Combtype<F,G>::type, dim> operator*(const G n) const;

    Mapping<F, dim>& operator+=(const Mapping<F, dim>& M);
    Mapping<F, dim>& scaleRows(const Point<F, dim>& s);
    Mapping<F, dim>& scaleCols(const Point<F, dim>& s);
    F determinant() const;
    Mapping<F, dim> inverse() const;
    Mapping<F, dim> adjugate() const;
    Mapping<F, dim> transpose() const;
    Mapping<F, dim> prodTranspose() const;
    Point<F, dim> column(const uint i) const;

      template<class H, class J, int dimy, int dimx>
      void operator()(const Point<H, dimy>& y, Point<J, dimx>& x) const;

      void zeros() { F* d = data_;
        for (uint i = 0; i < dim*dim; ++i) *d++ = 0; }
      
      template <class H>
      void mapTranspose(const Point<H, dim>& y, Point<H, dim>& x) const;
      
      void rank1Product(const Point<F, dim> x,
                        const Point<F, dim> y);
      
      std::ostream& info(std::ostream& os) const;
    private:
      F data_[dim*dim];
    };

  template<class F, int dim>
  template<class H>
  Point<typename Combtype<F,H>::type, dim> Mapping<F, dim>::operator*(const Point<H, dim>& p) const
  {
    Point<typename Combtype<F,H>::type, dim> res;
    this->operator()(p, res);
    return res;
  }

  template<class F, int dim>
  template<class H, class J, int dimy, int dimx>
  void Mapping<F, dim>::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  = data_;
      const H* yp = (const H*)y;
      H yy[dimy];
      for (int k = 0; k < dimy; k++)  yy[k] = yp[k];
      const H* dy = yy;
      for(int i = 0; i < (dim < dimx ? dim : dimx); ++dx, ++i) {
        dy = yy;
        *dx = 0;
        for(int j = 0; j < (dim < dimy ? dim : dimy); ++j) {
          *dx += d[j] * (*dy++);
        }
        d += dim;
      }
      if (dim < dimx) {
        for(int j = dim; j < (dimx < dimy ? dimx : dimy); ++j)
          *dx++ = *dy++;
        for(int k = 0; k < dimx - (dim < dimy ? dimy : dim); ++k)
          *dx++ = 0;
      }
      
    } // if (dimx == dimy && &x == &y)
    else {
      DEBUGL(MappingAppl_D, "x != y");
      J*       dx = (J*)x;
      const F* d  = data_;
      const H* dy = (const H*)y;
      for(int i = 0; i < (dim < dimx ? dim : dimx); ++dx, ++i) {
        dy = (const H*)y;
        *dx = 0;
        for(int j = 0; j < (dim < dimy ? dim : dimy); ++j) {
          *dx += d[j] * (*dy++);
        }
        DEBUGL(MappingAppl_D, "x = " << x);
        d += dim;
      }
      if (dim < dimx) {
        for(int j = dim; j < (dimx < dimy ? dimx : dimy); ++j)
          *dx++ = *dy++;
        for(int k = 0; k < dimx - (dim < dimy ? dimy : dim); ++k)
          *dx++ = 0;
      }
      
    } // else (dimx == dimy && &x == &y)
  }
  
  template<class F, int dim>
    template<class H>
    inline void Mapping<F, dim>::mapTranspose(const Point<H, dim>& y,
        Point<H, dim>& x) const 
    {
      H* dr = (H*)x;
      for(uint i = 0; i < dim; i++, dr++) {
        const H* dp = (const H*)y;
        const F* dm = data_ + i;
        *dr = 0;
        for(uint j = dim; j--; ++dp) {
          *dr += (*dm) * (*dp);
          dm += dim;
        }
      }
    }

  template<class F, int dim>
  std::ostream& operator<<(std::ostream& os, const Mapping<F, dim>& m) {
    return m.info(os);
  }

  // set the real type of a matrix to the real type of the components
  template<typename F, int dim>
  struct Realtype<Mapping<F,dim> > {
    typedef typename Realtype<F>::type type;
  };

  // set the data type of a matrix to the data type of the components
  template<typename F, int dim>
  struct Datatype<Mapping<F,dim> > {
    typedef typename Datatype<F>::type type;
  };

  // ******************************************************** 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<>
  template<class F, int dim>
  struct number<Point<F,dim> > {
    static inline std::string name() {
      return number<F>::name() + char('0'+dim) + 'd';
    }
  };

  template<>
  template<class F, int dim>
  struct number<Mapping<F,dim> > {
    static inline std::string name() {
      return "Map" + number<F>::name() + char('0'+dim) + 'd';
    }
  };

} // namespace concepts

namespace std {

  template<class F, int 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, int 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