TinyMatrixInline.cxx

#ifndef MONTJOIE_FILE_TINY_MATRIX_INLINE_CXX
 
#include "TinyMatrix.hxx"
 
#include "TinyMatrixExpressionInline.cxx"
 
namespace Seldon
{
  
  /***********************
   * TinyMatrix<General> *
   ***********************/
  
  
  template<class T, int m, int n>
  inline TinyMatrix<T, General, m, n>::TinyMatrix()
  {
    this->Zero();
  }
  
  
  template<class T, int m, int n>
  inline TinyMatrix<T, General, m, n>::TinyMatrix(int i)
  {
    this->Zero();
  }
 
 
  template<class T, int m, int n>
  inline TinyMatrix<T, General, m, n>::TinyMatrix(int i, int j)
  {
    this->Zero();
  }
 
 
  template<class T, int m, int n> template<class T0>
  inline TinyMatrix<T, General, m, n>::
  TinyMatrix(const TinyVector<TinyVector<T0, m>, n>& A)
  {
    TinyMatrixLoop<n>::Init(A, *this);
  }
  
 
  template<class T, int m, int n> template<class E>
  inline TinyMatrix<T, General, m, n>
  ::TinyMatrix(const TinyMatrixExpression<T, m, n, E>& A)
  {
    TinyMatrixLoop<m*n>::Copy(A, *this);
  }
  
  
  template<class T, int m, int n>
  inline int TinyMatrix<T, General, m, n>::GetM()
  {
    return m;
  }
  
  
  template<class T, int m, int n>
  inline int TinyMatrix<T, General, m, n>::GetN()
  {
    return n;
  }
  
  
  template<class T, int m, int n>
  inline int TinyMatrix<T, General, m, n>::GetSize()
  {
    return m*n;
  }
  
 
  template<class T, int m, int n>
  inline T* TinyMatrix<T, General, m, n>::GetData()
  {
    return data_;
  }  
 
 
  template<class T, int m, int n> template<class T1>
  inline TinyMatrix<T, General, m, n> & TinyMatrix<T, General, m, n>::
  operator *=(const T1 & a )
  {
    TinyMatrixLoop<m*n>::MltScal(a, *this);    
    return *this;
  }
  
  
  template<class T, int m, int n> template<class T1, class E>
  inline TinyMatrix<T, General, m, n>& TinyMatrix<T, General, m, n>::
  operator +=(const TinyMatrixExpression<T1, m, n, E> & B) 
  {
    TinyMatrixLoop<m*n>::AddCopy(B, *this);
    return *this;
  }
  
  
  template<class T, int m, int n> template<class T1, class E>
  inline TinyMatrix<T, General, m, n>& TinyMatrix<T, General, m, n>::
  operator -=(const TinyMatrixExpression<T1, m, n, E> & B) 
  {
    TinyMatrixLoop<m*n>::DiffCopy(B, *this);
    return *this;
  }
 
 
  template<class T, int m, int n>
  inline TinyMatrix<T, General, m, n> & TinyMatrix<T, General, m, n>::operator =(const T& x)
  {
    this->Fill(x);
    return *this;
  }
 
  
  template<class T, int m, int n> template<class T0, class E>
  inline TinyMatrix<T, General, m, n> & TinyMatrix<T, General, m, n>::
  operator =(const TinyMatrixExpression<T0, m, n, E>& A)
  {
    TinyMatrixLoop<m*n>::Copy(A, *this);
    return *this;
  }
  
  
  template<class T, int m, int n>
  inline T& TinyMatrix<T, General, m, n>::operator()(int a, int b)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(a, b, m, n, "TinyMatrix");
#endif
    
    return this->data_[a*n+b];
  }
  
  
  template<class T, int m, int n>
  inline const T& TinyMatrix<T, General, m, n>::operator()(int a, int b) const
  {
    
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(a, b, m, n, "TinyMatrix");
#endif
    
    return this->data_[a*n+b];
  }
 
 
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::Zero()
  {
    TinyMatrixLoop<m*n>::Zero(*this);
  }
  
 
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::SetIdentity()
  {
    
#ifdef SELDON_CHECK_BOUNDS
    if (m != n)
      throw WrongDim("TinyMatrix::SetIdentity()",
                     "The number of rows should be equal to the number of columns");
#endif
    
    TinyMatrixLoop<m*n>::SetIdentity(*this);
  }
  
 
  template<class T, int m, int n> template<class T0>
  inline void TinyMatrix<T, General, m, n>::SetDiagonal(const T0& alpha)
  {
    
#ifdef SELDON_CHECK_BOUNDS
    if (m != n)
      throw WrongDim("TinyMatrix::SetDiagonal(alpha)",
                     "The number of rows should be equal to the number of columns");
#endif
    
    TinyMatrixLoop<m*n>::SetDiagonal(*this, alpha);
  }
 
 
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::Fill()
  {
    TinyMatrixLoop<m*n>::Fill(*this);
  }
  
  
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::FillRand()
  {
    TinyMatrixLoop<m*n>::FillRand(*this);
  }
 
 
  template<class T, int m, int n> template<class T0>
  inline void TinyMatrix<T, General, m, n>::Fill(const T0& a)
  {
    TinyMatrixLoop<m*n>::Fill(*this, a);
  }
  
 
  template<class T, int m, int n>
  inline bool TinyMatrix<T, General, m, n>::IsZero() const
  {
    return TinyMatrixLoop<m*n>::IsZero(*this);
  }
  
 
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::Write(const string& file_name) const
  {
    ofstream out(file_name.data());
    Write(out);
    out.close();
  }
  
  
  template<class T, int m, int n>
  inline void TinyMatrix<T, General, m, n>::Write(ostream& out) const
  {
    int itmp = m;
    out.write(reinterpret_cast<char*>(&itmp), sizeof(int));
    itmp = n;
    out.write(reinterpret_cast<char*>(&itmp), sizeof(int));
    
    return TinyMatrixLoop<m>::Write(out, *this);    
  }
 
 
  /************************
   * TinyMatrixTripleLoop *
   ************************/
 
  
  template<int i, int j, int k> template<class T, int m>
  inline void TinyMatrixTripleLoop<i, j, k>::
  AddRow(TinyMatrix<T, General, m, m>& A, const T& val)
  {
    A(j, k-1) += val*A(i, k-1);
    TinyMatrixTripleLoop<i, j, k-1>::AddRow(A, val);
  }
  
 
  template<int j, int k, int m> template<class T, int p>
  inline void TinyMatrixTripleLoop<j, k, m>
  ::ModifyUpperCholesky(TinyMatrix<T, Symmetric, p, p>& A, T& invVal, T& vloc)
  {
    vloc -= A(k, m-1)*A(m-1, j);
    TinyMatrixTripleLoop<j, k, m-1>::ModifyUpperCholesky(A, invVal, vloc);    
  }
  
  
  /************************
   * TinyMatrixDoubleLoop *
   ************************/
 
  
  template<int p, int q> template<int m, int n, class T0, class T1>
  inline void TinyMatrixDoubleLoop<p, q>::
  Init(const TinyVector<T0, m>& x, TinyMatrix<T1, General, m, n>& A)
  {
    A(q-1, p-1) = x(q-1);
    TinyMatrixDoubleLoop<p, q-1>::Init(x, A);
  }
    
  
  template<int p, int q> template<int m, int n, class T, class E>
  inline void TinyMatrixDoubleLoop<p, q>::
  WriteText(ostream& out, const TinyMatrixExpression<T, m, n, E>& A)
  {
    out << A(m-p, n-q) << "  ";
    TinyMatrixDoubleLoop<p, q-1>::WriteText(out, A);
  }
  
  
  template<int p, int q> template<int m, int n, class T, class Prop>
  inline void TinyMatrixDoubleLoop<p, q>::
  Write(ostream& out, const TinyMatrix<T, Prop, m, n>& A)
  {    
    out.write(reinterpret_cast<const char*>(&A(m-p, n-q)), sizeof(T));
    
    TinyMatrixDoubleLoop<p, q-1>::Write(out, A);
  }
  
 
  template<int p, int q> template<int m, int n, class T, class E1, class E2>
  inline bool TinyMatrixDoubleLoop<p, q>
  ::IsEqual(const TinyMatrixExpression<T, m, n, E1>& A,
            const TinyMatrixExpression<T, m, n, E2>& B)
  {
    if (abs(A(p-1, q-1) - B(p-1, q-1)) > TinyVector<T, m>::threshold)
      return false;
 
    return TinyMatrixDoubleLoop<p, q-1>::IsEqual(A, B);
  }
  
 
  template<int i, int j> template<int m, int n, class T0, class E0,
                                  class T1, class E1, class T2>
  inline void TinyMatrixDoubleLoop<i, j>::
  Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
      const TinyVectorExpression<T1, n, E1>& x, T2& val)
  {
    val += A(i-1, j)*x(j);
    TinyMatrixDoubleLoop<i, j-1>::Mlt(A, x, val);
  }
  
 
  template<int p, int q>
  template<int m, int n, int k, class T0, class E0,
           class T1, class E1, class T2, class Prop2>
  inline void TinyMatrixDoubleLoop<p, q>::
  Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
      const TinyMatrixExpression<T1, n, k, E1>& B,
      TinyMatrix<T2, Prop2, m, k>& C)
  {
    C(TinyMatrixNode<Prop2, m, k, p-1 >::i,
      TinyMatrixNode<Prop2, m, k, p-1 >::j)
      += A(TinyMatrixNode<Prop2, m, k, p-1 >::i, q)
           *B(q, TinyMatrixNode<Prop2, m, k, p-1 >::j);
    
    TinyMatrixDoubleLoop<p, q-1>::Mlt(A, B, C);
  }
  
 
  template<int p, int q> 
  template<int m, int n, class T0, class T1, class E1, class T2, class E2, class T3>
  inline void TinyMatrixDoubleLoop<p, q>::
  Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
              const TinyVectorExpression<T2, n, E2>& y, TinyMatrix<T3, General, m, n>& A)
  {
    A(p-1, q-1) += alpha*x(p-1)*y(q-1);
    TinyMatrixDoubleLoop<p, q-1>::Rank1Update(alpha, x, y, A);
  }
 
 
  template<int p, int q> 
  template<int m, int n, class T1, class E1, class T2, class E2, class T3>
  inline void TinyMatrixDoubleLoop<p, q>::
  Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
              const TinyVectorExpression<T2, n, E2>& y, TinyMatrix<T3, General, m, n>& A)
  {
    A(p-1, q-1) = x(p-1)*y(q-1);
    TinyMatrixDoubleLoop<p, q-1>::Rank1Matrix(x, y, A);
  }
  
 
  template<int p, int q> 
  template<int m, class T0, class T1, class E1, class T3>
  inline void TinyMatrixDoubleLoop<p, q>::
  Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
              TinyMatrix<T3, Symmetric, m, m>& A)
  {
    A(p-1, q-1) += alpha*x(p-1)*x(q-1);
    TinyMatrixDoubleLoop<p, q-1>::Rank1Update(alpha, x, A);
  }
 
 
  template<int p, int q> 
  template<int m, class T1, class E1, class T3>
  inline void TinyMatrixDoubleLoop<p, q>::
  Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
              TinyMatrix<T3, Symmetric, m, m>& A)
  {
    A(p-1, q-1) = x(p-1)*x(q-1);
    TinyMatrixDoubleLoop<p, q-1>::Rank1Matrix(x, A);
  }
  
 
  template<int i, int j> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j>::GetMaximumColumn(TinyMatrix<T, General, m, m>& A,
                                                           int& jmax, T& val)
  {
    if (abs(A(j, i)) > abs(val))
      {
        jmax = j;
        val = A(j, i);
      }
    
    TinyMatrixDoubleLoop<i, (j-1)*(j-1 > i)>::GetMaximumColumn(A, jmax, val);
  }
  
  
  template<int i, int j> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j>::
  SwapRow(TinyMatrix<T, General, m, m>& A, int i2, T& val)
  {
    val = A(i, j-1);
    A(i, j-1) = A(i2, j-1);
    A(i2, j-1) = val;
    TinyMatrixDoubleLoop<i, j-1>::SwapRow(A, i2, val);
  }
 
  
  template<int i, int j> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j>::
  SwapColumn(TinyMatrix<T, General, m, m>& A, int i2, T& val)
  {
    val = A(j-1, i);
    A(j-1, i) = A(j-1, i2);
    A(j-1, i2) = val;
    TinyMatrixDoubleLoop<i, j-1>::SwapColumn(A, i2, val);
  }
 
  
  template<int i, int j> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j>::
  MltRow(TinyMatrix<T, General, m, m>& A, const T& coef)
  {
    A(i, j-1) *= coef;
    TinyMatrixDoubleLoop<i, j-1>::MltRow(A, coef);
  }
  
 
  template<int i, int j> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j>::
  PerformElimination(TinyMatrix<T, General, m, m>& A, const T& coef, T& val)
  {
    // val =  -a_ji / a_ii
    val = -A(j, i)*coef;
    
    // replacing Lj by Lj + val * Li
    TinyMatrixTripleLoop<i, j, m>::AddRow(A, val);
    
    // storing val in a_ij
    A(j, i) = val;
    
    TinyMatrixDoubleLoop<i, (j-1)*(j-1>i) >::PerformElimination(A, coef, val);
  }
  
 
  template<int i, int j1> template<class T, int m>
  inline void TinyMatrixDoubleLoop<i, j1>::
  PerformSolve(TinyMatrix<T, General, m, m>& A, T& val)
  {
    // j = m - 1 - j1 + i + 1
    // val =  -a_ij
    val = -A(i, m - 1 - j1 + i + 1);
    
    // replacing Li by Li + val * Lj
    TinyMatrixTripleLoop<m - 1 - j1 + i + 1, i, m>::AddRow(A, val);
    
    A(i, m - 1 - j1 + i + 1) = val*A(m - 1 - j1 + i + 1, m - 1 - j1 + i + 1);
    TinyMatrixDoubleLoop<i, (j1-1)*(j1-1>i) >::PerformSolve(A, val);
  }
  
 
  template<int j, int k> template<class T, int m>
  inline void TinyMatrixDoubleLoop<j, k>
  ::GetDiagonalCholesky(TinyMatrix<T, Symmetric, m, m>& A, T& val)
  {
    val -= A(k-1, j)*A(k-1, j);
    TinyMatrixDoubleLoop<j, k-1>::GetDiagonalCholesky(A, val);
  }
 
 
  template<int j, int k> template<class T, int m>
  inline void TinyMatrixDoubleLoop<j, k>
  ::ModifyUpperCholesky(TinyMatrix<T, Symmetric, m, m>& A, T& invVal, T& vloc)
  {
    vloc = A(j, m-k);
    TinyMatrixTripleLoop<j, m-k, j>::ModifyUpperCholesky(A, invVal, vloc);
    
    A(j, m-k) = vloc*invVal;
    TinyMatrixDoubleLoop<j, k-1>::ModifyUpperCholesky(A, invVal, vloc);
  }
  
 
  template<int i, int k> template<class T, class T2, int m>
  inline void TinyMatrixDoubleLoop<i, k>
  ::SolveCholesky(const class_SeldonNoTrans& trans,
                  const TinyMatrix<T, Symmetric, m, m>& A,
                  TinyVector<T2, m>& x)
  {
    x(i+k) -= A(i, i+k)*x(i);
    TinyMatrixDoubleLoop<i, k-1>::SolveCholesky(trans, A, x);
  }
  
 
  template<int i, int k> template<class T, class T2, int m>
  inline void TinyMatrixDoubleLoop<i, k>
  ::SolveCholesky(const class_SeldonTrans& trans,
                  const TinyMatrix<T, Symmetric, m, m>& A,
                  TinyVector<T2, m>& x, T2& val)
  {
    val -= A(i, i+k)*x(i+k);
    TinyMatrixDoubleLoop<i, k-1>::SolveCholesky(trans, A, x, val);
  }
  
 
  template<int i, int k> template<class T, class T2, int m>
  inline void TinyMatrixDoubleLoop<i, k>
  ::MltCholesky(const class_SeldonNoTrans& trans,
                const TinyMatrix<T, Symmetric, m, m>& A,
                TinyVector<T2, m>& x)
  {
    x(i+k) += A(i, i+k)*x(i);
    TinyMatrixDoubleLoop<i, k-1>::MltCholesky(trans, A, x);
  }
  
 
  template<int i, int k> template<class T, class T2, int m>
  inline void TinyMatrixDoubleLoop<i, k>
  ::MltCholesky(const class_SeldonTrans& trans,
                const TinyMatrix<T, Symmetric, m, m>& A,
                TinyVector<T2, m>& x, T2& val)
  {
    val += A(i, i+k)*x(i+k);
    TinyMatrixDoubleLoop<i, k-1>::MltCholesky(trans, A, x, val);
  }
  
  
  /******************
   * TinyMatrixLoop *
   ******************/
  
 
  template<int p> template<int m, int n, class T, class Prop>
  inline void TinyMatrixLoop<p>::Zero(TinyMatrix<T, Prop, m, n>& A)
  {
    FillZero(A.data_[p-1]);
    TinyMatrixLoop<p-1>::Zero(A);
  }
  
  
  template<int p> template<int m, int n, class T0, class T1>
  inline void TinyMatrixLoop<p>::
  Init(const TinyVector<TinyVector<T0, m>, n>& x, TinyMatrix<T1, General, m, n>& A)
  {
    TinyMatrixDoubleLoop<p, m>::Init(x(p-1), A);
    TinyMatrixLoop<p-1>::Init(x, A);
  }
 
  
  template<int k> template<int m, int n, class T0,
                           class Prop, class T1>
  inline void TinyMatrixLoop<k>::MltScal(const T0& alpha,
                                         TinyMatrix<T1, Prop, m, n>& A)
  {
    A.data_[k-1] *= alpha;
    TinyMatrixLoop<k-1>::MltScal(alpha, A);
  }
  
 
  template<int p> template<int m, int n, class T1, class E, class T0, class Prop>
  void TinyMatrixLoop<p>::Copy(const TinyMatrixExpression<T1, m, n, E>& x,
                               TinyMatrix<T0, Prop, m, n>& y)
  {
    y(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j)
      = x(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j);
 
    TinyMatrixLoop<p-1>::Copy(x, y);
  }
 
 
  template<int p> template<int m, int n, class T1, class E, class T0, class Prop>
  void TinyMatrixLoop<p>::AddCopy(const TinyMatrixExpression<T1, m, n, E>& x,
                                  TinyMatrix<T0, Prop, m, n>& y)
  {
    y(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j)
      += x(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j);
    
    TinyMatrixLoop<p-1>::AddCopy(x, y);
  }
 
 
  template<int p> template<int m, int n, class T1, class E, class T0, class Prop>
  void TinyMatrixLoop<p>::DiffCopy(const TinyMatrixExpression<T1, m, n, E>& x,
                                   TinyMatrix<T0, Prop, m, n>& y)
  {
    y(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j)
      -= x(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j);
    
    TinyMatrixLoop<p-1>::DiffCopy(x, y);
  }
    
 
  template<int p> template<int m, int n, class T, class Prop>
  inline void TinyMatrixLoop<p>::SetIdentity(TinyMatrix<T, Prop, m, n>& A)
  {
    SetComplexReal((TinyMatrixNode<Prop, m, n, p-1>::i == TinyMatrixNode<Prop, m, n, p-1>::j),
                   A(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j));
    
    return TinyMatrixLoop<p-1>::SetIdentity(A);
  }
  
 
  template<int p> template<int m, int n, class T, class Prop, class T0>
  inline void TinyMatrixLoop<p>::SetDiagonal(TinyMatrix<T, Prop, m, n>& A, const T0& alpha)
  {
    if (TinyMatrixNode<Prop, m, n, p-1>::i == TinyMatrixNode<Prop, m, n, p-1>::j)
      SetComplexReal(alpha, A(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j));
    else
      SetComplexZero(A(TinyMatrixNode<Prop, m, n, p-1>::i, TinyMatrixNode<Prop, m, n, p-1>::j));
    
    return TinyMatrixLoop<p-1>::SetDiagonal(A, alpha);
  }
 
 
  template<int p> template<int m, int n, class T, class Prop>
  inline void TinyMatrixLoop<p>::Fill(TinyMatrix<T, Prop, m, n>& A)
  {
    SetComplexReal(p-1, A.data_[p-1]);
    TinyMatrixLoop<p-1>::Fill(A);
  }
 
  
  template<int p> template<int m, int n, class T, class Prop>
  inline void TinyMatrixLoop<p>::FillRand(TinyMatrix<T, Prop, m, n>& A)
  {
    SetComplexReal(rand(), A.data_[p-1]);
    TinyMatrixLoop<p-1>::FillRand(A);
  }
 
 
  template<int p> template<int m, int n, class T, class Prop, class T0>
  inline void TinyMatrixLoop<p>::Fill(TinyMatrix<T, Prop, m, n>& A, const T0& alpha)
  {
    SetComplexReal(alpha, A.data_[p-1]);
    TinyMatrixLoop<p-1>::Fill(A, alpha);
  }
 
 
  template<int p> template<int m, int n, class Prop, class T>
  inline bool TinyMatrixLoop<p>::IsZero(const TinyMatrix<T, Prop, m, n>& A)
  {
    if (!IsComplexZero(A.data_[p-1]))
      return false;
    
    return TinyMatrixLoop<p-1>::IsZero(A);
  }
  
  
  template<int p> template<int m, int n, class T, class E>
  inline void TinyMatrixLoop<p>::
  WriteText(ostream& out, const TinyMatrixExpression<T, m, n, E>& A)
  {
    TinyMatrixDoubleLoop<p, n>::WriteText(out, A);
    out << '\n';
    TinyMatrixLoop<p-1>::WriteText(out, A);
  }
 
  
  template<int p> template<int m, int n, class T, class Prop>
  inline void TinyMatrixLoop<p>::Write(ostream& out, const TinyMatrix<T, Prop, m, n>& A)
  {
    TinyMatrixDoubleLoop<p, n>::Write(out, A);
    TinyMatrixLoop<p-1>::Write(out, A);
  }
 
 
  template<int p> template<int m, int n, class T, class E1, class E2>
  inline bool TinyMatrixLoop<p>
  ::IsEqual(const TinyMatrixExpression<T, m, n, E1>& A,
            const TinyMatrixExpression<T, m, n, E2>& B)
  {
    if (!TinyMatrixDoubleLoop<p, n>::IsEqual(A, B))
      return false;
    
    return TinyMatrixLoop<p-1>::IsEqual(A, B);
  }
 
 
  template<int i> template<int m, int n, class T0, class E0,
                           class T1, class E1, class T2>
  inline void TinyMatrixLoop<i>::
  Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
      const TinyVectorExpression<T1, n, E1>& x, TinyVector<T2, m>& y)
  {
    y(i-1) = A(i-1, 0)*x(0);
    TinyMatrixDoubleLoop<i, n-1>::Mlt(A, x, y(i-1));
    TinyMatrixLoop<i-1>::Mlt(A, x, y);
  }
  
 
  template<int i> template<int m, int n, class T0, class E0,
                           class T1, class E1, class T2>
  inline void TinyMatrixLoop<i>::
  MltAdd(const TinyMatrixExpression<T0, m, n, E0>& A,
         const TinyVectorExpression<T1, n, E1>& x, TinyVector<T2, m>& y)
  {
    y(i-1) += A(i-1, 0)*x(0);
    TinyMatrixDoubleLoop<i, n-1>::Mlt(A, x, y(i-1));
    TinyMatrixLoop<i-1>::MltAdd(A, x, y);
  }
  
 
  template<int i> template<int m, int n, class T0, class E0,
                           class T1, class E1, class T2, class T3>
  inline void TinyMatrixLoop<i>::
  MltAdd(const T3& alpha, const TinyMatrixExpression<T0, m, n, E0>& A,
         const TinyVectorExpression<T1, n, E1>& x, TinyVector<T2, m>& y)
  {
    T2 val = A(i-1, 0)*x(0);
    TinyMatrixDoubleLoop<i, n-1>::Mlt(A, x, val);
    y(i-1) += alpha*val;
    TinyMatrixLoop<i-1>::MltAdd(alpha, A, x, y);
  }
 
  
  template<int p> template<int m, int n, int k, class T0, class E0,
                           class T1, class E1, class T2, class Prop2>
  inline void TinyMatrixLoop<p>::
  Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
      const TinyMatrixExpression<T1, n, k, E1>& B,
      TinyMatrix<T2, Prop2, m, k>& C)
  {
    C(TinyMatrixNode<Prop2, m, k, p-1 >::i, TinyMatrixNode<Prop2, m, k, p-1 >::j)
      = A(TinyMatrixNode<Prop2, m, k, p-1 >::i, 0)*B(0, TinyMatrixNode<Prop2, m, k, p-1 >::j);
    
    TinyMatrixDoubleLoop<p, n-1>::Mlt(A, B, C);
    TinyMatrixLoop<p-1>::Mlt(A, B, C);
  }
  
 
  template<int p> template<int m, int n, class T0, class T1, class E1,
                           class T2, class E2, class T3>
  inline void TinyMatrixLoop<p>::
  Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
              const TinyVectorExpression<T2, n, E2>& y, TinyMatrix<T3, General, m, n>& A)
  {
    TinyMatrixDoubleLoop<p, n>::Rank1Update(alpha, x, y, A);
    TinyMatrixLoop<p-1>::Rank1Update(alpha, x, y, A);
  }
 
 
  template<int p> template<int m, int n, class T1, class E1,
                           class T2, class E2, class T3>
  inline void TinyMatrixLoop<p>::
  Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
              const TinyVectorExpression<T2, n, E2>& y, TinyMatrix<T3, General, m, n>& A)
  {
    TinyMatrixDoubleLoop<p, n>::Rank1Matrix(x, y, A);
    TinyMatrixLoop<p-1>::Rank1Matrix(x, y, A);
  }
  
 
  template<int p> template<int m, class T0, class T1, class E1, class T3>
  inline void TinyMatrixLoop<p>::
  Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
              TinyMatrix<T3, Symmetric, m, m>& A)
  {
    TinyMatrixDoubleLoop<p, p>::Rank1Update(alpha, x, A);
    TinyMatrixLoop<p-1>::Rank1Update(alpha, x, A);
  }
 
 
  template<int p> template<int m, class T1, class E1, class T3>
  inline void TinyMatrixLoop<p>::
  Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
              TinyMatrix<T3, Symmetric, m, m>& A)
  {
    TinyMatrixDoubleLoop<p, p>::Rank1Matrix(x, A);
    TinyMatrixLoop<p-1>::Rank1Matrix(x, A);
  }
 
  
  template<int p> template<int m, int n, class T1, class E1>
  inline void TinyMatrixLoop<p>::
  GetCol(const TinyMatrixExpression<T1, m, n, E1>& A, int k, TinyVector<T1, m>& x)
  {
    x(p-1) = A(p-1, k);
    TinyMatrixLoop<p-1>::GetCol(A, k ,x);
  }
  
  
  template<int p> template<int m, int n, class T0, class E0, class T1>
  inline void TinyMatrixLoop<p>::
  GetRow(const TinyMatrixExpression<T0, m, n, E0>& A, int k, TinyVector<T1, n>& x)
  {
    x(p-1) = A(k, p-1);
    TinyMatrixLoop<p-1>::GetRow(A, k ,x);
  }
  
  
  template<int p> template<int m, int n, class T1, class E1, class Prop>
  inline void TinyMatrixLoop<p>::
  SetCol(const TinyVectorExpression<T1, m, E1>& x, int k, TinyMatrix<T1, Prop, m, n>& A)
  {
    A(p-1, k) = x(p-1);
    TinyMatrixLoop<p-1>::SetCol(x, k, A);
  }
  
  
  template<int p> template<int m, int n, class T1, class E1, class Prop>
  inline void TinyMatrixLoop<p>::
  SetRow(const TinyVectorExpression<T1, n, E1>& x, int k, TinyMatrix<T1, Prop, m, n>& A)
  {
    A(k,p-1) = x(p-1);
    TinyMatrixLoop<p-1>::SetRow(x, k, A);
  }
  
 
  template<int p> template<int m, int n, class T, class Prop, class T0>
  inline void TinyMatrixLoop<p>::GetMaxAbs(const TinyMatrix<T, Prop, m, n>& A, T0& amax)
  {
    amax = max(amax, abs(A.data_[p-1]));
    TinyMatrixLoop<p-1>::GetMaxAbs(A, amax);
  }
  
 
  template<int i1> template<class T, int m>
  inline void TinyMatrixLoop<i1>::
  PivotGauss(TinyMatrix<T, General, m, m>& A, TinyVector<int, m>& pivot)
  {
    // finding the maximum coefficient in the column i
    // i = m -1 - i1
    int jmax = m -1 - i1; T val = A(m -1 - i1, m -1 - i1), coef;
    TinyMatrixDoubleLoop<m -1 - i1, m-1>::GetMaximumColumn(A, jmax, val);
    pivot(m -1 - i1) = jmax;
    
    // now swapping rows i and jmax
    if (m -1 - i1 != jmax)
      TinyMatrixDoubleLoop<m -1 - i1, m>::SwapRow(A, jmax, val);
    
    // loop over each row i+1 -> m, to make elimination (we solve by L as well)
    T one; SetComplexOne(one);
    coef = one/A(m -1 - i1, m -1 - i1);
    TinyMatrixDoubleLoop<m -1 - i1, m-1>::PerformElimination(A, coef, val);
    
    // multiplying row i with coef, and storing coef in a_ii
    TinyMatrixDoubleLoop<m - 1 - i1, m>::MltRow(A, coef);
    A(m -1 - i1, m -1 - i1) = coef;
    
    TinyMatrixLoop<i1-1>::PivotGauss(A, pivot);
  }
 
 
  template<int i1> template<class T, int m>
  inline void TinyMatrixLoop<i1>::SolveUpper(TinyMatrix<T, General, m, m>& A)
  {
    // i = m -1 - i1  
    T val;
    TinyMatrixDoubleLoop<m - 1 - i1, m-1>::PerformSolve(A, val);
    TinyMatrixLoop<i1-1>::SolveUpper(A);
  }
  
 
  template<int i> template<class T, int m>
  inline void TinyMatrixLoop<i>::
  PermuteColumn(TinyMatrix<T, General, m, m>& A, const TinyVector<int, m>& pivot)
  {
    T val;
    if (pivot(i-1) != i-1)
      TinyMatrixDoubleLoop<i-1, m>::SwapColumn(A, pivot(i-1), val);
    
    TinyMatrixLoop<i-1>::PermuteColumn(A, pivot);
  }
 
  
  template<int j> template<class T, int m>
  inline void TinyMatrixLoop<j>::GetCholesky(TinyMatrix<T, Symmetric, m, m>& A)
  {
    T val = A(m-j, m-j);
    TinyMatrixDoubleLoop<m-j, m-j>::GetDiagonalCholesky(A, val);
 
    val = sqrt(val);
    A(m-j, m-j) = val;
    T one; SetComplexOne(one);
    T invVal = one / val, vloc;
    
    TinyMatrixDoubleLoop<m-j, j-1>::ModifyUpperCholesky(A, invVal, vloc);
    
    TinyMatrixLoop<j-1>::GetCholesky(A);
  }
 
 
  template<int i> template<class T, class T2, int m>
  inline void TinyMatrixLoop<i>::SolveCholesky(const class_SeldonNoTrans& trans,
                                               const TinyMatrix<T, Symmetric, m, m>& A,
                                               TinyVector<T2, m>& x)
  {
    x(m-i) /= A(m-i, m-i);
    TinyMatrixDoubleLoop<m-i, i-1>::SolveCholesky(trans, A, x);
    TinyMatrixLoop<i-1>::SolveCholesky(trans, A, x);
  }
  
  
  template<int i> template<class T, class T2, int m>
  inline void TinyMatrixLoop<i>::SolveCholesky(const class_SeldonTrans& trans,
                                               const TinyMatrix<T, Symmetric, m, m>& A,
                                               TinyVector<T2, m>& x)
  {
    T2 val = x(i-1);
    TinyMatrixDoubleLoop<i-1, m-i>::SolveCholesky(trans, A, x, val);
    
    x(i-1) = val / A(i-1, i-1);
    TinyMatrixLoop<i-1>::SolveCholesky(trans, A, x);
  }
 
 
  template<int i> template<class T, class T2, int m>
  inline void TinyMatrixLoop<i>::MltCholesky(const class_SeldonNoTrans& trans,
                                             const TinyMatrix<T, Symmetric, m, m>& A,
                                             TinyVector<T2, m>& x)
  {
    TinyMatrixDoubleLoop<i-1, m-i>::MltCholesky(trans, A, x);
    
    x(i-1) *= A(i-1, i-1);
    TinyMatrixLoop<i-1>::MltCholesky(trans, A, x);
  }
 
 
  template<int i> template<class T, class T2, int m>
  inline void TinyMatrixLoop<i>::MltCholesky(const class_SeldonTrans& trans,
                                             const TinyMatrix<T, Symmetric, m, m>& A,
                                             TinyVector<T2, m>& x)
  {
    T val = x(m-i)*A(m-i, m-i);
    TinyMatrixDoubleLoop<m-i, i-1>::MltCholesky(trans, A, x, val);
    x(m-i) = val;
    TinyMatrixLoop<i-1>::MltCholesky(trans, A, x);
  }
  
  
  /*************************
   * TinyMatrix<Symmetric> *
   *************************/
  
  
  template<class T, int m>
  inline TinyMatrix<T, Symmetric, m, m>::TinyMatrix()
  {
    this->Zero();
  }
  
  
  template<class T, int m>
  inline TinyMatrix<T, Symmetric, m, m>::TinyMatrix(int i)
  {
    this->Zero();
  }
 
 
  template<class T, int m>
  inline TinyMatrix<T, Symmetric, m, m>::TinyMatrix(int i, int j)
  {
    this->Zero();
  }
  
 
  template<class T, int m> template<class E>
  inline TinyMatrix<T, Symmetric, m, m>
  ::TinyMatrix(const TinyMatrixExpression<T, m, m, E>& A)
  {
    TinyMatrixLoop<m*(m+1)/2>::Copy(A, *this);
  }
  
  
  template<class T, int m>
  inline int TinyMatrix<T, Symmetric, m, m>::GetM()
  {
    return m;
  }
  
  
  template<class T, int m>
  inline int TinyMatrix<T, Symmetric, m, m>::GetN()
  {
    return m;
  }
  
  
  template<class T, int m>
  inline int TinyMatrix<T, Symmetric, m, m>::GetSize()
  {
    return m*(m+1)/2;
  }
 
  
  template<class T, int m>
  inline T* TinyMatrix<T, Symmetric, m, m>::GetData()
  {
    return data_;
  }  
  
 
  template<class T, int m> template<class T1>
  inline TinyMatrix<T, Symmetric, m, m> & TinyMatrix<T, Symmetric, m, m>::
  operator *=(const T1& a )
  {
    TinyMatrixLoop<m*(m+1)/2>::MltScal(a, *this);
    
    return *this;
  }
  
  
  template<class T, int m> template<class T1, class E>
  inline TinyMatrix<T, Symmetric, m, m> & TinyMatrix<T, Symmetric, m, m>::
  operator +=(const TinyMatrixExpression<T1, m, m, E>& B )
  {
    TinyMatrixLoop<m*(m+1)/2>::AddCopy(B, *this);
    return *this;
  }
  
  
  template<class T, int m> template<class T1, class E>
  inline TinyMatrix<T, Symmetric, m, m> & TinyMatrix<T, Symmetric, m, m>::
  operator -=(const TinyMatrixExpression<T1, m, m, E>& B )
  {
    TinyMatrixLoop<m*(m+1)/2>::DiffCopy(B, *this);
    return *this;
  }
  
  
  template<class T, int m>
  inline TinyMatrix<T, Symmetric, m, m> & 
  TinyMatrix<T, Symmetric, m, m>::operator =(const T& x)
  {
    this->Fill(x);
    return *this;
  }
  
  
  template<class T, int m> template<class T0, class E>
  inline TinyMatrix<T, Symmetric, m, m> & TinyMatrix<T, Symmetric, m, m>::
  operator =(const TinyMatrixExpression<T0, m, m, E>& x)
  {
    TinyMatrixLoop<m*(m+1)/2>::Copy(x, *this);
    return *this;
  }
  
  
  template<class T, int m>
  inline T& TinyMatrix<T, Symmetric, m, m>::operator()(int i, int j)
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, m, m, "TinyMatrix");
#endif
    
    return this->data_[j > i ? i*m - (i*(i+1))/2 + j: j*m - (j*(j+1))/2 + i];
  }
  
  
  template<class T, int m>
  inline const T& TinyMatrix<T, Symmetric, m, m>::operator()(int i, int j) const
  {
 
#ifdef SELDON_CHECK_BOUNDS
    CheckBounds(i, j, m, m, "TinyMatrix");
#endif
    
    return this->data_[j > i ? i*m - (i*(i+1))/2 + j: j*m - (j*(j+1))/2 + i ];
  }
 
 
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::Zero()
  {
    TinyMatrixLoop<m*(m+1)/2>::Zero(*this);
  }
  
  
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::SetIdentity()
  {
    return TinyMatrixLoop<m*(m+1)/2>::SetIdentity(*this);
  }
 
 
  template<class T, int m> template<class T0>
  inline void TinyMatrix<T, Symmetric, m, m>::SetDiagonal(const T0& alpha)
  {
    return TinyMatrixLoop<m*(m+1)/2>::SetDiagonal(*this, alpha);
  }
  
  
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::Fill()
  {
    return TinyMatrixLoop<m*(m+1)/2>::Fill(*this);
  }
  
  
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::FillRand()
  {
    return TinyMatrixLoop<m*(m+1)/2>::FillRand(*this);
  }
  
  
  template<class T, int m> template<class T0>
  inline void TinyMatrix<T, Symmetric, m, m>::Fill(const T0& a)
  {
    return TinyMatrixLoop<m*(m+1)/2>::Fill(*this, a);
  }
  
 
  template<class T, int m>
  inline bool TinyMatrix<T, Symmetric, m, m>::IsZero() const
  {
    return TinyMatrixLoop<m*(m+1)/2>::IsZero(*this);
  }
  
 
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::WriteText(const string& file_name) const
  {
    ofstream out(file_name.data());
    out.precision(cout.precision());
    out << *this;
    out.close();
  }
 
 
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::Write(const string& file_name) const
  {
    ofstream out(file_name.data());
    Write(out);
    out.close();
  }
  
  
  template<class T, int m>
  inline void TinyMatrix<T, Symmetric, m, m>::Write(ostream& out) const
  {
    int itmp = m;
    out.write(reinterpret_cast<char*>(&itmp), sizeof(int));    
    out.write(reinterpret_cast<char*>(&itmp), sizeof(int));    
    
    return TinyMatrixLoop<m>::Write(out, *this);    
  }
 
 
  /*************
   * Operators *
   *************/
 
  
  template<class T, int m, int n, class E1, class E2>
  inline bool operator==(const TinyMatrixExpression<T, m, n, E1> & u,
                         const TinyMatrixExpression<T, m, n, E2> & v)
  {
    return TinyMatrixLoop<m>::IsEqual(u, v);
  }
  
  
  template<class T, int m, int n, class E1, class E2>
  inline bool operator!=(const TinyMatrixExpression<T, m, n, E1> & u,
                         const TinyMatrixExpression<T, m, n, E2> & v)
  {
    return !TinyMatrixLoop<m>::IsEqual(u, v);
  }
  
 
  template <class T, int m, int n, class E>
  inline ostream& operator <<(ostream& out, const TinyMatrixExpression<T, m, n, E> & A)
  {
    TinyMatrixLoop<m>::WriteText(out, A);
    return out;
  }
  
  
  /***************************
   * Matrix-vector functions *
   ***************************/
 
  
  template<class T, int m, int n, class E, class T1, class E1> inline TinyVector<T1, m>
  dot(const TinyMatrixExpression<T, m, n, E>& A,
      const TinyVectorExpression<T1, n, E1>& x)
  {
    TinyVector<T1, m> b;
    TinyMatrixLoop<m>::Mlt(A, x, b);
    
    return b;
  }
  
  
  template<class T, int m, int n, class E, class T1, int k, class E1>
  inline TinyMatrix<T, General, m, k>
  dot(const TinyMatrixExpression<T, m, n, E>& A, 
      const TinyMatrixExpression<T1, n, k, E1> & B)
  {
    TinyMatrix<T, General, m, k> C;
    TinyMatrixLoop<m*k>::Mlt(A, B, C);
    
    return C;
  }
 
  
  template<class T0, class E0, class T1, class E1, class T2, int m, int n>
  inline void Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
                  const TinyVectorExpression<T1, n, E1>& x, TinyVector<T2, m>& y)
  {
    TinyMatrixLoop<m>::Mlt(A, x, y);
  }
  
  
  template<class T1, class E1, class T2, class E2, class T3, int m, int n>
  inline void MltAdd(const TinyMatrixExpression<T1, m, n, E1>& A,
                     const TinyVectorExpression<T2, n, E2>& x, TinyVector<T3, m>& y)
  {
    TinyMatrixLoop<m>::MltAdd(A, x, y);
  }
 
  
  template<class T0, class T1, class E1, class T2, class E2, class T3, int m, int n>
  inline void MltAdd(const T0& alpha, const TinyMatrixExpression<T1, m, n, E1>& A,
                     const TinyVectorExpression<T2, n, E2>& x, TinyVector<T3, m>& y)
  {
    TinyMatrixLoop<m>::MltAdd(alpha, A, x, y);
  }
   
  
  template<class T0, class E0, class T1, class E1, class T2, int m, int n>
  inline void MltTrans(const TinyMatrixExpression<T0, m, n, E0>& A,
                       const TinyVectorExpression<T1, m, E1>& x, TinyVector<T2, n>& y)
  {
    TinyMatrixLoop<n>::Mlt(transpose(A), x, y);
  }
 
 
  template<class T0, class E0, class T1, class E1, class T2, int m, int n>
  inline void Mlt(const class_SeldonTrans&,
                  const TinyMatrixExpression<T0, m, n, E0>& A,
                  const TinyVectorExpression<T1, m, E1>& x, TinyVector<T2, n>& y)
  {
    TinyMatrixLoop<n>::Mlt(transpose(A), x, y);
  }
  
  
  template<int m, int n, class T0, class T1, class E1, class T2, class E2, class T3>
  inline void Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
                          const TinyVectorExpression<T2, n, E2>& y,
                          TinyMatrix<T3, General, m, n>& A)
  {
    TinyMatrixLoop<m>::Rank1Update(alpha, x, y, A);
  }
 
 
  template<int m, int n, class T1, class E1, class T2, class E2, class T3>
  inline void Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
                          const TinyVectorExpression<T2, n, E2>& y,
                          TinyMatrix<T3, General, m, n>& A)
  {
    TinyMatrixLoop<m>::Rank1Matrix(x, y, A);
  }
  
 
  template<int m, class T0, class T1, class E1, class T3>
  inline void Rank1Update(const T0& alpha, const TinyVectorExpression<T1, m, E1>& x,
                          TinyMatrix<T3, Symmetric, m, m>& A)
  {
    TinyMatrixLoop<m>::Rank1Update(alpha, x, A);
  }
 
 
  template<int m, class T1, class E1, class T3>
  inline void Rank1Matrix(const TinyVectorExpression<T1, m, E1>& x,
                          TinyMatrix<T3, Symmetric, m, m>& A)
  {
    TinyMatrixLoop<m>::Rank1Matrix(x, A);
  }
 
  
  template<int m, int n, class T1, class E1>
  inline void GetCol(const TinyMatrixExpression<T1, m, n, E1>& A, int k, TinyVector<T1, m>& x)
  {
    TinyMatrixLoop<m>::GetCol(A, k, x);
  }
  
  
  template<int m, int n, class T0, class E0, class T1>
  inline void GetRow(const TinyMatrixExpression<T0, m, n, E0>& A, int k, TinyVector<T1, n>& x)
  {
    TinyMatrixLoop<n>::GetRow(A, k, x);
  }
  
  
  template<int m, int n, class T1, class E1, class Prop>
  inline void SetCol(const TinyVectorExpression<T1, m, E1>& x, int k, TinyMatrix<T1, Prop, m, n>& A)
  {
    TinyMatrixLoop<m>::SetCol(x, k, A);
  }
  
  
  template<int m, int n, class T1, class E1, class Prop>
  inline void SetRow(const TinyVectorExpression<T1, n, E1>& x, int k, TinyMatrix<T1, Prop, m, n>& A)
  {
    TinyMatrixLoop<n>::SetRow(x, k, A);
  }
  
  
  /****************************
   * Matrix-Matrix operations *
   ****************************/
  
  
  template<class T, int m, int n, class E, class Prop>
  inline void Copy(const TinyMatrixExpression<T, m, n, E>& A, TinyMatrix<T, Prop, m, n>& B)
  {
    B = A;
  }
  
  
  template<class T0, class E0, class T1, class E1,
           class T2, class Prop2, int m, int n>
  inline void Add(const TinyMatrixExpression<T0, m, n, E0>& A,
                  const TinyMatrixExpression<T1, m, n, E1>& B,
                  TinyMatrix<T2, Prop2, m, n>& C)
  {
    C = A + B;
  }
  
  
  template<class T1, class E1, class T2,
           class Prop2, int m, int n>
  inline void Add(const T1& alpha, const TinyMatrixExpression<T1, m, n, E1>& A,
                  TinyMatrix<T2, Prop2, m, n>& B)
  {
    B += alpha*A;
  }
  
  
  template<class T0, class E0, class T1, class Prop1, int m, int n>
  inline void Add(const TinyMatrixExpression<T0, m, n, E0>& A, TinyMatrix<T1, Prop1, m, n>& B)
  {
    B += A;
  }
  
  
  template<class T0, class T1, class Prop, int m, int n>
  inline void Mlt(const T0& alpha, TinyMatrix<T1, Prop, m, n>& A)
  {
    A *= alpha;
  }
  
    
  template<class T0, class E0, class T1, class E1,
           class T2, class Prop2, int m, int n, int k>
  inline void Mlt(const TinyMatrixExpression<T0, m, n, E0>& A,
                  const TinyMatrixExpression<T1, n, k, E1>& B,
                  TinyMatrix<T2, Prop2, m, k>& C)
  {
    TinyMatrixLoop<TinyMatrix<T2, Prop2, m, k>::size_ >::Mlt(A, B, C);
  }
  
  
  template<class T0, class E0, class T1, class E1,
           class T2, class Prop2, int m, int n, int k>
  inline void MltTrans(const TinyMatrixExpression<T0, m, n, E0>& A,
                       const TinyMatrixExpression<T1, k, n, E1>& B,
                       TinyMatrix<T2, Prop2, m, k>& C)
  {
    TinyMatrixLoop<TinyMatrix<T2, Prop2, m, k>::size_ >::Mlt(A, transpose(B), C);
  }
 
 
  template<class T3, class T0, class E0, class T1, class E1,
           class T2, class Prop2, class T4, int m, int n, int k>
  inline void MltAdd(const T3& alpha, const class_SeldonTrans&,
                     const TinyMatrixExpression<T0, n, m, E0>& A,
                     const class_SeldonNoTrans&, const TinyMatrixExpression<T1, n, k, E1>& B,
                     const T4& beta, TinyMatrix<T2, Prop2, m, k>& C)
  {
    TinyMatrix<T2, Prop2, m, k> D;
    Mlt(transpose(A), B, D);
    C = beta*C + alpha*D;
  }
  
  
  template<class T, int m, int n, class E>
  inline void Transpose(const TinyMatrixExpression<T, m, n, E> & A,
                        TinyMatrix<T, General, n, m> & B)
  {
    B = transpose(A);
  }
  
  
  template<class T, int m>
  inline void Transpose(TinyMatrix<T, General, m, m>& B)
  {
    TinyMatrix<T, General, m, m> A(B);
    B = transpose(A);
  }
 
  
  template<class T, class Prop, int m, int n> inline
  typename ClassComplexType<T>::Treal MaxAbs(const TinyMatrix<T, Prop, m, n>& A)
  {
    typename ClassComplexType<T>::Treal amax(0);
    TinyMatrixLoop<TinyMatrix<T, Prop, m, n>::size_>::GetMaxAbs(A, amax);
    return amax;
  }
  
  
  /****************
   * 1x1 matrices *
   ****************/
  
  
  template<class T, class Prop>
  inline T Det(const TinyMatrix<T, Prop, 1, 1> & A)
  {
    return A(0, 0);
  }
  
  
  template<class T>
  inline void GetInverse(const TinyMatrix<T, General, 1, 1> & A,
                         TinyMatrix<T, General, 1, 1> & B)
  {
    T one; SetComplexOne(one);
    B(0,0) = one/A(0, 0);
  }
  
  
  template<class T>
  inline void GetInverse(const TinyMatrix<T, Symmetric, 1, 1> & A,
                         TinyMatrix<T, Symmetric, 1, 1> & B)
  {
    T one; SetComplexOne(one);
    B(0,0) = one/A(0, 0);
  }
 
  
  template<class T>
  inline void GetInverse(TinyMatrix<T, General, 1, 1> & B)
  {
    T one; SetComplexOne(one);
    B(0,0) = one/B(0, 0);
  }
 
  
  template<class T>
  inline void GetInverse(TinyMatrix<T, Symmetric, 1, 1> & B)
  {
    T one; SetComplexOne(one);
    B(0,0) = one/B(0, 0);
  }
 
 
  template<class T0, class T1>
  inline void GetNormalProjector(const TinyVector<T0, 1>& n, TinyMatrix<T1, Symmetric, 1, 1>& P)    
  {
    P(0,0) = n(0)*n(0);
  }
 
  
  template<class T0, class T1>
  inline void GetNormalProjector(const T0& n, TinyMatrix<T1, Symmetric, 1, 1>& P)    
  {
    P(0,0) = n*n;
  }
 
  
  /****************
   * 2x2 matrices *
   ****************/
  
  
  template<class T, class Prop>
  inline T Det(const TinyMatrix<T, Prop, 2, 2> & A)
  {
    return A(0,0)*A(1,1) - A(0,1)*A(1,0);
  }
 
  
  template<class T0, class T1>
  inline void GetTangentialProjector(const TinyVector<T0, 2>& n,
                                     TinyMatrix<T1, Symmetric, 2, 2>& P)
  {
    P(0, 0) = n(1)*n(1);
    P(0, 1) = -n(0)*n(1);
    P(1, 1) = n(0)*n(0);
  }
  
  
  template<class T0, class T1>
  inline void GetNormalProjector(const TinyVector<T0, 2>& n, TinyMatrix<T1, Symmetric, 2, 2>& P)
  {
    P(0, 0) = n(0)*n(0);
    P(0, 1) = n(0)*n(1);
    P(1, 1) = n(1)*n(1);
  }
 
 
  template<class T0, class T1>
  inline void GetNormalProjector(const TinyVector<T0, 2>& n, TinyMatrix<T1, General, 2, 2>& P)
  {
    P(0, 0) = n(0)*n(0);
    P(0, 1) = n(0)*n(1);
    P(1, 1) = n(1)*n(1);
    P(1, 0) = P(0, 1);
  }
  
  
  /**************************
   * Functions for any size *
   **************************/
  
  
  template<class T0, class T1, int m>
  inline void GetTangentialProjector(const TinyVector<T0, m>& n,
                                     TinyMatrix<T1, Symmetric, m, m>& P)
  {
    T0 one; SetComplexOne(one);
    P.SetIdentity();
    Rank1Update(-one, n, P);
  }
  
  
  template<class T0, class T1, int m>
  inline void GetNormalProjector(const TinyVector<T0, m>& n,
                                 TinyMatrix<T1, Symmetric, m, m>& P)
  {
    Rank1Matrix(n, P);
  }
 
    
  template<class T, class Prop, int m, int n>
  inline void FillZero(TinyMatrix<T, Prop, m, n>& X)
  {
    X.Zero();
  }
  
  
  template<class T, int m, int n>
  inline void DotProdCol(const TinyMatrix<T, General, m, n> &A,
                         const TinyMatrix<T, General, m, n>& B, TinyVector<T, n> &res)
  {    
    TinyVector<T, m> vectA, vectB;
    for (int p = 0; p < n ; p ++)
      {
        GetCol(A, p, vectA);
        GetCol(B, p, vectB);
        res(p) += DotProd(vectA, vectB);
      }
  }
  
} // end namespace
 
#define MONTJOIE_FILE_TINY_MATRIX_INLINE_CXX
#endif