//          Copyright Gunter Winkler 2004 - 2007.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

// authors: Gunter Winkler <guwi17 at gmx dot de>
//          Konstantin Kutzkow

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

#include <boost/random.hpp>

#include <boost/timer.hpp>

#include "scatter-gather.hpp"

#include <algorithm>
#include <iostream>

using namespace std;

template <class Vector>
void add_vectors(int size, int number, int interval){

  Vector v1(size);
  Vector v2(size);
  Vector tmp(size);
  Vector res(size);
    
  typedef std::vector< Vector > Matrix;
    
  Matrix M;
  int i = 0;
    
    
  boost::mt19937 rng(12345);          // produces randomness out of thin air
  boost::uniform_int<std::size_t> dist(0,size-1);
  boost::variate_generator<boost::mt19937&, boost::uniform_int<std::size_t> >
    gen(rng, dist);             // glues randomness with mapping

  for (int i = 0; i < 100; ++i){
    for (int j = 0; j < size; j += interval) {
      v2(gen()) += gen();
    }
    M.push_back(v2);
    v2.clear();
  }
    
    


    



  int cnt = 0;
    

  boost::timer t_scatter;
  working_row wr(size);

  cout << "Start scatter" << endl;

  t_scatter.restart();
  scatter(v1, wr);

  int k = 0;
    
  while (cnt < number){
    while (k < M.size()){
      sp_axpy(1, M[k], wr);
      ++k;
    }
    cnt ++;
  }
    
  gather(v1, wr);
  cout << "CPU time  = " << t_scatter.elapsed()/M[0].nnz() << endl;

  res = v1;
    
    
    
   

  cnt = 0;

  v1.clear();
    
  boost::timer t_dense;

  cout << "Start direct += " << endl;

  t_dense.restart();
  k = 0;

  while (cnt < number){
    while (k < M.size()){
      v1 += M[k];
      ++k;
    }
    cnt ++;
  }
    

  cout << "CPU time  = " << t_dense.elapsed()/M[0].nnz() << endl;
    
  cout << "Norm = " << boost::numeric::ublas::norm_1(v1 - res) << endl;

}

int main(int argc, char *argv[])
{
  using namespace std;

  int size = (int)4e4;
  int interval = 100;
  int number = (int)1e1;

  typedef boost::numeric::ublas::compressed_vector< double > comp_vector;
  typedef boost::numeric::ublas::coordinate_vector< double > coord_vector;
  typedef boost::numeric::ublas::mapped_vector< double > map_vector;
    
  cout << "Compressed vector" << endl;
  add_vectors< comp_vector >(size, number, interval);
    
  cout << "\nCoordinate vector" << endl;
  add_vectors< coord_vector >(size, number, interval);
    
  cout << "\nMapped vector" << endl;
  add_vectors< map_vector >(size, number, interval);
    
  return EXIT_SUCCESS;
}