【算法设计与分析基础】高斯消元

参考《算法导论》第七部分 算法问题选编中的第28章 矩阵运算。

#include<iostream>
using namespace std;

/**
 * 一些矩阵例子
 *//*
double Value[9][9]={
	{0,0,0,-85,0,59,0,0,72},
	{97,-137,122,81,97,55,0,132,0},
	{0,-69,145,70,82,0,66,0,88},
	{101,0,-103,76,0,0,-58,-143,0},
	{0,-100,0,-171,-93,0,89,0,0},
	{0,0,59,154,0,0,0,-50,0},
	{0,111,0,-161,-80,0,85,-173,0},
	{112,0,250,111,-167,208,96,0,116},
	{62,0,86,-111,-90,-228,132,-61,81}
},
b[9]={1,1,1,1,1,1,1,1,1};
*/
double Value[9][9]={
	{0,0,0,-85,0,59,0,0,72},
	{1,-137,122,81,97,55,0,132,0},
	{0,-69,145,70,82,0,66,0,88},
	{101,0,0,76,0,0,-58,-143,0},
	{0,-100,0,-171,-93,0,89,0,0},
	{0,0,59,154,0,0,0,-50,0},
	{0,111,0,-161,-80,0,85,0,0},
	{112,0,250,0,-167,208,96,0,116},
	{62,0,86,-111,0,0,132,-61,81}
},
b[9]={1,1,1,1,1,1,1,1,1};


double Example[4][4]={
	{2,3,1,5},
	{6,13,5,19},
	{2,19,10,23},
	{4,10,11,31},
};
double Examplee[3][3]={
	{1,2,0},
	{3,4,4},
	{5,6,3},
},
Ex[3]={3,7,8};

double Exampleee[4][4]={
	{2,0,2,0.6},
	{3,3,4,-2},
	{5,5,4,2},
	{-1,-2,3.4,-1},
},
Ex1[4]={1,1,1,1};

class Gauss{
private:
	double **A;//原矩阵
	double *x;//线性方程组的解
	double *y;//采用正向替换，对Ly=Pb求解y
	double *b;//线性方程组最右边一列的常数
	double **L;//LU分解的L数组
	double **U;//LU分解的U数组
	double **l;//L的逆矩阵
	double **u;//U的逆矩阵
	double **a;//A的逆矩阵 = l * u
	double **p;//根据数组P求出置换矩阵
	int *P;//行数排列
	int num;
public:
	Gauss(int n);//构造函数，构造一个 n x n 的矩阵
	~Gauss();//析构函数
	void setA(int x,int y,double v);//给数组A的[x][y]赋值为v
	void setB(int x,double v);//给数组B的[x]赋值为v
	
	void Add(int r1,int r2);//第r1行加上r2行
	void LU();// 计算 L 和 U 矩阵
	void LUP();// 计算出 L 和 U 矩阵 以及 排列行数P
	void solve_p();//计算 置换矩阵p 并 修改原始矩阵的顺序
	void LUP_SOLVE();//通过 L矩阵 和 U矩阵 解决线性方程组问题
	void U_();//计算 U 的逆矩阵
	void L_();//计算 L 的逆矩阵
	void A_();//计算 A 的逆矩阵
	void A_Next();//计算 原始矩阵 的逆矩阵 即 A的逆矩阵乘以置换矩阵

	void displayB();//输出 线性方程组 最右边一列常数
	void displayL();//输出 L 矩阵
	void displayU();//输出 U 矩阵
	void displayDiagonal();//输出对角线
	void displayx();//输出线性方程组的解
	void displayA();//输出原矩阵
	void displayy();//输出 y 矩阵
	void displayu();//输出 U 的逆矩阵
	void displayl();//输出 L 的逆矩阵
	void displaya();//输出 A 的逆矩阵
	void displayP();//输出 排列行数P
	void displayp();//输出 置换矩阵p
};

    Gauss::Gauss(int n){//构造函数，构造一个 n x n 的矩阵
		int i;
		num=n;
		A=new double*[n];
		for(i=-1;++i<n;){
			A[i]=new double[n];
			memset(A[i],0,n*sizeof(A[i]));
		}
		x=new double[n];
		b=new double[n];
	}
	Gauss::~Gauss(){//析构函数
	/*	int i;
		for(i=-1;++i<num;){
			delete[]A[i];
			delete[]L[i];
			delete[]U[i];
		}
		delete[]y;
		delete[]A;
		delete[]x;
		delete[]b;
		delete[]U;*/
	}
	void Gauss::setA(int x,int y,double v){//给数组A的[x][y]赋值为v
		if( (x<0 || x>=num) || (y<0 || y>=num)  )return;
		A[x][y]=v;
	}
	void Gauss::setB(int x,double v){//给数组B的[x]赋值为v
		b[x]=v;
	}

	void Gauss::Add(int r1,int r2){//r1行加上r2行
		int i=-1;
		for(;++i<num;){
			A[r1][i]+=A[r2][i];
		}//for
		b[r1]+=b[r2];
	}

	void Gauss::LU(){// 计算 L 和 U 矩阵
		//让L 和 U 成为 n x n 的矩阵
		L=new double*[num];
		U=new double*[num];
		int i=-1,j,k;
		for(;++i<num;){
			L[i]=new double[num];
			U[i]=new double[num];
			//L 和 U 全部初始化为0
			memset(L[i],0,num*sizeof(double));
			memset(U[i],0,num*sizeof(double));
			//L的对角线全是1
			L[i][i]=1;
		}
		for(k=-1;++k<num;){
			U[k][k]=A[k][k];
			for(i=k;++i<num;){
				L[i][k]=A[i][k]/U[k][k];
				U[k][i]=A[k][i];
			}
			for(i=k;++i<num;){
				for(j=k;++j<num;){
					A[i][j]=A[i][j]-L[i][k] * U[k][j];
				}
			}
		}	
	}

    void Gauss::LUP(){//计算出 L 和 U 矩阵 以及 排列行数P
		P=new int[num];
		double **_A=new double*[num];
		int i=-1,j,k,p,_k;
		double temp;
		/**
		 * 给P赋值为0 到 num-1
		 */
		for(;++i<num;){
			_A[i]=new double[num];
			P[i]=i;
			for(j=-1;++j<num;){
				_A[i][j]=A[i][j];
			}
		}//for
		for(k=-1;++k<num;){
			p=0;
			for(i=k;i<num;i++){
				if(abs(_A[i][k])>p){
					p=abs(_A[i][k]);
					_k=i;
				}
			}
			if(p==0){
				cout<<"singular matrix"<<endl;
				exit(0);
			}
			temp=P[k];
			P[k]=P[_k];
			P[_k]=temp;

			for(i=-1;++i<num;){
				temp=_A[k][i];
				_A[k][i]=_A[_k][i];
				_A[_k][i]=temp;
			}
			for(i=k;++i<num;){
				_A[i][k]/=_A[k][k];
				for(j=k;++j<num;){
					_A[i][j]-=_A[i][k]*_A[k][j];
				}
			}
		}
	}


	void Gauss::solve_p(){//计算 置换矩阵p 并 修改原始矩阵的顺序
		int i,j;
		double **k=new double*[num];
		p=new double*[num];
		for(i=-1;++i<num;){
			p[i]=new double[num];
			k[i]=new double[num];
			memset(p[i],0,sizeof(double)*num);
			for(j=-1;++j<num;){
				k[i][j]=A[i][j];
			}
		}
		for(i=-1;++i<num;){
			p[i][P[i]]=1;
		}
		for(i=-1;++i<num;){
			for(j=-1;++j<num;){
				A[i][j]=k[P[i]][j];
			}
		}
		
	}

    void Gauss::LUP_SOLVE(){//通过 L矩阵 和 U矩阵 解决线性方程组问题
		x=new double[num];
		y=new double[num];
		memset(x,0,sizeof(double)*num);
		memset(y,0,sizeof(double)*num);
		int i=-1,j;
		double sum;
		for(;++i<num;){
			sum=0;
			for(j=-1;++j<i;){
				sum+=L[i][j]*y[j];
			}
			y[i]=b[i]-sum;
		}
		for(i=num;--i>=0;){
			sum=0;
			for(j=num;--j>i;){
				sum+=U[i][j]*x[j];
			}
			x[i]=(y[i]-sum)/U[i][j];
		}
		
	}

     void Gauss::U_(){//计算 U 的逆矩阵
		int i,j,k;
		double sum;
		u=new double*[num];
		for(i=-1;++i<num;){
			u[i]=new double[num];
			memset(u[i],0,sizeof(double)*num);
		}
		for(i=-1;++i<num;){
			u[i][i]=1/U[i][i];
			for(k=i;--k>=0;){
				sum=0;
				for(j=k;++j<=i;){
					sum+=U[k][j]*u[j][i];
					u[k][i]=-sum/U[k][k];
				}
			}
		}
	}

    void Gauss::L_(){//计算 L 的逆矩阵
		int i,j,k;
		l=new double*[num];
		for(i=-1;++i<num;){
			l[i]=new double[num];
			memset(l[i],0,sizeof(double)*num);
		}
		for(i=-1;++i<num;){
			l[i][i]=1;//对角线依次赋值为1
			for(k=i;++k<num;){
				for(j=i-1;++j<=k-1;){
					l[k][i]=l[k][i]-L[k][j]*l[j][i];
				}
			}
		}
	}

    void Gauss::A_(){//计算 A 的逆矩阵
		int i,j,k;
		double sum;
		a=new double*[num];
		for(i=-1;++i<num;){
			a[i]=new double[num];
			memset(a[i],0,sizeof(double)*num);
		}
		for(i=-1;++i<num;){
			for(j=-1;++j<num;){
				sum=0;
				for(k=-1;++k<num;){
					sum+=u[i][k]*l[k][j];
				}
				a[i][j]=sum;
			}
		}
	}

    void Gauss::A_Next(){//计算 原始矩阵 的逆矩阵 即 A的逆矩阵乘以置换矩阵
		int i,j,k;
		double sum;
		double **aa=new double*[num];
		for(i=-1;++i<num;){
			aa[i]=new double[num];
			for(j=-1;++j<num;){
				aa[i][j]=a[i][j];
			}
		}
		for(i=-1;++i<num;){
			for(j=-1;++j<num;){
				for(j=-1;++j<num;){
					sum=0;
			    	for(k=-1;++k<num;){
					     sum+=aa[i][k]*p[k][j];
					}
				    a[i][j]=sum;
				}
			}
		}
	}

	void Gauss::displayB(){//输出 线性方程组 最右边一列常数
		int i=-1;
		for(;++i<num-1;){
			cout<<b[i]<<' ';
		}
		cout<<b[i]<<endl;
	}

    void Gauss::displayL(){//输出 L 矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<L[i][j]<<' ';
			}
			cout<<L[i][j]<<endl;
		}
		cout<<endl;
	}

	void Gauss::displayU(){//输出 U 矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<U[i][j]<<' ';
			}
			cout<<U[i][j]<<endl;
		}
		cout<<endl;
	}

	void Gauss::displayDiagonal(){//输出对角线
		int i=-1;
		for(;++i<num-1;){
			cout<<A[i][i]<<" ";
		}//for
		cout<<A[i][i]<<endl;
	}

	void Gauss::displayx(){//输出线性方程组的解
		int i=-1;
		for(;++i<num;){
			cout<<'x'<<i+1<<" = "<<x[i]<<endl;
		}//for
	}

	void Gauss::displayA(){//输出原矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<A[i][j]<<',';
			}
			cout<<A[i][j]<<endl;
		}
		cout<<endl;
	}


	void Gauss::displayy(){//输出 y 矩阵
		int i=-1;
		for(;++i<num;){
			cout<<'y'<<i+1<<" = "<<y[i]<<endl;
		}//for
	}
	void Gauss::displayu(){//输出 U 的逆矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<u[i][j]<<' ';
			}
			cout<<u[i][j]<<endl;
		}
		cout<<endl;
	}
	void Gauss::displayl(){//输出 L 的逆矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<l[i][j]<<' ';
			}
			cout<<l[i][j]<<endl;
		}
		cout<<endl;
	}
	void Gauss::displaya(){//输出 A 的逆矩阵
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<a[i][j]<<' ';
			}
			cout<<a[i][j]<<endl;
		}
		cout<<endl;
	}

	void Gauss::displayP(){//输出 排列行数P
		int i=-1;
		for(;++i<num;){
			cout<<'p'<<i+1<<" = "<<P[i]<<endl;
		}//for
	}

	void Gauss::displayp(){//输出 置换矩阵p
		int i=-1,j;
		for(;++i<num;){
			for(j=-1;++j<num-1;){
				cout<<p[i][j]<<' ';
			}
			cout<<p[i][j]<<endl;
		}
		cout<<endl;
	}

	
int main(){
	Gauss G(9);
	int i=-1,j;
	/**
	 * 初始化线性方程组的矩阵
	 */
	for(;++i<9;){
		G.setB(i,b[i]);
		for(j=-1;++j<9;){
			G.setA(i,j,Value[i][j]);
		}
	}
/*	Gauss G(3);
	int i=-1,j;
	for(;++i<3;){
		G.setB(i,Ex[i]);
		for(j=-1;++j<3;){
			G.setA(i,j,Examplee[i][j]);
		}
	}*/
/*	Gauss G(4);
	int i=-1,j;
	for(;++i<4;){
		G.setB(i,Ex1[i]);
		for(j=-1;++j<4;){
			G.setA(i,j,Exampleee[i][j]);
		}
	}*/
	/**
	 * 一旦主元选取为0，就会有除以0的问题，主元处于矩阵U的对角线上。
	 *
	 * 由于调用LU分解的时候有 除以0的现象
	 * 所以把矩阵中为 0 的数字 去除
	 * 方法是 这一行 加上 另一行
	 */
	/*	G.Add(0,1);
	    G.Add(0,3);
	    G.Add(1,0);
	    G.Add(2,1);
	    G.Add(3,2);
	    G.Add(4,3);
	    G.Add(5,4);
	    G.Add(6,5);
	    G.Add(7,6);
	    G.Add(8,7);*/
	/**
	 * 上面是靠初始行变换来创造调用 LU分解算法 的条件，这样的话无法求逆矩阵
	 * 这里是计算置换矩阵P来创造LU分解算法的条件
	 */
	cout<<"原始矩阵"<<endl;
	G.displayA();
	G.LUP();//计算重排列行数P
	cout<<"输出排列行数"<<endl;
	G.displayP();
	cout<<"输出置换矩阵P"<<endl;
	G.solve_p();//计算置换矩阵p
	G.displayp();


	cout<<endl;
	cout<<"b矩阵"<<endl;
	G.displayB();
	cout<<endl;
	cout<<"A矩阵"<<endl;
	G.displayA();
	G.LU();//计算LU矩阵
	G.LUP_SOLVE();//通过L矩阵和U矩阵计算线性方程组的解
	cout<<"L矩阵"<<endl;
	G.displayL();
	cout<<"U矩阵"<<endl;
	G.displayU();
	cout<<"y矩阵"<<endl;
	G.displayy();
	cout<<"线性方程组的解"<<endl;
	G.displayx();
	G.U_();//计算U的逆矩阵
	cout<<"U矩阵的逆矩阵"<<endl;
	G.displayu();
	G.L_();//计算L的逆矩阵
	cout<<"L矩阵的逆矩阵"<<endl;
	G.displayl();
	G.A_();//计算A的逆矩阵
	cout<<"A矩阵的逆矩阵"<<endl;
	G.displaya();
	G.A_Next();//计算原始矩阵的逆矩阵 = A的逆矩阵乘以置换矩阵p
	cout<<"原始矩阵的逆矩阵"<<endl;
	G.displaya();
	cout<<endl;
	return 0;
}

高斯消元LU分解PPT：http://download.csdn.net/detail/u013580497/8892729