全主元Gauss-Jordan消元法（Blitz++库）

全主元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]主元，又叫主元素，指用作除数的元素

 

本文地址：http://com.8s8s.com/it/it36711.htm