/* $B%"%=%7%"%F%#%S%F%#(B  */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define K 512 //$B0l<!%-%c%C%7%e%i%$%s(B
#define m 4096 //$B0l<!%-%c%C%7%e%5%$%:$NG\?t(B
#define m2 262144 //$BFs<!%-%c%C%7%e%5%$%:$NG\?t(B

void create_randam_array(int n, double *x){
  int i;

  for(i = 0; i < n; i++){
    x[i] = 123.45*(i+1);
  }
}

double cmp(int l, int J, double *x){
  int i, j, k;
  double s;
  
  s = 0.0;
  //  for(i = 0; i < n; i++)
    for(k = 0; k < K; k++)
      for(j = 0; j < J; j++)
	s += x[8*k + m*j] * x[8*k + m*j + l];

    return s;
}

double cmp0(int l, double *x){
  int k;
  double s;

  s = 0.0;
  for(k = 0; k < K; k++)
    s += x[8*k] * x[8*k + l];
  return s;
}

double perf(int n, int l, int J, double *x){
  int k, t;
  double s;
  
  t=clock();
  for(k = 0; k < 1000000000/n; k++){
    s = cmp(l, J, x);
  }
  return 2.0 * K * k * J / (clock() - t);
}

double perf0(int n, int l, double *x){
  int k, t;
  double s;

  t=clock();
  for(k = 0; k < 1000000000/n; k++){
    s = cmp0(l, x);
  }
  return 2.0 * K * k / (clock() - t);
}

int main(int argc, char *argv[]) {
  FILE *fp, *fp2, *fp3, *fp4;
  int num, start, end, l;
  double *data;
  double result;

  if (argc != 3){
    exit(1);
  }

  start= atoi(argv[1]);
  end = atoi(argv[2]);
  
  if((fp = fopen("cmp1.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }

  if((fp2 = fopen("cmp2.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }

  if((fp3 = fopen("cmp3.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }

  if((fp4 = fopen("cmp0.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }

  printf("l\tMflops\n");
  fprintf(fp, "l\tMflops\n");
  fprintf(fp2, "l\tMflops\n");
  fprintf(fp3, "l\tMflops\n");
  fprintf(fp4, "l\tMflops\n");

  /* n=1024$B$G8GDj(B */
  num = 1024;

  for(l=start; l<=end; l *= 2){
    data = (double*)malloc((8*K+l) * sizeof(double));
    create_randam_array((8*K+l), data);

    result = perf0(num, l, data);

    printf("%d\t%f\n", l, result);

    fprintf(fp4, "%d\t%f\n", l, result);

    free(data);
  }
  fclose(fp4);

  for(l=start; l<=end; l *= 2){
    data = (double*)malloc((8*K+1*m+l) * sizeof(double));
    create_randam_array((8*K+1*m+l), data);

    result = perf(num, l, 1, data);

    printf("%d\t%f\n", l, result);

    fprintf(fp, "%d\t%f\n", l, result);

    free(data);
  }
  fclose(fp);

  for(l=start; l<=end; l *= 2){
    data = (double*)malloc((8*K+2*m+l) * sizeof(double));
    create_randam_array((8*K+2*m+l), data);

    result = perf(num, l, 2, data);

    printf("%d\t%f\n", l, result);

    fprintf(fp2, "%d\t%f\n", l, result);

    free(data);
  }

  fclose(fp2);

  for(l=start; l<=end; l *= 2){
    data = (double*)malloc((8*K+3*m+l) * sizeof(double));
    create_randam_array((8*K+3*m+l), data);

    result = perf(num, l, 3, data);

    printf("%d\t%f\n", l, result);

    fprintf(fp3, "%d\t%f\n", l, result);

    free(data);
  }

  fclose(fp3);

}
