#include <stdio.h>
#include <math.h>
/* luinv.c : 1999  3 8 by Oguni */
#define N 3
main(){
  double A[N][N]={1., 4., 3., 2., 5., 4., 1., -3., -2.};
  double L[N][N], U[N][N];
  double b[N]={1.,4.,5.}, x[N];
  double tol=0.000001;
  int n=3, i, j;
  luinv(A,L,U,tol,n);
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("%f", L[i][j]); }
    printf("\n");
  }
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("%f", U[i][j]); }
    printf("\n");
  }
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("%f", A[i][j]); }
    printf("\n");
  }
}
/* function lu */
lu(A,L,U,tol,n)
double A[N][N], L[N][N], U[N][N], tol;
int n;
{
  double akk, aik;
  int i, j, k;
  for (k=0; k<n-1; k++){
    if (fabs(A[k][k])<tol) {
      printf("The pivot is too small at step # %d", k); }
    akk=1.0/A[k][k];
    for (i=k+1; i<n; i++){
      aik=-A[i][k]*akk;
      for (j=k+1; j<n; j++){ A[i][j]+=aik*A[k][j]; }
    }
    for (j=k+1; j<n; j++){ A[k][j]*=akk; }
  }
  for (i=0; i<n-1; i++){
    for (j=0; j<=i; j++){ L[i][j]=A[i][j]; }
    for (j=i+1; j<n; j++){ U[i][j]=A[i][j]; }
    U[i][i]=1.0;
  }
  U[n-1][n-1]=1.0;
  for (j=0; j<n; j++){ L[n-1][j]=A[n-1][j]; }
}
/* function luinv */
luinv(A, L, U, tol, n)
double A[N][N], L[N][N], U[N][N], tol;
int n;
{
  double r, s;
  int i, j, k;
  lu(A, L, U, tol, n);
  A[n-1][n-1]=1.0/L[n-1][n-1];
  for (k=n-2; k>=0; k--){
    for (i=k; i>=0; i--){
      r=0.0;  s=0.0;
      for (j=i+1; j<n; j++){
        r-=U[i][j]*A[j][k+1]; 
        s-=A[k+1][j]*L[j][i]; }
      A[i][k+1]=r;
      A[k+1][i]=s/L[i][i];
    }
    s=0.0;
    for (j=k+1; j<n; j++){
      s+=A[k][j]*L[j][k]; }
    A[k][k]=(1.0-s)/L[k][k];
  }
}