/**/ /////////////////////////////////////// /
//                                     //
//                                    //
//数值分析  直接三角分解法求方程组    //
//                                    //
//  在VC2005下编译通过，由于VC6.0的   //
//   不标准会有错误                   //
//                                    //
//   2006.5.30   v0.1                 //
/**/////////////////////////////////////////

#include "iostream"
using namespace std;

class Matrix
{
private:
 double** A;      //矩阵A
 double *b;       //向量b
public:
 int size;
 Matrix(int );
 ~Matrix();
friend double* Dooli(Matrix& );
 void Input();
 void Disp();

};

Matrix::Matrix(int x)
{
 size=x;

 //为向量b分配空间并初始化为0
 b=new double [x];
 for(int j=0;j<x;j++)
  b[j]=0;

 //为向量A分配空间并初始化为0
 A=new double* [x];
 for(int i=0;i<x;i++)
  A[i]=new double [x];
 for(int m=0;m<x;m++)
  for(int n=0;n<x;n++)
   A[m][n]=0;
}

Matrix::~Matrix()
{
    cout<<"正在析构中~~~~"<<endl;
    delete b;
    for(int i=0;i<size;i++)
        delete A[i];
    delete A;
}

void Matrix::Disp()
{
 for(int i=0;i<size;i++)
 {
  for(int j=0;j<size;j++)
   cout<<A[i][j]<<"  ";
  cout<<endl;
 }
}

void Matrix::Input()
{
 cout<<"请输入A:"<<endl;
 for(int i=0;i<size;i++)
  for(int j=0;j<size;j++){
   cout<<"第"<<i+1<<"行"<<"第"<<j+1<<"列:"<<endl;
  cin>>A[i][j];
  }
   cout<<"请输入b:"<<endl;
 for(int j=0;j<size;j++){
  cout<<"第"<<j+1<<"个:"<<endl;
   cin>>b[j];
 }
 
}
 
double* Dooli(Matrix& A)
{
 double *Xn=new double [A.size];
 Matrix L(A.size),U(A.size);

 //分别求得U，L的第一行与第一列
   for(int i=0;i<A.size;i++)
      U.A[0][i]=A.A[0][i];
   for(int j=1;j<A.size;j++)
      L.A[j][0]=A.A[j][0]/U.A[0][0];

//分别求得U，L的第r行，第r列
     double temp1=0,temp2=0;
 for(int r=1;r<A.size;r++){
     //U
     for(int i=r;i<A.size;i++){
         for(int k=0;k<r-1;k++)
            temp1=temp1+L.A[r][k]*U.A[k][i]; 
            U.A[r][i]=A.A[r][i]-temp1;
     }
     //L
     for(int i=r+1;i<A.size;i++){
          for(int k=0;k<r-1;k++)
            temp2=temp2+L.A[i][k]*U.A[k][r];
             L.A[i][r]=(A.A[i][r]-temp2)/U.A[r][r];
     }
 }
 cout<<"计算U得:"<<endl;
 U.Disp();
 cout<<"计算L的:"<<endl;
 L.Disp();
 
 double *Y=new double [A.size];

 Y[0]=A.b[0];
 for(int i=1;i<A.size;i++ ){
     double temp3=0;
     for(int k=0;k<i-1;k++)
         temp3=temp3+L.A[i][k]*Y[k];
     Y[i]=A.b[i]-temp3;
 }

 Xn[A.size-1]=Y[A.size-1]/U.A[A.size-1][A.size-1];
 for(int i=A.size-1;i>=0;i--){
     double temp4=0;
     for(int k=i+1;k<A.size;k++)
         temp4=temp4+U.A[i][k]*Xn[k];
     Xn[i]=(Y[i]-temp4)/U.A[i][i];
 }
 
return Xn;
}
 
int main()
{
 Matrix B(4);
 B.Input();
 double *X;
 X=Dooli(B);
 cout<<"~~~~解得:"<<endl;
 for(int i=0;i<B.size;i++)
     cout<<"X["<<i<<"]:"<<X[i]<<" ";
 cout<<endl<<"呵呵呵呵呵";
 return 0;
}

试验总结：将高斯消去法改写为紧凑形式，可以直接从矩阵 A 的元素的导计算 L ， U 元素的递推公式，而不需任何中间步骤，一旦实现了矩阵 A 的 U ， L 分解那么就等价于求解两个三角形方程组。

  注意： 编成语言中的数组以 ’0’ 为首元素，数组的一位偏移最容易出错；

             注意变量的作用域；