#ifndef _utl_SparseMatrixVectorOp_h_
#define _utl_SparseMatrixVectorOp_h_

/**
   \file
   Sparse Matrix-Vector operations

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

#include <numeric>

#include <boost/lambda/lambda.hpp>

#include <utl/SparseVector.h>
#include <utl/SparseMatrix.h>


namespace boost {

  namespace lambda {

    /// specialization for SMatrix<> * SVector<>
    template<typename T, typename U>
    class plain_return_type_2<
      arithmetic_action<multiply_action>,
      utl::SparseMatrix<T>,
      utl::SparseVector<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::SparseVector<res_type> type;
    };

  }

}


namespace utl {


  template<typename T, typename U>
  typename boost::lambda::return_type_2<
    boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
    SparseMatrix<T>,
    SparseVector<U>
  >::type
  operator*(const SparseMatrix<T>& m, const SparseVector<U>& v)
  {
    typedef typename boost::lambda::return_type_2<
      boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
      SparseMatrix<T>,
      SparseVector<U>
    >::type ResultType;

    ResultType result;

    for (typename SparseMatrix<T>::ConstIterator mIt = m.SparseBegin();
         mIt != m.SparseEnd(); ++mIt) {
      const unsigned int j = mIt.GetColumn();
      typename SparseVector<U>::ConstIterator aj = v.Find(j);
      if (aj != v.SparseEnd())
        result(mIt.GetRow()) += mIt() * aj();
    }

    return result;
  }


}


#endif