求解线性方程组之LU分解

//具有列选主元和从前到后代入的LU分解
#include <iostream>
#include <math.h>
#include <process.h>

using namespace std;

class lu
{
private:
 int i, j, k, n, *pivrow;
 double eps, pivot, sum, a1, a2, *b, **a;

public:
 void lu_input();
 void lu_decomposition();
 void lu_output();
 ~lu()
 {
  delete[] b;
  delete[] pivrow;
  for (i = 0; i < n; i++)
  {
   delete[] a[i];
  }
  delete[] a;
 }
};

void main()
{
 lu solution;
 solution.lu_input();
 solution.lu_decomposition();
 solution.lu_output();
}

void lu::lu_input()
{
 cout << "\n输入方程个数：";
 cin >> n;
 b = new double[n];
 pivrow = new int[n];
 a = new double*[n];
 for (i = 0; i < n; i++)
 {
  a[i] = new double[n];
 }
 for (i = 0; i < n; i++)
  for (j = 0; j < n; j++)
  {
   cout << "\n输入a[" << i << "][" << j << "] = ";
   cin >> a[i][j];
  }
 for (i = 0; i < n; i++)
 {
  cout << "\n输入b[" << i << "] = ";
  cin >> b[i];
 }
 cout << "\n输入最小主元素：";
 cin >> eps;
}

void lu::lu_decomposition()
{
//LU分解
 for (j = 0; j < n; j++)
 {
//计算上三角矩阵U
  if (j > 0)
  {
   for (i = 0; i < j; i++)
   {
    if (i > 0)
    {
     sum = 0.0;
     for (k = 0; k <= (i-1); k++)
     {
      sum += a[i][k]*a[k][j];
     }
    }
    else
    {
     sum = 0.0;
    }
    a[i][j] = a[i][j] - sum;
   }
  }
  pivot = 0.0;
//计算下三角矩阵L
  for (i = j; i < n; i++)
  {
   if (j > 0)
   {
    sum = 0.0;
    for (k = 0; k <= (j-1); k++)
    {
     sum += a[i][k] * a[k][j];
    }
   }
   else
   {
    sum = 0.0;
   }
   a[i][j] = a[i][j] - sum;
   if (fabs(a[i][j]) > pivot)
   {
    pivrow[j] = i;
    pivot = fabs(a[i][j]);
   }
  }
  if (pivot < eps)
  {
   cout << "\n矩阵可能奇异，失败..." << endl;
   exit(0);
  }
//列选主元
  if (pivrow[j] != j)
  {
   for (k = 0; k < n; k++)
   {
    a1 = a[pivrow[j]][k];
    a[pivrow[j]][k] = a[j][k];
    a[j][k] = a1;
   }
  }
//计算下三角矩阵L
  if (j != (n-1))
  {
   for (i = (j+1); i < n; i++)
   {
    a[i][j] /= a[j][j];
   }
  }
 }
//从前到后代入，先计算LY = B
 for (i = 0; i < n; i++)
 {
  sum = b[pivrow[i]];
  b[pivrow[i]] = b[i];
  if (i > 0)
  {
   for (j = 0; j < (i-1); j++)
   {
    sum -= a[i][j]*b[j];
   }
  }
  b[i] = sum;
 }
//再计算UX = Y
 for (i = (n-1); i >= 0; i--)
 {
  sum = b[i];
  if (i < (n-1))
  {
   for (j = (i+1); j < n; j++)
   {
    sum -= a[i][j]*b[j];
   }
  }
  b[i] = sum / a[i][i];
 }
}

void lu::lu_output()
{
 cout << "\n结果是：" << endl;
 for (i = 0; i < n; i++)
 {
  cout << "\nx[" << i << "] = " << b[i] << endl;
 }
}