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

lckboy

全主元Gauss-Jordan消元法

& O9 Z4 S4 _4 J( Q" E* S9 ~

6 E0 A3 \0 e. ]

& L7 } Q" Y) k9 [

0 Y G( [8 n, [# x: \) {' e

! u8 w, [. x" d- P

/ x/ F4 d3 F2 O

. t7 A Q( ^* m& Y+ u9 O5 r; E

* {. M- E" u0 u

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

. O0 y ?3 G& o$ d; ~/ Z

/ J; g' i+ Q& w% U8 O7 u

' d3 G$ D! l D5 A8 r; V

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

( A' e2 m+ }/ ^% ~( Q

8 ~( Z" {+ D* L* o5 r

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。 ; N0 S% R# \- g* J

' A0 {$ x2 G$ L/ X! u

" q; |+ C6 U! ^) u9 b

/ g4 {: Q# Y7 p @! I3 @; P1 r

Code： 0 J+ o/ z; i! p# c4 Z7 s

: y' _. X% V6 Q! c6 d' {7 ^: ~

. x1 a5 `+ E; w- K, C" A s* v

#include <blitz/array.h> 0 k* H8 G2 [$ K . K2 D5 D8 h8 t( f3 M8 k

a8 ]& e# k5 b$ P1 v

#include <cstdlib> 6 m3 X2 n I) m 7 {; n9 E, T# D o: o7 p; w+ ]: f! u

+ ^ i {) G1 R0 `: W: A" M

, L' Y l9 r0 n }% V# {! s5 \

#include <algorithm> ) o4 Q' s( @% @$ ]' f5 S

& T# g8 l7 m5 ~0 j5 R: d x# k1 m

#include <vector> 7 \; N3 r3 l/ @; b: w$ d* M

using namespace blitz; ( E H; ?- t2 M" L+ z" c, C) x% c

- }) e0 i8 q5 o7 |$ J4 B- @

- ?, f1 G; ]+ g3 m |$ H- o

) x2 T8 ^ g/ X) S' f. f7 S

" a+ J6 v! J! N4 U9 U# R1 k0 [: z

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) % J8 }# `$ |0 @5 g& W4 ~/ N2 q ; P( ]1 h; K2 @

: W* @# ?: S/ m2 z$ S1 t

{

6 N( @. e3 z2 Z& p" s; o Q* l

int n = A.rows(), m = b.cols(); % ]$ P- a4 N# L: I F& i" z \, s v/ y6 i2 c) M0 W

# }; U) e; n( a! L+ T

int irow, icol; : [! X6 i7 m4 [- t' |' y $ b7 d7 X5 U1 J' J+ G' S

vector<int> indexcol(n), indexrow(n), piv(n); p. w d' l9 }+ K% _ / y, a+ K8 l& w+ Q

+ R* V; m9 ~3 A, f3 M

for (int j=0; j<n; ++j) / I# y) o1 B0 w

/ U" i- V; J( X+ H% B. l

piv.at(j) = 0; " |' x6 l2 B# F2 C 0 b0 r1 F: X7 `7 Q1 w

6 y9 `7 o4 z) e: \- `

//寻找绝对值最大的元素作为主元 ! q* o* ^5 N. z) u+ R ( z3 c: s9 ]% V! \8 H9 q * J o; E. Q' L) d( ]

* a& n6 i5 R) x+ E, `9 N; M U+ L

for (int i=0; i<n; ++i) { # s3 E( S5 ~+ g/ S% Y & B: X1 X: P- H8 W) o8 R0 Q 2 ~: ?# i# W0 y9 D2 @, H

9 d4 p `# B. C2 u

# `6 F+ D8 v# v3 x3 Q, `

double big = 0.0; ) Y8 v! q ]! T6 X# K: ` * B. L4 V: |7 c+ V* N$ O

4 O- z. S- c% r- C/ X

for (int j=0; j<n; ++j) % d4 q0 H/ d, I) P8 R1 m : ], I$ o) c: t& M' m" i

+ N0 C; N& L7 x# }2 M

if (piv.at(j) != 1) ; s% N# j1 k: _8 @- i

for (int k=0; k<n; ++k) { 4 m9 p* d# u, l: q. p / q) p6 b( V' b" {1 o

. L" D( ~2 K3 [+ |! \7 [+ i$ ~# U

if (piv.at(k) == 0) { : t" L$ t: L- b! W: Y8 e ; z( B# q C/ ?, S: }6 y- Y1 R 2 ?+ B; T8 D' l7 u1 v3 D. Z# d

# t" b" a& y8 Q. S# ?

/ d$ P' D4 }0 d2 ~. x2 `

if (abs(A(j, k)) >= big) { 3 |5 \, Q2 k' ]* V; C

7 c8 |2 a/ M1 ]! d* P. u8 ^9 s

) s% q6 a" ? c c. ]$ X v

big = abs(A(j, k)); 2 \( Y1 ?- x; E2 l# ?2 p ; \( X# p$ e# \: J) U

, f( F7 L; a7 s) @

9 `4 N2 F0 b5 ]; ?. |- A" h4 A

irow = j; 6 Z4 u+ [8 R4 r* b* v O" Z" M ! C: N) S) z$ ^6 H$ B* L2 y

3 ^( ?0 d$ G( Z j% l' N

icol = k; % b* v7 r5 M, _/ W* |, r , c% {. Z0 m* d; b7 R% s

1 V j" g) Y3 p) O

3 ^5 \+ l& Y( r; E

if (irow == icol) break; ! U/ W* F; Z" }( f ^$ t' j. ?5 Q 5 A0 D# N/ P, J8 {1 f2 f$ N

. J M8 E) G2 n( s: l4 X

} ( ~. E2 B, Q( ]6 ^ 6 y" i$ N8 v1 o/ b* K) i

} 1 W4 R8 D; r6 Q2 L. a% p , f l, g- o2 j6 I5 C% H2 C0 w

0 h3 c$ J' R& z- y( O- |

} 5 i" M1 g+ O! T3 A; d - Q5 l4 R6 U+ V* |! Q

: H1 s% D. c# Z

# ~4 b9 n+ T8 D1 |

++piv.at(icol); + ]( M; t8 A$ u4 C9 i- o ' Z3 j9 j1 J0 W7 M ?/ L% R b! h% l3 ^

5 w- |' u6 Y' Z- [, ^- h' s9 @

+ W% r; G6 r$ v* E: G' h

9 A! e" A/ h" f9 l

//进行行交换，把主元放在对角线位置上，列进行假交换， ) d9 h5 N; V4 p3 T* x

, e. y T1 O4 v5 U3 F

( W1 w& J2 H0 U N* ?% Z

//使用向量indexrow和indexcol记录主元位置， 9 `0 H3 j5 Z/ X& u' j( c

8 Z* ]' u" d+ v: D, F7 H

7 Q- r. ~5 @8 C

//这样就可以得到最终次序是正确的解向量。 ( w- U0 S) }5 {, s+ k ; i0 c! v' w0 L; F

if (irow != icol) { # D5 Q; P. j1 ?: f8 L0 u ' p0 |+ J% o! D. ~* L. y6 @ 5 Y8 F9 B n" d! P8 R

for (int l=0; l<n; ++l) # _% K: ^# o+ r! V+ m2 w$ O- L! r. Y 3 _$ z! A, q8 i. M. A0 b& N

6 V$ m; t3 ^1 P# {% O9 Q7 S

swap(A(irow, l), A(icol, l)); - T/ P0 q6 B, V% ^ 8 ]( |' X. f# _, g9 ^/ a* {8 Z

) v$ y' c# P X; U% Q

1 s1 t, B% t/ a% H: |+ ^# q, m* P

for (int l=0; l<m; ++l) 7 `5 L) [- S( H% _. Y1 ~ 6 {% F% d7 X% B+ L% u

$ ]4 K; a5 |2 @3 p( Y+ L

swap(b(irow, l), b(icol, l)); 5 @9 o1 n, k. F: ^, n5 c

, ?, y/ m5 M7 ?& }3 D) ^. P$ A

} 0 ~; h: F+ A, T9 v" r/ K* k * N D! W# x% E6 a$ c( L- j + D. u8 w2 Q: k0 \7 B( G

b, {' Z, f \+ |8 H# r5 P% w9 Q

0 {; |* q$ | o; G+ b% T

) |+ C+ f, n% `2 A

1 |( K" B; K+ |) [

indexrow.at(i) = irow; 5 m8 I; s: i- I& P* X5 R' b

8 I' {3 B8 ?7 g3 L3 J+ Y/ }

indexcol.at(i) = icol; : d, y8 y4 k4 O& M ` / @. o+ @) D0 L- v( S+ h, p' Q* C6 P

( J# a/ H; E- _7 f

$ x+ ?2 C9 @3 F- f) C2 ? " T6 K' \& }( F' { % A' {8 T& _- g0 V& h/ Y

$ B6 M5 d; G# \8 ^9 M! v+ ]7 a

% E# P, E4 `* h! S& l

try { 0 J7 ^ p) I+ u# g& I5 k3 C ) M: W9 @* r+ ^0 Y/ @

8 o8 e S" e" r L; H/ X6 ^& V

double pivinv = 1.0 / A(icol, icol); - g y. Q' |% U V. ?3 w. v

for (int l=0; l<n; ++l) 4 J! J* [! `* c

- s1 ~ x R1 O) I* V# h5 E/ x, C

A(icol, l) *= pivinv;

( t+ o y2 A+ p! x9 O! G

7 h* z/ e2 |; h# D* C: A& x

for (int l=0; l<m; ++l) 3 ^' P5 V9 y+ o i9 c/ H. D8 f) f6 j 3 l5 H/ f0 E1 o! M0 Q7 p

b(icol, l) *= pivinv;

* x* K9 w# l4 l, C# d3 h0 o2 t

6 u+ a7 ? F# m/ S

//进行行约化 ' W/ W' k- J: r7 R& I

) A( o6 _5 a/ ^( y/ B* ?

7 J6 n* A: J3 H! c) B$ ^7 z: U

for (int ll=0; ll<n; ++ll) 8 @) R* b# a0 |4 z. h5 T

7 n$ z( |% }5 d% Y2 Y

if (ll != icol) { , d2 T4 e5 h8 Q7 s) j* @& d% [ 8 [& H/ s4 j# D$ t- {" e & j; e4 @; k" v M8 @" R$ W( _

3 d9 ]" y6 E9 X+ J' R

double dum = A(ll, icol); , r9 ?: Z1 U6 p9 o: q, J

! s% ?5 B- i; f; _4 X* k

9 p3 w; w' B3 T, B1 F$ }. A

' d- H" l) F, ~8 }& Q1 x# t2 L! z, T

6 \" L: Z- r/ V$ x2 S

for (int l=0; l<n; ++l) ( @* |- s& s' @' ?* l8 ^9 q * S$ u1 k, N( e l + l8 [3 _' Y, F3 J

% w) ?, k3 r1 V. L1 G' U9 ^

9 ~9 s6 M- R- N0 p3 U

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

for (int l=0; l<m; ++l) 9 M) Q; U0 a5 |' ^ ~& i% k

b(ll, l) -= b(icol, l)*dum; 9 P1 t$ B# s" ]; q# Z

" E9 x3 ~" O1 p! i8 s# m3 S

} + H1 U ?! a) O

1 r/ k" Q+ C$ V! c! c

} 5 [, ? H+ }3 x" N G* ]8 H* }" ?7 y" c % y2 I3 C) ]+ c7 B! a! D; Y

2 K" u5 L1 R3 T2 [( D

* {# N5 g* y3 M0 v7 W

catch (...) { & T6 n8 T6 S, d+ Z1 `# E v

' ~% l5 e. N. \' [ _" C/ K/ X

cerr << "Singular Matrix"; ; l( a& {3 e* f

/ Q/ ]9 }# v: B& q$ x0 p7 H/ z

} ( i2 H* I0 I5 I3 ? 6 l) L1 e ?" [8 t+ z7 U

L' x/ P& N) |# V+ k

" ^; ~: K' C5 h% }7 H

} , g* l' d2 N! n- Q r & q9 Z3 f% D3 M# a

) }* \ a. J4 B) n8 T9 v

} & f( e: ~" j z( C% I4 ^/ I0 E

+ S$ q( `6 v* o& O& j

+ Y, x' V* S0 ?% c" h$ A1 c

2 d; z' I3 n+ D1 P/ f$ \# {

% E3 v1 i: ]& p) ^& ^6 q

int main() 6 u" a: S) r- d ( U+ U4 d6 t8 \2 n2 A7 R+ \5 T

( t$ }2 Z- n8 t2 P: P! T

{ ; O0 n I# h( i, s! t* ^: @# n 1 W" C; }6 s! {) o2 z; o5 ?( Z / i- I2 a) [% n+ z

* S& b; `2 c s4 J6 @/ d' T

//测试矩阵

# d( z( B0 l- G Q! M/ b1 j

Array<double, 2> A(3,3), b(3,1); 9 V% H) g N# E! Q/ I5 m; S) c

9 S0 Q, k& z" i) g$ ~/ [, L

A = 10,-19,-2,

2 U8 O7 G. `* o

-20, 40, 1, Y& L7 W" b$ g) |; i, W$ v( P

: q6 t9 V; w9 \2 j1 I6 k5 G, o

1, 4, 5; 4 L' Q2 {0 I- Y

) ~3 P/ H1 ` H0 P( E

% z. s' }/ d3 W& D2 ~) F

- L5 c/ E8 [; [, C/ r" G+ _8 y

' b) R% h t4 N- o! {/ L

b = 3, ( _# @& P% F1 i V 1 w! A. @9 R$ G" ?% m* X ( H8 y$ R, s( |, s% I. ]/ r

9 h. @0 `( O& @; p5 ~6 Z }

4, 8 [6 S! t# ~6 S+ Q0 k/ `( {& |

. O+ Q% m$ n$ O1 @! c

5 ~9 a+ _3 t6 h& [, _8 w

5;

) ]" m, ]& e. L6 Y

~2 \! ?$ e0 ^ V$ \ " m0 b4 |$ ?( r' c8 ?) L: l

% I( X- }' Z' M ?8 i$ X6 z

" [' B7 ^6 ^, R

Gauss_Jordan(A, b); 9 ^5 A+ ?9 B- t9 c5 ^

/ H& [- M, a/ \: z' i U

) ?/ Y/ Z. I0 T( e$ y2 V

) h/ z' U0 l3 U* h( Z* V' d' O9 t

cout << "Solution = " << b <<endl; 9 N! c. C! O9 H. n' |9 H& W $ u6 X/ V2 Q" A0 n

} 7 \4 V. ?0 X* ~9 b* |" O6 D / p3 Q" i2 s& b! E

- d, L# U# q& E1 h4 e/ q

( O$ w/ Y+ f$ l- Y! A) G

1 B/ `% d( t( r, I% S

/ _; U1 `# _4 n- l2 t$ C! |

Result： 0 g9 g q1 s& l% E# R 1 C+ p3 l. C4 X2 M/ }

3 J. F* r, m% U3 Y+ u2 W7 I5 G) H$ m

' Z$ f& U/ ^4 F" k

5 _7 U* g5 a3 F1 [: M

Solution = 3 x 1 % _6 `+ i( F0 ~. v% L

* m3 u" z1 A/ E! H( O

- Q8 ]: b# y: x7 w/ p; F8 V. J+ A

[ 4.41637 : G; \1 F& f& h1 u, b

0 S7 | j' _0 \+ {

2.35231 * j1 l$ W3 p9 _ 3 ] X, ?( w3 Y2 y' d U ~- {5 H: L & H6 g0 b" Z3 m# O

-1.76512 ]

8 B/ ?+ i- O- c( T: j0 v3 z

& _5 P# r' B* b! T/ p% E1 k# O/ Z

r' N2 D$ d1 v3 ^) D! N/ Y. W

# t8 p- P6 R. T0 a

2 u% |7 `. x9 H7 Y6 ?# M

) a, X" G+ i8 }4 E1 m5 J1 c+ L

; z) x, G$ _5 ^! \

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

4 c% N" f2 y, J% T" f( I( w

* }( g. `' p$ q: z, L' `

% u% Z1 P( [3 S

( o$ ?0 O" q' }( {

( I% Q% o6 o6 M1 G1 ?

注释：[1]主元，又叫主元素，指用作除数的元素 / b9 s" H2 ~- j0 \8 C# U5 i

# l2 a! C: L( d' y" A" E3 }3 l

$ P0 B5 z6 X* L& H' b 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

路过。。。