#include<stdio.h>
#include<math.h>
#include<stdlib.h>

void LU(const int n, double a[][n], double L[][n], double U[][n])
{
   int i,j,k,l;

   for(i=0;i<n;i++)
    for(j=0;j<n;j++){   
     L[i][j] = 0.0;
     U[i][j] = 0.0;
    }
	   
//    for(j=0;j<n;j++){   
//     printf(" %lf %lf %lf \n",a[j][0],a[j][1],a[j][2]);
//    }
//     printf("\n");

/* The code below creates a perfect U matrix */   
	for(i=0;i<n;i++){
     	   L[i][i] = 1.0;
	  for(k=0;k<n;k++){
	   U[i][k] = a[i][k];
	   a[i][k] = a[i][k]/U[i][i];
	  }
	  
	  for(k=i+1;k<n;k++){
	   for(l=i+1;l<n;l++){
	    a[k][l] = a[k][l] - a[i][l]*a[k][i];
	   }

	   L[k][i] = a[k][i]/U[i][i];
	   
	  }
	  
	  for(k=i+1;k<n;k++){
	   for(l=0;l<i+1;l++){
	    a[k][l] = 0.0;
	   }
	  }
	}

//    for(j=0;j<n;j++){   
//     printf(" %lf %lf %lf \n",L[j][0],L[j][1],L[j][2]);
//    }

//     printf("\n");

//    for(j=0;j<n;j++){   
//     printf(" %lf %lf %lf \n",U[j][0],U[j][1],U[j][2]);
//    }

//    printf("\n");
}