#include <iostream.h>
#include <iomanip.h>
/*该程序实现将一个对称矩阵A分解成LDU
其中L为下三角矩阵;D为对角矩阵;U为上三角矩阵
该程序使用了二级指针,矩阵A的阶数可以由运行人员决定*/
class matrix
{
public:
matrix(); //初始化矩阵
~matrix();
void input(); //输入矩阵的元素
void output(); //输出矩阵的元素
void eye(); //把矩阵化为单位矩阵
void calldu(matrix &k); //将对称矩阵分解
private:
int m,n; //矩阵的行和列
double **p_matrix; //指向指针的指针,可以动态分配内存空间
};
void matrix::calldu (matrix &aa)
{
matrix dd,ll,uu;
dd.eye ();
ll.eye ();
uu.eye ();
for(int j=0; j<n; j++)
{
for(int i=0; i<j; i++)
{
uu.p_matrix [i][j]=aa.p_matrix [i][j];
for(int k=0; k<i; k++)
uu.p_matrix [i][j]=uu.p_matrix [i][j] -
dd.p_matrix [k][k] * ll.p_matrix [i][k] * uu.p_matrix [k][j];
uu.p_matrix [i][j] = uu.p_matrix [i][j] / dd.p_matrix [i][i];
}
dd.p_matrix [j][j] = aa.p_matrix [j][j];
for(int k=0; k<j; k++)
dd.p_matrix [j][j]=dd.p_matrix [j][j] -
dd.p_matrix [k][k] * ll.p_matrix [j][k] * uu.p_matrix [k][j];
for(int t=j+1; t<n; t++)
{
ll.p_matrix [t][j] = aa.p_matrix[t][j];
for(int k=0; k<j; k++)
ll.p_matrix [t][j] = ll.p_matrix [t][j] -
dd.p_matrix [k][k] * ll.p_matrix [t][k] * uu.p_matrix [k][j];
ll.p_matrix [t][j] = ll.p_matrix [t][j] / dd.p_matrix [j][j];
}
}
cout<<"the ll matrix is : \n";
ll.output ();
cout<<"the dd matrix is : \n";
dd.output ();
cout<<"the uu matrix is : \n";
uu.output ();
}
void matrix::eye ()
{
for(int i=0; i<m; i++)
for(int j=0; j<n; j++)
if(i==j)
p_matrix[i][j]=1;
else
p_matrix[i][j]=0;
}
matrix::matrix()
{
cout<<"input the col and row: \n";
cin>>m>>n;
p_matrix=new double*[m];
for(int i=0; i<m; i++)
p_matrix[i]=new double[n];
}
matrix::~matrix()
{
for(int i = 0; i < m; i++)
delete [] p_matrix[i];
delete p_matrix;
}
void matrix::input()
{
cout<<"input the element: \n";
for(int i=0; i<m; i++)
for(int j=0; j<n; j++)
cin>>p_matrix[i][j];
}
void matrix::output()
{
cout<<"the matrix is : \n";
for(int i=0; i<m; i++)
{
for(int j=0; j<n; j++)
cout<<setw(10)<<p_matrix[i][j];
cout<<endl;
}
}
void main()
{
matrix aa;
aa.input();
cout<<"the aa matrix is : \n";
aa.output();
cout<<"初始化LDU.\n";
aa.calldu (aa);
}