#ifndef SPARSEMULT_H
#define SPARSEMULT_H

#include <mtl/matrix.h>
#include <mtl/utils.h>



// mult -- sparse-sparse matrix multiply
//      -- sets out <- out + A*B
//      -- assumes matrices row-major for efficient implementation
template<class SpMat>
void mult(const SpMat &A, const SpMat &B, SpMat &out)
{
  typedef typename mtl::compressed1D<typename SpMat::value_type> SparseVec;
  typename mtl::rows_type<SpMat>::type Ar, Br;
  Ar = rows(A);
  Br = rows(B);

  // row for output -- need to re-use
  SparseVec temp_row;

  temp_row.reserve(10);

  typename SpMat::iterator i;
  for ( i = Ar.begin(); i != Ar.end(); ++i )
    {
      int row_num = (*i).begin().row();
      temp_row.clear();
      typename SpMat::OneD::iterator j;
      for ( j = (*i).begin(); j != (*i).end(); ++j )
        {
          int col_num = j.column();
          double val = *j;
          typename SpMat::Row Brow = Br[col_num];
          typename SpMat::Row::iterator k;
          for ( k = Brow.begin(); k != Brow.end(); ++k )
            temp_row[k.column()] = temp_row[k.column()] + val*(*k);
        }
      typename SparseVec::iterator p;
      for ( p = temp_row.begin(); p != temp_row.end(); ++p )
        out(row_num,p.index()) += *p;
    }
}

#endif // SPARSEMULT_H