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

void LU(const int n, double a[][n], double L[][n], double U[][n]);

void inverse(const int n, double A[][n], double InvA[][n])
{
	int i,j,k,l,m;
	double L[n][n],U[n][n];
	double C[n],z[n],x[n];

	LU(n,A,L,U);

	for(l=0;l<n;l++)
	{
	       	
	  for(m=0;m<n;m++)
	  {
	    
	    z[m] = 0.0;
	    x[m] = 0.0;
	    if(l==m){
	     C[m] = 1.0;
	    }
	    else{
	     C[m] = 0.0;
	    }
	  }

	  for(i=0;i<n;i++){
	   for(k=0;k<i;k++){
	     z[i] = z[i] - L[i][k]*z[k]/L[i][i];
	   }
	     z[i] = z[i] + C[i]/L[i][i]; 
	  }

	  for(i=n;i>0;i--){
	   for(k=n;k>i;k--){
	     x[i-1] = x[i-1] - U[i-1][k-1]*x[k-1]/U[i-1][i-1];
	   }
	     x[i-1] = x[i-1] + z[i-1]/U[i-1][i-1]; 
	  }

	  for(i=0;i<n;i++){
	    InvA[i][l] = x[i];
	  }
	  
	}
}