#include
#include
#include
#include
#include
#include "gauss_omp.h"

int NUM_THREADS = 2;
int n = 1000;
double* A;
double* x;

double LiToDouble(LARGE_INTEGER x)
{
 double result = ((double)x.HighPart) * 4.294967296E9 + (double)((x).LowPart);
 return result;
}
double GetTime() {
 LARGE_INTEGER lpFrequency, lpPerfomanceCount;
 QueryPerformanceFrequency (&lpFrequency);
 QueryPerformanceCounter (&lpPerfomanceCount);
 return LiToDouble(lpPerfomanceCount)/LiToDouble(lpFrequency);
}

int main(int argc, char** argv)
{
 int i = 0;  
 if (argc > 2)
 {
  n = atoi(argv[1]);
  NUM_THREADS = atoi(argv[2]);
 }
 A = new double[ n * (n + 1) ];
 x = new double[ n ];
 matrix_create(n, A);

 printf("%f",solve_gauss(n,A,x, NUM_THREADS)); 

 delete[]x;
 delete[]A;
 return 0;
}
void matrix_create(int n, double* A)
{
 int i,j;
 for (i = 0; i < n; i++)
 {
  for (j = 0; j < n+1; j++)
  {
   A[i * (n+1) + j] = ((double)(rand() % 1000))*0.001;
  }
 }
}
double solve_gauss(int n, double* matrix, double* x, int nt)
{ 
   double t1,t2;    
   double EPSILON = 10e-10;
   int i,j,k,p;
  
   t1 = GetTime();  

   for (i = 0; i < n; i++)
   {
      double pivot = matrix[i * (n+1) + i];
      if (fabs(pivot) > EPSILON)
      {
  for(j = i; j < n + 1; j++) matrix[i * (n+1) + j] /= pivot;
#pragma omp parallel for private (j) schedule(dynamic) num_threads(NUM_THREADS)
   for(k = i + 1; k < n; k++)
    {
     double koef = matrix[k * (n+1) + i];
                    for (j = 0; j < n+1; j++) matrix[k * (n+1) + j] -= matrix[i * (n+1) + j]*koef;
                }
      }  
   }
#pragma omp parallel for private (k,j) schedule(dynamic) num_threads(NUM_THREADS)
   for(i = n - 1; i > = 0; i--) // цикл по строкам   
    for(k = i-1; k > = 0; k--) // цикл по строкам
     for(j = n; j > = 0; j--) // цикл по столбцам
            matrix[k * (n+1) + j] -= matrix[i * (n+1) + j] * matrix[k * (n+1) + i];  

   t2 = GetTime();
   for (i = 0; i < n; i++) x[i] = matrix[i * (n+1) + n];

   return t2-t1;
}
//int select(int n, int row, double* matrix)

void print_matrix(int n, double* A)
{
 int i,j;
 for ( i = 0; i < n; i++ )
 {
  printf("| ");
  for ( j = 0; j < n + 1; j++ )
  {
   if (j < n)
   {
    printf("%.2f ", A[ i * (n+1) + j ]);
   } else
   {
    printf("| %.2f |\n", A[ i * (n+1) + j ]);
   }
  }      
 }
 printf("\n");
}
void print_x(int n, double *x)
{
 int i;
 printf("x = |");
 for ( i = 0; i < n; i++ ) printf("%.2f ", x[i]); 
 printf("|\n"); 
}