---

Matrix.H

Go to the documentation of this file.

#ifndef _MATRIX_H
#define _MATRIX_H

#include <Vector.H>

//#define MxM_VecOfVec Vector<Vector<T, M>, M>
#define MxN_VecOfVec Vector<Vector<T, M>, N>
#define NxM_VecOfVec Vector<Vector<T, N>, M>
#define NxL_VecOfVec Vector<Vector<T, N>, L>
#define MxL_VecOfVec Vector<Vector<T, M>, L>


namespace xchen
{
  template <class T, int M, int N>
  class Matrix
  {
  public:
      Matrix( )                                                                 { }
      Matrix(const MxN_VecOfVec& m) : vOfv(m)                                   { }

      MxN_VecOfVec& operator()()                                                { return vOfv; }       
      const MxN_VecOfVec& operator()() const                                    { return vOfv; }       

      T& operator() (int r, int c)                                              { return vOfv[c][r]; } 
      const T& operator() (int r, int c) const                                  { return vOfv[c][r]; } 
      void SetIdentity();
      void SetFromColumnMajorData(T* data);
      void GetColumnMajorData(T* data) const;


      Vector<T, M> operator()(const Vector<T, N>& x) const                      { return (*this) * x; }  
      Vector<T, M> operator*(const Vector<T, N>& x) const;  

      template <int L>
      Matrix<T,M,L> operator*(const Matrix<T,N,L>& B);                          

      void GetTranspose(Matrix<T,N,M>* t) const;
      void GetInverse(Matrix<T,M,M>* i) const;
      Matrix<T,N,M>* GetTranspose() const                                       { Matrix<T,N,M>* t = new Matrix<T,N,M>; GetTranspose(t); return t; }
      Matrix<T,M,M>* GetInverse() const                                         { Matrix<T,N,M>* i = new Matrix<T,N,M>; GetInverse(i); return t; }

      friend ostream &operator<<<>(ostream &os,const Matrix &m);

  private:
      MxN_VecOfVec vOfv;
  };



  /* matrix implementation */

  template <class T, int M, int N>
  inline void Matrix<T, M, N> :: SetFromColumnMajorData(T* data)
  {
    typename MxN_VecOfVec::iterator iter = vOfv.begin();
    for(; iter != vOfv.end(); iter++)
    {
      copy(data, data+M, iter->begin());
      data += M;
    }
  }

  template <class T, int M, int N>
  inline void Matrix<T, M, N> :: GetColumnMajorData(T* data) const
  {
    int i = 0;
    for(typename MxN_VecOfVec::const_iterator iter1 = vOfv.begin(); iter1 != vOfv.end(); iter1++)
      for(typename Vector<T,M>::const_iterator iter2 = iter1->begin(); iter2 != iter1->end(); iter2++)
      {
        data[i++] = *iter2;
      }
  }

  template <class T, int M, int N>
  inline void Matrix<T, M, N> :: SetIdentity()
  {
    fill(vOfv.begin(), vOfv.end(), Vector<T,M>(0.0));
    assert( M==N );
    for(int i=0; i<M; i++)
      vOfv[i][i] = 1.0;
  }

  template <class T, int M, int N>
  inline Vector<T, M>  Matrix<T, M, N> :: operator*(const Vector<T, N>& x) const
  {
    typename MxN_VecOfVec::const_iterator itr1 = vOfv.begin();
    typename Vector<T, N>::const_iterator itr2 = x.begin();
    Vector<T, M> ret(0.0);
    for(; itr1 != vOfv.end(); itr1++, itr2++)
    {
      ret += (*itr1) * (*itr2);
    }
    return ret;
  }

  template <class T, int M, int N>
  template <int L>
  inline  Matrix<T,M,L> Matrix<T, M, N> :: operator*(const Matrix<T,N,L>& B)
  {
    Vector<T,M> col_vecs[L];
    typename NxL_VecOfVec::const_iterator itr = B().begin();

    for(int i=0; itr != B().end(); itr++, i++)
    {
      col_vecs[i] = (*this) * (*itr);
    }
    return Matrix<T,M,L>(col_vecs);
  }

  template <class T, int M, int N>
  inline void Matrix<T, M, N> :: GetTranspose(Matrix<T,N,M>* t) const
  {
    typename MxN_VecOfVec::const_iterator iter1 = vOfv.begin();
    for(int col=0; iter1 != vOfv.end(); iter1++, col++)
    {
      typename Vector<T,M>::const_iterator iter2 = iter1->begin();
      for(int row=0; iter2 != iter1->end(); iter2++, row++)
      {
        (*t)(col, row) = *iter2;
      }
    }
  }

  template <class T, int M, int N>
  inline void Matrix<T,M,N> :: GetInverse(Matrix<T,M,M>* i) const
  {
    assert( M==N );
    cerr << "Matrix::GetInverse() not implemented yet.\n";
  }


  /* display column majors as columns as they should be */
  template <class T, int M, int N>
  inline ostream &operator<< (ostream &os,const Matrix<T,M,N >&m)
  {
    Matrix<T,N,M>* t = m.GetTranspose();
    copy((*t)().begin(), (*t)().end(), ostream_iterator< Vector<T,N> >(cout, "\n"));
    delete t;
    return os;
  }



  typedef Vector<Vector<double, 2>, 2> dVecOfVec22;
  typedef Vector<Vector<double, 2>, 3> dVecOfVec23;
  typedef Vector<Vector<double, 2>, 4> dVecOfVec24;
  typedef Vector<Vector<double, 2>, 5> dVecOfVec25;

  typedef Vector<Vector<double, 3>, 2> dVecOfVec32;
  typedef Vector<Vector<double, 3>, 3> dVecOfVec33;
  typedef Vector<Vector<double, 3>, 4> dVecOfVec34;
  typedef Vector<Vector<double, 3>, 5> dVecOfVec35;

  typedef Vector<Vector<double, 4>, 2> dVecOfVec42;
  typedef Vector<Vector<double, 4>, 3> dVecOfVec43;
  typedef Vector<Vector<double, 4>, 4> dVecOfVec44;
  typedef Vector<Vector<double, 4>, 5> dVecOfVec45;

  typedef Matrix<double, 2, 2> dMatrix22;
  typedef Matrix<double, 2, 2> dMatrix23;
  typedef Matrix<double, 3, 3> dMatrix24;
  typedef Matrix<double, 3, 3> dMatrix25;

  typedef Matrix<double, 3, 3> dMatrix32;
  typedef Matrix<double, 3, 3> dMatrix33;
  typedef Matrix<double, 3, 4> dMatrix34;
  typedef Matrix<double, 3, 5> dMatrix35;

  typedef Matrix<double, 4, 2> dMatrix42;
  typedef Matrix<double, 4, 3> dMatrix43;
  typedef Matrix<double, 4, 4> dMatrix44;
  typedef Matrix<double, 4, 5> dMatrix45;

  typedef dMatrix44 HMatrix;

}//end namespace xchen


#endif

---