//devin scott
//CS 407
//Numerical analysis #6


void do_gauss(int n, double a[10][10], double x[10], double b[10]);


int main()
{


	double A[10][10];	
	double X[10];
	double B[10];
	int x = 0, y = 0, n = 0;
	printf("Analysis #6 Devin Scott CS 407 \n");

	for (x = 0; x < 10; x++)
		for (y = 0; y < 10; y++)
		{
			A[x][y] = 0.0;
			B[x] = 0.0;
		}
	printf("please input the dimension n of your nxn matrix");	
	scanf("%d", &n);
	printf("please input a row of %d coeficients with spaces between them\n",n);
	
	for (x = 0; x < n; x++)
	{
		printf("row %d:  ",x+1);
		for(y = 0; y < n; y++)
		{
			//printf("row %d col %d  ",x,y);
			scanf("%lf",&A[x][y]);
			
		
		}
		
		
		
		printf("B %d:  ",x+1);
		scanf("%lf",&B[x]);
		for(y = 0; y < n; y++)
		{
		printf("%lf  ",A[x][y]);
		}
		printf("  |  %lf\n",B[x]);
	}
	
	printf("\n\n\n\n SOLVING MATRIX:\n");
	
	for (x = 0; x < n; x++)
	{
		for(y = 0; y < n; y++)
		{
		printf("%lf  ",A[x][y]);
		}
		printf("  |  %lf\n",B[x]);
	}
	
	do_gauss(n - 1, A, X, B);
	printf("\n\n\n");
   	printf("RESULTS\n");

	for (x = 1; x <n+1; x++)
	{
		printf("X[%d] = %lf\n",x,X[x-1]);
	}

return 0;



}

void do_gauss(int n, double a[10][10], double x[10], double b[10])

{

int l[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

//real array (aij)nxn, (li)n;
//real array allocate (si)n;
int i=0;
int j=0;
int k=0;

double r=0.0;
double rmax=0.0;
double smax=0.0;
double xmult=0.0;
double tempval =0.0;
double s[10];

double sum=0.0;


for (i = 0; i<=n; i++) 
{
	l[i] = i;
	smax = 0.0;
	
	for(j=0; j<=n; j++)
	{
		//smax = max(smax, |aij|)
		tempval = (double)a[i][j];
		if(tempval < 0)
		{
			tempval = -tempval;
		}
		if(tempval > smax)
		{
			smax = tempval;
		}
		
	}

	s[i] = smax;

}


for (k=0; k<=n-1; k++)
{
	rmax=0;

	for(i=k; i<=n; i++)
	{
		//r=|ali,k/sli|
		r = a[l[i]][k] / s[l[i]];
		if (r < 0)
		{
			r = -r;
		}
		
		if (r > rmax) 
		{
			rmax = r;
			j=i;
		}

	}

	//lj=lk;
    //lk=lj;
	tempval = l[j];
	
	l[j] = l[k];
	
	l[k] = tempval;
	
	for (i=(k+1); i<=n; i++)
	{
		xmult = a[l[i]][k]/a[l[k]][k];
		a[l[i]][k] = xmult;
		
		for (j=(k+1); j<=n; j++)
		{
		//ali,j = ali,j - (xmult)alk,j;
		
		a[l[i]][j] = a[l[i]][j] - xmult * a[l[k]][j];
		}
	}


}		

////second part of procedure...................

for (k=0; k<=(n-1); k++)
	{
		for (i=(k+1); i<=n; i++)
		{
			b[l[i]] = b[l[i]] - a[l[i]][k] * b[l[k]];
		}
	}



	x[n] = b[l[n]]/a[l[n]][n];

	for (i=(n-1); i >=0; i--)
	{
		sum = b[l[i]];
		
		for (j=(i+1); j<=n; j++)
		{
		sum = sum - a[l[i]][j]*x[j];
		}
		

		x[i] = sum/a[l[i]][i];
	}

}