// Checking speed of std::vector vs raw pointers for indirection
// Anders Logg 2013-01-09
//
// Results / seconds
//
//          raw     stl   raw_flattened  stl_flattened
//---------------------------------------------------
// g++      1.9    14.7   0.88           7.6
// g++ -O1  1.0    1.0    0.40           0.40
// g++ -O2  0.65   0.65   0.42           0.42
// g++ -O3  0.65   0.65   0.42           0.42

#include <iostream>
#include <vector>

#define N 100000000

void bench_raw()
{
  double** x;
  x = new double*[4];
  for (int i = 0; i < 4; i++)
    x[i] = new double[3];

  for (int i = 0; i < 4; i++)
    for (int j = 0; j < 3; j++)
      x[i][j] = 0.5;

  for (int n = 0; n < N; n++)
  {
    x[0][0] = x[1][0] - x[0][0];
    x[0][1] = x[2][0] - x[0][0];
    x[0][2] = x[3][0] - x[0][0];
    x[1][0] = x[1][1] - x[0][1];
    x[1][1] = x[2][1] - x[0][1];
    x[1][2] = x[3][1] - x[0][1];
    x[2][0] = x[1][2] - x[0][2];
    x[2][1] = x[2][2] - x[0][2];
    x[2][2] = x[3][2] - x[0][2];
    x[3][0] = x[0][2] - x[0][2];
    x[3][1] = x[1][1] - x[0][2];
    x[3][2] = x[2][0] - x[0][2];
  }

  for (int i = 0; i < 4; i++)
    for (int j = 0; j < 3; j++)
      std::cout << x[i][j] << " ";
  std::cout << std::endl;
}

void bench_raw_flattened()
{
  double x[12];

  for (int i = 0; i < 12; i++)
      x[i] = 0.5;

  for (int n = 0; n < N; n++)
  {
    x[0]  = x[4] - x[0];
    x[1]  = x[6] - x[0];
    x[2]  = x[9] - x[0];
    x[4]  = x[5] - x[1];
    x[5]  = x[7] - x[1];
    x[6]  = x[10] - x[1];
    x[6]  = x[6] - x[2];
    x[7]  = x[5] - x[2];
    x[8]  = x[11] - x[2];
    x[9]  = x[2] - x[2];
    x[10] = x[3] - x[2];
    x[11] = x[6] - x[2];
  }

  for (int i = 0; i < 12; i++)
    std::cout << x[i] << " ";
  std::cout << std::endl;
}

void bench_stl()
{
  std::vector<std::vector<double> > x(4, std::vector<double>(3, 0.5));

  for (int n = 0; n < N; n++)
  {
    x[0][0] = x[1][0] - x[0][0];
    x[0][1] = x[2][0] - x[0][0];
    x[0][2] = x[3][0] - x[0][0];
    x[1][0] = x[1][1] - x[0][1];
    x[1][1] = x[2][1] - x[0][1];
    x[1][2] = x[3][1] - x[0][1];
    x[2][0] = x[1][2] - x[0][2];
    x[2][1] = x[2][2] - x[0][2];
    x[2][2] = x[3][2] - x[0][2];
    x[3][0] = x[0][2] - x[0][2];
    x[3][1] = x[1][1] - x[0][2];
    x[3][2] = x[2][0] - x[0][2];
  }

  for (int i = 0; i < 4; i++)
    for (int j = 0; j < 3; j++)
      std::cout << x[i][j] << " ";
  std::cout << std::endl;
}

void bench_stl_flattened()
{
  std::vector<double> x(12);

  for (int n = 0; n < N; n++)
  {
    x[0]  = x[4] - x[0];
    x[1]  = x[6] - x[0];
    x[2]  = x[9] - x[0];
    x[4]  = x[5] - x[1];
    x[5]  = x[7] - x[1];
    x[6]  = x[10] - x[1];
    x[6]  = x[6] - x[2];
    x[7]  = x[5] - x[2];
    x[8]  = x[11] - x[2];
    x[9]  = x[2] - x[2];
    x[10] = x[3] - x[2];
    x[11] = x[6] - x[2];
  }

  for (int i = 0; i < 12; i++)
    std::cout << x[i] << " ";
  std::cout << std::endl;
}

int main(int argc, char* argv[])
{
  if (argc != 2)
  {
    std::cout << "Missing argument: 0 = raw 1 = stl 2 = raw_flattened 3 = stl_flattened" << std::endl;
    return 1;
  }

  if (argv[1][0] == '0')
    bench_raw();
  else if (argv[1][0] == '1')
    bench_stl();
  else if (argv[1][0] == '2')
    bench_raw_flattened();
  else if (argv[1][0] == '3')
    bench_stl_flattened();

  return 0;
}