#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

const int kdim=20;

void crout(int kdim1, int n, double a[kdim][kdim])
{
	double eps = 1.e-8;
	
	for(int j=1;j<n;j++)
		a[0][j] = a[0][j] / a[0][0];
	
	for(int k=1;k<n;k++)
	{
		for(int i=k;i<n;i++)
		{	
			double sum=0.e0;
			for(int j=0;j<=k-1;j++)
				sum+=a[i][j]*a[j][k];
			a[i][k]-=sum;
		}
		for(int j=k+1;j<n;j++)
		{
			double sum=0.e0;
			for(int i=0;i<=k-1;i++)
				sum+=a[k][i]*a[i][j];
			a[k][j]=(a[k][j]-sum)/a[k][k];

		}
	}	
	
}

int main(int argc, char** argv)
{
	int n, kdim1 = 20;
	ifstream in("matrix.dat");
	in>>n;
	
	double a[kdim][kdim];
	double b[kdim];
	double aux[kdim];
	double sol[kdim];

	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
			in>>a[i][j];
		in>>b[i];
	}
	in.close();
	crout(kdim1,n,a);
	
	ofstream out("res_lu.dat");
	out<<" L-matrix: "<<endl;
	for(int k1=0;k1<n;k1++)
	{
		for(int i=0;i<=k1;i++)
			out<<a[k1][i]<<"\t";
		out<<endl;
	}
	out<<" U-matrix: (UNIT diagonal omitted)"<<endl;
	for(int k1=0;k1<n;k1++)
	{
		for(int i=k1+1;i<n;i++)
			out<<a[k1][i]<<"\t";
		out<<endl;
	}
	out.close();

	aux[0]=b[0]/a[0][0];
	for(int i=1;i<n;i++)
	{	
		double sum=0.e0;
		for(int j=0;j<=i-1;j++)
		{
			sum+=a[i][j]*aux[j];
		}
		aux[i]=(b[i]-sum)/a[i][i];
	}
	sol[n-1]=aux[n-1];
	for(int i=n-2;i>=0;i--)
	{
		double sum = 0.e0;
		for(int j=n-1;j>=0;j--)
			sum+=a[i][j]*sol[j];
		sol[i]=aux[i]-sum;
	}

	ofstream out1("solution.dat");
	for(int i=0;i<n;i++)
		out1<<i<<"\t"<<sol[i]<<endl;
	out1.close();
}