/**
   \file
   Static Matrix and Vector operations

   \author Darko Veberic
   \version $Id$
   \date 12 Jun 2008
*/

#ifndef _utl_SMatrixSVector_h_
#define _utl_SMatrixSVector_h_

#include <boost/lambda/lambda.hpp>

#include <utl/SVector.h>
#include <utl/SMatrix.h>


namespace boost {

  namespace lambda {

    /// specialization for SMatrix<> * SVector<>
    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::SVector<m, U>
    > {
    private:
      // delegate to the type of T * U
      typedef typename
        return_type_2<
          arithmetic_action<multiply_action>,
          T,
          U
        >::type res_type;

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

  }

}


namespace utl {

  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>,
    utl::SMatrix<n, m, T>,
    utl::SVector<m, U>
  >::type
  operator*(const SMatrix<n, m, T>& a, const SVector<m, U>& x)
  {
    typedef typename boost::lambda::return_type_2<
      boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
      utl::SMatrix<n, m, T>,
      utl::SVector<m, U>
    >::type VectorResultType;

    VectorResultType y;

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

    return y;
  }


  template<class Stream, unsigned int n, unsigned int m, typename T>
  Stream&
  operator<<(Stream& s, const SMatrix<n, m, T>& a)
  {
    const SVector<n, SVector<m, T> > v(a);

    return s << v;
  }

}


#endif