#include
#include
#include
using namespace std;

class Matrix
{
public:
 Matrix();
 ~Matrix();
 
 void Crout();
 void Doolittle();

private:
 int n;
 double **Matrix_data;

};

Matrix::Matrix()
{
 ifstream fin("fin.txt",ios::in);
 fin >> n;
 cout << n << endl;

 Matrix_data=new(std::nothrow) double*[n];
 for(int i = 0; i < n;i++)
 {
  Matrix_data[i] = new(std::nothrow) double[n];

 }

 cout << "The Matrix =" << endl;
 for(int i = 0; i < n; i++)
 {
  for(int j = 0; j < n; j++)
  {
   fin >> Matrix_data[i][j];
   
   cout << setw(10) << setiosflags(ios::fixed) <<setprecision(4) <<Matrix_data[i][j];
  }
  cout << endl;
 }

 

 fin.close();
}

void Matrix::Crout()
{
 double **C_L;
 double **C_U;
 C_L = new(std::nothrow) double*[n];
 C_U = new(std::nothrow) double*[n];

 for(int i = 0; i < n; i++)
 {
  C_L[i] = new(std::nothrow) double[n];
  C_U[i] = new(std::nothrow) double[n];

 }

 for(int i = 0; i < n; i++)
 {
  for(int j = 0; j < n; j++)
  {
   C_L[i][j] = 0;
   C_U[i][j] = 0;
  }
 }

 double temp;
 
 for(int i = 0; i < n; i++)
 {
  C_U[i][i] = 1;
  C_L[i][0] = Matrix_data[i][0];
 }

 for(int i = 1; i < n; i++)
 {
  C_U[0][i] = Matrix_data[0][i] / C_L[0][0];
 }

 for(int k = 1; k < n; k++)
 {
  for(int i = k; i < n; i++)
  {
   temp = 0;
   for(int r = 0; r < k; r++)
   {
    temp = temp + C_L[i][r] * C_U[r][k];
   }
   C_L[i][k] = Matrix_data[i][k] - temp;
  }

  for(int j = k+1; j < n; j++)
  {
   temp = 0;
   for(int r = 0; r < k; r++)
   {
    temp = temp + C_L[k][r] * C_U[r][j];
   }
   C_U[k][j] = (Matrix_data[k][j] - temp) / C_L[k][k];
  }

 }

 ofstream fout("fout.txt",ios::out);
 fout << "Crout分解：" << endl;
 fout << "Crout分解的L=" << endl;
 for(int i = 0; i < n; i++)
 {
  for(int j = 0; j < n; j++)
  {
   fout << setw(10) << setiosflags(ios::fixed) <<setprecision(4) <<C_L[i][j];
   
  }
  fout << endl;
 }

 fout << endl;
 fout << "Crout分解的U=" << endl;
 for(int i = 0; i < n; i++)
 {
  for(int j = 0; j < n; j++)
  {
   fout << setw(10) << setiosflags(ios::fixed) <<setprecision(4) <<C_U[i][j];
  }
  fout << endl;
 }

 fout.close();
}

void Matrix::Doolittle()
{
 double **D_L;
 double **D_U;
 D_L = new(std::nothrow) double*[n];
 D_U = new(std::nothrow) double*[n];

 for(int i =0; i < n; i++)
 {
  D_L[i] = new(std::nothrow) double[n];
  D_U[i] = new(std::nothrow) double[n];
 }

 for(int i = 0; i < n; i++)
 {
  for(int j = 0; j < n; j++)
  {
   D_L[i][j] = 0;
   D_U[i][j] = 0;
  }
 }

 double temp;
 for(int i = 0; i < n; i++)
 {
  D_U[0][i] = Matrix_data[0][i];
  D_L[i][0] = Matrix_data[i][0] / D_U[0][0];
 }

 
 for(int k = 1; k < n; k++)
 {
   D_L[k][k] = 1;

  for(int j = k; j < n; j++)
  {
   temp = 0;
   for(int m = 0; m < k; m++)
   {
    temp = temp + D_L[k][m] * D_U[m][j];
   }
   D_U[k][j] = Matrix_data[k][j] - temp;

  }

  for(int i = k+1; i < n; i++)
  {
   temp = 0;
   for(int m = 0; m < k; m++)
   {
    temp = temp + D_L[i][m] * D_U[m][k];
   }
   D_L[i][k] = (Matrix_data[i][k] - temp) / D_U[k][k];
  }
 }

  ofstream fout("fout.txt",ios::out|ios::app);
  fout << endl;
  fout << endl;
  fout << "Doolittle分解：" << endl;
  fout << "Doolittle分解的L=" << endl;
  for(int i = 0; i < n; i++)
  {
   for(int j = 0; j < n; j++)
   {
    fout << setw(10) << setiosflags(ios::fixed) <<setprecision(4) <<D_L[i][j];
   }
   fout << endl;

  }
  fout << endl;
  fout << "Doolittle分解的U=" << endl;
  for(int i = 0; i < n; i++)
  {
   for(int j = 0; j < n; j++)
   {
    fout << setw(10) << setiosflags(ios::fixed) <<setprecision(4) <<D_U[i][j];
   }
   fout << endl;
  }
 
 
 fout.close();

}

Matrix::~Matrix()
{
 for(int i = 0; i < n; i++)
 {
  delete []Matrix_data[i];
 
 }
 delete[]Matrix_data;
}

int main()
{
 Matrix A;
 
 A.Crout();
 A.Doolittle();

 return 0;
}