/Users/jab/Documents/Teaching/02585/GEL2_and_demos/GEL/src/CGLA/ArithMatFloat.h

#ifndef __CGLA_ARITHMATFLOAT_H__
#define __CGLA_ARITHMATFLOAT_H__

#include <vector>
#include <iostream>
#include <numeric>

#include "CGLA.h"


namespace CGLA 
{

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    class ArithMatFloat
    { 
    public:

      typedef HVT HVectorType;

      typedef VVT VVectorType;

      typedef typename HVT::ScalarType ScalarType;
                
    protected:

      HVT data[ROWS];

    protected:

      ArithMatFloat() 
        {
#ifndef NDEBUG
          std::fill_n(data, ROWS, HVT(CGLA_INIT_VALUE));
#endif
        }

      explicit ArithMatFloat(ScalarType x)
        {
          std::fill_n(data, ROWS, HVT(x));
        }

      explicit ArithMatFloat(HVT _a)
        {
          std::fill_n(data, ROWS, _a);
        }

      ArithMatFloat(HVT _a, HVT _b)
        {
          assert(ROWS==2);
          data[0] = _a;
          data[1] = _b;
        }

      ArithMatFloat(HVT _a, HVT _b, HVT _c)
        {
          assert(ROWS==3);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
        }

      ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)
        {
          assert(ROWS==4);
          data[0] = _a;
          data[1] = _b;
          data[2] = _c;
          data[3] = _d;
        }
                
    public:

      static unsigned int get_v_dim() {return VVT::get_dim();}

      static unsigned int get_h_dim() {return HVT::get_dim();}


      const ScalarType* get() const 
        {
          return data[0].get();
        }

      ScalarType* get()
        {
          return data[0].get();
        }

      //----------------------------------------------------------------------
      // index operators
      //----------------------------------------------------------------------

      const HVT& operator [] ( unsigned int i ) const
        {
          assert(i<ROWS);
          return data[i];
        }

      HVT& operator [] ( unsigned int i ) 
        {
          assert(i<ROWS);
          return data[i];
        }

      //----------------------------------------------------------------------

      bool operator==(const MT& v) const 
        {
          return std::inner_product(data, &data[ROWS], &v[0], true,
                                    std::logical_and<bool>(), std::equal_to<HVT>());
        }

      bool operator!=(const MT& v) const 
        {
          return !(*this==v);
        }

      //----------------------------------------------------------------------

      const MT operator * (ScalarType k) const
        {
          MT v_new;
          std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::multiplies<HVT>(), k));
          return v_new;
        }

      const MT operator / (ScalarType k) const
        {
          MT v_new;
          std::transform(data, &data[ROWS], &v_new[0], std::bind2nd(std::divides<HVT>(), k));
          return v_new;      
        }

      const MT& operator *=(ScalarType k) 
        {
          std::transform(data, &data[ROWS], data, std::bind2nd(std::multiplies<HVT>(), k));
          return static_cast<const MT&>(*this);
        }

      const MT& operator /=(ScalarType k) 
        { 
          std::transform(data, &data[ROWS], data, std::bind2nd(std::divides<HVT>(), k));
          return static_cast<const MT&>(*this);
        }

      //----------------------------------------------------------------------

      const MT operator + (const MT& m1) const
        {
          MT v_new;
          std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::plus<HVT>());
          return v_new;
        }

      const MT operator - (const MT& m1) const
        {
          MT v_new;
          std::transform(data, &data[ROWS], &m1[0], &v_new[0], std::minus<HVT>());
          return v_new;
        }

      const MT& operator +=(const MT& v) 
        {
          std::transform(data, &data[ROWS], &v[0], data, std::plus<HVT>());
          return static_cast<const MT&>(*this);
        }

      const MT& operator -=(const MT& v) 
        {
          std::transform(data, &data[ROWS], &v[0], data, std::minus<HVT>());
          return static_cast<const MT&>(*this);
        }

      //----------------------------------------------------------------------

      const MT operator - () const
        {
          MT v_new;
          std::transform(data, &data[ROWS], &v_new[0], std::negate<HVT>());
          return v_new;
        }
    };

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v) 
    {
      VVT v2;
      for(unsigned int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);
      return v2;
    }


#ifndef WIN32

  template <class VVT, class HVT, 
    class HV1T, class VV2T,
    class MT1, class MT2, class MT,
    unsigned int ROWS1, unsigned int ROWS2>
    inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,
                    const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,
                    ArithMatFloat<VVT,HVT,MT,ROWS1>& m)
    {
      unsigned int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
      for(unsigned int i=0;i<ROWS1;i++)
        for(unsigned int j=0;j<cols;j++)
          {
            m[i][j] = 0;
            for(unsigned int k=0;k<ROWS2;k++)
              m[i][j] += m1[i][k] * m2[k][j]; 
          }
    }


  template <class VVT, class HVT, class M1T, class M2T, unsigned int ROWS, unsigned int COLS>
    inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,
                          ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)
    {
      for(unsigned int i=0;i<M2T::get_v_dim();++i)
        for(unsigned int j=0;j<M2T::get_h_dim();++j)
          m_new[i][j] = m[j][i];
    }

#else

  //----------------- win32 -------------------------------
  // Visual studio is not good at deducing the args. to these template functions.
  // This means that you can call the two functions below with 
  // matrices of wrong dimension.

  template <class M1, class M2, class M>
    inline void mul(const M1& m1, const M2& m2, M& m)
    {
      unsigned int cols = M::get_h_dim();
      unsigned int rows1 = M1::get_v_dim();
      unsigned int rows2 = M2::get_v_dim();

      for(unsigned int i=0;i<rows1;++i)
        for(unsigned int j=0;j<cols;++j)
          {
            m[i][j] = 0;
            for(unsigned int k=0;k<rows2;++k)
              m[i][j] += m1[i][k] * m2[k][j];
          }
    }


  template <class M1, class M2>
    inline void transpose(const M1& m1, M2& m2)
    {
      for(unsigned int i=0;i<M2::get_v_dim();++i)
        for(unsigned int j=0;j<M2::get_h_dim();++j)
          m2[i][j] = m1[j][i];
    }

#endif

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    void outer_product(const VVT& a, const HVT& b, 
                       ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      unsigned int R = VVT::get_dim();
      unsigned int C = HVT::get_dim();
      for(unsigned int i=0;i<R;++i)
        for(unsigned int j=0;j<C;++j)
          {
            m[i][j] = a[i] * b[j];
          }
    }

  template <class VVT, class HVT, class MT, int ROWS, class BinOp>
    void outer_product(const VVT& a, const HVT& b, 
                       ArithMatFloat<VVT,HVT,MT,ROWS>& m, BinOp op)
    {
      int R = VVT::get_dim();
      int C = HVT::get_dim();
      for(int i=0;i<R;++i)
        for(int j=0;j<C;++j)
          {
            m[i][j] = op(a[i], b[j]);
          }
    }

  template <class M1, class M2>
    void copy_matrix(const M1& inmat, M2& outmat)
    {
      const unsigned int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
      const unsigned int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
      for(unsigned int i=0;i<R;++i)
        for(unsigned int j=0;j<C;++j)
          outmat[i][j] = inmat[i][j];
    }

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline std::ostream& 
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      os << "[\n";
      for(unsigned int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
      os << "]\n";
      return os;
    }

  template <class VVT, class HVT, class MT, unsigned int ROWS>
    inline std::istream& operator>>(std::istream&is, 
                                    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      for(unsigned int i=0;i<ROWS;i++) is>>m[i];
      return is;
    }
}
#endif