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

lckboy

全主元Gauss-Jordan消元法

$ m5 k$ g. i* C. T0 F

1 y4 L& f6 r/ d% P0 ]. U0 v

* C! g) H f' Y" c

6 U% \. o9 W% c6 j$ c

Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法，该方法的优点是稳定，而全主元[1]法可能更加稳定，同时，它也存在着弱点——求解速度可能较慢。 ' o3 q$ t% T! c* [1 H4 _5 z. ]1 E

" [4 P7 K6 g/ p9 j6 B3 x9 j

: R U& E4 x, Z4 m) H# [( }2 S

( _" k8 G% C( a0 `* x

) Z& C" }4 ^8 ^: f; \4 X' i

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

( w! U& F9 a' I8 [

$ @5 I6 h5 d1 W6 a+ y0 e

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。

( `" E$ f1 E ^- _$ t- I" a; q

/ b7 ]8 @- n5 Z g. b0 y1 A& c

* v% ?* { Z) P

Code： * H9 D$ |/ x9 _9 t+ Q N

7 `- a7 K) J+ Z, a5 l4 m- n0 v

4 r; _* i& p- C) f/ R9 I' X

#include <blitz/array.h> 8 Y! v' A2 P+ k( T- w: W* p

& Y8 {* J, B0 N/ L, i

2 R1 ?9 |( Z& H3 w& f

#include <cstdlib> ' V/ ]0 K$ \/ S! s 2 ?# [7 M$ d9 J* e+ b - o, F3 [& |4 g% b, Q: h

: j3 ~+ Q' R: \

#include <algorithm> ( {' |: O7 B9 X- z

3 c/ o) }/ f( O& [. L" D% D; |

#include <vector> 3 O1 D6 {, X g( g, D% u 7 z' L/ F" H) N8 d W1 U& J* a q, H

; L, G; e/ I+ C* i/ p

using namespace blitz; ' ~. f6 p+ Y" C. ]9 c* ~ + \. M+ T. {) ~* v( ^3 L6 t& d2 u

6 N$ I" h3 q# p6 k, l1 d

X! @# r* H8 r9 x; o( n. y7 M

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) / p) l9 V7 z4 G" }+ H" g% C% Y; C

! ]/ }% Z' `( d9 m p

{ $ k) n4 \! @1 I ) N* M2 o, k; \+ p3 ?( c( X; | 9 b( ?! K+ A( G( n9 S# R8 t3 I; I

7 @% m9 ]- E% i: d% t! t

6 C+ @) D7 j/ o- v, w1 ^4 g

int n = A.rows(), m = b.cols();

* [3 [2 k! R6 |3 g

int irow, icol; * C6 F2 c0 q: K# }- a! m

vector<int> indexcol(n), indexrow(n), piv(n);

2 A' p$ i. Q7 [4 r! D# x1 \8 Q R

0 P! q! m- C+ F! G# T/ a

- w# s$ B! V$ W: v+ ~9 B8 J: F5 [

" O) }) L0 n; q7 o; y4 _1 W/ i! W

for (int j=0; j<n; ++j) 2 d- K4 e! z; f: A" X

8 \/ D0 m. h' T K! o

$ U( e4 V; f1 |4 n) ?9 A" l

piv.at(j) = 0; # n5 v+ m! i; ?* p4 w

( C6 |0 Y! k# J

; {% {, r T9 F7 d1 I) U

5 ~, I% z4 k( [

7 u) T% ?' F+ |' l3 k4 {

9 w; V$ ?3 K+ ~, L3 \

//寻找绝对值最大的元素作为主元 " U1 ~) H) a' H3 o( S# n * p- e- i0 v( F5 y- a & }1 h9 L. U0 E5 C

. l8 k1 w8 c9 h' e6 u/ v

8 A! y3 w' X8 e( }, T L, C

for (int i=0; i<n; ++i) { , [% T4 O2 v9 w3 r) Z

# C y3 C0 y2 h% `' F- A

double big = 0.0; 2 u& ~( s8 H/ S8 C& W ) |* E( w2 b" h1 j6 I

( u/ w: b. X$ u, u1 N9 Q$ t) b

8 B2 j( }+ ?4 n* j' @+ R/ s! ~

, C7 O7 R% B. {, t

for (int j=0; j<n; ++j) + j9 r5 p+ }; k 1 @1 x& U% T" o/ a" g, |. c k2 H

if (piv.at(j) != 1) 1 D4 y5 n8 `4 A/ Q9 q - u& `- B1 C% H 0 r4 I7 G2 a- ~3 s

* u+ e+ }* e- B

for (int k=0; k<n; ++k) { * ?) n9 h0 R9 W# @ 4 `/ u' S, e# O1 `

7 {2 b3 R! D2 _ O0 d, N, L. ?

if (piv.at(k) == 0) { 6 n0 q q8 S& J6 N 9 z1 Y" \4 }( ?7 a& ~1 Q+ J

# U0 s8 i$ z' p; w# s

if (abs(A(j, k)) >= big) { . \, l$ n3 h' Q, x/ f

* [# q- H0 i2 x( K- w* |# f6 f! j

big = abs(A(j, k)); + X' K- P3 g( A1 ~1 P% u( l1 t % { K. w7 h' C2 ^( ]9 O0 n# ~

irow = j; , n" V' f+ h; [ ?7 {( ^! L 9 I) v) F9 c% X5 h' t/ {

; k0 a; r; O* N2 g# D

icol = k; 2 J6 g2 G, {# Y1 w9 V& N, J

. r# D: Q& A. F. N. ]

1 _% p* ^( Q' [5 ^: N+ ^: _

if (irow == icol) break; 4 j3 M' l6 x$ L! d 2 X0 H" p& B+ s4 h

) s2 {/ h1 S0 H' A" R# C

2 o$ e( H, H1 \# ?' N9 a Z

} * @ ]1 |4 [ K; _7 F. x5 b3 I

} 1 j/ R0 G/ |- g0 t) E

: R% j! S% F0 a# H; V7 D d2 r

} ; h! L8 q4 k. O- d ! k$ ^. D% ]2 k8 L6 m

' S: W* }# x& J2 m8 e& r$ R; q" U( v

& p, O% {" N. M2 C4 J

7 k4 o5 Q: ~0 }: c' |

++piv.at(icol); 8 c: E t8 c* ]- `4 w

& N& d, A/ @) i* i8 n( b5 [

9 o0 S. A1 ~* o0 S& P0 ^. m2 K

, S) S. W! Z0 u

//进行行交换，把主元放在对角线位置上，列进行假交换， 2 k4 I3 P6 H$ b& j

//使用向量indexrow和indexcol记录主元位置， % R6 N- a0 T$ @! N- d% A. u7 ~ * H( @: z6 t1 i- w# u6 K

//这样就可以得到最终次序是正确的解向量。 6 d3 ]3 U5 |- _8 C$ I" B! s3 Y$ I a* `' ]# l! `$ X. {5 F1 a

if (irow != icol) { 6 \9 o6 T5 D% E: W% n. ~' h9 |: z 1 f0 D* P$ R$ m! K

for (int l=0; l<n; ++l) 0 A# D; D& X( p' ?- n" K( i

" I% [6 A J6 S6 d. ~1 |

: B2 f! p4 x7 i4 Z; f3 i

swap(A(irow, l), A(icol, l)); : M7 D5 i9 W/ h3 ?3 L

( u6 O. f |- Q

% l( P; j: M% @, U

# \( S4 l4 q7 ?8 w- S0 m

for (int l=0; l<m; ++l) # a& d+ p- u+ t, H# z9 L6 s 6 [' \' ?; j! ^& x- l+ d! R

$ g' R1 y4 X# H _: p4 S

swap(b(irow, l), b(icol, l)); % z2 {4 G3 ^" m# R# `; V f $ k9 h3 V" U, ? s

# P B. x- F) f; Y9 W- N" K9 B

} 7 M' _0 K/ @3 v Z- g5 s/ {" t0 z B9 u; U8 Z: i4 t3 I$ A & H% g8 J# ^. s0 ^

/ b' ]/ p6 {& W- D I

( d4 I- B* a# A3 I

& @1 y7 W, Z+ z+ o$ L" s

& H8 p2 `5 _ k! a* b9 R

1 w. B/ Y/ s% G

indexrow.at(i) = irow;

indexcol.at(i) = icol; / w* |" [8 O4 C$ f ) Z* B' j4 Z& q+ w/ }' a5 P, M

6 P1 N( ]7 f7 I % D+ m) H" W% J- r; r+ d, E* a! k5 ^ # d% K& K# c' N" M! q* F- E

) O2 @, R! i/ r) O. P

6 j3 R( C' t- f- p& {" Y

try { & w1 Z/ t0 h& S : F% E8 t0 K/ Q0 g : ]* \. Z3 a4 n) H

1 n+ Q- S3 h2 _, Y, z1 U

double pivinv = 1.0 / A(icol, icol); 0 D" K8 k3 {3 z7 b, ~% U : `, T+ J, }2 J K2 u/ | M; o3 q! O& X& G' u! O2 a: I

; X: C5 F, p5 w, i$ n+ E

. v' [- ?$ t6 |' c

for (int l=0; l<n; ++l) 5 y) A. p! Y( f+ N$ s! ` , \- ~8 Q) k+ `1 ]- D

0 V3 u" J& G% l9 J8 U) v6 X1 ]4 @

A(icol, l) *= pivinv; ' ^) C5 v; Z% ^% U6 h' t 3 `. x* p7 H0 S/ p3 B+ E 2 m$ o7 a# \) L! c* a) i) z

7 X9 d: R* a/ M+ [* Y G1 _

# }& R1 O5 V+ D( l; L% e- j

for (int l=0; l<m; ++l) ! P/ X* U! \- A" t

b(icol, l) *= pivinv; * r6 c( Y b9 |; P' |8 y. ` 9 `: B! x- c& m0 Q9 Z

" v' h# K: k. F6 N

//进行行约化 8 r$ d/ b( ?( \/ H' f$ n* j, s) w 1 t7 T$ p. F! d8 `

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

$ k" ?- N" d1 t0 ~8 k8 O8 [: W

5 m( o. h( I2 {0 h" f

if (ll != icol) { ) z7 |' D# B9 @, I& R M ! }5 o4 E! E) O3 l* p

8 j6 c! \; U! b+ g1 \) |

double dum = A(ll, icol); H, D2 m \4 a) l; t

$ x+ r! M, K$ ~% B

for (int l=0; l<n; ++l) * p0 X6 N9 c, U" c5 ^( l

' i3 D' O# l9 z

A(ll, l) -= A(icol, l)*dum; $ g) f P$ m2 g

% a& b9 ]2 y) m; |9 j% Z' v

" m* A8 j& K: B8 x1 J

for (int l=0; l<m; ++l) - L/ R0 G& E# h0 ` 0 u; _! P; X6 |' v1 |

b(ll, l) -= b(icol, l)*dum; ( v# `- H D! k J! K 3 |' M5 j. i' z0 C7 j# H$ Y

}

- x4 X. D, b9 I, W, _, t7 q

6 F6 \9 J. Z/ k1 O2 }- B2 ^6 C

} . s* W3 U7 e' p# L( W: U9 Z! F 3 G& F0 Q0 C: S* Z9 s

7 B3 q4 x/ h0 W+ A* y

catch (...) { * |1 ^; l, J/ R5 P+ q * J# |) F# h6 Z, b- A5 Y7 Z

$ O& t5 S4 w( n- N

) d2 }. T$ ]- a' k( _& M

cerr << "Singular Matrix"; 6 v& t8 A% o- F" @ 8 _- L% ^! s$ s! D; }0 f

: @0 U f a0 O1 h# [

} : v9 u2 Z& l; L" y0 B$ `

. ~" P. I' _0 Q4 E- q) d/ J" b

} Y8 R% I! i- U) y/ L6 C9 K

( @" s/ |% N8 t3 b5 n* o

* [- X3 G9 X! {* p- C: q: Y

} 8 C4 m% ?4 r) M. X7 |# U+ P7 t " x# V7 r% C, O. v I. |

# m# ~3 W, }) J" d4 T

' E B. ?: |8 T$ i; h, W! }" X

6 l6 P. T! k& g6 q

int main() # H1 T6 r. P/ D* P( h) ~% i7 C! b

8 I+ k. W& j, C

{ . P8 s1 Z8 D$ J& F* v1 M 0 F% d! Z7 W7 p, N L

- s/ X; [% x; @5 @5 R

//测试矩阵

Array<double, 2> A(3,3), b(3,1); ! R" V# Z, R9 s6 l) f : Q0 f1 L# S2 Z4 N$ B1 { / c: e% `- Y- M! F& x

; Y4 _5 _1 D* n& d

A = 10,-19,-2, * J2 E1 w Z6 ?5 W, z% L) @ ) ]2 U- M, S8 f

* s3 |4 f0 U" D3 i

1 ?& J1 @. P3 j

-20, 40, 1, 8 P( `% H9 {) d- ]5 G. @: g * v A$ i. O. l( q s / G4 i' A6 |' x. ^) u3 U

, m m. t3 X% f/ L# O0 a( m" M

" F( ]9 T% F4 z' [

1, 4, 5; V; I5 N$ D- a3 C# z! j9 Z! k2 K 7 h! ^* V$ h9 W: Z' o3 D4 x ( v# [+ k; N1 h* D7 a" T

b = 3, * Y/ y. D4 ]8 r 9 y' B3 D g( e & s; V5 z7 H* ]& r

" p5 c3 f6 h. e

* q- s: G( Q* R/ b1 J" x

4, 3 I4 U, V! r" z# H

/ }6 D+ e9 a: V% D$ R0 i0 o

& e3 v# K& U7 u% _/ D

5; $ H6 I. B8 p; J3 T/ _( z

9 ]+ m1 ^; b2 ^4 Q! Y( S; m

6 c& u8 Q: }: M5 c' F) A+ O 6 V6 P& V' X# V+ K% j0 X

Gauss_Jordan(A, b); . W* I2 Z0 I& G; [& d

* V% @9 S( m$ D! M# P

' I* K3 D8 ?# q4 |

cout << "Solution = " << b <<endl; 5 E$ y9 l5 x2 l; c ( U% a$ ^$ @3 D: u% x

}

) C- T. Q( m9 }6 K

* c9 n: P5 [* s' m6 v

1 `. H# ^9 i: |5 d. r) H3 ?

Result：

# I& \, S5 J }3 v w7 S

# T+ j( o8 T1 a! p

Solution = 3 x 1 1 R1 s1 p# {3 I+ [; p, u) e, q f ( J7 b" \' s: ^5 X+ b" d: ~

[ 4.41637 . `- [) f' K% L! ^( L 5 z; `: S) X& m; @

8 z) g; o. l6 O; n3 z4 ?3 _% ^

2.35231 % M1 K, X6 Q0 o7 o: x2 U+ \ % w* u4 l( B' w! z

* X6 j' O) `, M" {, ^$ q. O

-1.76512 ]

) A# l: X: d3 N8 [' o: m" E0 s- a

A( x. H& f( J( B: f# ^

% S7 ?- U4 w: L

( A1 a( I' a& Q+ E8 h# y

- x% s8 i0 `( m- F. y

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

; v0 Z" _. K1 f. _: E3 _

5 c% M, ]4 O- n: w% c* d

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

7 u. c$ R/ c$ V( @6 _: q

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

MichaeLonger

赞赞。。。

MichaeLonger

楼主的帖子怎么样？赶紧试试这里的快速回复给楼主点评论吧。。。。。

MichaeLonger

学习学习。。。

MichaeLonger

看看。。。

MichaeLonger

路过。。。

zqyzixin

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

后青春期的诗

数值计算强烈支持Fortran

loooog12

数值计算强烈支持Fortran

ilikenba

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