#include<iostream.h>
#include<math.h>
void Doolittle(double A[][4],double L[][4],double U[][4],int n)
{  double s,t;
   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(r=0;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 zhijiejiefa(double L[][4],double U[][4],double b[],double x[],int n)
{  int i,k; 
   double s,t;
   double y[4];
   
   y[0]=b[0];
   for(i=1;i<n;i++)
   {  s=0.0;
      for(k=0;k<i;k++)
      {  s=s+L[i][k]*y[k];}
      y[i]=b[i]-s;
   }

   x[n-1]=y[n-1]/U[n-1][n-1];
   for(i=n-2;i>=0;i--)
   {  t=0.0;
      for(k=i+1;k<=n-1;k++)
      {  t=t+U[i][k]*x[k];}
      x[i]=(y[i]-t)/U[i][i];
   }
}
      
   
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];
   double x[4],b[4]={5,3,17,-13};
   
   int n=4,r,i;
   Doolittle(A,L,U,4);

   cout<<"L矩陣為:"<<endl;
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  cout<<L[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
   
   cout<<endl<<"U矩陣為:"<<endl;
   for(r=0;r<n;r++)
   {  for(i=0;i<n;i++)
      {  cout<<U[r][i]<<" "<<" "<<" ";
      }
      cout<<endl;
   }
  
  zhijiejiefa(L,U,b,x,4);
  cout<<endl<<"方程組的解為:"<<endl;
  for(i=0;i<n;i++)
  cout<<x[i]<<endl;
 
}