#ifndef FILE_NGBLA_EXPR
#define FILE_NGBLA_EXPR

/**************************************************************************/
/* File:   expr.hpp                                                       */
/* Author: Joachim Schoeberl                                              */
/* Date:   01. Jan. 02                                                    */
/**************************************************************************/



template <int H, int W, typename T> class Mat;
template <int S, typename T> class Vec;


/*
  Matrix expression templates
 */


#ifdef GNU_PRE3

template <class T, int num>
inline T & Access (Vec<num, T> & vec, int nr)
{
  return vec(nr);
}
template <class T, int num>
inline T Access (const Vec<num, T> & vec, int nr) 
{
  return vec(nr);
}

#else

template <typename TVEC> 
inline typename TVEC::TELEM & Access (TVEC & vec, int nr)
{
  return vec(nr);
}

template <typename TVEC> 
inline typename TVEC::TELEM Access (const TVEC & vec, int nr)
{
  return vec(nr);
}

#endif // GNU_PRE3



inline double & Access (double & vec, int nr)
{
  return vec;
}

inline double Access (const double & vec, int nr)
{
  return vec;
}

inline Complex & Access (Complex & vec, int nr)
{
  return vec;
}

inline Complex Access (const Complex & vec, int nr)
{
  return vec;
}








/**
   Trait to access small matrix dimensions
 */
template <class T>
class mat_traits
{
public:
  /// matrix element
  typedef typename T::TELEM TELEM;
  /// field of matrix element
  typedef typename T::TSCAL TSCAL;
  /// type of column vector
  typedef typename T::TV_COL TV_COL;
  /// type of row vector
  typedef typename T::TV_ROW TV_ROW;
  /// matrix height
  enum { HEIGHT = T::HEIGHT };
  /// matrix with
  enum { WIDTH  = T::WIDTH  };
};


// some compiler thinks it needs Mat<int>...
template <>
class mat_traits<Mat<1,1,int> >
{
public:
   typedef int TELEM;
   typedef int TSCAL;
   typedef int TV_COL;
   typedef int TV_ROW;
   enum { HEIGHT = 1 };
   enum { WIDTH = 1 };
};

template <>
class mat_traits<int>
{
public:
  typedef int TELEM;
  typedef int TSCAL;
  typedef int TV_COL;
  typedef int TV_ROW;
  enum { HEIGHT = 1 };
  enum { WIDTH = 1 };
};


template <>
class mat_traits<double>
{
public:
  typedef double TELEM;
  typedef double TSCAL;
  typedef double TV_COL;
  typedef double TV_ROW;
  enum { HEIGHT = 1 };
  enum { WIDTH = 1 };
};

template <>
class mat_traits<Complex>
{
public:
  typedef Complex TELEM;
  typedef Complex TSCAL;
  typedef Complex TV_COL;
  typedef Complex TV_ROW;
  enum { HEIGHT = 1 };
  enum { WIDTH = 1 };
};

template <>
class mat_traits<const double>
{
public:
  typedef double TELEM;
  typedef double TSCAL;
  typedef double TV_COL;
  typedef double TV_ROW;
  enum { HEIGHT = 1 };
  enum { WIDTH = 1 };
};

template <>
class mat_traits<const Complex>
{
public:
  typedef Complex TELEM;
  typedef Complex TSCAL;
  typedef Complex TV_COL;
  typedef Complex TV_ROW;
  enum { HEIGHT = 1 };
  enum { WIDTH = 1 };
};


/// matrix type from column and row vectors
template <typename TV_COL, typename TV_ROW>
class mat_from_vecs
{
  enum { HEIGHT = mat_traits<TV_COL>::HEIGHT };
  enum { WIDTH  = mat_traits<TV_ROW>::HEIGHT };
  typedef typename mat_from_vecs<typename TV_COL::TELEM, typename TV_ROW::TELEM>::TMAT TELEM;
  typedef Mat<HEIGHT,WIDTH,TELEM> TMAT;
};

template <> class mat_from_vecs<double,double> { typedef double TMAT; };
template <> class mat_from_vecs<double,Complex> { typedef Complex TMAT; };
template <> class mat_from_vecs<Complex,double> { typedef Complex TMAT; };
template <> class mat_from_vecs<Complex,Complex> { typedef Complex TMAT; };



/// matrix type of product
template <typename TA, typename TB>
class mat_prod_type
{
public:
  enum { HEIGHT = mat_traits<TA>::HEIGHT };
  enum { WIDTH  = mat_traits<TB>::WIDTH };
  typedef typename mat_prod_type<typename TA::TELEM, 
				 typename TB::TELEM>::TMAT TELEM;
  typedef Mat<HEIGHT,WIDTH,TELEM> TMAT;
};

template <> class mat_prod_type<double,double> { public: typedef double TMAT; };
template <> class mat_prod_type<double,Complex> { public: typedef Complex TMAT; };
template <> class mat_prod_type<Complex,double> { public: typedef Complex TMAT; };
template <> class mat_prod_type<Complex,Complex> { public: typedef Complex TMAT; };




/// matrix type of sum (important for double+Complex)
template <typename TA, typename TB>
class mat_sum_type
{
  enum { HEIGHT = mat_traits<TA>::HEIGHT };
  enum { WIDTH  = mat_traits<TA>::WIDTH };
  typedef typename mat_sum_type<typename TA::TELEM, 
				typename TB::TELEM>::TMAT TELEM;
  typedef Mat<HEIGHT,WIDTH,TELEM> TMAT;
};

template <> class mat_sum_type<double,double> { public: typedef double TMAT; };
template <> class mat_sum_type<double,Complex> { public: typedef Complex TMAT; };
template <> class mat_sum_type<Complex,double> { public: typedef Complex TMAT; };
template <> class mat_sum_type<Complex,Complex> { public: typedef Complex TMAT; };








/// matrix type of scale
template <typename TM, typename TS>
class mat_scale_type
{
public:
  enum { HEIGHT = mat_traits<TM>::HEIGHT };
  enum { WIDTH  = mat_traits<TM>::WIDTH };
  typedef typename mat_scale_type<typename TM::TELEM, TS>::TMAT TELEM;
  typedef Mat<HEIGHT,WIDTH,TELEM> TMAT;
};

template <> class mat_scale_type<double,double> { public: typedef double TMAT; };
template <> class mat_scale_type<double,Complex> { public: typedef Complex TMAT; };
template <> class mat_scale_type<Complex,double> { public: typedef Complex TMAT; };
template <> class mat_scale_type<Complex,Complex> { public: typedef Complex TMAT; };







/// Height of matrix
template <class TM> 
inline int Height (const TM & m)
{
  return m.Height();
}

/// Width of matrix
template <class TM> 
inline int Width (const TM & m)
{
  return m.Width();
}

template <> inline int Height<double> (const double&) { return 1; }
template <> inline int Height<Complex> (const Complex&) { return 1; }
template <> inline int Width<double> (const double&) { return 1; }
template <> inline int Width<Complex> (const Complex&) { return 1; }


/// Complex to double assignment called
class Complex2RealException : public Exception
{
public:
  Complex2RealException()
    : Exception("Assignment of Complex 2 Real") { ; }
};

/// converts Complex to type T. Throw exception if illegal
template<class T>
inline T ReduceComplex (Complex x)
{
  return x;
}

template<> 
inline double ReduceComplex<double> (Complex x)
{
  throw Complex2RealException();
}


/// converts type T to double. Throw exception if illegal
template<class T>
inline double ReduceToReal (T x)
{
  return x;
}

template<>
inline double ReduceToReal (Complex x)
{
  throw Complex2RealException();
}





/// matrices do not fit for matrix-matrix operation
class MatrixNotFittingException : public Exception
{
public:
  MatrixNotFittingException (const string & where, int h1, int w1,
			     int h2, int w2)
    : Exception ("") 
  {
    stringstream str;
    str << where << " mat1 = " << h1 << " x " << w1 << ", mat2 = " << h2 << " x " << w2 << endl;
    Append (str.str());
  }
};



/**
  Expr is the base class for all matrix template expressions.
  Barton and Nackman Trick for template polymorphism, function Spec.

  provides Height and Width of matrix.
  IsLinear allows linear matrix element access.
 */

template <typename T>
class Expr 
{
public:
  Expr () { ; }

  typedef T TConv;  
  T & Spec() { return static_cast<T&> (*this); }
  const T & Spec() const { return static_cast<const T&> (*this); }

  int Height() const { return Spec().T::Height(); }
  int Width() const { return Spec().T::Width(); }

  void Dump (ostream & ost) const { Spec().Dump(ost); }
};



template <typename T>
class SymExpr : public Expr<SymExpr<T> >
{
  const T a;
public:
  typedef typename T::TELEM TELEM;
  typedef typename mat_traits<TELEM>::TSCAL TSCAL;

  SymExpr (const T & aa) : a(aa) { ; }

  TELEM operator() (int i) const { return a(i); }
  TELEM operator() (int i, int j) const { return a(i,j); }
  int Height() const { return a.Height(); }
  int Width() const { return a.Width(); }
  enum { IS_LINEAR = T::IS_LINEAR };
  void Dump (ostream & ost) const
  { ost << "Sym ("; a.Dump(ost); ost << ")"; }
};


template <typename T>
inline SymExpr<T> Symmetric (const Expr<T> & a)
{
  return SymExpr<T> (static_cast<const T&> (a));
}



/**
   Reference to matrix expression.
   T is matrix class.
 */
template <typename T>
class RefMatExpr : public Expr<RefMatExpr<T> >
{
protected:
  const T & data;
public:
  typedef typename T::TELEM TELEM;
  typedef typename T::TSCAL TSCAL;

  RefMatExpr (const T & d) : data(d) { ; }

  TELEM operator() (int i) const  { return data(i); }
  TELEM operator() (int i, int j) const  { return data(i, j); }

  int Height() const { return data.T::Height(); }
  int Width() const { return data.T::Width(); }

  enum { IS_LINEAR = T::IS_LINEAR };
  void Dump (ostream & ost) const { ost << "Matrix-Ref(" << data << ")"; }
};


/**
   The base class for matrices.
 */
template <class T>
class MatExpr : public Expr<T>
{
public:
  typedef RefMatExpr<T> TConv;

  int Height() const { return Spec().T::Height(); }
  int Width() const { return Spec().T::Width(); }

  MatExpr () : Expr<T>() { ; }

  T & Spec() { return static_cast<T&> (*this); }
  const T & Spec() const { return static_cast<const T&> (*this); }

  enum { IS_LINEAR = 1 };


  template<typename TB>
  T & operator= (const Expr<TB> & v)
  {
#ifdef CHECK_RANGE
    if (Height() != v.Height() || Width() != v.Width())
      {
	throw MatrixNotFittingException ("operator=", 
					 Height(), Width(),
					 v.Height(), v.Width());
      }
#endif
    /*
    cout << "operator=, type expr = " << typeid(v).name() << endl
	 << ", size = " << sizeof(Expr<TB>) << endl;
    */


    if (TB::IS_LINEAR)
      {
	int hw = Expr<T>::Height() * Expr<T>::Width();
	for (int i = 0; i < hw; i++)
	  Spec()(i) = v.Spec()(i);
      }
    else
      {
	int k = 0;
	int h = Expr<T>::Height();
	int w = Expr<T>::Width();
	// TB::TELEM * pelem = & Spec()(0);

	for (int i = 0; i < h; i++)
	  for (int j = 0; j < w; j++)
	    {
	      Spec()(k) = v.Spec()(i,j);
	      // pelem[k] = v.Spec()(i,j);
	      k++;
	    }
      }
    return Spec();
  }


  template<typename TB>
  T & operator+= (const Expr<TB> & v)
  {
#ifdef CHECK_RANGE
    if (Height() != v.Height() || Width() != v.Width())
      throw MatrixNotFittingException ("operator+=", 
				       Height(), Width(),
				       v.Height(), v.Width());
#endif

    if (TB::IS_LINEAR)
      {
	int hw = Expr<T>::Height() * Expr<T>::Width();
	for (int i = 0; i < hw; i++)
	  Spec()(i) += v.Spec()(i);
      }
    else
      {
	int k = 0;
	int h = Height();
	int w = Width();
	for (int i = 0; i < h; i++)
	  for (int j = 0; j < w; j++)
	    {
	      Spec()(k) += v.Spec()(i,j);
	      k++;
	    }
      }
    return Spec();
  }



  template<typename TB>
  T & operator+= (const Expr<SymExpr<TB> > & v)
  {
#ifdef CHECK_RANGE
    if (Height() != v.Height() || Width() != v.Width())
      throw MatrixNotFittingException ("operator+=", 
				       Height(), Width(),
				       v.Height(), v.Width());
#endif
    int h = Height();
    for (int i = 0; i < h; i++)
      {
	for (int j = 0; j < i; j++)
	  {
	    double val = v.Spec()(i,j);
	    Spec()(i,j) += val;
	    Spec()(j,i) += val;
	  }
	Spec()(i,i) += v.Spec()(i,i);
      }
    return Spec();
  }



  T & operator+= (double scal)
    {
      int hw = Height() * Width();
      for (int i = 0; i < hw; i++)
	Spec()(i) += scal;
    }
  

  template<typename TB>
  T & operator-= (const Expr<TB> & v)
  {

#ifdef CHECK_RANGE
    if (Height() != v.Height() || Width() != v.Width())
      throw MatrixNotFittingException ("operator-=", 
				       Height(), Width(),
				       v.Height(), v.Width());
#endif

    if (TB::IS_LINEAR)
      {
	int hw = Height() * Width();
	for (int i = 0; i < hw; i++)
	  Spec()(i) -= v.Spec()(i);
      }
    else
      {
	int k = 0;
	int h = Height();
	int w = Width();

	for (int i = 0; i < h; i++)
	  for (int j = 0; j < w; j++)
	    {
	      Spec()(k) -= v.Spec()(i,j);
	      k++;
	    }
      }
    return Spec();
  }


  template <class SCAL2>
  T & operator*= (const SCAL2 & s)
  {
    int hw = Height() * Width();
    for (int i = 0; i < hw; i++)
      Spec()(i) *= s;
    return Spec();
  }

  template <class SCAL2>
  T & operator/= (const SCAL2 & s)
  {
    return (*this) *= (1/s);
  }

};





/* *************************** SumExpr **************************** */

/**
   Sum of 2 matrix expressions
 */


template <class TA, class TB> 
class SumExpr : public Expr<SumExpr<TA,TB> >
{
  const TA a;
  const TB b;
public:
  typedef typename mat_sum_type<typename TA::TELEM,
				typename TB::TELEM>::TMAT TELEM;
  typedef typename mat_traits<TELEM>::TSCAL TSCAL;

  SumExpr (const TA & aa, const TB & ab) : a(aa), b(ab) 
  {
    ;
    //    cout << "SUM, TA = " << typeid(TA).name() << endl
    //	 << ", TB = " << typeid(TB).name() << endl;
  }

  TELEM operator() (int i) const { return a(i)+b(i); }
  TELEM operator() (int i, int j) const { return a(i,j)+b(i,j); }
  int Height() const { return a.Height(); }
  int Width() const { return a.Width(); }

  enum { IS_LINEAR = TA::IS_LINEAR && TB::IS_LINEAR };
  void Dump (ostream & ost) const
  { ost << "("; a.Dump(ost); ost << ") + ("; b.Dump(ost); ost << ")"; }
};

template <typename TA, typename TB>
inline SumExpr<typename TA::TConv, typename TB::TConv>
operator+ (const Expr<TA> & a, const Expr<TB> & b)
{
  return SumExpr<typename TA::TConv, typename TB::TConv> (static_cast <const TA&> (a), static_cast <const TB&> (b));
}



/* *************************** SubExpr **************************** */


/**
   Matrix-expr minus Matrix-expr
 */

template <class TA, class TB> 
class SubExpr : public Expr<SubExpr<TA,TB> >
{
  const TA a;
  const TB b;
public:
  typedef typename TA::TELEM TELEM;
  typedef typename TA::TSCAL TSCAL;

  SubExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }

  TELEM operator() (int i) const { return a(i)-b(i); }
  TELEM operator() (int i, int j) const { return a(i,j)-b(i,j); }
  int Height() const { return a.Height(); }
  int Width() const { return a.Width(); }
  // static bool IsLinear() { return TA::IsLinear() && TB::IsLinear(); }
  enum { IS_LINEAR = TA::IS_LINEAR && TB::IS_LINEAR };
};


template <typename TA, typename TB>
inline SubExpr<typename TA::TConv, typename TB::TConv>
operator- (const Expr<TA> & a, const Expr<TB> & b)
{
  return SubExpr<typename TA::TConv, typename TB::TConv> 
    (static_cast <const TA&> (a), static_cast <const TB&> (b));
}







/* *************************** MinusExpr **************************** */


/**
   minus Matrix-expr
 */

template <class TA>
class MinusExpr : public Expr<MinusExpr<TA> >
{
  const TA a;
public:
  typedef typename TA::TELEM TELEM;
  typedef typename TA::TSCAL TSCAL;

  MinusExpr (const TA & aa) : a(aa) { ; }

  TELEM operator() (int i) const { return -a(i); }
  TELEM operator() (int i, int j) const { return -a(i,j); }
  int Height() const { return a.Height(); }
  int Width() const { return a.Width(); }

  enum { IS_LINEAR = TA::IS_LINEAR };
};



template <typename TA>
inline MinusExpr<typename TA::TConv>
operator- (const Expr<TA> & a)
{
  return MinusExpr<typename TA::TConv> (static_cast<const TA&> (a) );
}



/* *************************** ScaleExpr **************************** */


/**
   Scalar times Matrix-expr
 */
template <class TA, class TS> 
class ScaleExpr : public Expr<ScaleExpr<TA,TS> >
{
  TA a;
  TS s;
public:
  typedef typename mat_scale_type<typename TA::TELEM, TS>::TMAT TELEM;
  typedef typename mat_traits<TELEM>::TSCAL TSCAL;
  enum { IS_LINEAR = TA::IS_LINEAR };


  ScaleExpr (const TA & aa, TS as) : a(aa), s(as) { ; }

  TELEM operator() (int i) const { return s * a(i); }
  TELEM operator() (int i, int j) const { return s * a(i,j); }

  int Height() const { return a.Height(); }
  int Width() const { return a.Width(); }
  //  static bool IsLinear() { return TA::IsLinear(); }
  void Dump (ostream & ost) const
  { ost << "Scale, s=" << s << " * "; a.Dump(ost);  }
};

template <typename TA>
inline ScaleExpr<typename TA::TConv, double> 
operator* (double b, const Expr<TA> & a)
{
  return ScaleExpr<typename TA::TConv, double> (static_cast<const TA&> (a), b );
}

template <typename TA>
inline ScaleExpr<TA, Complex> 
operator* (Complex b, const Expr<TA> & a)
{
  return ScaleExpr<TA, Complex> (static_cast<const TA&> (a), b);
}


/* ************************* MultExpr ************************* */


/**
   Matrix-expr timex Matrix-expr
 */
template <class TA, class TB> class MultExpr : public Expr<MultExpr<TA,TB> >
{
  const TA a;
  const TB b;
public:
  typedef typename mat_prod_type<typename TA::TELEM, 
				 typename TB::TELEM>::TMAT TELEM;
  typedef typename mat_traits<TELEM>::TSCAL TSCAL;

  MultExpr (const TA & aa, const TB & ab) : a(aa), b(ab) { ; }

  TELEM operator() (int i) const { return TELEM(TSCAL(0)); }
  TELEM operator() (int i, int j) const
  { 
    int wa = a.Width();
    TELEM sum;
    if (wa >= 1)
      {
	sum = a(i,0) * b(0,j);
	for (int k = 1; k < wa; k++)
	  sum += a(i,k) * b(k, j);
      }
    else
      sum = TSCAL(0);
    
    return sum;
  }

  int Height() const { return a.Height(); }
  int Width() const { return b.Width(); }
  enum { IS_LINEAR = 0 };
};


template <typename TA, typename TB>
inline MultExpr<typename TA::TConv, typename TB::TConv>
operator* (const Expr<TA> & a, const Expr<TB> & b)
{
  return MultExpr<typename TA::TConv, typename TB::TConv> 
    (static_cast <const TA&> (a), static_cast <const TB&> (b));
}


/* ************************** Trans *************************** */


/**
   Transpose of Matrix-expr
 */
template <class TA> class TransExpr : public Expr<TransExpr<TA> >
{
  const TA a;
public:
  typedef typename TA::TELEM TELEM;

  TransExpr (const TA & aa) : a(aa) { ; }

  int Height() const { return a.Width(); }
  int Width() const { return a.Height(); }

  TELEM operator() (int i, int j) const { return Trans (a(j,i)); }
  TELEM operator() (int i) const { return 0; }

  enum { IS_LINEAR = 0 };
};


/// Transpose 
template <typename TA>
inline TransExpr<typename TA::TConv>
Trans (const Expr<TA> & a)
{
  return TransExpr<typename TA::TConv> (static_cast <const TA&> (a));
}

inline double Trans (double a) { return a; }
inline Complex Trans (Complex a) { return a; }





/* ************************* InnerProduct ********************** */


/**
   Inner product
 */
inline double InnerProduct ( const double& a, const double& b )
{
  return a * b;
}

inline Complex InnerProduct ( const Complex& a, const Complex b)
{
  return a * b;
}

template <class TA, class TB>
inline typename TA::TSCAL
InnerProduct (const MatExpr<TA> & a, const MatExpr<TB> & b)
{
  typedef typename TA::TSCAL TSCAL;
  TSCAL sum = 0;
  for (int i = 0; i < a.Height(); i++)
    sum += InnerProduct (a.Spec()(i), b.Spec()(i));
  return sum;
}


/* **************************** Trace **************************** */

template <class TA>
inline typename TA::TELEM
Trace (const Expr<TA> & a)
{
  typedef typename TA::TELEM TELEM;
  TELEM sum = 0;
  for (int i = 0; i < Height(a); i++)
    sum += a.Spec()(i,i);
  return sum;
}

/* **************************** L2Norm **************************** */

/// Euklidean norm squared

inline double L2Norm2 (double v)
{
  return v*v;
}

inline double L2Norm2 (Complex v)
{
  return v.real()*v.real()+v.imag()*v.imag();
}

template <class TA>
inline double L2Norm2 (const MatExpr<TA> & v)
{
  typedef typename TA::TSCAL TSCAL;
  double sum = 0;
  for (int i = 0; i < v.Height(); i++)
    sum += L2Norm2 (v.Spec()(i));  // REval
  return sum;
}

template <class TA>
inline double L2Norm (const MatExpr<TA> & v)
{
  return sqrt (L2Norm2(v));
}

/* *************************** Output ****************************** */



/// Print matrix-expr
template<typename T>
std::ostream & operator<< (std::ostream & s, const Expr<T> & v)
{
  for (int i = 0; i < v.Height(); i++)
    {
      for (int j = 0 ; j < v.Width(); j++)
	s << " " << setw(7) << v.Spec()(i,j);
      s << endl;
    }
  return s;
}



/* **************************** Inverse *************************** */

template <typename T> class FlatMatrix;
template <typename T> class Matrix;
template <int H, int W, typename T> class Mat;

/// Calculate inverse. Gauss elimination with row pivoting
template <class T, class T2>
extern void CalcInverse (const FlatMatrix<T> m, FlatMatrix<T2> inv);


inline void CalcInverse (const double & m, double & inv)
{
  inv = 1 / m;
}

inline void CalcInverse (const Complex & m, Complex & inv)
{
  inv = 1.0 / m;
}

template <int H, int W, typename T, typename TINV>
inline void CalcInverse (const Mat<H,W,T> & m, TINV & inv)
{
  switch (H)
    {
    case 1:
      {
	inv(0,0) = 1.0 / m(0,0);
	return;
      }
    case 2:
      {
	T idet = 1.0 / (m(0,0) * m(1,1) - m(0,1) * m(1,0));
	inv(0,0) = idet * m(1,1);
	inv(0,1) = -idet * m(0,1);
	inv(1,0) = -idet * m(1,0);
	inv(1,1) = idet * m(0,0);
	return;
      }
    case 3:
      {
	T det = m(0) * m(4) * m(8)
	  + m(1) * m(5) * m(6)
	  + m(2) * m(3) * m(7)
	  - m(0) * m(5) * m(7)
	  - m(1) * m(3) * m(8)
	  - m(2) * m(4) * m(6);
	det = 1.0 / det;

	inv(0,0) =  det * (m(4) * m(8) - m(5) * m(7));
	inv(0,1) = -det * (m(1) * m(8) - m(2) * m(7));
	inv(0,2) =  det * (m(1) * m(5) - m(2) * m(4));
	
	inv(1,0) = -det * (m(3) * m(8) - m(5) * m(6));
	inv(1,1) =  det * (m(0) * m(8) - m(2) * m(6));
	inv(1,2) = -det * (m(0) * m(5) - m(2) * m(3));
	
	inv(2,0) =  det * (m(3) * m(7) - m(4) * m(6));
	inv(2,1) = -det * (m(0) * m(7) - m(1) * m(6));
	inv(2,2) =  det * (m(0) * m(4) - m(1) * m(3));
	return;
      }
    default:
      {
	FlatMatrix<T> fm(m);
	FlatMatrix<T> finv(inv);
	CalcInverse (fm, finv);
      }
    }
}


template <int H, int W, typename T>
inline Mat<H,W,T> Inv (const Mat<H,W,T> & m)
{
  Mat<H,W,T> inv;
  CalcInverse (m, inv);
  return inv;
}

template <typename T>
inline Matrix<T> Inv (const FlatMatrix<T> & m)
{
  Matrix<T> inv(m.Height(),m.Height());
  CalcInverse (m, inv);
  return inv;
}



/* ********************** Determinant *********************** */

template <class T>
inline typename T::TSCAL
Det (const MatExpr<T> & m)
{
  const T & sm = m.Spec();
  switch (sm.Height())
    {
    case 1: 
      {
	return sm(0,0); 
	break;
      }
    case 2:
      {
	return ( sm(0,0)*sm(1,1) - sm(0,1)*sm(1,0) ); 
	break;
      }
    case 3:
      {
	return 
	  sm(0) * sm(4) * sm(8)
	  + sm(1) * sm(5) * sm(6)
	  + sm(2) * sm(3) * sm(7)
	  - sm(0) * sm(5) * sm(7)
	  - sm(1) * sm(3) * sm(8)
	  - sm(2) * sm(4) * sm(6);
	break;
      }
    default:
      {
	cerr << "general det not implemented" << endl;
      }
    }
  return 0;
}


#endif