#include <iostream> 
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/operation.hpp>

using namespace boost::numeric::ublas;

  typedef float index_type;
  template<class T, class IT>
  class si_allocator : public std::allocator<T>
  {
  public:
    typedef IT size_type;
    typedef IT difference_type;
  };
  
  typedef basic_row_major<index_type, index_type> si_row_major;

  typedef coordinate_matrix<
  	float, si_row_major, 0,
  	unbounded_array<index_type>,
  	unbounded_array<float> > si_coordinate_matrix;


int main(int argc, char** argv)
{
  index_type nrows = 3;
  index_type ncols = 4;
  si_coordinate_matrix m(nrows, ncols);
  vector<float> v(ncols), r1(nrows);
  
  for (float pctNZ = 0.0; pctNZ < 1.01; pctNZ += 0.1) {

    // Initialize matrices
    // Values don't matter except for % NZ
    m.clear ();
    for (index_type i = 0; i < nrows; ++i)
      for (index_type j = 0; j < ncols; ++j)
         if (float(rand()%100)/100.0 < pctNZ)
         	m(i,j) = 1;

    // Initialize RHS std::vectors
    // Values don't matter
    for (index_type i = 0; i < ncols; ++i)
      v[i] = i;
   }

   m= prod(m,m);
}