// *LU decomposition & determinant at O(n^3)* by Krzysztof Bosak, 1998.12.05

#include <iostream.h>
#include <assert.h>
#include <memory.h>

#define OUTPUT

#define TRANSTABSEMI(x, y) (tab[(x)+indexy[(y)]])
// semi ==> possinle rows swapped, no columns swapping
#define TRANSTABFULL(x, y) (tab[indexx[(x)]+indexy[(y)]])
// full ==> possible rows and columns swapping

void testout(const int size, const double*const tab, const int*const indexy, const int*const indexx)
{
	for(int x, y=0; y<size; y++)
	{
		for(x=0; x<size; x++)
			cout<<TRANSTABFULL(x, y)<<'\t';
		cout<<endl;
	}
}

void paluout(const int size, const double*const src, const double*const tab, const int*const indexx, const int*const indexy)
{
	int x, y;
	cout.precision(4);
	cout<<"--- swapped A matrix == right side of equation ---\n";
	for(y=0; y<size; y++)
	{
		for(x=0; x<size; x++)
			cout<<src[indexx[x]+indexy[y]]<<'\t';
		if(y<size-1)
			cout<<endl;
	}
	cout<<"P1==[";
	for(x=0; x<size; x++)
		cout<<(indexy[x]/size+1);
	cout<<"]\tP2==[";
	for(x=0; x<size; x++)
		cout<<(indexx[x]+1);
	cout<<"]\n";
	cout<<"--- P1*A*P2 == L*U - result ----------------------\n";
	for(y=0; y<size; y++)			//P1
	{
		for(x=0; x<size; x++)
		{
			if(x*size==indexy[y])
				cout<<"1\t";
			else
				cout<<"0\t";
		}
		cout<<endl;
	}
	cout<<"*\n";					//A
	for(y=0; y<size; y++)
	{
		for(x=0; x<size; x++)
			cout<<src[x+y*size]<<'\t';
		cout<<endl;
	}
	cout<<"*\n";
	for(y=0; y<size; y++)			//P2
	{
		for(x=0; x<size; x++)
		{
			if(y==indexx[x])
				cout<<"1\t";
			else
				cout<<"0\t";
		}
		cout<<endl;
	}
	cout<<"=\n";
	for(y=0; y<size; y++)
	{
		for(x=0; x<size; x++)
		{
			if(x>y)
				cout<<"0\t";
			if(x==y)
				cout<<"1\t";
			if(x<y)
				cout<<tab[indexx[x]+indexy[y]]<<'\t';
		}
		cout<<endl;
	}
	cout<<"*\n";
	for(y=0; y<size; y++)
	{
		for(x=0; x<size; x++)
		{
			if(x<y)
				cout<<"0\t";
			if(x>=y)
				cout<<tab[indexx[x]+indexy[y]]<<'\t';
		}
		cout<<endl;
	}
}

template<class typea, class typeb>
void swap(typea &aa, typeb &bb)
{
	typea temp=aa;
	aa=bb;
	bb=temp;
}

double DET(const int size, const double*const src)
{

	if(size==1)
		return *src;
	const int max=size*size;
	bool sign=false;
	double*const tab=new double[max];
	memcpy(tab, src, max*sizeof(double));
	int*const indexy=new int[size<<1], *const indexx=indexy+size;
	int x, mul;
	for(x=0, mul=0; x<size; x++, mul+=size)
	{		
		indexx[x]=x;
		indexy[x]=mul;
	}
	for(x=0; x<=size-2; x++)
	{
		////////////////Axial factor semi-selection below
		double maxabs=0;
		int ybest=0xFFFF;
		for(int ytest=x; ytest<size; ytest++)
		{
			double value=TRANSTABSEMI(x, ytest);
			if(value<0)
				value=-value;
			if(value>maxabs)
			{
				maxabs=value;
				ybest=ytest;
			}
		}
		if(maxabs==0.0)
			return 0;
		if(x!=ybest)
			swap(indexy[x], indexy[ybest]);
		if(int(x)!=ybest)
			sign=!sign;
		//////////////comment this region to avoid axial factor semi-selection
		for(int y=x+1; y<size; y++)
		{
			TRANSTABSEMI(x, y)/=TRANSTABSEMI(x, x);
			for(int xx=x+1; xx<size; xx++)
				TRANSTABSEMI(xx, y)-=TRANSTABSEMI(x, y)*TRANSTABSEMI(xx, x);
		}
	}
	double result=TRANSTABSEMI(0, 0);
	for(x=1; x<size; x++)
		result*=TRANSTABSEMI(x, x);
	if(sign)
		return -result;
	return result;
}

double* PALU(const int size, const double*const src, const double*const vec=NULL)
{
	const int max=size*size;
	double*const tab=new double[max];
	memcpy(tab, src, max*sizeof(double));
	int*const indexy=new int[size<<1], *const indexx=indexy+size;
	int x, mul;
	for(x=0, mul=0; x<size; x++, mul+=size)
	{		
		indexx[x]=x;
		indexy[x]=mul;
	}
	for(x=0; x<=size-2; x++)
	{
		////////////////Axial factor full selection below
		/*int xbest=0xFFFF, ybest=0xFFFF;
		double maxabs=0;
		for(int ytest=x; ytest<size; ytest++)
		{
			for(int xtest=x; xtest<size; xtest++)
			{
				double value=TRANSTABFULL(xtest, ytest);
				if(value<0)
					value=-value;
				if(value>maxabs)
				{
					maxabs=value;
					xbest=xtest;
					ybest=ytest;
				}
			}
		}
		assert(maxabs>0);
		if(x!=ybest)
			swap(indexy[x], indexy[ybest]);//row swapping (~P1)
		if(x!=xbest)
			swap(indexx[x], indexx[xbest]);//column swapping (~P2)*/
		//////////////comment this region to avoid axial factor full selection
		for(int y=x+1; y<size; y++)
		{
			TRANSTABFULL(x, y)/=TRANSTABFULL(x, x);
			for(int xx=x+1; xx<size; xx++)
				TRANSTABFULL(xx, y)-=TRANSTABFULL(x, y)*TRANSTABFULL(xx, x);
		}
		//testout(size, tab, indexy, indexx);
		//cout<<endl;
	}
#ifdef OUTPUT
	paluout(size, src, tab, indexx, indexy);
#endif
	if(vec)
	{
		double *semi=new double[size<<1], *result=semi+size;
		int x, y;
		for(y=0; y<size; y++)
		{
			semi[y]=vec[y];//[indexy[y]/size];
			for(x=0; x<y; x++)
				semi[y]-=TRANSTABFULL(x, y)*semi[x];
		}
		cout<<"Vector y: Ly=b          (y=Ux)\n ";
		int i;
		for(i=0; i<size; i++)
			cout<<semi[i]<<endl;

		/*for(y=size-1; y>=0; y--)
		{
			result[y]=semi[y];
			for(x=y+1; x<size; x++)
				result[y]-=TRANSTABFULL(x, y)*result[x];
			result[y]/=TRANSTABFULL(y, y);
		}
		cout<<"Vector x: Ux=y          (y=Ux)\n ";
		for(i=0; i<size; i++)
			cout<<result[i]<<endl;*/
	}
	return tab;
}

int main()
{
	const double tab[]={1, 1, 1, 2, 5, 6, 7 ,8, 9}, vec[]={4, 2, 3};
	const int size=3;
	PALU(size, tab, vec);
	//cout<<"===========\n|A| = "<<DET(size, tab)<<endl;
	return 0;
}