//metoda eliminacji Gaussa_Jordana
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <omp.h>

int main(int argc, char *argv[])
{
	double **A, *b, *x;
	double c;
	double t1, t2;
	long int i, j, k, l, n;
	
	n = atol(argv[1]);
	A = (double**)malloc((n + 1) * sizeof(double));
	b = (double*)malloc((n + 1) * sizeof(double));
	x = (double*)malloc((n + 1) * sizeof(double));
	for(i = 1; i <= n; i++)
		A[i] = (double*)malloc((n + 1) * sizeof(double*));
	
	for(i = 1; i <= n; i++)
	{
		for(j = 1; j <= n; j++)
			A[j][i] = (double)(i + j) / 1000.0;
		
		b[i] = (double)(i) / 1000.0;
	}
	
	
	A[1][1] = 1;
	A[1][2] = 2;
	A[2][1] = 3;
	A[2][2] = 4;
	b[1] = 5;
	b[2] = 6;
	
t1 = omp_get_wtime();

	//algorytm eliminacji Gaussa - etap eliminacji zmiennych
	//pętla wierszy eliminujących
	for(k = 1; k <= n; k++)
	{
		//pętla modyfikacji wierszy poniżej k-tego wiersza
    #pragma omp parallel for shared(A,b,n,k) private(i,j,c)
		for(i = 1; i <= n; i++)
		{
			if(i==k) continue;
			c = A[i][k] / A[k][k];
			
			//pętla kolumn w danym wierszu
      
			for(j = k; j <= n; j++)
				A[i][j] = A[i][j] - c * A[k][j];
				
			b[i] = b[i] - c * b[k];
		}	
	}
	for(l=1;l<=n;l++)
  x[l]=b[l]/A[l][l];
	//postępowanie odwrotne
	x[n] = b[n] / A[n][n];
	for(k = n - 1; k >= 1; k--)
	{
		x[k] = b[k];
		
		for(j = k + 1; j <= n; j++)
			x[k] = x[k] - A[k][j] * x[j];
			
		x[k] = x[k] / A[k][k];
	}
t2 = omp_get_wtime();
		
	
	printf("Czas obliczeń = %f\n", t2 - t1);
	printf("Rozw.: x[1] = %f, x[2] = %f\n", x[1], x[2]);
	
	
	return 0;
}