#include "iostream.h"

#define SIZE 21

//A*X=Z, A=L*U, U*X=Y, L*Y=Z
double A[SIZE][SIZE];
double Z[SIZE];
double X[SIZE];
double Y[SIZE];

void InitLU(int maxsize)
{
	for(int i=1;i<=maxsize;i++)
	{
		for(int k=1;k<=maxsize;k++)
		{
			A[i][k] = 0;
		}
		Z[i] = 0;
		X[i] = 0;
	}
	return;
}

//Move max column to the top, work from a[step][step] to a[maxsize][step]
void FindMax(int maxsize, int step)
{
	double max=0;
	double temp=0;
	int row=step;
	for(int i=step;i<=maxsize;i++)
	{
		if(A[i][step]<0) temp=-1.0*A[i][step];
		else temp=A[i][step];

		if(temp>max)
		{
			max = temp;
			row = i;
		}
	}

	if(row!=step)
	{
		for(int k=1;k<=maxsize;k++)
		{
			temp = A[row][k];
			A[row][k] = A[step][k];
			A[step][k] = temp;
		}
		
		temp = Z[row];
		Z[row] = Z[step];
		Z[step] = temp;
	}
	return;
}


// LU Decomposition from A[1][1] to A[maxsize][maxsize]
int LUDecomposition(int maxsize)
{
	int step = 1;
	for(; step<maxsize; step++)
	{
		FindMax(maxsize, step);
		for(int j=step+1; j<=maxsize; j++)
		{
			A[step][j] /= A[step][step];
		}

		for(int i=step+1;i<=maxsize;i++)
			for(int j=step+1;j<=maxsize;j++)
				A[i][j] -= A[i][step]*A[step][j];
	}

	for(int i=1;i<=maxsize;i++)
	{
		for(int k=1; k<=i-1; k++)
			Y[i] += A[i][k]*Y[k];
		Y[i] = (Z[i]-Y[i])/A[i][i];
	}

	for(int ii=maxsize;ii>=1;ii--)
	{
		for(int k=ii+1; k<=maxsize; k++)
			X[ii] +=A[ii][k]*X[k];
		X[ii] = Y[ii] - X[ii];
	}

	return 1;
}