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

lckboy

全主元Gauss-Jordan消元法 L8 ~9 H- q1 G/ U* ~8 o8 Z

" R. ^! O$ i M1 ^, w

' k- n- H' `4 ?. ?' S

& {0 \; `$ q! h, b: r$ Q

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

4 a: V1 N3 U, H2 l9 d- D' e: p

5 j6 h, X4 j. X1 i" W

2 ^/ Y! A( H( P1 @, k! ?

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

+ N( I9 P; E6 o' t& E* L p

6 U. a9 r l0 T5 H

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

L D, `+ G& X

+ [7 _( p) B; P5 }. ^+ n7 c+ q0 |

8 P/ C) i+ l+ d$ R, k( H. ^4 v# ~

Code：

# k9 S4 ~9 q# C0 M

/ o8 _/ [$ a3 w9 u3 Z

$ g+ d3 w0 ]* n' e! Q

#include <blitz/array.h> 4 u* ^" s0 E& L! G4 t # i1 ^- M" \4 P& ^! c$ o

# R' ^1 n( t- h; P" n5 V% d

#include <cstdlib> : h, ~; s* |& s1 s0 {

+ f! X1 L) u2 P3 j

#include <algorithm> ( |8 G( i$ L7 c/ A9 W 8 G7 k: r4 F" P* L

, f. t8 I6 x8 a- J# J

#include <vector> - [+ | V( \$ U' q6 \% _ 2 f9 n3 r* r4 w& B w

6 D% c/ V6 `( Q& Z/ E9 }- c* Y4 X

& i* a+ c" _, Q. F3 Y

using namespace blitz; 2 K! z2 O+ x' `, i! P* U, Y6 x

0 p7 M/ ]$ m' ?/ b0 j( E/ d

. u# @; _$ q' ?4 L5 P y

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) ) u& l/ x! @ ] o: _ 5 [. U' z2 `# s, ]9 g

! J7 x+ v3 B% o F) c( s/ N0 v( K

{ {1 @9 o3 G! I

int n = A.rows(), m = b.cols(); 8 ]/ z3 U5 s$ }0 d( ?: ] # T& f( D# L7 A# a# m1 z ) G8 ~) C# r; C+ v

- j/ [9 J! D& p; {; z

int irow, icol; 6 c `4 N7 U; N! g% E# A% G 2 p) @: [ u, b" ~

7 w1 S6 Y6 h* F( H6 O

vector<int> indexcol(n), indexrow(n), piv(n); ) h8 W- `% e: r; c; S) h" n. f

, U" E* k& l: V1 \ t

/ k5 `$ p1 ?& n& @: S& j

: F& R# v9 n# t9 F$ J. x" p3 E

. R) v$ Q; F+ B6 e' O- z4 K

! w R+ o9 I9 o; H; O3 y

for (int j=0; j<n; ++j) ) {7 M' k R+ x) i! O 0 `, k, l' ^+ t- s" G, K. Z: j3 o

5 i" W7 t% Y4 q7 n( h( T4 _1 z

piv.at(j) = 0; . l1 p" F, F$ n( G

, i1 W0 J4 Z6 H' u# n

6 a, K/ x& Q- X1 j% C( P

( ?! j2 n2 N9 I" I3 j

//寻找绝对值最大的元素作为主元 ; f+ c# e2 A# ^9 `9 _( m# q6 Z

0 H0 A/ G1 Q$ ~0 l; M* J8 G& J& ]

for (int i=0; i<n; ++i) { 8 P$ Q2 j' G* H# \ 1 K8 D9 r% d, ]

6 ^. K0 r3 |& Q

double big = 0.0; + n+ b6 n" j+ n' f ( ~- ]0 u b& r0 B/ l( K9 {2 i

, D! O1 W% I; K% \1 ^! D

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

! n; |) j, z2 i4 E: \

- s `8 G: b# d6 c( {

if (piv.at(j) != 1) 7 u- o) q+ n3 `( v1 d: g/ {/ p* N6 F , e' b$ ^6 a. q1 s- H

for (int k=0; k<n; ++k) { , X5 N: {" `/ L: d$ x6 P& ], P " G9 Y2 N3 Y' m1 `! s

% _; |3 z( A! Y) [3 z! K

if (piv.at(k) == 0) { " s( W4 j. s! H" \& X2 k. V

, }! u5 V! T; x6 I! X

1 W: A* Z1 {& u( W: o# ~- g

if (abs(A(j, k)) >= big) { ; q, C, }8 H" ? P2 F1 _ * ^1 a; v; A; k

big = abs(A(j, k)); 2 g" b6 J* {3 H( ?; i* l : `. n7 }6 a& I" y4 i( z

irow = j; % `, z" p7 f& ]" c9 G: i

icol = k; & \0 Z) | \! o& P1 H 0 D8 k/ |. E& P+ e8 A8 ] # v k) l! z2 w; `; ^% X

if (irow == icol) break; 1 m$ J* l, y e, @% X " X/ d2 |0 J6 G" C3 x9 i2 C

% X5 `; }5 i) [1 `9 ]6 G, [; [" m

} ) @+ H' U' n+ |# I0 A: X + O& ?. u6 s1 _& K1 O0 F

; j' b* b H, W3 F5 X3 l

}

} : `( w, O! Y3 Q2 |) }; z3 x0 p 3 m' h2 N9 B. z8 d; b0 i

; `* z/ X7 k8 m$ Y, @; S

1 ?) K. T; b) ~% ~) d- |1 |2 p% M

++piv.at(icol); 3 Q( F$ @! M4 {& t% u+ V8 N $ ?9 s3 A( P6 m% v! q) k L! t, y7 Z6 X # f x2 o, ^4 Q3 @6 T/ s

- W& ^5 E5 H. L9 J: h5 l' [. G

. \! M9 i+ Y' v( L 7 w# Y% n4 \- _% Q3 c. R

: ~4 F3 V. _4 E# z: o( \

//进行行交换，把主元放在对角线位置上，列进行假交换， 7 a7 _( Y- X7 p" e + W- r7 r# h i6 |( B1 x# E

; v2 O* c, ~$ S- f$ s% W

# b& ?& u" y' F5 Q' x

//使用向量indexrow和indexcol记录主元位置， ! k9 l7 [6 `. L& s" Q0 K

//这样就可以得到最终次序是正确的解向量。 + X, B" n, ?. r3 G * z) w' F* Y$ G7 Q

if (irow != icol) { ( a9 W& u, L6 g" r7 b

+ ~& z* E- X& Y `

for (int l=0; l<n; ++l) 3 m2 k: r6 b3 `0 U0 ^

( i9 z: q4 \6 T9 B/ U. @) R9 C1 [* K

) w% g( L) w, G9 O1 C& ~3 b/ [6 U

swap(A(irow, l), A(icol, l)); & G7 q2 o0 w1 n* ~4 k7 {+ N 2 W3 _% ~, s0 x1 c L- B/ i* L6 S

2 T& D/ }9 T; V4 |7 i& X" ]. \

1 Y/ s0 S, @# ]9 g' I

0 x, `& m. M: P! K

for (int l=0; l<m; ++l) 3 _1 _1 i1 u- _7 W. B+ B% S

" L& N- F1 {! j M! o5 t' D9 r

swap(b(irow, l), b(icol, l)); $ J$ r# ]+ d8 }5 h. d9 Z T - C3 D/ ~; D7 h+ \1 M3 L

( `7 ]1 G( S7 E( b7 R$ m2 H

} + A" W1 ?8 @4 T5 `, b7 Z

0 o z# S' L# C

indexrow.at(i) = irow; / U# p0 u; c9 F* h/ H+ s2 V7 Q 3 i" y, B1 R4 h

8 {& C7 u2 v0 h5 `; n) h

+ I- ]4 Y& S: ~ m' ~3 F

indexcol.at(i) = icol; , F& p3 z7 i$ W6 P) ?( q1 p

try { 3 m/ n- ]/ x7 V/ T1 i- c 6 m9 ]6 X: ]! ?9 k

5 f5 ~2 n; d2 w& n

+ }( k9 P% f% B. n+ c9 N

double pivinv = 1.0 / A(icol, icol); ! @* o1 Y2 ~+ t$ D

/ S3 q, k' v/ d

( @% _3 l- m. ]# W/ j

for (int l=0; l<n; ++l) 2 P$ t3 i& m \' j/ A ) Z! Z. Q1 ~7 @' L2 S

( {! D5 H- R4 b% D' I

A(icol, l) *= pivinv; 9 D8 P8 \% Y" z1 @6 u! o1 Q. f

7 z, g) x) z( D2 v. d) X1 c

for (int l=0; l<m; ++l) + S: n) ?' ?# h. X3 w R% D 8 x% t5 v' p6 ?+ S 2 C+ y3 O( W9 T3 y

" A q) l9 y5 H+ f

b(icol, l) *= pivinv; ! |* x1 Y9 {; s B& F. l) `; A ' I5 d( S# D# \+ \! \% s

' U8 r5 Q! m( f

2 B6 ?1 {$ U! |; @7 ^3 q

4 r- |, S/ D% `

//进行行约化 $ \( d0 _/ |* Q. A% v8 P' V8 v- T

for (int ll=0; ll<n; ++ll) [: {8 ^8 h* p+ l+ @# {8 ] ^" l8 t9 s' B2 L6 z

7 P- h8 Q( \9 j K# Q: f

% J* w3 S- N1 |; B9 r) L& M

if (ll != icol) { 0 x, f1 Q7 n6 `5 V& {. r 1 E( R$ C3 \, S4 t

/ c& O# U% J' K1 ^1 J+ [

double dum = A(ll, icol); 0 d9 L+ ]2 C: U6 i" _: P' L4 |0 p! @ 2 D8 R- A1 V: @ C* x. Y/ A/ }: q) \ 9 D. Y; e& L- V( E: B6 y

4 Q0 E% V* N* F3 Y& n* z; f. c

) l" i2 d8 v) l; [ V7 u. n4 h! H+ p! s

for (int l=0; l<n; ++l) 8 i) Z5 I* l, }0 W9 u& i" b: P

% R; T9 W6 Z& I

A(ll, l) -= A(icol, l)*dum; 1 c8 z2 E. `8 r8 z5 ]. P, N5 P& [8 W " `# |7 o# t) n6 n- V

5 r j5 ^& }) o. Q+ w

for (int l=0; l<m; ++l) - J/ |. t7 ^9 _) {1 c; n/ r . W9 m" H+ y& I1 |7 a" `5 a

. D1 {' U; R$ }3 g; p/ G

6 B* v6 b" [9 T; O( T! R

b(ll, l) -= b(icol, l)*dum;

: {, ^! Z2 y+ E& ]& e! ?9 M

' e. R- O6 H4 x

} 0 Q9 G; v7 H7 M* s# X4 | 8 ? ^ N v- H' |. z; Q ; n1 D0 m3 e \2 U- I* c& x

2 T& D; V# a+ ?+ @/ J+ k% C

} ( @" Q. {2 [, y, p 9 D4 x# Z8 m6 i" Q( r7 C! A

3 S6 x% o" u3 Z

9 F# ^" t; K/ s9 A% q1 t

catch (...) { / r9 Z8 K! f" m% C1 p 2 l. ^5 Y' s3 t9 V* t/ M' E + a ?5 d# ?% w. B2 m

cerr << "Singular Matrix"; ' G$ A7 V7 N3 `. A- D' R9 i) w/ Z

} # P' a8 Y& |. a7 E 9 E: e# ?& M& R/ D# B

+ C* \; |* [. k* o- O3 S ~, d, \7 a

6 q0 Y6 V$ e4 ?- k; G5 F% l

} " G6 ^/ x+ q& {* k4 N+ K 6 S* p! M9 c; n- z5 a( q

! v1 ]6 M6 P9 c( s; j: R$ n6 v

} - U. F- q6 w2 H & z/ h' D8 l" b2 p' M0 C

" j% Z( J2 w7 u7 a4 A

0 G0 x1 U$ y& m* s( p. r

8 Q* f. \; M3 i( N5 c

0 k( C4 u1 L: F+ |& [- J) L7 W8 |: W

int main() * P1 R( O1 f3 h! t- q , A+ N3 P4 ]3 H/ K! R5 P' ]$ m/ b w

{ % ?: V' S3 f2 W- w- f 0 r9 T/ P* S* B# B3 F) a4 W

( I( J* u6 L) D3 A7 u) L3 d

//测试矩阵 8 p% I% E5 y: ]7 E4 v

0 W5 `+ i' ]2 a- @$ | Z9 D

Array<double, 2> A(3,3), b(3,1); / ^1 M; R% y4 n6 R, G5 [, u

6 i9 S2 b8 j, \- o+ G

A = 10,-19,-2, 2 k. u% U3 F7 T! j8 o4 X

j D) E8 {" o6 v

-20, 40, 1,

5 t9 f( p" l" b4 u

- G- X: H3 @2 y/ G

1, 4, 5; ( e- `7 p( t; j! U" H1 G$ U8 W

8 [0 |' g# ]( m: q& Z) _1 ^3 o- @1 Y" n

, A: @6 }" W9 u" x8 l2 {& ?

: E3 X# v$ v/ }- V! q- o

% q. J- c/ `9 o* K- Z

b = 3, # S$ w* M( P$ X9 }+ F 5 R) o3 |1 L: s& b

4, ) j' g8 m- C, Q+ Q- T

5; v u3 \' H: {6 x+ t " E5 E* U" p8 e

& W% u# K+ A/ q8 R" ?1 H

) k3 V$ D% l/ C- O ' }+ S% V$ `3 I V ) v. K& ]! Z: j

J) z" f" x4 N; O6 u* E( z8 R; t

8 i+ |- I. w& F9 a0 s! l& r

Gauss_Jordan(A, b);

9 L; k7 K% R; Q, f/ F% L7 ^. H y# c

; Y7 J5 {- U' O6 p2 ~( v5 c- R! y

q- @. }- B1 Y- W# S

cout << "Solution = " << b <<endl; 6 J$ D& [$ z9 S* F; G + ^# O, v) L# Y" R( p

9 g* ^9 D. S/ Z- q

} 4 F5 z0 S$ Z O! a8 k % L$ e" P8 h: |9 H! k e

" @6 h! q$ F6 \5 V3 O

( }$ C% l, Y/ i

Result： 2 P7 F* t' c U9 P 9 P0 p( @) {* ?8 V4 k 5 o N6 L: ?4 D2 _

) j- i6 R+ x7 _! A9 j1 n5 N1 X

' A6 ?4 u8 J. @/ x5 |

Solution = 3 x 1 * r: D' @! M0 l0 d) O5 B3 s% H * Y/ T6 Y1 |# X

9 c3 I' R$ {$ Q

[ 4.41637 6 B* h- ?- |; B6 T5 `' O - s [) K; z: q ( w( I% o5 ~9 \. q2 N

9 Y7 Z9 H: \) L: U( W

2.35231 5 p' O% d& c7 b2 `5 g. {

' }" G+ M; Q7 Q+ ~* L/ Q

, \) V9 D5 ]9 F8 A9 | X9 f

-1.76512 ]

( N3 [, M/ o7 y: b0 h: w: [

+ a8 J6 k8 O- G1 i W& u& O

* Q' R/ D# P1 W( a, ?

, l3 x7 p6 |( p3 c5 v- A) E

0 I" \) R% f9 l4 O

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

' C- F' I) D/ J$ G" H j

' h- F( u( e3 J: ]% ^- T9 n" y

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

2 ^4 o+ }5 T: N- G: |

0 M4 `% n K' {+ p

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

路过。。。