#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <string.h>

void calc_mkl(int n, double* A, double* w);
void calc_magma(int n, double* A, double *w);
double max_diff(int N, double *array1, double *array2);
void genRandomMatrix(int n, double *A);

int main( int argc, char** argv){

  int i,j,k;
 
  //size of main matrix
  int N = 8999;
  //size of small matrix
  int n;
 if (argc != 1){
   if (strcmp("-N", argv[1])==0) 
     n = atoi(argv[2]);
 }else{
   n = 8999;
 }
 //Read files
  printf("\n\nReading input file \n");
  
  FILE *fileptr;
  double *arrayd = new double[N*N];
  //fileptr=fopen("/exports/work/physics_epcc/genomic_prediction/genpred_benchmark_data/V.bin","rb");
  //fileptr=fopen("/exports/work/scratch/lcebaman/ACTA_REF_OUT/CPU_Input_A.bin","rb");
  //fread(arrayd, sizeof(double), N*N, fileptr);
  //fclose(fileptr);
  genRandomMatrix( N, arrayd);
  printf("Read succesful\n");

  double *A_magma = new double[n*n];
  double *A_mkl   = new double[n*n];
  
  for( i = 0; i < n; i++ )
    for( j = 0; j < n; j++ ){
      A_magma[i*n+j] = arrayd[i*N+j];
      A_mkl[i*n+j] = arrayd[i*N+j];
    }

  //eigenvalues 
  double *w_mkl   = new double[n];
  double *w_magma = new double[n];

  //calculate eigenvalues and eigenvectors using MKL
  printf("Calculating eigenvalues with MKL\n");
  calc_mkl(n, A_mkl, w_mkl);

  //calculate eigenvalues and eigenvectors using MAGMA
  printf("Calculating eigenvalues with MAGMA\n");
  calc_magma(n, A_magma, w_magma);

  
  for( i =0; i< n; i++){
    if (fabs(w_mkl[i] - w_magma[i]) > 0.0005){
      printf("%lf\n",w_mkl[i]);
      printf("%lf\n",w_magma[i]);
      exit(0);
    }
  }
  printf("\n");
  printf("Eigenvalues: SUCESS!!!!\n");
  
  printf("Calculating max difference eigenvectors ...\n");
  double d = max_diff(n, A_mkl, A_magma);
   
  if( d > 0.0005){
    printf("ERROR: max diff = %lf\n",d);
  }else{
    printf("Eigenvectors: SUCCESS: max diff = %lf\n",d);
  }
  
}

//calculates the max difference between the two input arrays
double max_diff(int N, double *array1, double *array2)
{
  double tmpdiff, diff;
  int i,j,k;
  diff=0.0;

  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      tmpdiff = fabs(array1[i*N+j] - array2[i*N+j]);
      if (tmpdiff > diff) diff = tmpdiff;
      
      }
    }
  return diff;
}
