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

#ifndef __CGLA_ARITHSQMATFLOAT_H__
#define __CGLA_ARITHSQMATFLOAT_H__

#include "ArithMatFloat.h"

namespace CGLA 
{

  template <class VT, class MT, unsigned int ROWS>
    class ArithSqMatFloat: public ArithMatFloat<VT,VT,MT,ROWS> 
    { 
    public:

      typedef VT VectorType;

      typedef typename VT::ScalarType ScalarType;

    protected:

      ArithSqMatFloat() {}

      explicit ArithSqMatFloat(ScalarType _a):
        ArithMatFloat<VT,VT,MT,ROWS>(_a) {}

      ArithSqMatFloat(VT _a, VT _b): 
        ArithMatFloat<VT,VT,MT,ROWS>(_a,_b) {}

      ArithSqMatFloat(VT _a, VT _b, VT _c): 
        ArithMatFloat<VT,VT,MT,ROWS>(_a,_b,_c) {}

      ArithSqMatFloat(VT _a, VT _b, VT _c, VT _d): 
        ArithMatFloat<VT,VT,MT,ROWS>(_a,_b,_c,_d) {}
                
    public:

      const MT& operator*=(const MT& m2)
        {
          (*this) = (*this) * m2;
          return static_cast<const MT&>(*this);
        }
                
      const MT& operator *=(ScalarType k) 
        {
          return ArithMatFloat<VT,VT,MT,ROWS>::operator*=(k);
        }

      void identity()
        {
          for(unsigned int i=0;i<ROWS;++i)
            {
              for(unsigned int j=0;j<ROWS;++j)
                (*this)[i][j] = ScalarType(0);
              (*this)[i][i] = ScalarType(1);
            }
        }

                
    };

  template <class VT, class MT, unsigned int ROWS>
    inline MT operator*(const ArithSqMatFloat<VT,MT,ROWS>& m1,
                        const ArithSqMatFloat<VT,MT,ROWS>& m2) 
    {
      MT n;
      for(unsigned int i=0;i<ROWS;i++)
        for(unsigned int j=0;j<ROWS;j++)
          {
            n[i][j] = 0;
            for(unsigned int k=0;k<ROWS;k++)
              n[i][j] += m1[i][k] * m2[k][j]; 
          }
      return n;
    }

  template <class VT, class MT, unsigned int ROWS>
    inline MT transpose(const ArithSqMatFloat<VT,MT,ROWS>& m) 
    {
      MT m_new;
      for(unsigned int i=0;i<MT::get_v_dim();i++)
        for(unsigned int j=0;j<MT::get_h_dim();j++)
          m_new[i][j] = m[j][i];
      return m_new;
    }

  template <class VT, class MT, unsigned int ROWS>
    inline typename MT::ScalarType trace(const ArithSqMatFloat<VT,MT,ROWS>& M)
    {
      typename ArithSqMatFloat<VT,MT,ROWS>::ScalarType s=0;
      for(unsigned int i=0;i<ROWS;i++)
        s += M[i][i];
      return s;
    }

}
#endif