#include <iostream>
#include <boost/timer.hpp>

// Compute y=Ax

using namespace std;

const unsigned M= 5000, N= M;

 

int main(int argc, char** argv) 
{
  double *a= new double[M*N], *x= new double[N], *y= new double[M];

  fill(&a[0], &a[M*N], 2.0);
  fill(&x[0], &x[N], 3.0);
  
  boost::timer t1;
  for(unsigned k= 0; k< 5; ++k)
  for (unsigned i= 0; i < M; ++i) {
    double tmp= 0.0;
    for (unsigned j= 0; j < N; ++j)
      tmp+= a[i*N+j] * x[j];
    y[i]= tmp;
  }
  cout << "Over rows first: time = " << t1.elapsed() << "\n";

  boost::timer t2;
  for(unsigned k= 0; k< 5; ++k)
  for (unsigned i= 0; i < M; ++i) 
    y[i]= 0.0;
  for(unsigned k= 0; k< 5; ++k)
  for (unsigned j= 0; j < N; ++j) {
    double tmp= x[j];
    for (unsigned i= 0; i < M; ++i)
      y[i]+= a[i*N+j] * tmp;
  }      
  cout << "Over columns first: time = " << t2.elapsed() << "\n";

//   boost::timer t3; 
//   for (unsigned i= 0; i < M; ++i) {
//     y[i]= 0.0;
//     for (unsigned j= 0; j < N; ++j)
//       y[i]+= a[i+j*M] * x[j];
//   }
//   cout << "Over rows first: time = " << t3.elapsed() << "\n";

//   boost::timer t4;
//   for (unsigned i= 0; i < M; ++i) 
//     y[i]= 0.0;
//   for (unsigned j= 0; j < N; ++j)
//     for (unsigned i= 0; i < M; ++i)
//       y[i]+= a[i+j*M] * x[j];      
//   cout << "Over columns first: time = " << t4.elapsed() << "\n";



  return 0;
}