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

lckboy

全主元Gauss-Jordan消元法 & _' l3 N! }7 f- J; w

( Y6 F" _5 V; m) U3 f: d

: X6 g: u' n8 R, L+ `

9 n ^* Q2 _7 @" j- L2 Y: x7 ~( `

6 D B" }+ U6 }+ ^) V& O s

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

* l% Y1 O$ H, ^; j7 M5 K

( p& o# \! g. l# [

9 x6 [* A8 z `5 D) u; \) ]

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

0 K$ g! g9 F" z: N% D" P/ s$ o

; y( e; V3 X/ C0 Z2 y1 o

4 z2 L, o( r- K" g

# V) M( m Z2 f o: J$ V0 ^$ `

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

' o# L( l( |) g& }

; w9 n0 T& m/ c0 g/ }

Code： & C+ f7 S) n6 V8 J. }+ L+ }

+ Y# P" L1 x; K0 ^- e5 l

1 s* z6 u, a7 M+ r; ]$ ?5 x

# }: u1 m5 K2 Q0 Y+ `

5 a# Q" }( V7 N# j

#include <blitz/array.h>

7 l# e7 A9 b9 ?, c8 y; X' i% H

#include <cstdlib> - ^, C. B: S j' z, D# h & j! l# F4 ]; L- x1 z ( Y& d' S1 b$ T0 [

#include <algorithm> 7 @& t4 ?7 e: w) E0 U% b1 J6 @ . R& H1 e" V4 K8 b: z- K( U

: }6 p, y7 i* a) c$ T8 @

#include <vector> ; q- x \+ u5 d! z8 y( Z $ ^% C4 _( C) @1 a0 P5 F

1 e' F( @3 Z" B, R

' p3 d. i6 B# [8 a, Z

using namespace blitz; & V8 x. t, k' w& R- \# h e0 l7 W7 m& N

" r# O) n. q: E8 q% U5 }3 j$ v3 F

. U0 l9 H% |, A) N& t

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) & g6 V+ m( C9 `4 K3 r- @: p

. \1 Z( K* X( e3 I

{

/ r5 @# Q: _6 J; c2 J

int n = A.rows(), m = b.cols(); ) O* ]# E7 j* C& ^( O. l) c8 L# a

; v3 [, [. J9 i( L, s

int irow, icol; 3 X' O$ U! f6 C* |: c9 } / X3 d- C- ~/ V* H2 b, m8 T# H

vector<int> indexcol(n), indexrow(n), piv(n); + C6 l7 I% A( i* k9 r % ^$ B$ z, k" c- B8 A& N

' H* a/ P; O* E( }1 w

7 s% S7 A& q( z# V- i0 Q

for (int j=0; j<n; ++j) 3 e& P+ N9 C* e( u, g6 u$ p

2 }4 c( |, a; ` {7 w9 }$ a$ U1 X8 V

piv.at(j) = 0; 4 |( N2 w' D$ R: @ 5 E- \- X2 N% p6 [* @ ! Q5 ?7 v S' S- o- e1 V: `* n) w" T5 x- i

% G5 z8 M$ k$ ?9 v2 s0 U

$ k/ ]" N3 ?6 U1 Y+ ?$ @" R3 m ) b7 _- {( K$ p2 l8 a1 X$ @ 9 K/ ] g7 B, N, J8 F

% c' w% m8 w C$ Q/ a2 t& @

//寻找绝对值最大的元素作为主元 9 Z" @1 g4 F; Q! l) a) W {0 I" r* a6 \8 b9 d " B2 @. T3 N M' d0 D

* J. d3 e5 w5 L: b( v$ V- R. v

for (int i=0; i<n; ++i) { . j% y( v( j: o/ s+ d8 c' \* r. H5 n+ d . I* f, R- ?2 V% u8 l7 b8 x2 m% Z

double big = 0.0; ( Q$ c' A# l' r( H5 e9 ^ b : t/ v6 {+ q( T9 X5 p1 F ! T& e) }. ?9 W: P5 x+ V/ W: q

4 y% ]. D9 o6 {7 l3 k/ b

; F0 h* a. F. Z# `

- a0 D2 B. h. B% w, G

) d' L9 G; r0 }9 L9 j. k* i% D$ W9 E

for (int j=0; j<n; ++j) 0 b+ {; d( P/ ^; S

if (piv.at(j) != 1) , m7 e5 o& j6 o- a( P& F# G9 h a+ c 9 N, b8 _! ?; ~* W1 e& C

for (int k=0; k<n; ++k) { 0 m; f$ G) A9 U 1 E- j$ K1 T! D% d8 ?8 w

7 S2 r) i' x" h+ V5 F3 W$ o

8 X# \/ J! |9 Q- m' I1 E! j

if (piv.at(k) == 0) { ! e( s& W( p6 ]' s$ d 7 h& b8 Y' L6 D

if (abs(A(j, k)) >= big) {

L" u2 \' b# F, k* q* d

7 w; H$ a4 A9 ?8 ~( r& H5 i

big = abs(A(j, k)); 7 m9 X; F/ u+ k0 }, A6 @! ` 6 D# {: _6 U4 s8 L4 C1 M+ X8 \+ E; h' d

irow = j; / b, \3 W& F+ j3 {8 e ) l6 Q' W/ Q2 S: _8 @& H- P

3 e ~4 I3 C$ Z# |8 a q& ~

icol = k;

9 X$ Y% F6 F) y0 _7 e2 P" S

" p' |1 m8 L0 u* p" _' s& x: U, W

if (irow == icol) break; 2 s5 y' i" q& Y6 \1 x* I

8 D" v! ]0 V; D/ V

" f: K, p# _6 o0 J* v

} 2 b- e9 w8 Y0 ^7 F 6 M7 S7 x5 ]# d) N4 o

) f; w* R" ~2 _% p: g

! x4 \* y4 b' S! w. J& T( K, q/ {

} ; s$ K) X/ R+ u! X3 F7 s! ?' \ $ _! {- i6 A! j9 C0 G' }. A 0 d |( N; B. q2 s( I

5 m0 \% y* }% P* x& ~( K5 [

} 1 U4 ?1 q" B- D, ~

; i+ @* E2 }* t, v, ] |

0 _- N& a0 S$ d- a9 z

++piv.at(icol); ; {! C! N4 M! Y% p3 @4 W2 t0 j 0 x7 v3 A2 s$ c' c B1 }* j. F# B 1 {( b3 y6 d4 s0 d

" t. N v+ ~8 X, V0 r# M U

3 L1 I& v1 t4 h+ e1 G* R7 L

2 w1 O0 `" ?) G# D

//进行行交换，把主元放在对角线位置上，列进行假交换， 2 ` X2 [7 Z+ {# ~# N/ o/ {

+ i' K" J3 V) F3 K4 m5 }

//使用向量indexrow和indexcol记录主元位置， " j( U5 g9 P* c6 h: v9 B- r% } : l M |8 I9 \- D8 d

: [8 E) R+ C8 Q

//这样就可以得到最终次序是正确的解向量。 " V; J5 ~8 c5 R. i; u, a 3 k3 h3 p: v2 H2 z! w' n

6 i( r1 Y0 V+ T

if (irow != icol) { & Z- {7 l7 E! c! {, J # [/ f, ^5 o- H4 s/ B* d

/ l9 @- ^: _, R3 X# |

W$ y4 B' Q" i% s* k/ X

for (int l=0; l<n; ++l) + u# S" J+ u( w" O

: a$ X4 A1 K9 @+ A

swap(A(irow, l), A(icol, l)); # l% y) v9 L9 }. @) B

% L3 X2 z& O# {0 C( R

) V$ P2 O4 I0 }4 m4 o8 E- `5 S

, i, ^. f) v/ l1 P" y, x q" [

for (int l=0; l<m; ++l) 0 q$ b9 C5 ~+ u" X, ^- j

swap(b(irow, l), b(icol, l)); : @( Q9 V3 `, z9 A+ n

) u i/ x# {$ p% Z5 Z" W0 |( ]. O

} 1 Z: L/ c9 ^; m/ f9 |

& B$ z; f U6 z1 j. q$ B3 r# W% ~

1 ?& S- c! [! [. K& s2 B

, L0 \& a$ d" J! [1 t8 E

4 `2 I& l) `7 o% @9 b2 i

indexrow.at(i) = irow; 3 o) U7 J4 q8 h- W7 \8 m6 d1 [. Z

2 ~% L. D+ ~9 p' b: e0 U

* H3 u4 v) L4 A9 @; o

indexcol.at(i) = icol; - |$ f, n+ B9 Y. x ' T& l$ D* z: l i! A

, g7 m% B2 X1 F& ]7 H

% j8 }0 j0 u" v) L4 T- W

% z% \: w0 M# j; q! B8 l; @- A8 _

try { 5 F6 C6 x$ b( Q0 \ M7 M 1 Q6 l( j/ c& Q: t- k ~3 X) h6 |# h ( w: z+ P! |( G; }. d5 c% E

2 U4 K' b q& c% O& ~ a

; R3 D& N: l1 V0 ?

double pivinv = 1.0 / A(icol, icol); 5 u" L2 K/ v5 n5 I

, g6 U% l; j3 a' e( W

# x! d' B) E& f) G6 z) G" v# a

for (int l=0; l<n; ++l) 5 x4 j- R/ K& K( P! E 5 W e; H. H5 h q& u) {+ u) G

2 z3 l% C- X- |3 K; B' H

A(icol, l) *= pivinv;

4 _7 v/ B* \: z$ D1 O

for (int l=0; l<m; ++l) 8 o p$ M1 }% M- ~! z5 i

s" e: \7 P: \5 ]

b(icol, l) *= pivinv; # W& ]0 F. p$ ]8 \% c

. R8 j6 ~/ f$ k* f+ U8 w

! q5 ~( G, s2 f5 X( s" Q9 z/ j: _

I7 I; E' d% y7 S+ Q; q0 \

( b8 x. ]: K! s0 q

//进行行约化 ; o2 d) M7 e9 W ]% Y# {! f6 ~

* B- U" a+ a6 U( J

for (int ll=0; ll<n; ++ll) 0 m4 X) N# \9 } s! P# s/ \ ( r0 y& t. g6 o5 ^* x$ ~ ) ^7 n4 s2 B; c

3 A' J* V# s& f) Z. }

if (ll != icol) {

% v% x* n7 C1 P h

double dum = A(ll, icol);

3 Y) |2 h% t0 J T$ g* b

9 j* j& g2 W% J' J8 R) U8 I

for (int l=0; l<n; ++l) - ] y1 k# u# {% ]( O6 P + b0 r8 C6 J8 y% {

2 |" g( Y* c0 [3 ]

A(ll, l) -= A(icol, l)*dum; 5 A W8 C' z! D- U/ L% x2 `

, t O0 u3 x8 ?1 b0 M: _

7 V' w4 t/ `4 ]9 B) r

for (int l=0; l<m; ++l) ) x7 ]9 u+ N5 F, u+ D 2 ]- ?$ Q6 Y4 n" K4 Q/ |

! O1 v& J) ^1 g4 V

b(ll, l) -= b(icol, l)*dum; * D+ j8 ~0 {4 ^/ Y2 b

} % v, e0 M: S6 G& \ ~+ b

: |' s* D8 c6 p% Y% ?& Q" W- r1 l

( ] L" G3 M1 h- L I1 R7 {

} 4 d2 M+ k& F6 k* r. N 0 \' o2 k3 ~- c1 q5 i8 p ( Y0 K4 u9 F2 T) n' [

0 v8 i5 a0 C; G; s& V, ?0 V/ h0 v

catch (...) { % r Q, U; v* i. V 9 M9 p- o+ e8 X) Z

* \' x( |: | v& ` ]

cerr << "Singular Matrix";

/ W; S# U+ o1 Y

3 E# y7 M$ \$ `' u# s2 }3 J/ h

} % l% C6 | J3 S5 g * k9 Q! p- \% o

, H, S" r) _, W

8 @( ~; T/ l6 y4 A0 Y; ?7 q

} 4 C% ]- F0 f# E5 q) J% h* F3 A7 D! L y( J/ D) f' Y; W

* |6 L/ A0 n2 u/ x3 ~5 E$ n

# O0 W, Q3 }; f( B# l" z8 s1 j

} % V+ ?( `7 s- z5 l: O: Q( ?. w, y& H 6 g" g! ~% J5 [0 M

; i- k$ o+ y d8 d+ w

int main() 6 Q+ p$ \. ^& O- J8 V9 v # M. g( c: R- l% E* X' y7 n

; e; t- s& F; F+ F4 w. d

{ . l+ o/ r: _. O; D, V0 O

//测试矩阵 & H" \8 C9 M& O! i

7 _2 i1 c! V% W& d

Array<double, 2> A(3,3), b(3,1); / V$ i; b& U- c4 T0 Z5 k2 ?

4 s) Z, E: b. q9 o, Z

8 X k( P9 h X6 D. h/ H8 U8 H$ n" v

A = 10,-19,-2, # N2 E' L9 J2 N( e Y7 r( m: Z) i; b ; ?8 C: c% m. x& r $ e" w' U8 L, {6 I, |

1 b# E# I) }' U$ D3 d1 Q, R

-20, 40, 1, ( n" n, b1 ^3 n( P$ j1 m, | % d/ Q1 X/ T: o! d g) o5 `! _ % z V& s2 w; o) v) |) e5 _

( w4 q4 a2 P% p2 M0 S% }6 s( V

/ T' j$ w- l, M/ }" @+ J

1, 4, 5; 7 F8 X; `0 N9 b8 Z& r / Y Z4 [- e) T0 J" d$ v6 k

* g, X$ @& @/ j8 ~% s) W, h

5 s. d4 N* q8 T {. S

8 \' l* G. R& v* O1 d4 }

b = 3, 3 d6 T i8 p+ A t+ M3 D" ]8 U# {2 m 3 v/ z3 P" J: _3 L0 B

' O* @" r7 d8 G2 q9 u- p

4, ( f4 `: j* I9 H9 {

5; 7 [' ?, l- d/ ^' g4 f& i

0 b+ h& g! V* M, x9 {$ I$ x- T. S

; T; \9 T0 K" v4 i

; i% x. x( a D& `( y9 e

$ H% ]3 A# ?9 {3 O

' a- k8 h1 e0 \! {6 F

Gauss_Jordan(A, b); ( ]/ n4 `3 M' P1 T, ~3 X

cout << "Solution = " << b <<endl; , Z) g) n b2 L6 g% S( F: R ' s. U* ]% ]# Q$ \% q: ~

; B$ L2 p) y# W* J+ @: h

} / U# `' p. O R

. b" l1 U9 X; w Y0 C8 @5 n

Result： 6 s1 |, m' [0 a' q8 F, s7 J3 V * W) {3 F) m- P 5 w9 a+ s- _; t8 u# {

" T. C) R+ L* U# l0 Q

Solution = 3 x 1 - M i3 L \' m ; L& J3 ^# D7 z7 p" N) h

: l6 C R6 P0 n5 A

[ 4.41637 ( V# F+ }! c( |0 |* b8 B 3 V( }: @" V0 z1 [" r A9 M

! q; ~7 }/ G* K2 x. ]2 Z/ \8 w) a- G" O

2.35231 4 Y; q1 y" R* I2 t0 f* J! P * H" y5 z7 @& ?7 M. E9 e4 L

\5 {$ ^; s) u- E; |

-1.76512 ]

$ R1 l: ^' E+ w6 C2 M8 N0 H6 Z

3 `( L X1 J+ F( t

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

0 a% ^9 o3 z' U5 _# k+ L2 g+ F0 s

* T4 [/ V+ Y9 g: b

6 z* n- |3 A/ [. N: A3 c8 [2 _" o

2 g% r$ W+ [( Y$ y

' z5 L0 U$ p* R8 Y9 a7 d) \

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

+ W2 u3 e5 M) @

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

路过。。。