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

lckboy

全主元Gauss-Jordan消元法 : p. X, j" T% P7 o- V0 n: l

) y6 I- b5 A# ?3 q i; ~5 P' | {1 ?- J

& D* f$ l' y: g* z3 I

( [+ g* a- ]: T( o

( V7 s$ i" p8 w/ k, l* o: `9 }

1 o8 i5 Y: a. ~) L/ g4 b5 q

3 f- N& D: J0 V" c& p+ Q, |

8 [- n4 S+ N. M% d! T6 T2 [ |9 v

Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法，该方法的优点是稳定，而全主元[1]法可能更加稳定，同时，它也存在着弱点——求解速度可能较慢。 + g0 r4 f$ N% R$ H* N, s/ I

9 v! n. n- R, _1 v- W& U

e. m. s& V/ b" u

9 M0 h3 W# i0 y

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

. Z6 H$ R8 f0 y- I

/ d3 I! p) j: N

( p/ u! M2 t; w: u, R

7 H5 X8 b$ ]1 B. \/ Y

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

; x( L1 T7 {, Z

Code： 4 I# h5 o" r. U7 Y& w( n# L

$ M- T/ L/ X2 i! H' Z

" V# _+ l' r- i Z7 `& D' \, V

9 W3 m S5 i. h; }* T. k) A* I

#include <blitz/array.h> + R5 S/ q. z2 q# S% T7 k! e

. k) N E) h+ v' F! t9 [; m

& Z9 a5 F; b% ]- i* F6 c: y

#include <cstdlib> ; S! ^1 q- ]+ X7 w7 N$ C

+ {1 E a1 R! O0 |9 o: c0 }# G

#include <algorithm> 4 N' i/ \ v3 V' s" g3 Y ) ?+ n; `, g" h9 _' M 2 S3 T. C5 p$ c! f, j2 T

8 ^7 J& s8 e) d) s5 e. \

#include <vector> 3 K6 M9 Z# r+ K" M; Y $ U* [& k. Q; ]4 T) M0 Q: h/ _9 r

using namespace blitz; ) N* N) H7 @, a) ^ 7 l: |- Z; ]( _3 P; _4 \ V) @

# v. v: A+ T/ c2 f7 g# @' B" q

! E" v+ N$ w; ]8 ]9 s

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) ' l: p F% M- t; z

. {" S* e/ R+ R% E

( f! P: U/ N+ |5 [

{ 1 r( u( L* q) w$ V2 h* q/ W6 S* x

9 e* d/ V, p: O# Q, T! A7 q8 |

int n = A.rows(), m = b.cols(); 2 c# B6 A. F' B ! J! r7 X k- S/ d% O

% ~: F5 d) l+ C$ \" u2 Z

int irow, icol; 1 b! I/ V3 [6 R4 N

7 [: a% f2 d9 K5 b6 I3 H

vector<int> indexcol(n), indexrow(n), piv(n); 5 ^$ S- ]8 b ~! h 4 N& w+ N' B; N& I2 Z

$ A9 Q% Z! Y8 t: I9 D

5 w h! [) w0 [8 |4 C' ~

4 G0 L* b3 I; M

for (int j=0; j<n; ++j) 7 i- K. ~3 `$ v e5 G) {, }% @

piv.at(j) = 0; 1 J9 @$ N8 c9 U

% `' B9 G5 N4 i, {

0 w4 H4 t1 p9 p( d9 T" E' Y7 |

//寻找绝对值最大的元素作为主元 ) U* D6 l/ e S& S5 ` " @& f2 U) z) ]; o4 z' x

5 e" k7 p9 e( `8 f( N- M

. k4 w5 D5 d/ P; e4 B

for (int i=0; i<n; ++i) { 7 H! I/ q; J/ }" q9 V( m g+ ]+ b ) G- F9 d% `, n& ?* z

* B4 v5 f, [- I: K. E

double big = 0.0; M3 _) q$ G' F+ d ( c2 Y, L0 ?" Y) e+ {! M

$ b! i9 w; f9 k! a/ S( R2 [! P

% g& Y" F1 K Q& \, f

: j' d: G6 V& C

& J, d+ w* H; {) B, I, @! o

for (int j=0; j<n; ++j) ' K: l: j9 t7 y4 q$ @ 7 m& m* e$ `0 M0 ` 0 H' {4 }/ Z0 t8 ~4 W

7 A+ p! ^, Y7 o5 Q

if (piv.at(j) != 1) 6 y @2 W2 g9 F 5 C+ u7 H: |) U0 W& A6 u/ L0 |( E

; v" K j, P1 z6 m- B6 ]6 U

for (int k=0; k<n; ++k) { # w: p+ v9 z% `+ R" B' B$ X 4 A( k; v3 B0 B4 A# k6 g( L 9 `# y' }7 `# _

if (piv.at(k) == 0) { ; d% o" D, [6 h$ A( h* b

$ A' l/ F$ G+ n& R: g; T3 w. T, y

8 E& i$ N$ N. Y

if (abs(A(j, k)) >= big) { & W5 c+ e, G1 ` n3 Q, w

( a3 k9 S! q# K7 A

big = abs(A(j, k)); * g& D' d$ r O' s5 T. I! q$ Y

irow = j; 1 `* h0 }) l! t' C2 j 6 P) H, z+ r, r 8 V& N$ d8 d5 `

$ o! N# {6 z/ H; G

icol = k; $ H8 N1 X: ~5 L6 u9 O1 R0 A ; r }* C' d7 d1 O$ `

( }) D' F& V7 F

" K: N5 I B- K3 ^2 T" w% q" l

if (irow == icol) break; 6 a. g: v! e7 w5 G/ ~ # o0 _' l, u2 e8 [/ R

9 d3 |: y0 E3 c @. g: B- [

} 2 Y) z8 F* d! U% G 1 S2 Q3 q" i% e

% i% Y E( Z% z$ p

} . k; [, S- B/ N9 ?$ U2 g* X

1 v5 ?) D5 X2 ]: `

5 \. C5 H6 e: n1 O3 L; j4 S9 a# f

}

8 w# i, J) X p! Q+ v# e

& G8 }% g$ ^% b' Y. M+ o; b5 E

' ]* k# Z& _9 p1 [7 M

++piv.at(icol);

; ?. K3 E7 h+ @& ^

; x! U& u. n) q: m; C

' G( }/ b6 ~$ d: `6 p+ s

& l1 l* v7 P8 [

//进行行交换，把主元放在对角线位置上，列进行假交换， " F+ t) T; e0 A. c W6 O/ }/ U

5 s: ?: Z( K3 x; B/ T+ f9 w5 t, {9 C

) i P7 y1 a# _; ^

//使用向量indexrow和indexcol记录主元位置， 6 m0 X) N' d3 d 7 M) K( L% N4 `: ` . Z) S% \# e: J

D' V8 W @7 s) O* H" H0 u" n1 W6 [

//这样就可以得到最终次序是正确的解向量。 ! I1 K) P9 y; `( L+ Z6 e ' E$ O2 E" A1 M a9 k6 P5 c, T$ m- D

! J8 v- T% v! C4 _) t

if (irow != icol) { 1 @3 [ o9 J+ G* R! [7 m, | 3 n9 u( S# L1 |5 a" N0 a

5 k, M2 R5 g2 t

for (int l=0; l<n; ++l) 1 h1 L. s P$ R0 v, d

swap(A(irow, l), A(icol, l)); ! R3 ?: W; O+ T! w t

3 p8 M. @+ @* n1 g6 ~1 V# A

8 B) X9 p8 L, X) V3 |, ?7 ~

! c; u- u% a5 j- N8 \+ ?6 n- Y% s

2 }( L2 P6 X% B* r2 I! p

for (int l=0; l<m; ++l) 3 c% b/ \. ?/ s* d) o. Q1 F + Y/ }; b3 N) u

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

1 k5 r( G# F) Z3 Z

} # }5 @% z1 d- g - N- k4 ~3 j; H3 B, v" G+ S

8 }5 f2 E5 K: |9 m$ I/ N

indexrow.at(i) = irow; % A o& x, i7 p ! b A+ @$ H- k3 ~9 N

6 Q* O5 x/ B+ F

indexcol.at(i) = icol; G* C6 S: t; W" t+ G2 a 4 ~" _! a( L0 W! d4 Q ' j( Y# b: y/ T0 L; R" W

( R) |" D5 y! f* h

4 A/ m# N. ?1 H. v$ C

; l3 d; D0 i9 U 5 Y% C I5 t$ F

6 h" c; v+ \+ s

9 f; S9 ^! U& x+ [! O8 h

try { 0 t7 o+ }! {5 W) q+ |. J5 O

; L6 T, |) h# P& p# g* j- n# O' w

+ Z2 F9 {& @1 ~

double pivinv = 1.0 / A(icol, icol); / ]+ {$ i+ Z9 X1 `: a" ]& O8 ~5 x q5 V; l, J# `9 o9 M+ F3 r

5 g) o# O% ?3 o4 b' R7 h

for (int l=0; l<n; ++l) ]1 z$ x. ~0 b" g: C, l$ s9 W

' l6 b: g I& P5 w& j

! D/ q. B. b' l1 @4 x

A(icol, l) *= pivinv;

( H- Q# Y0 m: E7 _

" C5 H l- q7 J; m9 o

for (int l=0; l<m; ++l) $ _% D" E4 W/ v W$ r1 G% _ 8 w( u! [! A M4 T0 V

b(icol, l) *= pivinv;

8 Y5 A- c6 F3 _4 |" e. }

2 A& {' E- n g1 O7 m2 L

* c1 |' Y$ r0 ?3 @

//进行行约化 ( g q y& l" V; {8 J8 ] ! Q3 J" T$ i8 D: I5 w) x

1 C8 J7 e; ^" r1 z; ]) n$ w5 l

for (int ll=0; ll<n; ++ll) * n/ q% M- `0 b! _: x ( X5 l( ^6 X% a3 v" R( W, U* v

4 |) {, J0 P( J" S, l/ q9 w

O( @ x) G; X& {5 C9 k

if (ll != icol) { , N' d4 l5 B6 }- e6 \4 x3 [ C 8 S4 `& ~' l5 s0 P4 ~ 5 B1 d1 i4 K$ w) c: }6 L9 y

n# a% {* _0 ?/ v* v5 c# q

% z" y, g* X$ _' s& Q

double dum = A(ll, icol); 5 m8 D0 J s9 s3 }& b5 a

4 h r+ l7 [& b

/ X0 r- S+ q" P$ Z

; c! _% Q) `% F- l

for (int l=0; l<n; ++l) . h8 e# J) g$ G: _4 s6 a; J " C) s* D( R t) V, T3 P1 ]

/ o/ V2 Q. O& P2 A" X

8 U' E( H! ?7 O

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

: S1 I3 Y5 Z# |6 Q6 ]* q

& S8 ^9 }" Q4 K! ?# t8 L

for (int l=0; l<m; ++l) , t0 S' U( l8 k P8 H3 g8 X; [& U

b(ll, l) -= b(icol, l)*dum; + x& @# }; W5 r; J. |

! q% n5 P! |+ V0 f7 s: E

; Z+ Z# A, Y6 T+ y* s8 k2 V

} # K& }' V+ V5 `0 \: b- y5 j" A * k" G& M6 |- k9 o X! s& D, e$ M" X8 O

* {4 H; K k8 s* `- c8 s

8 O6 E, H x F$ X3 G; B

} 0 n( S6 ?( g1 @" W, L . n5 t0 T) u. E! d, D" H. y0 J/ z3 Y 8 U. ~) ?. l+ s4 A, i) p) E

; M/ O* W2 r4 Z: R4 b

catch (...) { 5 U( e" A# x7 m

( ]; j! z& \" s7 P* D- T

cerr << "Singular Matrix"; " j0 K* I. l/ ?( I( Y. l8 s * V1 C2 X/ l8 W! b6 A

4 ]" p- o% W% [. k' v

3 I9 C( f% B- ?$ @; |' _( P

} 3 I$ F, W. Q0 u" g

} # Y1 ^7 s- u' j. l9 x! V 5 ?3 P* X+ D% ~6 ^1 T& F0 }

} ' y) v/ {6 k8 u- ]. ~% H

7 S! @& h/ @! [" s( R2 }# m9 D8 S7 R$ o

% F2 D+ L# Y2 ]. d8 z

2 r" M+ C B- w7 n n& S

int main() : B2 E1 m8 w% U / Z) g/ a- Z. V3 c! _. @" C

7 Z# ?. U2 B' p, z

{ 5 ^$ M! ~- U: C G, ^. _ " a4 u6 [: _/ _; _) E8 {$ @+ b + ?; e; P7 a- n& n- u

//测试矩阵 3 S# {0 F2 p* Z

# G# `/ N0 o6 O

9 | N' T) Y7 D

Array<double, 2> A(3,3), b(3,1); 6 x( K+ f5 \( S6 V' b" ]1 ]

! l6 k* A+ I4 ?+ j. F& s ?* T# b- k

9 @2 k& [9 o8 n/ v+ q+ K

A = 10,-19,-2, 5 m% g: s: G, u0 l$ N6 X7 r - Y e! g& Q7 { `! e , A7 n4 D- G7 T0 R8 f5 a+ g$ O. p

, F; X6 i8 U( M+ ]" \' _% B

-20, 40, 1, 8 M2 d( e: {% v( C 7 }7 V# n& m. o) H2 v- e4 J9 z + `1 @& f3 v6 `/ P# u$ k

& [3 {) F. ]8 s& F

1, 4, 5; ; V3 T7 j! w8 a' Y1 \- w$ k % T* p, d l; C# z( X/ {; F& v

; m# r5 U! o9 o9 d% u! M1 v0 r

2 ~ G2 j p# ~* z w) d

b = 3, * b% e: N# N! f! `% a0 |: W+ o8 h; [ , X/ {& A$ V, H+ T k

. f8 @2 C1 L1 {$ w/ {

4, ) M8 ^! ]1 E8 x4 P

+ ?/ R8 @/ p8 d9 E

; A; t3 B% l7 {2 t8 R: ~4 Y) I& o

5; - V& [6 Q* h4 F; c; d 8 ?4 {5 a. s$ A/ h! y( r

& e+ [: e, ?7 G# D

% \7 I* |/ T2 M$ x$ C' T

; q0 i( D, G: i* G/ E

( L% m# a$ r) v, ?

- P/ B: D+ y' N* O: d, O5 `

Gauss_Jordan(A, b); 6 T2 }5 ^0 g1 b8 b1 M# M

P; e( R! t0 o- M1 J$ n: |' X/ v

( S" i7 M9 K+ b% l

, R) [; ^3 A: {2 H

cout << "Solution = " << b <<endl; + M1 w( ]6 f/ J* R) e- A 6 b7 t. H" x: U2 b, n% D% P 7 F$ v, r' u' v Q- W& x

' t- ]: m/ f% x& ?! ^& c. z

} T1 [* _2 m' ~6 [! X6 l 8 b. r& q8 N1 o" @6 Y0 \

" j" E+ T0 D x" {4 ^! r5 T

' I7 c" N4 @) s9 I3 w# q; y7 w' v

Result： - I; i3 d' ^9 p0 m4 Y 1 B( B% ^+ W2 ]) W1 t

/ @- Q7 Q- Y8 ~/ B4 U @; B. r" ], D

* Q$ F- b5 l d s2 I' Y

Solution = 3 x 1 - j! p" A" ?* ?7 t

[ 4.41637

, w$ ^1 U% ?$ T! C6 h$ d x

, L5 f4 U% _8 l8 u- g, s5 n

2.35231 - E9 P! ?9 d' H' X , Y: `+ N9 ]- S

8 @0 U' [9 [3 x, _. n

-1.76512 ] . J* K% H! Z: Z: U$ q

: E7 K4 q& g$ D2 x( |0 j1 c s3 ]

0 y+ q, H# H1 o2 j

% g4 u0 I9 G0 H& D3 G2 o* r

5 S, l( s. p" Q2 P- @/ b- D$ U

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 , | h8 c( P1 a; r7 A: g# `" s

+ M5 }* t( L; Q

4 Y2 c# F1 ]9 D

7 {! K" Q8 m: O# @1 y. }! W. t

) p A' A5 \8 H% ~: t

4 S) q8 O7 F$ g* O6 |- x' A

注释：[1]主元，又叫主元素，指用作除数的元素 7 P4 z& d4 h+ m( D3 j3 |) S' i5 O2 \

5 j+ O, W$ o2 F: s5 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

路过。。。