sanjiao.c

矩阵直接三角分解法的c++语言实现

#include<iostream.h>

class GS2
{
public:
  GS2 (int i,int j,double e)               //构造矩阵
  {
   x=i;
   y=j;
         eps=e;
   a=new double[x*(x+y)];
  }
  void build();
  void return_x();                   //回代求 x[i]
  bool judge (int);              //判断是否|P|<=EPS
  bool Doolittle ();              //doolittle直接三角分解法的主体程序
  void changcemain (int);       //选主元
  void show();                  //打印x
  void print();
  ~GS2()
  {
   delete []a;
  }
protected:
 double *a,eps;
    int x,y;
};

void GS2::build()
{
 for(int i=0;i<x;i++)
  for(int j=0;j<x+y;j++)
  {
   cout<<"x["<<i+1<<"]["<<j+1<<"]=";

   cin>>a[j+(x+y)*i];
  }
}
//
void GS2::print()
{
 for(int i=0;i<x;i++)
  for(int j=0;j<x+y;j++)
  {
   cout<<"x["<<i+1<<"]["<<j+1<<"]="

   <<a[j+(x+y)*i]<<endl;
  }
}
//

bool GS2::Doolittle()
{
 double temp;
    

 for(int r=0;r<x-1;r++)
 {  
 
  if(r==0)
  {
   changcemain(r);
        if (!judge(r))
               return false;
        else
     for(int i=r+1;i<x;i++)
       a[r+(x+y)*i]/=a[r+(x+y)*r];
  }  
        else
  {
      for(int i=r;i<x;i++)
   {
    temp=0;
     for(int k=0;k<r;k++)
      temp+=a[k+(x+y)*i]*a[r+(x+y)*k];
             
     a[r+(x+y)*i]-=temp;
   }
  
   changcemain(r);
   
   if (!judge(r))
           return false;

   for( i=r+1;i<x;i++)
   {
    a[r+(x+y)*i]/=a[r+(x+y)*r];
   }
   for(i=r+1;i<x+y;i++)
   {
    temp=0;
    for(int k=0;k<r;k++)
    {
     temp+=a[k+(x+y)*r]*a[i+(x+y)*k];
    }
         a[i+(x+y)*r]-=temp;
   }
  }  //end else

 }

 for(int i=x-1;i<x+y;i++)
 {
  temp=0;
  for(int k=0;k<x-1;k++)
  {
   temp+=a[k+(x+y)*(x-1)]*a[i+(x+y)*k];
  }
  a[i+(x+y)*(x-1)]-=temp;
 }
     
    if (!judge(x-1))
           return false;
    
 return true;
}

void GS2::return_x()
{
 double temp;
 for(int j=0;j<y;j++)
 {

    for(int i=x-1;i>=0;i--)
 {
  temp=0;
  for(int k=x-1;k>i;k--)
  {
          temp+=a[k+(x+y)*i]*a[(x+j)+(x+y)*k];
  }
    a[(x+j)+(x+y)*i]-=temp;
 
 a[(x+j)+(x+y)*i]/=a[i+(x+y)*i];
 }
 }
}

void GS2::show()
{
  for(int j=0;j<y;j++)
  {
   cout<<endl;
   for(int i=0;i<x;i++)
   {
    cout<<"x["<<j+1<<"]["<<i+1<<"]="<<a[(x+j)+(x+y)*i]<<'\t';
   }
  }
}

void GS2::changcemain(int r)
{
 double temp=a[r+(x+y)*r];
    int t=r;
 for(int i=r+1;i<x;i++)
 {
  if( a[r+(x+y)*i] > temp )
  {  
   temp=a[r+(x+y)*i];
       t=i;
  }
 }
 if (t!=r)
  for(i=0;i<x+y;i++)
  {
   temp=a[i+(x+y)*r];
   a[i+(x+y)*r]=a[i+(x+y)*t];
   a[i+(x+y)*t]=temp;
  }
}

bool GS2::judge(int r)
{ 
 if (a[r+(x+y)*r]>0)
   return ((a[r+(x+y)*r]>eps )?true:false);
 else
      return (((-a[r+(x+y)*r])>eps )?true:false);
}
   
int main()
{
 int n,m;
 double e;
 cout<<"please input n , m and EPS :";
 cin>>n>>m>>e;

 GS2 gs2(n,m,e);
 gs2.build();
   
 if(!gs2.Doolittle())
 {
  cout<<"此方程无解!"<<endl;
  return 1;
 }
 gs2.print();
 gs2.return_x();
 gs2.show();
 return 0;
}