全主元Gauss-Jordan消元法
Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法，该方法的优点是稳定，而全主元[1]法可能更加稳定，同时，它也存在着弱点——求解速度可能较慢。
Gauss-Jordan消元法主要特点是通过交换任意的行列，把矩阵A约化为单位矩阵，约化完成后，方程组右端项向量b即为解向量。我们知道，选择绝对值最大的元素作为主元是很好的办法，所以，全主元法目的就是在一个大的范围里面寻找主元，以达到比较高的精度。
下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。
Code：
#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 n = A.rows(), m = b.cols();
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)
swap(A(irow, l), A(icol, l));
for (int l=0; l<m; ++l)
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; l<n; ++l)
A(icol, l) *= pivinv;
for (int l=0; l<m; ++l)
b(icol, l) *= pivinv;
//进行行约化
for (int ll=0; ll<n; ++ll)
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; l<m; ++l)
b(ll, l) -= b(icol, l)*dum;
}
}
catch (...) {
cerr << "Singular Matrix";
}
}
}
int main()
{
//测试矩阵
Array<double, 2> A(3,3), b(3,1);
A = 10,-19,-2,
-20, 40, 1,
1, 4, 5;
b = 3,
4,
5;
Gauss_Jordan(A, b);
cout << "Solution = " << b <<endl;
}
Result：
Solution = 3 x 1
[ 4.41637
2.35231
-1.76512 ]
从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。
注释：[1]主元，又叫主元素，指用作除数的元素