#include "Matrix.h"

void LUDecomposition(struct matrix A, struct matrix L, struct matrix U){
	int dim,i,j,k;
	float pivotvalue, mult;
	
	/* Checking Dimensions */
	if(A.m!=A.n){
		fprintf(stderr, "Non square (%ix%i) input matrix",A.m,A.n);
		};
	dim = A.m;
	if(L.m!=dim || L.n!=dim){
		fprintf(stderr, "Innapropriately dimensioned L ");
		};
	if(U.m!=dim || U.n!=dim){
		fprintf(stderr, "Innapropriately dimensioned U ");
		};
	/* Initializing */
	for(i=0;i<dim;i++){
		for(j=0;j<dim;j++){
			L.mat[i][j]=0;
			U.mat[i][j]=A.mat[i][j];
		};
	};
	/* LU Decomp: Note currently recomputes all the zero elements */	
	for(i=0;i<dim;i++){
		L.mat[i][i]=1;
		pivotvalue =U.mat[i][i];
		for(j=i+1;j<dim;j++){
			mult = U.mat[j][i]/pivotvalue;
			L.mat[j][i] = mult;
			/* should be k=i+1 -> to not compute the zeros */
			for(k=0;k<dim;k++)
				U.mat[j][k] -= mult*(U.mat[i][k]);
			};
		};
	};