我听取了我们老师的想法，把原书中的命名改为以意义为中心的，不简化的命名方式。原本函数名为gaussj，现在改为gaussJordan

Gauss-Jordan消去法解的是这个线性方程组集合：（我完全不会LaTex呢，而且鉴于复制贴贴的方便，或者转载分享之后会出现不能显示的问题，我会尽量少用Latex公式的）

[tex]A*(X_0 \cup X_1 \cup X_2 \cup Y) = (b_0 \cup b_1 \cup b_2 \cup 1)[/tex]

如果上面一行的公式显示不出，请看这行：A * (X0 || X1 || X2 || Y) = (b0 || b1 || b2 || 1)

其中A是一个n*n的非奇异矩阵，b是n*m的矩阵。

整个式子相当于拆分成4个子式子：

A * X0 = b0

A * X1 = b1

A * X2 = b2

A * Y = 1

算法的运算思想是：

施行初等变换将A变为单位矩阵，则

Xi = A的逆矩阵 * b

Y = A的逆矩阵

算法的步骤：

（1）对于每一列k：

（2）选主元pivot element（绝对值最大的那个）（为了减小舍入误差），记为a[i][j]

（3）如果选出的主元不在对角线上，我们把它换到对角线上

（4）如果主元在对角线上，做消元。

（5）最后再把交换过的元素交换回来以得到正确的A的逆矩阵。

/**
*@file gaussJordan.h
*/

#ifndef GAUSSJORDAN_H_INCLUDED
#define GAUSSJORDAN_H_INCLUDED

#include "nr3.h"

/**
*@brief Solve linear equations A*X = b using Gauss-Jordan elimination.
*@param	a[0..n-1][0..n-1]		Input as Matrix A in A * x = b, Output its matrix inverse
*@param	b[0..n-1][0..m-1]		Input as Matrix b in A * x = b, Output the corresponding set of solution vectors
*Linear equation solution by Gauss-Jordan elimination, equation example as follow,
*	A * (X0 || X1 || X2 || Y) = (b0 || b1 || b2 || 1):
*
*	(	a00 a01 a02 a03	)		(	(x00)		(x01)		(x02)		(y00 y01 y02 y03)	)
*	(	a10 a11 a12 a13	)		(	(x10)		(x11)		(x12)		(y10 y11 y12 y13)	)
*	(	a20 a21 a22 a23	)	*	(	(x20)	||	(x21)	||	(x22)	||	(y20 y21 y22 y23)	)
*	(	a30 a31 a32 a33	)		(	(x30)		(x31)		(x32)		(y30 y31 y32 y33)	)
*
*		(	(b00)		(b01)		(b02)		(1 0 0 1)	)
*		(	(b10)		(b11)		(b12)		(0 1 0 0)	)
*	=	(	(b20)	||	(b21)	||	(b22)	||	(0 0 1 0)	)
*		(	(b30)		(b31)		(b32)		(0 0 0 1)	)
*
*@attention A and Y should be square matrixes; b and X should be the column-vector matrix.
*This function simultaneously solves the linear sets:
*	A * x0 = b0
*	A * x1 = b1
*	A * x2 = b2
*	A * Y = 1
*
*Input matrix:
*	a[0..n-1][0..n-1], b[0..n-1][0..m-1] is input containing
*	the m right-hand side vectors
*Output matrix:
*	a is replaced by its matrix inverse
*	b is replaced by the corresponding set of solution vectors
*/
void gaussJordan(MatDoub_IO &a, MatDoub_IO &b){
	Int i, icol, irow, j, k ,l , ll, n = a.nrows(), m = b.ncols();
	Doub big, dum, pivinv;
	VecInt indxc(n), indxr(n), ipiv(n, 0); //These integer arrays are used for bookkeeping on the pivoting.
	for(i = 0; i < n; ++i){
		//Main loop over the columns to be reduced
		big = 0.0;
		for(j = 0; j < n; ++j){
			//Outer loop of the search for a pivot element.
			if(ipiv[j] != 1){
				for(k = 0; k < n; ++k){
					if(ipiv[k] == 0){
						if(abs(a[j][k]) >= big){
							big = abs(a[j][k]);
							irow = j;
							icol = k;
						}//if
					}//if
				}//for
			}//if
		} //for
		++(ipiv[icol]);
		/*
		*We now have the pivot element, so we interchange rows, if needed, to put the pivot element on the diagonal.
		*The columns are not physically interchanged, only reabeled: indxc[i], the column of the (i+1)-th pivot element,
		*is the (i+1)-th column that is reduced, while indxr[i] is the row in which that pivot element was originally located.
		*If indxr[i] != indxc[i], there is an implied column interchange. With this form of bookkeeping, the solution b's will
		*end up in the correct order, and the inverse matrix will be scrambled by columns.
		*/
		if(irow != icol){
			for(l = 0; l < n; ++l)
				SWAP(a[irow][l], a[icol][l]);
			for(l = 0; l < m; ++l)
				SWAP(b[irow][l], b[icol][l]);
		} //if
		indxr[i] = irow; //We are now readey to divide the pivot row by the pivot element, located at irow and icol.
		indxc[i] = icol;
		if(a[icol][icol] == 0.0)
			NRthrow("gaussj: Singular Matrix");
		pivinv = 1.0/a[icol][icol];
		a[icol][icol] = 1.0;
		for(l = 0; l < n; ++l)
			a[icol][l] *= pivinv;
		for(l = 0; l < m; ++l)
			b[icol][l] *= pivinv;
		for(ll = 0; ll < n; ++ll){
			//Next, we reduce the rows
			if(ll != icol){
				//except for the pivot one, of course.
				dum = a[ll][icol];
				a[ll][icol] = 0.0;
				for(l = 0; l < n; ++l)
					a[ll][l] -= a[icol][l] * dum;
				for(l = 0; l < m; ++l)
					b[ll][l] -= b[icol][l] * dum;
			} //if
		} //for
	} //for
	/*
	*This is the end of the main loop over columns of the reduction. It only remains to unscramble
	*the solution in view of the column interchanges. We do this by interchanging pairs of columns
	*in the reverse order that the permutation was built up.
	*/
	for(l = n-1; l >= 0; --l){
		if(indxr[l] != indxc[l])
			for(k = 0; k < n; ++k)
				SWAP(a[k][indxr[l]], a[k][indxc[l]]);
	} //for... And we are done.
} //gaussj

/**
*@brief Overloaded version with no right-hand sides. Replaces a by its inverse.
*@param	a[0..n-1][0..n-1]	Input as Matrix A in A * x = 0, Output its matrix inverse
*/
void gaussJordan(MatDoub_IO &a){
	MatDoub_IO b(a.nrows(), 1, 0.0); //Dummy vector with zero columns.
	gaussJordan(a,b);
}

#endif // GAUSSJ_H_INCLUDED

大家如果想要做测试的话，可以试试以下代码：

/**
*@file gaussJordan_test.cpp
*/
#include <iostream>
#include "nr3.h"
#include "gaussJordan.h"
using namespace std;

int main(){
	MatDoub_IO a,b;
	int rows, cols;
	//INPUT DATA:
	cout << "Please input rows, cols and the matrix data for A: " << endl;
	cin >> rows >> cols;
	a.resize(rows, cols);
	readMatrix_console(a);

	cout << "Please input rows, cols and the matrix data for B: " << endl;
	cin >> rows >> cols;
	b.resize(rows, cols);
	readMatrix_console(b);
	//Do the Linear Solve
	gaussJordan(a, b);
	//PRINT
	cout << "Gives the Inverse of A: " << endl;
	printMatrix_console(a);
	cout << "Gives the solution matrix: " << endl;
	printMatrix_console(b);
	return 0;
}

• Relative Posts
• Global Relative Posts