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

lckboy

全主元Gauss-Jordan消元法 * e% S1 g1 H$ c# _" D

. ~4 L/ ^5 e5 a; u3 |6 }5 R1 _

/ u! p7 t: f" k

$ {* O, \/ G1 G$ R: |

- h8 g4 T! T+ y8 O2 J$ i

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

9 F; M1 P& Y2 P. s6 f" f

. T- ~+ r) c# T% H* `0 a

; J" S- m8 \/ y* g

' q6 q8 }6 L7 W+ R \3 D5 F

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

2 W* [7 F" g- k$ R

, q! @* V/ m2 {$ |5 C* F7 |, \

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

* N7 K7 H( u+ P) O

% X( s/ c' e5 v \: [% k

7 g3 o: |" u R G( A# ~, n

Code：

: q9 x( I+ k1 W6 p( e

" c& _7 z9 d; T/ I3 ^* R

#include <blitz/array.h> $ g' R) P3 F; M8 n) w8 ? ) ?; r2 O& w; m* t8 N & q, Q: x) s4 M( u4 g% [' x

; ~+ n1 Y2 w3 W% a

#include <cstdlib> 8 r/ r" \" k3 _. K: n7 Y! R $ L, y% e, [. e) v; n

7 j; o2 c% K8 l' r

#include <algorithm> : c t0 [% f$ h/ ]2 Z. ~2 G ( g6 O4 X% B9 ^* y

* L& G4 Q9 g% o' q

- P b# U8 R6 {3 P

#include <vector> G7 P. ~% J+ G8 S! C

3 m# m2 l% U& L- Z# g, V

using namespace blitz; ! b' L1 A5 y3 O0 e1 X 3 k8 i) [0 X% O: k* _1 p; S

3 O; W. w+ T9 Q$ q/ n3 `6 p

' }6 N0 n/ c( n2 ]

# U4 H, i* _" M5 a( R( b

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) 7 g1 o7 Z1 A' u! B ]( S

' R1 ~+ T* x/ n# r- w

% U- _% t2 O7 l) a: V

{ " d+ ~( p4 {4 |" T4 d/ k

+ ]1 w( a# S. T

int n = A.rows(), m = b.cols(); 9 i |2 d' f! p2 Z+ j # O* I. ?9 l' }- K$ Q5 H

int irow, icol; ! o7 h1 `1 z1 P$ r9 K

& F' W( ]: E! [" U4 M

) X; G9 J \( H* n1 N8 T% o

vector<int> indexcol(n), indexrow(n), piv(n); + r# J* c7 w' t; X

' @( k- _% h, L* @

% f @2 I6 E4 s) C* g

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

' @$ R- ]% s% H7 L

piv.at(j) = 0; 7 b2 n% v% v9 W+ j s5 }/ k! i+ e2 H( f# d# v- K z; z

" E# `. `9 m0 x; e 0 ~4 S/ D( [5 c; k

//寻找绝对值最大的元素作为主元 , T* f1 ]& e1 Z# ? _ P' V * T2 l/ g+ v* ?5 n' S 4 R1 k# t' `9 F0 ^+ Z7 Y

; w. ]% Y8 R6 F8 d3 g

for (int i=0; i<n; ++i) { 3 t0 d* c, f9 o- D

9 [: c1 ], X6 c5 ?# p& V

double big = 0.0; + x4 d: G" d4 a) ` A

6 J7 s. w$ g: d7 R, p+ K

5 L/ P m, ~) x4 n X1 f

& c7 n' ^. j4 i, z/ @$ ~

' a* E$ C: M* `7 f. b' A; m

for (int j=0; j<n; ++j) & l; G1 {0 s6 R) O

: _. H3 U6 S2 ?# S: X% k2 j

$ G4 X* @, c& I, {1 K

if (piv.at(j) != 1) / t. p" j, D% F$ ?. O

. s0 P; P. H Z& B; X; a

for (int k=0; k<n; ++k) { & ?0 j6 d' R# E% \8 v) d# y / {2 B2 ?, S3 I* X1 J

if (piv.at(k) == 0) {

% m# K+ b* l4 s' R: A a& F, Y" L

- A" {8 y3 ~/ w$ ^ O( f

if (abs(A(j, k)) >= big) { 4 k& c8 p0 Y: c; P* X6 R; R

0 C9 q9 Z0 k& Q& h. c5 S7 {0 u

# C; F% k! F7 w* }' |

big = abs(A(j, k)); : c, P* ]8 W& m) E! e1 n" i 2 {, W0 q z5 y( c

irow = j; ) m4 P1 F% |( v) [ A! n' J4 _

icol = k; $ {8 {3 y' h9 ^5 l! u9 |3 v7 z

7 B/ b8 E! P3 V9 ]+ ]

if (irow == icol) break; ! a1 h' l7 S: @* H

7 g! n: h9 Z" ]6 `5 Y5 b& K5 l

} : C6 ^" T H1 ^* z2 \ ) ~; ?, Z1 M7 }

} d X; f2 W1 ?2 S4 @

}

7 o1 Y% g* ?; u9 f0 K7 o* P

' y1 j2 K8 u7 U, a4 X

++piv.at(icol); 4 `: P2 J; d( t w. `; E7 l$ n - p6 G2 w+ s8 N$ A2 B- K

) |# p: }5 H7 G ) Y4 o }$ M; A1 D# V& O5 g

1 r L+ ]' |* g- D. f/ `# ]

2 I: A% j& T2 G% J

//进行行交换，把主元放在对角线位置上，列进行假交换， . y/ r% ] G2 j. u* a- p ; ~8 F0 u3 p) M8 N( E7 C4 u $ Y) a- c% z5 F. B: K

. z% }* V# d7 @+ g+ M4 k) A

" ?: U' O0 g% a

//使用向量indexrow和indexcol记录主元位置， 4 Z: _. {" B$ M / R& ? K3 {3 y) z2 Y5 t1 S7 K

; W: g; `" z$ F# Z" c( c% b

//这样就可以得到最终次序是正确的解向量。

if (irow != icol) { 7 G5 v8 |$ T; Q% F) n- o# C , P/ Y+ j1 x; g: \( J

2 L9 z7 D6 [: k' `( J

, \, c/ U# u% s2 j" k' m9 ~

for (int l=0; l<n; ++l) ! Z: l1 R \7 @2 ` % _! @/ b* D6 F2 q, N$ S

' y t4 x% ~/ \5 y9 G8 @2 \

swap(A(irow, l), A(icol, l));

) \0 `, s0 Z+ I2 W4 b

$ s# ^* j+ E; G5 ~5 w

for (int l=0; l<m; ++l)

4 ]9 S. n* K/ a+ |% y+ ~

swap(b(irow, l), b(icol, l)); ) Y; J& L) s4 ~0 z1 c $ V H! ?% B6 m% L' g5 k+ a

, U" U/ x9 @; L) f/ r i

}

6 D: K* m' ~! K' h6 @$ E

- @2 W' q: H+ k7 A h

! a& w$ j6 Q! H# i

9 }* N0 a3 D1 S7 m7 Z

indexrow.at(i) = irow; 3 ]/ B# o: a3 C" Q5 d

0 n( M5 V# H/ Y# K( @# J+ Z

indexcol.at(i) = icol; 3 p& s v' O) c- M% @/ J

7 j2 K$ n" z8 b; Z: q* z% X

; ] V( F) }( }# \+ C* C

6 G2 m) `! e* f, B4 l

try { / j6 z# M7 d0 p1 P

' f" [+ L6 s$ U4 H5 m% Q' k" S

double pivinv = 1.0 / A(icol, icol);

; |4 W% `- ]$ q# j1 [

+ Z r+ {5 Q9 A7 x; u4 {

; k& S3 `/ Z0 M7 G; q/ i" o; Q

+ c4 J0 L; \# s1 S H

( q, r& y. s* h% z! |! R' P; X3 {

for (int l=0; l<n; ++l) # J4 b7 p9 }5 L5 l4 z5 G. s " J3 \. I& L/ q4 v% ?0 z/ n8 B Q

( r4 `, K, _- ^

A(icol, l) *= pivinv;

for (int l=0; l<m; ++l) : i* s% F# \; z C+ u4 y 5 g- f$ S9 m' B& n

b(icol, l) *= pivinv; , {; s+ T& p6 ?2 m' N

) \: Y) @% K2 z9 ~; |

: x. t G7 p4 M2 t

//进行行约化

for (int ll=0; ll<n; ++ll) 4 I( K8 @) v9 J% u; d2 @ 1 d% c, A4 A/ K

if (ll != icol) { ) I8 t0 R' D, o' R, H0 N

4 u% G' n8 I( R. z' d$ X, [4 `

double dum = A(ll, icol); - T2 e0 B" S# Q4 a. W8 w! `+ U+ [

6 `1 r( `0 K h; \2 I. ]' j

# F2 Z- U% N: S+ J

$ h. I% R7 z) V" D4 Z

for (int l=0; l<n; ++l) 5 B! k. x1 q* Y3 C8 G1 e

A(ll, l) -= A(icol, l)*dum; # k! s9 F* ]9 `2 j: B

; i& N9 L/ H# [/ ?0 |

for (int l=0; l<m; ++l) ; W' s& Z5 G1 U! L

! V, `4 a) F6 K$ C

, _5 J& m- H( x Q9 w

b(ll, l) -= b(icol, l)*dum; $ r4 _% ~; E/ J: m3 g U% x) j+ t5 c/ q, ]& A

}

4 K0 F0 E( C( T8 G, e2 P

} " P& p" A3 C1 R& ]4 A ; {- k8 E' x9 I2 j ) e/ O$ x( L3 U: S

catch (...) { 3 a. v: V% V6 Q" ` @ B

) r! {9 O/ U/ A0 k; M$ d

cerr << "Singular Matrix"; ( _) [+ n8 ^3 K5 C7 J5 |0 F) F# E# ^ V& H$ ~" @+ \( j7 ^3 i$ _

3 D5 \( ^; H1 w7 n7 h8 x

} ' X, {% {. x: \( e! e7 s

& y, `9 I8 m6 T7 Y" Y/ \; i

} $ o- Z A4 c$ m, U' { 5 A! Z; `' l+ m1 O" S! `

0 ]% f0 f- b% c3 H( I

- A6 }- A/ m) X* u+ [

} 0 J6 b/ e1 e q* V I0 y / [* z- I& Q: r! t/ B8 {5 M

; E" l3 D$ z; C- |9 }

+ x0 w* A9 u" J9 I) I0 Q0 }# g% L k

int main()

{ # J) X1 I$ y6 N. p3 ^! Q

8 L/ M! R; k+ e/ P8 F3 G# K

8 W5 O/ A* m1 y( W

//测试矩阵 % Z {& h; b, k) o! k2 Q . Z0 d- h8 E5 `; m

; Z2 @. q& p( t. ~4 J

Array<double, 2> A(3,3), b(3,1); * X9 p3 `9 S# Y+ x$ C$ S / L* a# H0 Z& r: j d; Q. Z 7 I# ?( A# T: s' p3 ^8 s6 j

: V, n" P3 t8 V) c

A = 10,-19,-2, 4 c" }; z0 k- M, X2 ^

9 b9 ]7 D# p; L m! E2 F

-20, 40, 1, " B% w4 \6 I$ k+ k- v9 u

0 C0 q; a" H& ~# q/ @

1, 4, 5; + }. B v- W* @. X5 \0 S " H3 _9 l8 ~5 b3 ]2 q

; T: J6 _ Q$ K4 F/ R/ U; s6 O Y

i( o+ h* p0 b# p8 k3 e6 ~$ u8 }; \- g% a

4 s. D. o8 h8 u }' A/ Y

+ x9 m2 l2 A: B, o2 d e

b = 3, " M$ c8 j6 b5 T! r

* ~* T/ U' c; U" F: V

3 x H# u4 n/ h" ?! Z

4, 9 E. M- |5 t) U2 S& ]1 Z ( ]8 V* O; _/ X$ p

) z, b( d/ _! s( j' p% d- I# K3 E

5 G% _4 s* B* U D' I7 D( m7 O. i

5;

" e9 o: V$ [. X% Z' @) s 7 i* E2 e9 H5 B: S% |, g

Gauss_Jordan(A, b); $ O9 p" U9 N( Z % A1 X! \) h- K5 O+ D4 `

+ K4 O4 s% A" t! o% v! T; U

cout << "Solution = " << b <<endl; 6 l v& g, t- S. ] 4 X4 N3 P: w4 p- U2 m4 Y

5 z8 L# } r, v$ n& A

& u6 r) q( w0 I C

} % t, t" v! J6 L3 a+ D 0 F# ~4 O( J3 S# F9 T& y3 G @" _7 l9 b

3 i! O$ c A2 n( B: D

6 _3 `( E# ^5 R9 C& Y: t

Result： * ] S; m5 E* j/ y E% B0 X! l

7 d/ F5 x0 G \9 h; Y

Solution = 3 x 1 8 p6 i! H# y" M9 ?0 b$ D& z ! i8 G2 F2 d. k; O. j( d" t9 {1 C

0 n1 y+ W9 T: D/ X& ]: l2 X2 z& W% H7 G) H

+ C5 S. ]5 G% {

[ 4.41637

5 ?0 _$ Y5 v8 G: b% ~1 \5 B

2.35231 # T8 H; v3 C/ O2 p( v) p5 L6 I

-1.76512 ] " G6 l$ e6 {( B4 }; R

; m0 v( U4 m# }: G# I4 q, I, w% C

$ t7 n$ @$ R2 j6 R! C

/ R3 d, i$ l* o8 h; [# V2 R' b

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

. @ {$ A/ J [0 w2 k1 e

' l# G) z) J }

: }- |0 D5 o* F- a2 E% Z

1 j: {5 o4 A7 L# y) f9 P" s

: A4 m+ T6 |1 v# t8 N% _) I( B

" E5 b8 p0 L8 _/ ?

5 h8 x6 ?. D2 F s( H9 W

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

6 M, }% L) \' m" J) [ R6 |4 H# w

[此贴子已经被作者于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

路过。。。