#include "Matrix.h"

struct matrix* XLUDecomposition(struct matrix A){
	int dim,i,j,k;
	struct matrix L, U;
	struct matrix LU[2];
	
	/* Checking Dimensions */
	if(A.m!=A.n){
		fprintf(stderr, "Non square (%ix%i) matrix",A.m,A.n);
		L.m=0;L.n=0;L.mat=NULL;
		U.m=0;U.n=0;U.mat=NULL;
		LU[1]=L; LU[2]=U;
		return LU;
		};
	/* Initializing */
	dim = A.m;
	L = InitializeMatrix(dim,dim);
	U = InitializeMatrix(dim,dim);
	/* LU Decomp: Note currently recomputes all the zero elements */	
	for(i=0;i<dim;i++){
		L.mat[i][i]=1;
		for(j=i+1;j<dim;j++){
			L.mat[i][j] = (U.mat[i][j])/(U.mat[i][j]);
			/* should be k=i+1 -> to not compute the zeros */
			for(k=0;k<dim;k++)
				U.mat[i][k] -= (L.mat[i][j])*(U.mat[i][j]);
			};
		};
	/* Building the Array LU of matrix elements */
	LU[0]=L; LU[1]=U;
	return LU;
	};