#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>

typedef boost::numeric::ublas::mapped_vector_of_mapped_vector<double> SPARSE_MATRIX;
typedef boost::numeric::ublas::mapped_vector_of_mapped_vector<double, boost::numeric::ublas::column_major> SPARSE_MATRIX2;
typedef boost::numeric::ublas::vector<double> DENSE_VECTOR;

template <class V>
void print(const boost::numeric::ublas::vector_expression<V>& v) {
  std::cout << "size: " << v().size() << std::endl;
  typename V::const_iterator i = v().begin();
  for (; i != v().end(); i++) {
	std::cout << i.index() << ":" << (*i) << "  ";
  }
  std::cout << std::endl;
}

template <class V>
void print(const boost::numeric::ublas::matrix_expression<V>& v) {
  std::cout << "size: " << v().size1() << "x" << v().size2() << std::endl;
  typename V::const_iterator1 i = v().begin1();
  for (; i != v().end1(); ++i) {
    std::cout << " row " << i.index1() << ": ";
	typename V::const_iterator2 j = i.begin();
	for (; j != i.end(); ++j) {
	  std::cout << j.index1() << "," << j.index2() << ":" << (*j) << "  ";
	}
  }
  std::cout << std::endl;
}

template <class V>
void print2(boost::numeric::ublas::matrix_expression<V>& v) {
  std::cout << "size: " << v().size1() << "x" << v().size2() << std::endl;
  typename V::iterator1 i = v().begin1();
  for (; i != v().end1(); ++i) {
	typename V::iterator2 j = i.begin();
	for (; j != i.end(); ++j) {
	  std::cout << j.index1() << "," << j.index2() << ":" << (*j) << "  ";
	}
  }
  std::cout << std::endl;
}

int main(int /* argc */ , char * /* argv */[])
{

  // check row major

  SPARSE_MATRIX A(10,10);
  DENSE_VECTOR  x(10);
  DENSE_VECTOR  y(10);

  std::fill(x.begin(), x.end(), 1.0);
  std::fill(y.begin(), y.end(), 0.0);
  A.insert_element(7,8,1.0);
  A.insert_element(7,9,2.0);
  A.insert_element(9,8,3.0);
  A.insert_element(9,9,4.0);

  std::cout << "matrix A  " << A << "\n";
  std::cout << "vector x  " << x << "\n";

  std::cout << "prod(A,x) " << prod(A,x) << "\n";

  std::cout << "matrix A   "; print( A );
  std::cout << "matrix A^T ";print( trans(A) );
  std::cout << "vector Ax  ";print( prod(A,x) );

  y = prod(A,x);
  std::cout << "prod(A,x) " << y << "\n";

  // check column major
  SPARSE_MATRIX2 B(10,10);
  B.insert_element(7,8,1.0);
  B.insert_element(7,9,2.0);
  B.insert_element(9,8,3.0);
  B.insert_element(9,9,4.0);

  std::cout << "matrix B   "; print( B );
  std::cout << "matrix B^T "; print( trans(B) );
  std::cout << "vector Bx  "; print( prod(B,x) );

  y = prod(B,x);
  std::cout << "prod(B,x) " << y << "\n";

  // check mutable iterators
  print2( A );
  print2( B );

  return EXIT_SUCCESS;
}