/****************************************************************************************************
* Program 3.3: LU-decomposition Method, based on Guass elimination, for systems of linear equations.*
*				Coded by Dilip Datta                                                                *
*                                                                                                   *
*****************************************************************************************************/

# include <stdio.h>
# include <math.h>
# include <stdlib.h>
# include <malloc.h>

/********************************************************************/

int main()
{
	int i, j, n;
	double **a, *b, *x;

	void LU_decompose_gauss(int n, double **a);
	void LU_solution(int n, double **a, double *b, double *x);

	printf("\nNumber of equations: ");										/* Step-1 */
	scanf("%d", &n);

	x = (double *)malloc( (n+1)*sizeof(double) );
	b = (double *)malloc( (n+1)*sizeof(double) );
	a = (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++)
		{
			printf("a[%d][%d] = ", i, j);
			scanf("%lf", &a[i][j]);
		}
		printf("b[%d] = ", i);
		scanf("%lf", &b[i]);
	}

	LU_decompose_gauss(n, a);
	LU_solution(n, a, b, x);

	printf("\n\n");
	printf("Solution of the given system of linear equations\n");
	printf("------------------------------------------------\n");
	printf("x = (");
	for(i = 1; i <= n-1; i++)	printf("%lf, ", x[i]);
	printf("%lf)\n\n", x[n]);

	return(0);
}
/********************************************************************/

void LU_decompose_gauss(int n, double **a)
/* Sub-routine for LU-decomposition based on Guass elimination method */
{
	int i, j, k;
	int p;
	double factor;
	double big, dummy;

	// LU decompision through forward elimination with partial pivoting
	for(k = 1; k <= n-1; k++)
	{
		// Partial (column) pivoting
		p = k;
		big = fabs(a[k][k]);
		for(i = k+1; i <= n; i++)
		{
			if(big < fabs(a[i][k]))
			{
				big = fabs(a[i][k]);
				p = i;
			}
		}

		if(p != k)
		{
			for(j = k; j <= n; j++)
			{
				dummy = a[p][j];
				a[p][j] = a[k][j];
				a[k][j] = dummy;
			}
		}

		// Forward elimination
		for(i = k+1; i <= n; i++)
		{
			factor = a[i][k]/a[k][k];
			a[i][k] = factor;										// L matrix
			for(j = k+1; j <= n; j++)	a[i][j] -= a[k][j]*factor;	// U matrix
		}
	}

	return;
}
/********************************************************************/

void LU_solution(int n, double **a, double *b, double *x)
/* Sub-routine for solutions of systems of linear equation using L & U matrices */
{
	int i, j;
	double sum;

	// Obtain d (in b only) vector by forward substitution of L matrix
	for(i = 2; i <= n; i++)
	{
		sum = b[i];
		for(j = 1; j <= i-1; j++)	sum -= a[i][j]*b[j];
		b[i] = sum;
	}

	// Obtain x by back subsitution
	x[n] = b[n]/a[n][n];
	for(i = n-1; i >= 1; i--)
	{
		sum = 0;
		for(j = i+1; j <= n; j++)	sum += a[i][j]*x[j];
		x[i] = (b[i]-sum)/a[i][i];
	}

	return;
}
/********************************************************************/