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

using namespace std;

const unsigned M= 20, N= M;


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

  fill(&a[0], &a[M*N], 2.0);
  fill(&b[0], &b[M*N], 3.0);
  
  boost::timer t1; 
  for(unsigned k= 0; k< 3000000; ++k)
  for (unsigned i= 0; i < M; ++i)
    for (unsigned j= 0; j < N; ++j)
      c[i*N+j]= a[i*N+j] + b[i*N+j];
  cout << "Over rows first: time = " << t1.elapsed() << "\n";

//   boost::timer t2;
//   for(unsigned k= 0; k< 10; ++k)
//   for (unsigned j= 0; j < N; ++j)
//     for (unsigned i= 0; i < M; ++i)
//       c[i*N+j]= a[i*N+j] + b[i*N+j];
//   cout << "Over columns first: time = " << t2.elapsed() << "\n";

//   boost::timer t3; 
//   for(unsigned k= 0; k< 10; ++k)
//   for (unsigned i= 0; i < M; ++i)
//     for (unsigned j= 0; j < N; ++j)
//       c[i+j*M]= a[i+j*M] + b[i+j*M];
//   cout << "Over rows first: time = " << t3.elapsed() << "\n";

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

  boost::timer t5; 
  for(unsigned k= 0; k< 3000000; ++k)
    for (unsigned i= 0; i < M*N; ++i)
      c[i]= a[i] + b[i];
  cout << "Vector: time = " << t5.elapsed() << "\n";

  return 0;
}