wokr3.cpp

系数矩阵为对称正定的线性方程组的求解程序(LDL~分解法)

/**************************************************************************************/
                                         
                                                                                     
                 //编制系数矩阵为对称正定的线性方程组的求解程序(LDL~分解法)     
                                                                                      
/***************************************************************************************/
#include<iostream.h>

/*由于本程序出现矩阵或者是对称矩阵，或者是只用到矩阵的下三角元素，故本程序全部用一位数组存储矩阵*/

int withkeyboard_A(int nn,double *A)//从键盘获得系数矩阵A，按照a[i][j]=A[i*(i-1)/2+j-1]的对照关系
                                    //存入数组A中，A为nn=n*(n-1)/2+n-1数组
{    	
	cout<<"输入A矩阵下三角元素:"<<endl;
	for (int a=0;a<nn;a++)
		{		 		
			cin>>A[a];
		}    	
	return 0;
}
int withkeyboard_D(int n,double *D)//从键盘获得右端项d，存入数组D中，D为n数组
{   
	cout<<"输入D矩阵:"<<endl;		
	for (int b=0;b<n;b++)
		{		 	
			cin>>D[b];
		}
	return 0;
}
void LdLT(double *A,int n,int nn,double *xD,double *L)//由系数矩阵A计算LDL~求解法中的
                                                    //L矩阵和D矩阵（为区别于右端项，D矩阵记为*xD）
{
   double *AA=new double[nn];//AA矩阵为算法中的a~[i][j],按照a~[i][j]=AA[i*(i-1)/2+j-1]的对照关系
                             //存入数组AA中，AA为nn=n*(n-1)/2+n-1数组
   
   for (int c=1;c<=n;c++)
   {
	   int d=c*(c-1)/2+c-1;//给下三角矩阵L赋值1
	   L[d]=1;	   
   }
    
   double dof=A[0];        
   xD[0]=dof;             //给D矩阵第一个元素赋值         
   double res1=0;
   double res2=0;
   for(int j=2;j<=n;j++)
   {
	   for(int k=1;k<=j-1;k++)
	   {          
		  res1=0;///恶心死我了，此处注意重新初始化，在循环进入之前。否则会有莫名其妙的结果
		            ///结果老是不对，还以为是下标想糊涂了，改下标改了数小时之后不见效才发现
		            //是此处的小错误导致结果的错误。		            
		  for(int m=1;m<=k-1;m++)
		  {			  
			   res1=res1+AA[j*(j-1)/2+m-1]*L[k*(k-1)/2+m-1];
		  }
		  int jk=j*(j-1)/2+k-1;          
		  AA[jk]=A[jk]-res1;   //               //迭代式1
		  L[jk]=AA[jk]/xD[k-1];//               //迭代式2
          	  
	   }
	   res2=0;      ///同上一个res数，进入for以前要初始化
	   for(int l=1;l<=j-1;l++)
		{			  
		res2=res2+AA[(j-1)*j/2+l-1]*L[(j-1)*j/2+l-1];
		}	
	   xD[j-1]=A[j*(j-1)/2+j-1]-res2;//         //迭代式3
	                        ////运行三个迭代式，求解xD和L矩阵的各个元素
   }
   
   cout<<"下三角矩阵L元素为:  ";//输出由系数矩阵A计算LDL~求解法中的L矩阵和D矩阵*xD
   for(int q2=0;q2<nn;q2++)
	{
	   cout<<L[q2]<<"   ";
	}
	cout<<endl;
	cout<<"LDLT中的D矩阵元素为:  ";
	for(int q3=0;q3<n;q3++)
	{
	   cout<<xD[q3]<<"   ";
	}
	cout<<endl;
}//获得xD和L矩阵
void getYfromLDLT_equto_D(double *L,double *D,int n,double *Y)//由已知的右端项D和求得的L矩阵
                                                              //以及式子LY=D利用迭代法求解矩阵Y
{
  
  Y[0]=D[0]/L[0];
  double  res3=0;
  for (int i=1;i<n;i++)
  {
	  res3=0;///注意事项同以上res，初值化
	  for (int j=0;j<i;j++)
	  {		  
		  res3=res3+Y[j]*L[i*(i+1)/2+j];
	  }
	  Y[i]=(D[i]-res3)/L[i*(i+1)/2+i];///利用迭代法求解矩阵Y
  }
  
  cout<<"Y矩阵元素为:  ";//输出由已知的右端项D和求得的L矩阵以及式子LY=D求得的矩阵Y
  for(int a=0;a<n;a++)////////////////
   {
	   cout<<Y[a]<<"  ";
   }
  cout<<endl;
}
void result(double *Y,double *xD,double *L,int n,double *resul)
                       //由以上函数获得的xD，Y以及L,求解方程
					   //L~X=Z=Y/xD;求得的结果存放于数组*resul中，并输出结果
{
	double *Z=new double[n];

	for (int i=0;i<n;i++)
	{
		Z[i]=Y[i]/xD[i];
	}

    cout<<"Z矩阵元素为:  ";
    for(int a=0;a<n;a++)////////////////
    {
	   cout<<Z[a]<<" ";
    }
    cout<<endl;

	resul[n-1]=Z[n-1]/L[n*(n-1)/2+n-1];
	double res4=0;
    for (int ii=n-1;ii>=0;ii--)
	{
	  res4=0;
	  for (int j=n-1;j>=ii;j--)
	  {		  
		  res4=res4+resul[j]*L[j*(j+1)/2+ii-1];
	  }
	  resul[ii-1]=(Z[ii-1]-res4)/L[ii*(ii-1)/2+ii-1];//利用迭代法求解矩阵*resul
  }
	
	cout<<"最终解为:  ";
    for (int p=0;p<n;p++)/////////////////////
	{		
		cout<<resul[p]<<"  ";
	}
    cout<<endl;
}
int main()
{
	int n,m=1;	
	cout<<"输入A对称矩阵行列数n: ";
	cin>>n;
	int nn=n*(n+1)/2;
	double *A=new double[nn];
	double *D=new double[n];
	double *xD=new double[n];
    double *L=new double[nn];
	double *resul=new double[n];
	double *Y=new double[n];

    withkeyboard_A(nn,A);
	LdLT(A,n,nn,xD,L);
    while(m==1)// 是否继续计算另一组右端项控制
	{
	withkeyboard_D(n,D);	
	getYfromLDLT_equto_D(L,D,n,Y);
    result(Y,xD,L,n,resul);
    cout<<"是否计算继续另一组右端项？:(ture=1 false=0)";
	cin>>m;
	}

	/*delete []A;
	delete []D;
	delete []xD;
	delete []L;
	delete []resul;
	delete []Y;*/

	return 0;
}