#include<iostream.h>
#include<math.h>
void Doolittle(double A[][4],double L[][4],double U[][4],int n)
{  double s=0.0,t=0.0;
   int i,r,h,k;
   for(i=0;i<n;i++)
   {  for(r=0;r<n;r++)
      {  if(i==r)
         L[i][r]=1;
         if(i<r)
         L[i][r]=0;
         else U[i][r]=0;
      }
   }
   
   for(i=0;i<n;i++)
   {  U[0][i]=A[0][i];}
   for(i=1;i<n;i++)   
   {  L[i][0]=A[i][0]/U[0][0];}
   for(r=1;r<n;r++)
   {  for(i=r;i<n;i++)
      {  s=0.0;
         t=0.0;
         for(k=0;k<r;k++)
         {  s=s+L[r][k]*U[k][i];}
         U[r][i]=A[r][i]-s;
         h=i+1;
         for(k=0;k<r;k++)
         {  t=t+L[h][k]*U[k][r];}
         L[h][r]=(A[h][r]-t)/U[r][r];

      }
    }

}

void main()
{  double A[4][4]={{1,2,3,-1},{2,-1,9,-7},{-3,4,-3,19},{4,-2,6,-21}};
   double L[4][4],U[4][4];
   int n=4,r,i;
   Doolittle(A,L,U,4);
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  cout<<L[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
   cout<<endl;
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  cout<<U[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
}