// -*- c++ -*- //
/** \file main.cpp  \brief main routines */
/***************************************************************************
 -   begin                : 2003-09-05
 -   copyright            : (C) 2003 by Gunter Winkler
 -   email                : guwi17@gmx.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include <iostream>
#include <algorithm>

#include <boost/timer.hpp>

#include <boost/random.hpp>
#include <boost/random/mersenne_twister.hpp>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/vector_sparse.hpp>
#include <boost/numeric/ublas/vector_of_vector.hpp>
#include <boost/numeric/ublas/vector_expression.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/numeric/ublas/lu.hpp>


#ifndef NDEBUG
const size_t default_size = 100;
#else
const size_t default_size = 1000;
#endif

const size_t subsize = 10;

using namespace boost::numeric::ublas ;

using boost::mt19937 ;

using std::cout;
using std::endl;

template <class MAT>
void test_matrix_fill_plus_assign(MAT& A)
{
  typedef typename MAT::size_type  size_type;

  boost::numeric::ublas::vector<size_type>    index(subsize);

  mt19937   mtrand;

  size_type  size = A.size1();
  size_type count = size * 5;

  for (size_type n=0; n<count; ++n)	{
	for (size_type i=0; i<subsize; ++i) {
	  size_type  r = mtrand();
	  index(i) = size_type( (r * double(size)) / double(mtrand.max()) );
	}
	for (size_type i=0; i<subsize; ++i) {
	  for (size_type j=0; j<subsize; ++j) {
		A(index(i),index(j)) += (1.0);
	  }
	}
  }

  double sum = 0.0;
  for (size_type i=0; i<subsize; ++i)
    sum += A(i,0);
  cout << "sum(column(A,0)) " << sum << endl;
}

template <class MAT>
void test_matrix_fill_insert(MAT& A)
{
  typedef typename MAT::size_type  size_type;

  boost::numeric::ublas::vector<size_type>    index(subsize);

  mt19937   mtrand;

  size_type  size = A.size1();
  size_type count = size * 5;

  for (size_type n=0; n<count; ++n)	{
	for (size_type i=0; i<subsize; ++i) {
	  size_type  r = mtrand();
	  index(i) = size_type( (r * double(size)) / double(mtrand.max()) );
	}
	for (size_type i=0; i<subsize; ++i) {
	  for (size_type j=0; j<subsize; ++j) {
		A.insert_element(index(i),index(j),1.0);
	  }
	}
  }

  double sum = 0.0;
  for (size_type i=0; i<subsize; ++i)
    sum += A(i,0);
  cout << "sum(column(A,0)) " << sum << endl;
}

template <class CONTAINER>
class op_reserve
{
private: 
  size_t size;
public:
  op_reserve() : size(0) {}
  op_reserve(const size_t n) : size(n) {}
  void operator()(CONTAINER & c) const
  {
    c.reserve(size);
  }
};

template<class CONTAINER>
op_reserve<CONTAINER>
do_reserve(const size_t n, const CONTAINER& /* c */)
{
  return op_reserve<CONTAINER>(n);
}

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

  typedef size_t  size_type;

  size_type size = default_size;

  if (argc > 1)
	size = ::atoi (argv [1]);

  cout << "size " << size << endl;

  {
	coordinate_matrix<double, row_major> A(size,size,2*subsize*size);
	boost::timer t;
	t.restart();
	test_matrix_fill_insert(A);
	cout << "test_matrix_fill_insert coordinate_matrix: " << t.elapsed() << endl;
	compressed_matrix<double, row_major> B(size,size,A.nnz());
	t.restart();
    B.assign(A);
	cout << "assign from coordinate_matrix: " << t.elapsed() << endl;
  }

  {
	generalized_vector_of_vector< double,
      row_major, vector<coordinate_vector<double> > > A(size,size);
	boost::timer t;
	t.restart();
#ifdef RESERVE
    std::for_each( A.data().begin(), A.data().end(), do_reserve(2*subsize+1, A.data()[0]) );
#endif
	test_matrix_fill_insert(A);
	cout << "test_matrix_fill_insert vector_of_coordinate_vector: " << t.elapsed();
    cout << "  nnz: " << A.nnz() << " nnz capacity: " << A.nnz_capacity() << "\n";
	compressed_matrix<double, row_major> B(size,size,A.nnz());
	t.restart();
    B.assign(A);
	cout << "assign from vector_of_coordinate_vector: " << t.elapsed() << endl;
  }

#ifdef SLOW
  {
	compressed_matrix<double, row_major> A(size,size,2*subsize*size);
	boost::timer t;
	t.restart();
	test_matrix_fill_plus_assign(A);
	cout << "test_matrix_fill_plus_assign compressed_matrix: " << t.elapsed() << endl;
  }
#endif

  {
	generalized_vector_of_vector< double,
	  row_major, vector<compressed_vector<double> > > A(size,size);
	boost::timer t;
	t.restart();
#ifdef RESERVE
    std::for_each( A.data().begin(), A.data().end(), do_reserve(2*subsize+1, A.data()[0]) );
#endif
	test_matrix_fill_plus_assign(A);
	cout << "test_matrix_fill_plus_assign vector_of_compressed_vector: " << t.elapsed() ;
    cout << "  nnz: " << A.nnz() << " nnz capacity: " << A.nnz_capacity() << "\n";
	compressed_matrix<double, row_major> B(size,size,A.nnz());
	t.restart();
    B.assign(A);
	cout << "assign from vector_of_compressed_vector: " << t.elapsed() << endl;
  }

#ifdef SLOW
  {
	sparse_matrix<double, row_major> A(size,size);
	boost::timer t;
	t.restart();
	test_matrix_fill_plus_assign(A);
	cout << "test_matrix_fill_plus_assign sparse_matrix: " << t.elapsed() << endl;
  }
#endif

  return EXIT_SUCCESS;

}