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

lckboy

全主元Gauss-Jordan消元法

7 d0 P" O" A; Y3 B- a7 H3 L

$ n$ Q( B7 {/ ?6 g* e1 k d

6 l8 i* ^ r! M" U- t4 [# q

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

, E( W1 k+ G& s1 h" G+ K

' p y. ~( b) i5 x" c

9 {. K" ]/ h7 `) ]

8 ~7 H% b/ `' }* a! l

Gauss-Jordan消元法主要特点是通过交换任意的行列，把矩阵A约化为单位矩阵，约化完成后，方程组右端项向量b即为解向量。我们知道，选择绝对值最大的元素作为主元是很好的办法，所以，全主元法目的就是在一个大的范围里面寻找主元，以达到比较高的精度。

; c* }! x$ W7 j0 f2 i

) O" s$ N$ I' S8 A: d' H- F9 S. J+ M: b

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。 ( G. M5 Y/ _7 |0 ?" C f8 L8 }8 Q" N S

" n; W6 n- l9 P8 R! Z2 M* D

Code： 8 z: }4 q( N% t2 V! M. k0 L6 E/ A : L0 O% U1 X: l4 u$ t: Q% }0 Z 4 z ~5 X5 f& D# v* L, e

* l% W3 r: K" f; J5 A3 e l, N

; ~9 z- C5 f5 o; ~+ M; j) U

4 g2 [: b$ x' k0 n* J7 l- E! H

% k1 v# _5 z- O7 c V/ `1 ~

#include <blitz/array.h> 2 O. T8 A% D) }7 f; M& |9 H6 l / V, r( A( ?! ^6 R

/ y) ?0 E0 e6 z

#include <cstdlib> 5 s0 V2 r; q- w4 K1 R , }6 }$ Y" [/ y/ O8 S

. ~+ {+ W( r0 z0 U4 N8 o

1 Z% G4 l. s8 c: g

#include <algorithm> , w' v, W3 O' A& X

#include <vector> : O9 m. \3 _" R- Z/ s . a3 j4 p, w9 e% d3 B5 U0 m4 j

; }$ ~9 ]$ X) H6 W

, ]1 [0 k* j( h8 J) {' r

using namespace blitz; % |- w+ f: }, \( U6 M2 y) S0 D 9 B) j2 l# r- |/ D" Q ) f2 `. H+ o9 G# E8 T

1 L$ `9 {/ o' h5 `# ?+ H8 b

' s/ v( ~& a- o8 z. R( Y

, A Q1 Z4 f' ^9 t, l$ z2 _

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) : E9 R ]: z( E# Z" S4 W2 e( G* v- J " n" y5 }" k- A% ]& ]2 F 9 q( e- v0 K |4 I

7 x" D1 E: |8 S; i1 A2 e

{ 6 Z7 \7 {. X9 N: X

int n = A.rows(), m = b.cols(); 4 h2 |: F! |6 K r( w# {- D / [; d7 I6 ~- D, R/ Q% `

int irow, icol; # r4 ?0 n7 a& f. v; O

( z& z* A4 [- U8 p q

vector<int> indexcol(n), indexrow(n), piv(n); : }: T4 [- l }* y. n) o3 ~" F% h) N , G, e+ Y9 G7 n 2 C; g& s4 S) R8 f

for (int j=0; j<n; ++j)

9 ?- K3 y8 n9 [

piv.at(j) = 0; * P. v$ ]8 j& q3 S4 Z+ g ( }% y& G* n T- K( ~6 s

( M8 Z' ~- _4 q2 [( \9 d, Y$ O, R d

9 j: u8 c: [5 H( |9 o% Q% H3 M

//寻找绝对值最大的元素作为主元 $ s+ P/ y& P y9 k( H2 Q- a9 r . g1 G& ~4 w" t% T8 b5 E7 G6 f ! Q. k: d: D- R( z" X. ]0 `1 |3 C

7 n- ]7 p8 Q$ S

for (int i=0; i<n; ++i) { / v3 }2 Q7 e$ ~# }+ D# r

9 G6 M* l5 ^; S C: Y! r& U

, @3 E% }. C1 L$ x- ?

double big = 0.0;

3 @& \: z% ~. E7 F8 n2 u x/ k

for (int j=0; j<n; ++j)

; ? y$ V& F0 U; [ @8 E

if (piv.at(j) != 1) ' Z& g/ h8 N2 l ] A - P, ~. a2 |3 C: n1 v2 S/ f7 R9 _

6 [3 V y2 c3 I4 j+ k& r

! E( f: R) C, I. p! ^/ ]

for (int k=0; k<n; ++k) { * E$ M/ h+ e: j% ]2 e' x' S; r

if (piv.at(k) == 0) { 5 p4 ^* I5 |' m5 U; @

( {$ @- h* K) h

. ^, i+ R# l! @3 h3 d3 l

if (abs(A(j, k)) >= big) { ' ]7 c: f# m$ }3 G . l; U, C% [& @, m8 E! e n

, K6 I' R$ p- F3 ]+ H

big = abs(A(j, k)); # l) v1 a3 ^$ ]+ n

) K) J5 E7 G: c$ _3 F. @

irow = j; % O- Q# w$ K7 e& b $ O. [4 |7 |8 g' R

icol = k; # O& [5 h/ Y7 y: V! G. k* W ! }( ]+ C4 b7 E- }' }. ?

' g' K3 ^$ f1 P' J

if (irow == icol) break; ! W; X. P: g z( I9 M7 ^1 F

* f! I& A( G" m( y5 k. c5 ~

} - v; L) Q/ [5 {6 s * C3 O) b: }% m1 N$ a0 L [

, w* m* M/ c7 g* [+ a; b

} ; u0 H' X5 J$ j. @1 H

/ R1 e' R4 M% V. c$ N& r

} y: p4 m' _) m: r' W

, O' E$ V" ^2 `/ y/ x

! L- i/ v) |: }& I

++piv.at(icol); 0 V! N; {0 v, M# p; s6 g

0 _. `4 ~/ i9 Y3 V. u1 W" f7 n

) G! I' C2 y3 g {. j3 L$ F

//进行行交换，把主元放在对角线位置上，列进行假交换， , X3 s8 K3 j' L" |* V ( n h% b" h1 D0 t1 i) l ( C" O2 S, `" k- B$ P

6 \2 }% p( l Z0 M+ y

/ z0 J- M6 Z* \( h3 I

//使用向量indexrow和indexcol记录主元位置， ; Q& l2 D) r( Q+ Q) v

//这样就可以得到最终次序是正确的解向量。 ' n! P, f9 M1 [4 p" i& p 6 Q: S! n7 b0 V* x$ J! _: \7 _. Z

& v$ l0 a3 l) b# h+ z

if (irow != icol) {

- ]( }3 K& A$ N0 I5 m

for (int l=0; l<n; ++l) + `; x s$ j( m$ | 8 ?$ S! x$ l/ J7 f' Q/ |

, S. O8 I& t# |) [

swap(A(irow, l), A(icol, l)); & P+ a* u0 A; ^

for (int l=0; l<m; ++l) % G+ M/ i* ^3 x) {) u2 I2 V2 \ - F$ V$ J- g0 E9 J2 G

" X+ P& F7 X9 `1 J- g3 U

6 g+ [$ x) K9 ]# ?9 o' }+ Q

swap(b(irow, l), b(icol, l)); ) q! n* X% O4 i/ S; \8 {& _

} . y$ C- ]2 x5 m/ u8 n1 ^

1 m9 z2 l8 x* B1 W7 m

, i/ Z4 |8 d. N. F- p, \

+ @) j! f$ V2 y; _ ~9 I3 \3 L

+ O6 z A1 A4 e3 R( q' |

indexrow.at(i) = irow; ) E% v6 W; \ S6 x4 t

; L% Y4 S! [! p: t

1 u0 ~0 e# v: B8 C8 O( D; }4 G

indexcol.at(i) = icol;

& t" O- |4 X$ {4 x 3 I( [9 d m- Z1 _" @1 u2 U! u

F, @* i j: w; }0 c6 u i3 q

( k) Z3 G; M( Y; C

try { 0 r- ~3 }7 u$ j+ W& w( M& B

/ K( W5 s3 {/ a7 U" x. ~5 w/ Q

1 @8 a% P# Q! ~& M" Z

double pivinv = 1.0 / A(icol, icol); / d+ d( X; ^. t- @) f6 r ) S* ?0 T; x( @4 Q

. O9 w) q, g7 i- E1 W) g, i5 m! S

# [, y+ [7 q! C- M* l7 l

! ]; P3 Q, |. G/ s* D) [

for (int l=0; l<n; ++l) 3 Z7 L( u8 Z! x) s4 X: [( o : W$ f$ a5 Q) y: P; j

4 ?0 B* X# p, U/ i

A(icol, l) *= pivinv; + c; U* Z! W. }8 u2 s : z% j9 c+ y4 M! H% u( w: j- M- j

# u7 e `4 ^2 u2 H# `. f u1 |

for (int l=0; l<m; ++l) # f6 V. j: ?1 t 9 A% n {& v& B

4 p( q# O# o2 F; C

b(icol, l) *= pivinv; # C. [3 t! y) z ! v0 K; l5 a% t y% a6 S% s5 l# b, u F

: e" T3 }" p/ q8 C3 X/ [$ m/ n L* Z

" c( e0 j& H* L, R5 k4 d

) P3 T; e/ j0 g- m! s* d

* h8 U; a3 |% G( ~) f

//进行行约化 # u; Q) [7 q, C9 i( v; q$ ~

% p3 E3 B! f3 ^7 a

( A. u/ ?" A. I# q: x8 V+ f

for (int ll=0; ll<n; ++ll) , }2 g. b+ z% b

8 {$ p! E! ]: b1 ?

+ \. r& T% Y. f4 t1 t4 F

if (ll != icol) { 5 E( t7 z) Y: P/ C6 `% {9 u, h

+ A! X0 N! v- e8 S E! b8 R

double dum = A(ll, icol);

. w* Y9 b" K, ]3 @) H0 t3 X. u

; D( Y- _' [, ^! ` d0 C$ m

) b: {& k6 C& \. Y2 D. H

for (int l=0; l<n; ++l) - Z. w) K3 ^6 x0 v

6 ~+ F1 ~/ S, J. g- @

A(ll, l) -= A(icol, l)*dum; - j% \9 ~" v0 v

* t. S7 W2 l6 T: ~% C+ w

for (int l=0; l<m; ++l) ; P+ i& k# K4 {8 z" e3 L5 L

* l! k4 q; m: d: ]) W

b(ll, l) -= b(icol, l)*dum; 7 \9 ]' S8 x, d& C& ` 6 s6 T1 D$ M& x; ?0 o1 y

* g8 d, J# @2 j6 x' k0 f! D. V4 L

} + e) }5 R3 f4 {( S/ C1 E7 N: P

} 5 n# F' f6 ?" R6 T+ q& K

! A2 v9 W$ s- E( B3 R7 D

catch (...) {

4 ?- g4 J6 x2 e7 j5 C1 p d

cerr << "Singular Matrix"; 3 P# \3 ]; t4 O . t* Q3 x4 Q& b- u % N, |% @+ V: A$ `. V' @

. b# F" v% y7 I' f4 s1 a+ b

} , a6 E) C4 m* H% L' @/ q- S, Z$ ~6 D : Q0 u* D; N$ }5 M" Z

7 Y. j$ x6 i2 _0 o7 {

} 0 m0 a. R7 _, I4 O) N4 ^

& @2 O: A9 p7 r# q4 {

" g1 Q: F: ^0 t/ d

} 1 M; f+ z8 o, |! x) W/ @# F! b

. ]/ ^7 W- J% D1 l

int main() / k: m4 l- w T1 @: r; v

{ 9 W- t: x# s) E + J( ~# [7 [0 C- ?! _! m# U

6 f. N. n' r+ X- }3 |

//测试矩阵 % g6 a1 a g" J P I

2 s4 u" D: p5 g5 ?( Y

Array<double, 2> A(3,3), b(3,1);

, `8 Q8 D; U- O3 O8 Q

% g" [* H5 M. f/ W

A = 10,-19,-2, 4 ^: M( { p O' f

-20, 40, 1,

8 p# a R5 g+ G, @6 v

' F4 Q G1 J( ^ ~7 }/ b i

1, 4, 5;

) w$ t9 A* [; h$ _" t7 f6 p: ]

5 q( [2 g+ F# p( e; S, u" A; B

2 G0 y/ u5 [$ G. G" A# R; e

b = 3, ; ~* o! j7 r4 K) f0 P t; B, v : ]+ V' x9 l c# X, |+ B

3 \0 l8 F% e' N. q

- J5 j/ j+ n, L! w" H& q! W( T

4, ! L+ Y8 N/ `" C# Q; D8 M! P" O 9 C$ h5 x! j9 b( @' m* N+ X+ V$ n ) C+ u* D+ a7 d

2 O3 Z3 d1 t2 u4 y

/ G" r: C+ b& y3 G# b% ]

5; i9 \* v4 S# _" _1 ]) H& J* v

* k; _% S1 H+ h0 R& A* M

: q* K- V3 W# |

& K- r$ g* J# T8 q7 Y! n

6 @ e9 I& `1 U, E) f+ U

- v% I' |1 v) e! E

Gauss_Jordan(A, b); . Z. \8 d6 y, l6 i

1 P k4 B' _% e5 B/ Z+ k/ U

* U8 x4 D2 s% \

! M i2 @) U: Z

cout << "Solution = " << b <<endl; ) ~$ A" P' g1 b2 L4 y: O" a1 `* _

" O$ D( Y- C2 f4 Y3 t

} 1 e ^) \) U4 Q7 L. U. c1 x2 c6 ~ . t+ E+ T( N% u

Result： / j3 L: u, o4 Z3 x+ ^$ h5 ?

8 f3 l2 z% x5 {9 z: a' k( u: Q

# o4 h! y- x1 d) j5 K

: N! ?$ [, V* B

Solution = 3 x 1 & [7 Z; j. K/ \( m& R

: i0 A, A" ?' _4 ^ E

[ 4.41637 6 c: e" _! X$ T

( Z# I/ w/ [+ j7 y# g

# ?" i1 q) z( O; d* _

2.35231 5 t5 k$ L! m, j3 k1 S7 w % @( x3 F" S r

-1.76512 ] ( S, u9 k0 ]5 u P7 O, _3 P3 ~, m, v0 n& z4 X

4 w; ^# T- m2 R$ e

$ v9 v, l. H9 @1 j6 I% L. Y

8 X1 {, L) D' a2 U9 m; z# N

0 H8 ]" I: A* W# Z8 N; L. H

- z, W$ l: G5 }4 M% D4 W( O6 j

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

* ]2 }0 `5 S0 G. H3 b: k

/ w6 k& e8 g+ w( k

$ t& N: b& G1 |% P. @# y

7 y+ `4 Z! x2 F7 H) w8 h1 E$ k/ m

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

2 _4 D* Q1 I' B: y$ M5 A5 a3 `' d

! u& N2 Q0 Z e

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

ilikenba

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

lckboy

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

ilikenba

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

lckboy

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

ilikenba

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

loooog12

数值计算强烈支持Fortran

后青春期的诗

数值计算强烈支持Fortran

zqyzixin

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

MichaeLonger

路过。。。