/* This demonstrates how data cache misses can affect the performance of an 
application. We show how the time/counts for a simple matrix multiplication
algorithm dramatically reduce when we employ a strip mining optimization. */

#include <stdio.h>

#include <Profile/Profiler.h>
#include <unistd.h>
#define SIZE 512
#define CACHE 64

double A[SIZE][SIZE], B[SIZE][SIZE], C[SIZE][SIZE];

double multiply(void)
{
  int i, j, k, n, m;
  int vl, sz, strip;
  TAU_PROFILE("multiply", "void (void)", TAU_USER);
  TAU_PROFILE_TIMER(t1,"regular", "", TAU_USER);
  TAU_PROFILE_TIMER(strip_timer,"stripmine", "", TAU_USER);


  for (n = 0; n < SIZE; n++)
    for (m = 0; m < SIZE; m++)
     {
      A[n][m] = B[n][m] = n + m ;
      C[n][m] = 0;
     }
  TAU_PROFILE_START(t1);
  sleep(2);
  for (i = 0; i < SIZE; i ++)
  { 
    for (j = 0; j < SIZE; j++)
    {
      for (k = 0; k < SIZE; k++)
  	C[i][j] += A[i][k] * B[k][j];
    }
  }
  TAU_PROFILE_STOP(t1);

  /* Now we employ the strip mining optimization */

  for(n = 0; n < SIZE; n++)
    for(m = 0; m < SIZE; m++)
      C[n][m] = 0; 
  
  TAU_PROFILE_START(strip_timer);
  sleep(1);
  for(i=0; i < SIZE; i++)
    for(k=0; k < SIZE; k++)
      for(sz = 0; sz < SIZE; sz+=CACHE)
      {
        //vl = min(SIZE-sz, CACHE);
  	vl = (SIZE - sz < CACHE ? SIZE - sz : CACHE); 
        for(strip = sz; strip < sz+vl; strip++)
          C[i][strip] += A[i][k]*B[k][strip];
      }
  TAU_PROFILE_STOP(strip_timer);

 
  return C[SIZE-10][SIZE-10];
  // So KCC doesn't optimize this loop away.
}
       
       


int checkResults() {
  const char **inFuncs;
  /* The first dimension is functions, and the second dimension is counters */
  double **counterExclusiveValues;
  double **counterInclusiveValues;
  int *numOfCalls;
  int *numOfSubRoutines;
  const char **counterNames;
  int numOfCouns;

  int numFunctions;
  const char **functionList;

  TAU_GET_FUNC_NAMES(functionList, numFunctions);


  int numCounters;
  const char ** counterList;

  TAU_GET_COUNTER_NAMES(counterList, numCounters);


  int WALL = -1;
  int VIRT = -1;
  int L1_DCM = -1;

  for(int j=0;j<numCounters;j++){
    if (strcmp("P_WALL_CLOCK_TIME",counterList[j]) == 0) {
      WALL = j;
    }
    if (strcmp("P_VIRTUAL_TIME",counterList[j]) == 0) {
      VIRT = j;
    }
    if (strcmp("PAPI_L1_DCM",counterList[j]) == 0) {
      L1_DCM = j;
    }
  }

  if (WALL == -1 || VIRT == -1 || L1_DCM == -1) {
    printf ("Error: Couldn't find P_WALL_CLOCK_TIME, P_VIRTUAL_TIME and PAPI_L1_DCM among the counters\n");
    return -1;
  }


  double regularWALL=-1, regularVIRT, regularL1=-1;
  double stripWALL=-1, stripVIRT, stripL1=-1;

  TAU_GET_FUNC_VALS(functionList, numFunctions,
		    counterExclusiveValues,
		    counterInclusiveValues,
		    numOfCalls,
		    numOfSubRoutines,
		    counterNames,
		    numOfCouns);


  for (int i=0; i < numFunctions; i++) {
    if (strcmp("regular",functionList[i]) == 0) {
      regularWALL = counterExclusiveValues[i][WALL];
      regularVIRT = counterExclusiveValues[i][VIRT];
      regularL1 = counterExclusiveValues[i][L1_DCM];
    }

    if (strcmp("stripmine",functionList[i]) == 0) {
      stripWALL = counterExclusiveValues[i][WALL];
      stripVIRT = counterExclusiveValues[i][VIRT];
      stripL1 = counterExclusiveValues[i][L1_DCM];
    }
  }

  printf ("Regular: P_WALL_CLOCK_TIME = %G\n", regularWALL);
  printf ("Regular: P_VIRTUAL_TIME = %G\n", regularVIRT);
  printf ("Regular: PAPI_L1_DCM = %G\n", regularL1);
  printf ("Strip: P_WALL_CLOCK_TIME = %G\n", stripWALL);
  printf ("Strip: P_VIRTUAL_TIME = %G\n", stripVIRT);
  printf ("Strip: PAPI_L1_DCM = %G\n", stripL1);


  int error = 0;

  if (regularWALL < 0 || regularL1 < 0 || stripWALL < 0 || stripL1 < 0) {
    printf ("Error, Couldn't find regular and stripmine timers!\n");
    error = 1;
  }

  if (stripL1 > regularL1) {
    printf ("Uh oh, strip mining optimization increased L1_DCM!!!\n");
    error = 1;
  }

  // OptixICC7 usually fails if we only check that strip > regular)
  if (stripVIRT > (regularVIRT * 1.5)) {
    printf ("Hmm, stripmining took longer, please investigate!\n");
    error = 1;
  }

  if (regularVIRT > regularWALL) {
    printf ("Hmm, regular virtual time is more than wall clock!\n");
    error = 1;
  }

  if (stripVIRT > stripWALL) {
    printf ("Hmm, strip virtual time is more than wall clock!\n");
    error = 1;
  }

  if (error == 0) {
    printf ("\nMultiple Counters and PAPI seem to be working!\n");
  }

  return error;
}


int main(int argc, char **argv) {
  TAU_PROFILE("main()", "int (int, char **)", TAU_DEFAULT);
  TAU_PROFILE_SET_NODE(0);

  multiply();

  //  return checkResults();
  checkResults();
  return 0;
}