[原创]全主元Gauss-Jordan消元法（Blitz++库）

lckboy 发表于 2004-6-3 22:11

全主元Gauss-Jordan消元法










 Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法，该方法的优点是稳定，而全主元法可能更加稳定，同时，它也存在着弱点——求解速度可能较慢。





 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 = " <

}






Result：






Solution = 3 x 1



[ 4.41637



 2.35231



 -1.76512 ]









 从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。








注释：主元，又叫主元素，指用作除数的元素



[此贴子已经被作者于2004-6-3 22:15:49编辑过]

ilikenba 发表于 2004-6-3 22:28

消元法相当于在一个多面体的上，遍历各个边去寻找，所以很慢的！

lckboy 发表于 2004-6-3 22:32

嗯，就是慢，不过精度还算可以，用了blitz++库，发挥C++到极点了，现在应该比Fortran编写的要快的

ilikenba 发表于 2004-6-3 22:51

不会吧，Frotran和C++的速度一样，很难分出上下的！

lckboy 发表于 2004-6-3 23:01

如果C++不用模板，Frotran是比C++快的，尤其在数值算法上，但Blitz++库就针对科学技术开发的，非常的快～～上千条方程的方程组很快就可以算好，当然还要使用编译器的优化

ilikenba 发表于 2004-6-29 10:34

Intel出了Fortran 8了，据说性能又提高了20%！

loooog12 发表于 2010-7-27 13:28

数值计算强烈支持Fortran

后青春期的诗 发表于 2012-2-5 14:23

数值计算强烈支持Fortran

zqyzixin 发表于 2012-8-1 02:18

我继续顶你！太好的帖子了支持

MichaeLonger 发表于 2014-6-30 18:17

路过。。。

页: [1] 2

数学建模社区-数学中国's Archiver

[原创]全主元Gauss-Jordan消元法（Blitz++库）