Gauss-Jordan消元法主要特点是通过交换任意的行列,把矩阵A约化为单位矩阵,约化完成后,方程组右端项向量b即为解向量。我们知道,选择绝对值最大的元素作为主元是很好的办法,所以,全主元法目的就是在一个大的范围里面寻找主元,以达到比较高的精度。"
下面,我使用Blitz++的Array库,编写全主元Gauss-Jordan消元法。
#include <blitz/array.h>#include <cstdlib>#include <algorithm>
#include <vector>
using namespace blitz;void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b)int irow, icol;vector<int> indexcol(n), indexrow(n), piv(n); for (int j=0; j<n; ++j)
piv.at(j) = 0; for (int i=0; i<n; ++i) { double big = 0.0;for (int j=0; j<n; ++j)if (piv.at(j) != 1)
for (int k=0; k<n; ++k) {if (piv.at(k) == 0) {
if (abs(A(j, k)) >= big) {big = abs(A(j, k));irow = j;icol = k;if (irow == icol) break;}$} }
{++piv.at(icol);//进行行交换,把主元放在对角线位置上,列进行假交换//使用向量indexrow和indexcol记录主元位置//这样就可以得到最终次序是正确的解向量。if (irow != icol) {for (int l=0; l<n; ++l) for (int l=0; l<m;swap(b(irow, l), b(icol, l));}indexrow.at(i) = irow;indexcol.at(i) = icol;try {'double pivinv = 1.0 / A(icol, icol);for (int l=0; for (int l=0;for (int ll=0;if (ll != icol) {double dum = A(ll, icol);;for (int l=0; l<n; ++l) ;
A(ll, l) -= A(icol, l)*dum;for (int l=0;b(ll, l) -= b(icol, l)*dum;
catch (...) {
}0}1int main()
{ //测试矩阵Array<double, 2> b = 3,Gauss_Jordan(A, b);}+Result:从代码的过程可以看出,矩阵A的逆在A中逐步构造,最终矩阵A演变成单位矩阵,解向量X也逐步替代右端项向量,且使用同一存储空间。
注释:[1]主元,又叫主元素,指用作除数的元素


雷达卡


京公网安备 11010802022788号







