// @(#)root/smatrix:$Id$
// Authors: T. Glebe, L. Moneta    2005  

#ifndef ROOT_Math_MatrixFunctions
#define ROOT_Math_MatrixFunctions
// ********************************************************************
//
// source:
//
// type:      source code
//
// created:   20. Mar 2001
//
// author:    Thorsten Glebe
//            HERA-B Collaboration
//            Max-Planck-Institut fuer Kernphysik
//            Saupfercheckweg 1
//            69117 Heidelberg
//            Germany
//            E-mail: T.Glebe@mpi-hd.mpg.de
//
// Description: Functions/Operators special to Matrix
//
// changes:
// 20 Mar 2001 (TG) creation
// 20 Mar 2001 (TG) Added Matrix * Vector multiplication
// 21 Mar 2001 (TG) Added transpose, product
// 11 Apr 2001 (TG) transpose() speed improvment by removing rows(), cols()
//                  through static members of Matrix and Expr class
//
// ********************************************************************

//doxygen tag
/**
   @defgroup MatrixFunctions Matrix Template Functions 
   @ingroup SMatrixGroup

   These function apply to matrices (and also Matrix expression) and can return a 
   matrix expression of a particular defined type, like in the matrix multiplication or 
   a vector, like in the matrix-vector product or a scalar like in the Similarity vector-matrix product.
*/

#ifndef ROOT_Math_BinaryOpPolicy
#include "Math/BinaryOpPolicy.h"
#endif

namespace ROOT { 

  namespace Math { 



#ifdef XXX
//==============================================================================
// SMatrix * SVector
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
SVector<T,D1> operator*(const SMatrix<T,D1,D2,R>& rhs, const SVector<T,D2>& lhs)
{
  SVector<T,D1> tmp;
  for(unsigned int i=0; i<D1; ++i) {
    const unsigned int rpos = i*D2;
    for(unsigned int j=0; j<D2; ++j) {
      tmp[i] += rhs.apply(rpos+j) * lhs.apply(j);
    }
  }
  return tmp;
}
#endif


// matrix-vector product: 
// use apply(i) function for matrices. Tested  (11/05/06) with using (i,j) but 
// performances are slightly worse (not clear  why)

//==============================================================================
// meta_row_dot
//==============================================================================
template <unsigned int I>
struct meta_row_dot {
  template <class A, class B>
  static inline typename A::value_type f(const A& lhs, const B& rhs,
                                         const unsigned int offset) {
    return lhs.apply(offset+I) * rhs.apply(I) + meta_row_dot<I-1>::f(lhs,rhs,offset);
  }
};


//==============================================================================
// meta_row_dot<0>
//==============================================================================
template <>
struct meta_row_dot<0> {
  template <class A, class B>
  static inline typename A::value_type f(const A& lhs, const B& rhs,
                                         const unsigned int offset) {
    return lhs.apply(offset) * rhs.apply(0);
  }
};

//==============================================================================
// VectorMatrixRowOp
//==============================================================================
template <class Matrix, class Vector, unsigned int D2>
class VectorMatrixRowOp {
public:

  typedef typename Vector::value_type T; 

  ///
  VectorMatrixRowOp(const Matrix& lhs, const Vector& rhs) :
    lhs_(lhs), rhs_(rhs) {}

  ///
  ~VectorMatrixRowOp() {}

  /// calc \f$ \sum_{j} a_{ij} * v_j \f$
  inline typename Matrix::value_type apply(unsigned int i) const {
    return meta_row_dot<D2-1>::f(lhs_, rhs_, i*D2);
  }

   // check if passed pointer is in use  
   // check only the vector since this is a vector expression
  inline bool IsInUse (const T * p) const { 
    return rhs_.IsInUse(p);
  }


protected:
  const Matrix& lhs_;
  const Vector& rhs_;
};


//==============================================================================
// meta_col_dot
//==============================================================================
template <unsigned int I>
struct meta_col_dot {
  template <class Matrix, class Vector>
  static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
                                              const unsigned int offset) {
    return lhs.apply(Matrix::kCols*I+offset) * rhs.apply(I) + 
           meta_col_dot<I-1>::f(lhs,rhs,offset);
  }
};


//==============================================================================
// meta_col_dot<0>
//==============================================================================
template <>
struct meta_col_dot<0> {
  template <class Matrix, class Vector>
  static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
                                              const unsigned int offset) {
    return lhs.apply(offset) * rhs.apply(0);
  }
};

//==============================================================================
// VectorMatrixColOp
//==============================================================================
/**
   Class for Vector-Matrix multiplication

   @ingroup Expression
 */
template <class Vector, class Matrix, unsigned int D1>
class VectorMatrixColOp {
public:

  typedef typename Vector::value_type T; 
  ///
  VectorMatrixColOp(const Vector& lhs, const Matrix& rhs) :
    lhs_(lhs), rhs_(rhs) {}

  ///
  ~VectorMatrixColOp() {}

  /// calc \f$ \sum_{j} a_{ij} * v_j \f$
  inline typename Matrix::value_type apply(unsigned int i) const {
    return meta_col_dot<D1-1>::f(rhs_, lhs_, i);
  }

   // check if passed pointer is in use  
   // check only the vector since this is a vector expression   
   inline bool IsInUse (const T * p) const { 
    return lhs_.IsInUse(p);
  }


protected:
  const Vector&    lhs_;
  const Matrix&    rhs_;
};

/**
   Matrix *  Vector multiplication   \f$ a(i) = \sum_{j} M(i,j) * b(j) \f$
   returning a vector expression 

   @ingroup MatrixFunctions
 */
//==============================================================================
// operator*: SMatrix * SVector
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2>, T, D1>
 operator*(const SMatrix<T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
  typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: SMatrix * Expr<A,T,D2>
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>, VecExpr<A,T,D2>, D2>, T, D1>
 operator*(const SMatrix<T,D1,D2,R>& lhs, const VecExpr<A,T,D2>& rhs) {
  typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,VecExpr<A,T,D2>, D2> VMOp;
  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: Expr<A,T,D1,D2> * SVector
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, SVector<T,D2>, D2>, T, D1>
 operator*(const Expr<A,T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
  typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: Expr<A,T,D1,D2> * VecExpr<B,T,D2>
//==============================================================================
template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, VecExpr<B,T,D2>, D2>, T, D1>
 operator*(const Expr<A,T,D1,D2,R>& lhs, const VecExpr<B,T,D2>& rhs) {
  typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,VecExpr<B,T,D2>, D2> VMOp;
  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: SVector * SMatrix
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
 operator*(const SVector<T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
  typedef VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: SVector * Expr<A,T,D1,D2>
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1>, T, D2>
 operator*(const SVector<T,D1>& lhs, const Expr<A,T,D1,D2,R>& rhs) {
  typedef VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1> VMOp;
  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: VecExpr<A,T,D1> * SMatrix
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
 operator*(const VecExpr<A,T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
  typedef VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
}

//==============================================================================
// operator*: VecExpr<A,T,D1> * Expr<B,T,D1,D2>
//==============================================================================
template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1>, T, D2>
 operator*(const VecExpr<A,T,D1>& lhs, const Expr<B,T,D1,D2,R>& rhs) {
  typedef VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1> VMOp;
  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
}

//==============================================================================
// meta_matrix_dot
//==============================================================================
template <unsigned int I>
struct meta_matrix_dot {

  template <class MatrixA, class MatrixB>
  static inline typename MatrixA::value_type f(const MatrixA& lhs, 
                                               const MatrixB& rhs,
                                               const unsigned int offset) {
    return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols + I) *
           rhs.apply(MatrixB::kCols*I + offset%MatrixB::kCols) + 
           meta_matrix_dot<I-1>::f(lhs,rhs,offset);
  }

  // multiplication using i and j indeces
  template <class MatrixA, class MatrixB>
  static inline typename MatrixA::value_type g(const MatrixA& lhs, 
                                               const MatrixB& rhs,
                                               unsigned int i, 
                                               unsigned int j) {
    return lhs(i, I) * rhs(I , j) + 
           meta_matrix_dot<I-1>::g(lhs,rhs,i,j);
  }
};


//==============================================================================
// meta_matrix_dot<0>
//==============================================================================
template <>
struct meta_matrix_dot<0> {

  template <class MatrixA, class MatrixB>
  static inline typename MatrixA::value_type f(const MatrixA& lhs, 
                                               const MatrixB& rhs,
                                               const unsigned int offset) {
    return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols) *
           rhs.apply(offset%MatrixB::kCols);
  }

  // multiplication using i and j 
  template <class MatrixA, class MatrixB>
  static inline typename MatrixA::value_type g(const MatrixA& lhs, 
                                               const MatrixB& rhs,
                                               unsigned int i, unsigned int j) {
    return lhs(i,0) * rhs(0,j);
  }

};

//==============================================================================
// MatrixMulOp
//==============================================================================
/**
   Class for Matrix-Matrix multiplication 

   @ingroup Expression
 */
template <class MatrixA, class MatrixB, class T, unsigned int D>
class MatrixMulOp {
public:
  ///
  MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
    lhs_(lhs), rhs_(rhs) {}

  ///
  ~MatrixMulOp() {}

  /// calc \f$\sum_{j} a_{ik} * b_{kj}\f$
  inline T apply(unsigned int i) const {
    return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
  }

  inline T operator() (unsigned int i, unsigned j) const {
    return meta_matrix_dot<D-1>::g(lhs_, rhs_, i, j);
  }

  inline bool IsInUse (const T * p) const { 
    return lhs_.IsInUse(p) || rhs_.IsInUse(p);
  }


protected:
  const MatrixA&    lhs_;
  const MatrixB&    rhs_;
};


/**
   Matrix *  Matrix multiplication , \f$ C(i,j) = \sum_{k} A(i,k) * B(k,j)\f$ 
   returning a matrix expression 

   @ingroup MatrixFunctions
 */
//==============================================================================
// operator* (SMatrix * SMatrix, binary)
//==============================================================================
template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, T,D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,
    typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//==============================================================================
// operator* (SMatrix * Expr, binary)
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,
    typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//==============================================================================
// operator* (Expr * SMatrix, binary)
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,
    typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//==============================================================================
// operator* (Expr * Expr, binary)
//==============================================================================
template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}



#ifdef XXX
//==============================================================================
// MatrixMulOp
//==============================================================================
template <class MatrixA, class MatrixB, unsigned int D>
class MatrixMulOp {
public:
  ///
  MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
    lhs_(lhs), rhs_(rhs) {}

  ///
  ~MatrixMulOp() {}

  /// calc $\sum_{j} a_{ik} * b_{kj}$
  inline typename MatrixA::value_type apply(unsigned int i) const {
    return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
  }

protected:
  const MatrixA&    lhs_;
  const MatrixB&    rhs_;
};


//==============================================================================
// operator* (SMatrix * SMatrix, binary)
//==============================================================================
template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//==============================================================================
// operator* (SMatrix * Expr, binary)
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//==============================================================================
// operator* (Expr * SMatrix, binary)
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}

//=============================================================================
// operator* (Expr * Expr, binary)
//=============================================================================
template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
 operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D> MatMulOp;
  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
}
#endif

//==============================================================================
// TransposeOp
//==============================================================================
/**
   Class for Transpose Operations

   @ingroup Expression
 */
template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
class TransposeOp {
public:
  ///
  TransposeOp( const Matrix& rhs) :
    rhs_(rhs) {}

  ///
  ~TransposeOp() {}

  ///
  inline T apply(unsigned int i) const {
    return rhs_.apply( (i%D1)*D2 + i/D1);
  }
  inline T operator() (unsigned int i, unsigned j) const {
    return rhs_( j, i);
  }

  inline bool IsInUse (const T * p) const { 
    return rhs_.IsInUse(p); 
  }

protected:
  const Matrix& rhs_;
};


/**
   Matrix Transpose   B(i,j) = A(j,i) 
   returning a matrix expression

   @ingroup MatrixFunctions
 */
//==============================================================================
// transpose
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline Expr<TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
 Transpose(const SMatrix<T,D1,D2, R>& rhs) {
  typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;

  return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
}

//==============================================================================
// transpose
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
 Transpose(const Expr<A,T,D1,D2,R>& rhs) {
  typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;

  return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
}


#ifdef ENABLE_TEMPORARIES_TRANSPOSE 
// sometimes is faster to create a temp, not clear why

//==============================================================================
// transpose
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
 Transpose(const SMatrix<T,D1,D2, R>& rhs) {
  typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;

  return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType> 
    ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
}

//==============================================================================
// transpose
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline  SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
 Transpose(const Expr<A,T,D1,D2,R>& rhs) {
  typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;

  return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType> 
    ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
}

#endif 


#ifdef OLD
//==============================================================================
// product: SMatrix/SVector calculate v^T * A * v
//==============================================================================
template <class T, unsigned int D, class R>
inline T Product(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: SVector/SMatrix calculate v^T * A * v
//==============================================================================
template <class T, unsigned int D, class R>
inline T Product(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: SMatrix/Expr calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Product(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/SMatrix calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Product(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: SVector/Expr calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Product(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: Expr/SVector calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Product(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/Expr calculate v^T * A * v
//==============================================================================
template <class A, class B, class T, unsigned int D, class R>
inline T Product(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/Expr calculate v^T * A * v
//==============================================================================
template <class A, class B, class T, unsigned int D, class R>
inline T Product(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}
#endif

/**
   Similarity Vector - Matrix Product:  v^T * A * v
   returning a scalar value of type T   \f$ s = \sum_{i,j} v(i) * A(i,j) * v(j)\f$ 

   @ingroup MatrixFunctions
 */

//==============================================================================
// product: SMatrix/SVector calculate v^T * A * v
//==============================================================================
template <class T, unsigned int D, class R>
inline T Similarity(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: SVector/SMatrix calculate v^T * A * v
//==============================================================================
template <class T, unsigned int D, class R>
inline T Similarity(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: SMatrix/Expr calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Similarity(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/SMatrix calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Similarity(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: SVector/Expr calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Similarity(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}

//==============================================================================
// product: Expr/SVector calculate v^T * A * v
//==============================================================================
template <class A, class T, unsigned int D, class R>
inline T Similarity(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/Expr calculate v^T * A * v
//==============================================================================
template <class A, class B, class T, unsigned int D, class R>
inline T Similarity(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
  return Dot(rhs, lhs * rhs);
}

//==============================================================================
// product: Expr/Expr calculate v^T * A * v
//==============================================================================
template <class A, class B, class T, unsigned int D, class R>
inline T Similarity(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
  return Dot(lhs, rhs * lhs);
}


/**
   Similarity Matrix Product :  B = U * A * U^T for A symmetric 
   returning a symmetric matrix expression: 
   \f$ B(i,j) = \sum_{k,l} U(i,k) * A(k,l) * U(j,l) \f$ 

   @ingroup MatrixFunctions
 */
//==============================================================================
// product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
// return matrix will be nrows M x nrows M
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const SMatrix<T,D1,D2,R>& lhs, const SMatrix<T,D2,D2,MatRepSym<T,D2> >& rhs) {
  SMatrix<T,D1,D2, MatRepStd<T,D1,D2> > tmp = lhs * rhs;
  typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
  SMatrixSym mret; 
  AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  ); 
  return mret; 
}

//==============================================================================
// product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
// return matrix will be nrowsM x nrows M
// M is a matrix expression
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const Expr<A,T,D1,D2,R>& lhs, const SMatrix<T,D2,D2,MatRepSym<T,D2> >& rhs) {
  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = lhs * rhs;
  typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
  SMatrixSym mret; 
  AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  ); 
  return mret; 
}

#ifdef XXX
    // not needed (
//==============================================================================
// product: SMatrix/SMatrix calculate M * A * M where A and M are symmetric matrices
// return matrix will be nrows M x nrows M
//==============================================================================
template <class T, unsigned int D1>
inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const SMatrix<T,D1,D1,MatRepSym<T,D1> >& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
  SMatrix<T,D1,D1, MatRepStd<T,D1,D1> > tmp = lhs * rhs;
  typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
  SMatrixSym mret; 
  AssignSym::Evaluate(mret,  tmp * lhs  ); 
  return mret; 
}
#endif


/**
   Transpose Similarity Matrix Product :  B = U^T * A * U for A symmetric
   returning a symmetric matrix expression: \f$ B(i,j) = \sum_{k,l} U(k,i) * A(k,l) * U(l,j) \f$ 

   @ingroup MatrixFunctions
 */
//==============================================================================
// product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
// return matrix will be ncolsM x ncols M
//==============================================================================
template <class T, unsigned int D1, unsigned int D2, class R>
inline SMatrix<T,D2,D2,MatRepSym<T,D2> > SimilarityT(const SMatrix<T,D1,D2,R>& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
  typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym; 
  SMatrixSym mret; 
  AssignSym::Evaluate(mret,  Transpose(lhs) * tmp ); 
  return mret; 
}

//==============================================================================
// product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
// return matrix will be ncolsM x ncols M
// M is a matrix expression
//==============================================================================
template <class A, class T, unsigned int D1, unsigned int D2, class R>
inline SMatrix<T,D2,D2,MatRepSym<T,D2> > SimilarityT(const Expr<A,T,D1,D2,R>& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
  typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym; 
  SMatrixSym mret; 
  AssignSym::Evaluate(mret,  Transpose(lhs) * tmp ); 
  return mret; 
}





// //==============================================================================
// // Mult * (Expr * Expr, binary) with a symmetric result
// // the operation is done only for half
// //==============================================================================
// template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
// inline Expr<MatrixMulOp<Expr<A,T,D,D,MatRepSym<T,D> >, Expr<B,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
//  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
//   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
//   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
// }



//==============================================================================
// TensorMulOp
//==============================================================================
/**
   Class for Tensor Multiplication (outer product) of two vectors   
   giving a matrix 

   @ingroup Expression
 */
template <class Vector1, class Vector2>
class TensorMulOp {
public:
  ///
  TensorMulOp( const Vector1 & lhs, const Vector2 & rhs) :
    lhs_(lhs),
    rhs_(rhs) {}

  ///
  ~TensorMulOp() {}

  /// Vector2::kSize is the number of columns in the resulting matrix
  inline typename Vector1::value_type apply(unsigned int i) const {
    return lhs_.apply( i/ Vector2::kSize) * rhs_.apply( i % Vector2::kSize );
  }
  inline typename Vector1::value_type operator() (unsigned int i, unsigned j) const {
    return lhs_.apply(i) * rhs_.apply(j);
  }

  inline bool IsInUse (const typename Vector1::value_type * ) const { 
    return false; 
  }


protected:

  const Vector1 & lhs_;
  const Vector2 & rhs_;

};



/**
   Tensor Vector Product : M(i,j) = v(i) * v(j) 
   returning a matrix expression 

   @ingroup VectFunction
 */

#ifndef _WIN32

    // Tensor Prod (use default MatRepStd for the returned expression
    // cannot make a symmetric matrix
//==============================================================================
// TensorProd (SVector x SVector)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline Expr<TensorMulOp<SVector<T,D1>, SVector<T,D2>  >, T, D1, D2 >
 TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
  typedef TensorMulOp<SVector<T,D1>, SVector<T,D2> > TVMulOp;
  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
}

//==============================================================================
// TensorProd (VecExpr x SVector)
//==============================================================================
 template <class T, unsigned int D1, unsigned int D2, class A>
inline Expr<TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2>  >, T, D1, D2 >
 TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
  typedef TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> > TVMulOp;
  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
}

//==============================================================================
// TensorProd (SVector x VecExpr)
//==============================================================================
 template <class T, unsigned int D1, unsigned int D2, class A>
inline Expr<TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2>  >, T, D1, D2 >
 TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
  typedef TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> > TVMulOp;
  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
}


//==============================================================================
// TensorProd (VecExpr x VecExpr)
//==============================================================================
 template <class T, unsigned int D1, unsigned int D2, class A, class B>
inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2>  >, T, D1, D2 >
 TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
  typedef TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> > TVMulOp;
  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
}

#endif
#ifdef _WIN32 
/// case of WINDOWS - problem using Expression (  C1001: INTERNAL COMPILER ERROR )

//==============================================================================
// TensorProd (SVector x SVector)
//==============================================================================
template <class T, unsigned int D1, unsigned int D2>
inline SMatrix<T,D1,D2>  TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
  SMatrix<T,D1,D2> tmp;
  for (unsigned int i=0; i< D1; ++i) 
    for (unsigned int j=0; j< D2; ++j) {
      tmp(i,j) = lhs[i]*rhs[j];
    }

  return tmp;
}
//==============================================================================
// TensorProd (VecExpr x SVector)
//==============================================================================
 template <class T, unsigned int D1, unsigned int D2, class A>
inline SMatrix<T,D1,D2>  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
  SMatrix<T,D1,D2> tmp;
  for (unsigned int i=0; i< D1; ++i)
    for (unsigned int j=0; j< D2; ++j)
      tmp(i,j) = lhs.apply(i) * rhs.apply(j);

  return tmp;
}
//==============================================================================
// TensorProd (SVector x VecExpr)
//==============================================================================
 template <class T, unsigned int D1, unsigned int D2, class A>
inline SMatrix<T,D1,D2> TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
  SMatrix<T,D1,D2> tmp;
  for (unsigned int i=0; i< D1; ++i)
    for (unsigned int j=0; j< D2; ++j)
      tmp(i,j) = lhs.apply(i) * rhs.apply(j);

  return tmp;
}

//==============================================================================
// TensorProd (VecExpr x VecExpr)
//==============================================================================

 template <class T, unsigned int D1, unsigned int D2, class A, class B>
inline SMatrix<T,D1,D2  > TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
  SMatrix<T,D1,D2> tmp;
  for (unsigned int i=0; i< D1; ++i)
    for (unsigned int j=0; j< D2; ++j)
      tmp(i,j) = lhs.apply(i) * rhs.apply(j);

  return tmp;
}


#endif

// solving a positive defined symmetric linear system using Choleski decompositions
// matrix will be decomposed and the returned vector will be overwritten in vec
// If the user wants to pass const objects need to copy the matrices
// It will work only for symmetric matrices
template <class T, unsigned int D>
bool SolveChol( SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  SVector<T, D> & vec ) { 
   CholeskyDecomp<T, D> decomp(mat);
   return decomp.Solve(vec); 
}

/// same function as before but not overwriting the matrix and returning a copy of the vector
/// (this is the slow version)
template <class T, unsigned int D>
SVector<T,D> SolveChol( const SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  const SVector<T, D> & vec, int & ifail  ) { 
   SMatrix<T, D, D, MatRepSym<T, D> > atmp(mat); 
   SVector<T,D> vret(vec); 
   bool ok = SolveChol( atmp, vret); 
   ifail =  (ok) ? 0 : -1;
   return vret; 
}



  }  // namespace Math

}  // namespace ROOT
          

#endif  /* ROOT_Math_MatrixFunctions */