#ifndef _utl_SMatrix_h_
#define _utl_SMatrix_h_

#include <boost/lambda/lambda.hpp>

#include <utl/ListAssignmentProxy.h>


namespace utl {

  /**
    \class SMatrix SMatrix.h "utl/SMatrix.h"
    \brief Static (small and dense) matrix class

    Power of unrolled loops.

    \author Darko Veberic
    \version $Id$
    \date 11 Jun 2008
    \ingroup math
  */

  template<unsigned int n, unsigned int m, typename T = double>
  class SMatrix {
  public:
    typedef T Type;

    static const unsigned int NRows = n;
    static const unsigned int NColumns = m;
    static const unsigned int Size = n * m;

    static unsigned int GetNRows() { return n; }

    static unsigned int GetNColumns() { return m; }

    static unsigned int GetSize() { return Size; }

    SMatrix() { }

    explicit SMatrix(const T& init) { Clear(init); }

    template<class Matrix>
    explicit SMatrix(const Matrix& a) { Assign(a); }

    template<typename Matrix>
    void
    Assign(const Matrix& a)
    {
      for (unsigned int i = 0; i < n; ++i)
        for (unsigned int j = 0; j < m; ++j)
          fMatrix[i][j] = a[i][j];
    }

    ListMatrixAssignmentProxy<SMatrix, m>
    operator=(const T& e)
    {
      fMatrix[0][0] = e;
      return ListMatrixAssignmentProxy<SMatrix, m>(*this);
    }

    template<typename Matrix>
    SMatrix& operator=(const Matrix& a) { Assign(a); return *this; }

    void
    Clear(const T& init = T())
    {
      T* a = &fMatrix[0][0];
      *a = init;
      for (unsigned int i = 1; i < Size; ++i)
        *(++a) = init;
    }

    T* operator[](const unsigned int i) { return fMatrix[i]; }

    const T* operator[](const unsigned int i) const { return fMatrix[i]; }

    template<unsigned int i, unsigned int j>
    T& Get() { return fMatrix[i][j]; }

    template<unsigned int i, unsigned int j>
    const T& Get() const { return fMatrix[i][j]; }

    template<typename U>
    SMatrix&
    operator*=(const U& fact)
    {
      T* a = &fMatrix[0][0];
      *a *= fact;
      for (unsigned int i = 1; i < Size; ++i)
        *(++a) *= fact;
      return *this;
    }

    template<typename U>
    SMatrix&
    operator/=(const U& div)
    {
      T* a = &fMatrix[0][0];
      *a /= div;
      for (unsigned int i = 1; i < Size; ++i)
        *(++a) /= div;
      return *this;
    }

    template<typename U>
    SMatrix&
    operator+=(const SMatrix<n, m, U>& mat)
    {
      T* a = &fMatrix[0][0];
      const T* b = &mat.fMatrix[0][0];
      *a += *b;
      for (unsigned int i = 1; i < Size; ++i)
        *(++a) += *(++b);
      return *this;
    }

    template<typename U>
    SMatrix&
    operator-=(const SMatrix<n, m, U>& mat)
    {
      T* a = &fMatrix[0][0];
      const T* b = &mat.fMatrix[0][0];
      *a -= *b;
      for (unsigned int i = 1; i < Size; ++i)
        *(++a) -= *(++b);
      return *this;
    }

    template<unsigned int nn, unsigned int mm, typename U>
    bool
    operator==(const SMatrix<nn, mm, U>& a)
      const
    {
      if (n != nn || m != mm)
        return false;
      for (unsigned int i = 0; i < n; ++i)
        for (unsigned int j = 0; j < m; ++j)
          if (fMatrix[i][j] != a[i][j])
            return false;
      return true;
    }

    template<unsigned int nn, unsigned int mm, typename U>
    bool operator!=(const SMatrix<nn, mm, U>& a) const
    { return !operator==(a); }

  private:
    T fMatrix[n][m];
  };

  template<unsigned int n, unsigned int m, typename T>
  const unsigned int SMatrix<n, m, T>::NRows;

  template<unsigned int n, unsigned int m, typename T>
  const unsigned int SMatrix<n, m, T>::NColumns;

  template<unsigned int n, unsigned int m, typename T>
  const unsigned int SMatrix<n, m, T>::Size;


  template<unsigned int n, unsigned int m, typename T, typename U>
  inline
  typename boost::lambda::return_type_2<
    boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
    SMatrix<n, m, T>,
    SMatrix<m, n, U>
  >::type
  operator*(const SMatrix<n, m, T>& a, const SMatrix<m, n, U>& b)
  {
    typedef typename boost::lambda::return_type_2<
      boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
      SMatrix<n, m, T>,
      SMatrix<m, n, U>
    >::type ProductMatrixType;

    ProductMatrixType r;

    for (unsigned int i = 0; i < n; ++i)
      for (unsigned int j = 0; j < n; ++j) {
        r[i][j] = a[i][0] * b[0][j];
        for (unsigned int k = 1; k < m; ++k)
          r[i][j] += a[i][k] * b[k][j];
      }

    return r;
  }

}


namespace boost {

  namespace lambda {

    /// specialization for custom class SMatrix<>
    template<unsigned int n, unsigned int m, typename T, typename U>
    class plain_return_type_2<
      arithmetic_action<multiply_action>,
      utl::SMatrix<n, m, T>,
      utl::SMatrix<m, n, U>
    > {
    private:
      // find type of T * U
      typedef typename
        return_type_2<
          arithmetic_action<multiply_action>,
          T,
          U
        >::type res_type;

    public:
      typedef utl::SMatrix<n, n, res_type> type;
    };

  }

}


#include <utl/SMatrixSVectorOp.h>


#endif