// -*- 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 <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;
const size_t default_matrix_size = 5000;
#else
const size_t default_size = 1000;
const size_t default_matrix_size = 50000;
#endif


using namespace boost::numeric::ublas ;

using boost::mt19937 ;

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

template <class MAT>
void fill_matrix(MAT& A, const size_t N)
{
  mt19937   mtrand(273645);

  size_t  size1 = A.size1();
  size_t  size2 = A.size2();

  for (size_t n=0; (n<10*N) && (A.nnz() < N); ++n)	{
	size_t  r1 = mtrand();
	size_t  i = size_t( (r1 * double(size1)) / double(mtrand.max()) );
	size_t  r2 = mtrand();
	size_t  j = size_t( (r2 * double(size2)) / double(mtrand.max()) );
	A(i,j) = 1.0;
  }

}

int main(int argc, char *argv[])
{
  typedef size_t  size_type;

  size_type size = default_size;
  size_type matrix_size = default_matrix_size;

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

  if (argc > 2) {
	matrix_size = ::atoi (argv [1]);
	size = ::atoi (argv [2]);
  }

  cout << matrix_size << "\t" << size << flush;

  boost::timer t;

  mapped_matrix<double>     mt(matrix_size, matrix_size);
  compressed_matrix<double> ms1(matrix_size, matrix_size),	ms2(matrix_size, matrix_size);

  fill_matrix(mt, size);
  fill_matrix(ms2, size);

  cout << "\t" << mt.nnz() ;

  t.restart();
  compressed_matrix<double> ms3(-ms1); // O(N^2) time
  cout << "\t" << t.elapsed() << flush;

  t.restart();
  ms1 = mt ; //O(N^2) time
  cout << "\t" << t.elapsed() << flush;

  t.restart();
  mt = prod(ms1, trans(ms2)); // seems O(N^2) ,I didn't have time to check it individually
  cout << "\t" << t.elapsed() << "\t" << mt.nnz() << flush;

  t.restart();
  ms3 = prod(ms1, trans(ms2)); // seems O(N^2) ,I didn't have time to check it individually
  cout << "\t" << t.elapsed() << "\t" << mt.nnz() << flush;

  cout << std::endl;

  return EXIT_SUCCESS;
}


/*

build:

g++ -I $HOME/include -o sparse_copy sparse_copy.cpp -g -DNDEBUG -O3

execute:

for ms in 100 1000 2000 4000 8000 16000; do 
  echo -n "$ms " ; 
  for s in 100200 400 800 1600 3200 6400 12800; do 
    echo -n "$s " ; 
    ./sparse_copy $ms $s >>sparse_copy.log ; 
  done ; 
  echo ; 
done

*/