/* ラインアドレス */

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

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

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

double norm(int n, int m, double *x){
  int i;
  double s;
  
  s = 0.0;
  for(i = 0; i < n; i++){
    s += x[i*m] * x[i*m];
  }
  return s;
}

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

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

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

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

  if((fp2 = fopen("stride-line64.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }
  
  if((fp3 = fopen("stride-line512.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }
  
  if((fp4 = fopen("stride-line-8.dat","w")) == NULL){
    printf("file open error\n");
    exit(1);
  }

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

  /* n=64で固定 */
  num = 64;

  for(m=start; m<=end; m*=2){

    data = (double*)malloc(num*(m-1) * sizeof(double));
    data2 = (double*)malloc(num*m * sizeof(double));
    data3 = (double*)malloc(num*(m+1) * sizeof(double));
    create_randam_array(num*(m-1), data);
    create_randam_array(num*m, data2);
    create_randam_array(num*(m+1), data3);

    result = perf(num, m-1, data);
    result2 = perf(num, m, data2);
    result3 = perf(num, m+1, data3);

    printf("%d\t%f\n", m-1, result);
    printf("%d\t%f\n", m, result2);
    printf("%d\t%f\n", m+1, result3);
    fprintf(fp, "%d\t%f\n", m-1, result);
    fprintf(fp, "%d\t%f\n", m, result2);
    fprintf(fp, "%d\t%f\n", m+1, result3);

    free(data);
    free(data2);
    free(data3);

  }

  /* 64の倍数 */
  for(m=64; m<=3000; m+=64){

    data = (double*)malloc(num*m * sizeof(double));
    create_randam_array(num*m, data);

    result = perf(num, m, data);

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

    free(data);

  }

  /* 512の倍数 */
  for(m=512; m<=30000; m+=512){

    data = (double*)malloc(num*m * sizeof(double));
    create_randam_array(num*m, data);

    result = perf(num, m, data);

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

    free(data);

  }

  /* 8の倍数 */
  for(m=8; m<=320; m+=8){

    data = (double*)malloc(num*m * sizeof(double));
    create_randam_array(num*m, data);

    result = perf(num, m, data);

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

    free(data);

  }

  fclose(fp);
  fclose(fp2);
  fclose(fp3);
  fclose(fp4);

}