ldlt分解.cpp

本人的计算方法课程作业

//ldlt.cpp--矩阵的LDLT分解
#include<iostream>
#include<cmath>
using namespace std;
const int ROW=20;
void print_array(int x[ROW][ROW]);
void print_d(int d[ROW]);
void main()
{
	int a[ROW][ROW],l[ROW][ROW],d[ROW],d1[ROW][ROW];
	int i,j,k,p;
    for(i=1;i<=ROW;i++)          //求A矩阵的各元素值
		for(j=1;j<=ROW;j++)
		{
			if(j<=i) a[i][j]=j;
            else a[i][j]=i;
		}
    cout<<"在进行分解之前，先输出矩阵A:"<<endl;
	print_array(a);              //输出A矩阵
	for(i=1;i<=ROW;i++)
		l[i][i]=1;               //L矩阵对角线上的元素均为1
		d[1]=a[1][1];            //D矩阵对角线上的D[1][1]为1
	for(k=1;k<=ROW;k++)
	{
		{   
			int t1=0;            
		    for(p=1;p<k;p++)
	           t1+=l[k][p]*d[p]*l[k][p];
	        d[k]=a[k][k]-t1;     //求D矩阵对角线上的其他元素值
		}
	    if(d[k]==0) cout<<"Failure!"<<endl;    //若D阵对角线上的为零，输出错误信息
		for(i=k+1;i<=ROW;i++)
		{
			int t2=0;
			for(p=1;p<k;p++)
				t2+=l[i][p]*d[p]*l[k][p];
			l[i][k]=(a[i][k]-t2)/d[k];    //求L矩阵下三角部分的元素值
		}
		for(i=1;i<k;i++)
			l[i][k]=0;                    //L矩阵上三角部分的元素值均为0
	}
	cout<<"通过分解，可得矩阵L为:"<<endl;
	print_array(l);                       //输出L矩阵
    print_d(d);                           //输出D矩阵对角线上的元素
	cout<<"在得到对角线上的元素后，可得矩阵D:"<<endl;
	for(i=1;i<=ROW;i++)
		for(j=1;j<=ROW;j++)               //求D矩阵其他位置的元素值
		{
			if(j==i) d1[i][j]=d[i];
		    else d1[i][j]=0;
		}
		print_array(d1);                  //输出D矩阵
}

//print_array()用于输出矩阵
void print_array(int x[ROW][ROW])
{
	int i,j;
	for(i=1;i<=ROW;i++)
	{
		for(j=1;j<=ROW;j++)
			cout<<x[i][j]<<" ";
		cout<<endl;
	}
}

//print_d()输出矩阵对角线上的元素
void print_d(int d[ROW])
{
	int i;
	cout<<"矩阵D对角线上的元素为:"<<endl;
	for(i=1;i<=ROW;i++)
	cout<<d[i]<<"  ";
	cout<<endl;
}