基于 MATRIX 类的矩阵分解和方程组求解

前一篇日志主要是一份 MATRIX 的类说明书，经过扩展之后现在可以可以对矩阵进行几种常见的分解运算，可以用来求解线性方程组。

例程1：

 

#include "MATRIX.h"

using namespace std;
int main(int argc,char* argv[])
{

	double a[]={16,4,8,4,4,10,8,4,8,8,12,10,4,4,10,12};
	MATRIX A(4,4,a);
	MATRIX P(4,4);//elementary transformal matrix
	MATRIX B=LUDecomposition(A);
	MATRIX C=Column_PivotLU(A,P);
	MATRIX D=Cholesky(A);//Cholesky Decomposition
	MATRIX L,U;

	cout<<"LU Decomposition:"<<endl;
	U=B.Upper();//get upper triangular matrix
	L=B.LowerI();//get lower triangular matrix
	U.PrintMatrix(7,4,true);//width=7,precision=4,fixed point=true
	cout<<endl;
	L.PrintMatrix(7,4,true);
	cout<<endl;

	cout<<"LU Decomposition under coloumn-pivot rule:"<<endl;
	U=C.Upper();
	L=C.LowerI();
	U.PrintMatrix(7,4,true);
	cout<<endl;
	L.PrintMatrix(7,4,true);
	cout<<endl;

	cout<<"Cholesky Decomposition:"<<endl;
	L=D.Lower();
	U=L;
	U.trans().PrintMatrix(7,4,true);
	cout<<endl;
	L.PrintMatrix(7,4,true);
	cout<<endl;

	double barray[]={3,2,0,5};
	MATRIX x,y;
	MATRIX b(4,1,barray);
	y=SolveL(L,b);
	y.PrintMatrix(4);
	cout<<endl;
	x=SolveU(U,y);
	x.PrintMatrix(4);

	return 0;
}

结果如下：

LU Decomposition:
16.0000  4.0000  8.0000  4.0000
 0.0000  9.0000  6.0000  3.0000
 0.0000  0.0000  4.0000  6.0000
 0.0000  0.0000  0.0000  1.0000

 1.0000  0.0000  0.0000  0.0000
 0.2500  1.0000  0.0000  0.0000
 0.5000  0.6667  1.0000  0.0000
 0.2500  0.3333  1.5000  1.0000

LU Decomposition under coloumn-pivot rule:
16.0000  4.0000  8.0000  4.0000
 0.0000  9.0000  6.0000  3.0000
 0.0000  0.0000  6.0000 10.0000
 0.0000  0.0000  0.0000 -0.6667

 1.0000  0.0000  0.0000  0.0000
 0.2500  1.0000  0.0000  0.0000
 0.2500  0.3333  1.0000  0.0000
 0.5000  0.6667  0.6667  1.0000

Cholesky Decomposition:
 4.0000  1.0000  2.0000  1.0000
 0.0000  3.0000  2.0000  1.0000
 0.0000  0.0000  2.0000  3.0000
 0.0000  0.0000  0.0000  1.0000

 4.0000  0.0000  0.0000  0.0000
 1.0000  3.0000  0.0000  0.0000
 2.0000  2.0000  2.0000  0.0000
 1.0000  1.0000  3.0000  1.0000

0.75
0.42
-1.17
7.33

2.79
5.42
-11.58
7.33

Process returned 0 (0x0)   execution time : 0.214 s
Press any key to continue.

 

例程2：

#include "MATRIX.h"

using namespace std;
int main(int argc,char* argv[])
{
	double aArray[]={2,-1,3,4,2,5,2,1,2};
	double bArray[]={1,4,5};
	MATRIX A(3,3,aArray);
	MATRIX L,U,y,x,b(3,1,bArray);

	cout<<"A:"<<endl;
	A.PrintMatrix(3);
	cout<<endl<<"b:"<<endl;
	b.PrintMatrix(3);
	L=LUDecomposition(A);
	U=L.Upper();
	L=L.LowerI();//When using Cholesky method,L=L.Lower()

	y=SolveL(L,b);
	x=SolveU(U,y);

	cout<<endl<<"The answer is:"<<endl;
	x.PrintMatrix(3);

	return 0;
}

 

结果如下：

A:
2.00 -1.00 3.00
4.00 2.00 5.00
2.00 1.00 2.00

b:
1.00
4.00
5.00

The answer is:
9.00
-1.00
-6.00

Process returned 0 (0x0)   execution time : 0.188 s
Press any key to continue.

 

新的 MATRIX 类已经上传到资源，免费开源，欢迎修改完善。