#include <petscsys.h>
#include <petsctime.h>

extern int BlastCache(void);
extern int test1(void);
extern int test2(void);

int main(int argc, char **argv)
{
  PetscCall(PetscInitialize(&argc, &argv, 0, 0));
  PetscCall(test1());
  PetscCall(test2());
  PetscCall(PetscFinalize());
  return 0;
}

int test1(void)
{
  PetscLogDouble t1, t2;
  double         value;
  int            i, ierr, *z, *zi, intval;
  PetscScalar   *x, *y;
  PetscRandom    r;

  PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
  PetscCall(PetscRandomSetFromOptions(r));
  PetscCall(PetscMalloc1(20000, &x));
  PetscCall(PetscMalloc1(20000, &y));

  PetscCall(PetscMalloc1(2000, &z));
  PetscCall(PetscMalloc1(2000, &zi));

  /* Take care of paging effects */
  PetscCall(PetscTime(&t1));

  /* Form the random set of integers */
  for (i = 0; i < 2000; i++) {
    PetscCall(PetscRandomGetValue(r, &value));
    intval = (int)(value * 20000.0);
    z[i]   = intval;
  }

  for (i = 0; i < 2000; i++) {
    PetscCall(PetscRandomGetValue(r, &value));
    intval = (int)(value * 20000.0);
    zi[i]  = intval;
  }
  /* fprintf(stdout,"Done setup\n"); */

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[i] = y[i];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[i]", (t2 - t1) / 2000.0);

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 500; i += 4) {
    x[i]     = y[z[i]];
    x[1 + i] = y[z[1 + i]];
    x[2 + i] = y[z[2 + i]];
    x[3 + i] = y[z[3 + i]];
  }
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[idx[i]] - unroll 4", (t2 - t1) / 2000.0);

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[i] = y[z[i]];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[idx[i]]", (t2 - t1) / 2000.0);

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 1000; i += 2) {
    x[i]     = y[z[i]];
    x[1 + i] = y[z[1 + i]];
  }
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[idx[i]] - unroll 2", (t2 - t1) / 2000.0);

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[z[i]] = y[i];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[z[i]] = y[i]", (t2 - t1) / 2000.0);

  PetscCall(BlastCache());

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[z[i]] = y[zi[i]];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[z[i]] = y[zi[i]]", (t2 - t1) / 2000.0);

  PetscCall(PetscArraycpy(x, y, 10));
  PetscCall(PetscArraycpy(z, zi, 10));
  PetscCall(PetscFree(z));
  PetscCall(PetscFree(zi));
  PetscCall(PetscFree(x));
  PetscCall(PetscFree(y));
  PetscCall(PetscRandomDestroy(&r));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int test2(void)
{
  PetscLogDouble t1, t2;
  double         value;
  int            i, ierr, z[20000], zi[20000], intval, tmp;
  PetscScalar    x[20000], y[20000];
  PetscRandom    r;

  PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
  PetscCall(PetscRandomSetFromOptions(r));

  /* Take care of paging effects */
  PetscCall(PetscTime(&t1));

  for (i = 0; i < 20000; i++) {
    x[i]  = i;
    y[i]  = i;
    z[i]  = i;
    zi[i] = i;
  }

  /* Form the random set of integers */
  for (i = 0; i < 20000; i++) {
    PetscCall(PetscRandomGetValue(r, &value));
    intval    = (int)(value * 20000.0);
    tmp       = z[i];
    z[i]      = z[intval];
    z[intval] = tmp;
  }

  for (i = 0; i < 20000; i++) {
    PetscCall(PetscRandomGetValue(r, &value));
    intval     = (int)(value * 20000.0);
    tmp        = zi[i];
    zi[i]      = zi[intval];
    zi[intval] = tmp;
  }
  /* fprintf(stdout,"Done setup\n"); */

  /* PetscCall(BlastCache()); */

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[i] = y[i];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[i]", (t2 - t1) / 2000.0);

  /* PetscCall(BlastCache()); */

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) y[i] = x[z[i]];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[i] = y[idx[i]]", (t2 - t1) / 2000.0);

  /* PetscCall(BlastCache()); */

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) x[z[i]] = y[i];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[z[i]] = y[i]", (t2 - t1) / 2000.0);

  /* PetscCall(BlastCache()); */

  PetscCall(PetscTime(&t1));
  for (i = 0; i < 2000; i++) y[z[i]] = x[zi[i]];
  PetscCall(PetscTime(&t2));
  fprintf(stdout, "%-27s : %e sec\n", "x[z[i]] = y[zi[i]]", (t2 - t1) / 2000.0);

  PetscCall(PetscRandomDestroy(&r));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int BlastCache(void)
{
  int          i, ierr, n = 1000000;
  PetscScalar *x, *y, *z, *a, *b;

  PetscCall(PetscMalloc1(5 * n, &x));
  y = x + n;
  z = y + n;
  a = z + n;
  b = a + n;

  for (i = 0; i < n; i++) {
    a[i] = (PetscScalar)i;
    y[i] = (PetscScalar)i;
    z[i] = (PetscScalar)i;
    b[i] = (PetscScalar)i;
    x[i] = (PetscScalar)i;
  }

  for (i = 0; i < n; i++) a[i] = 3.0 * x[i] + 2.0 * y[i] + 3.3 * z[i] - 25. * b[i];
  for (i = 0; i < n; i++) b[i] = 3.0 * x[i] + 2.0 * y[i] + 3.3 * a[i] - 25. * b[i];
  for (i = 0; i < n; i++) z[i] = 3.0 * x[i] + 2.0 * y[i] + 3.3 * a[i] - 25. * b[i];
  PetscCall(PetscFree(x));
  PetscFunctionReturn(PETSC_SUCCESS);
}