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

lckboy

全主元Gauss-Jordan消元法 : }& \5 R3 e5 n

" z. \/ {* d' Y8 {

- Q+ H2 i$ k9 b

7 w/ e2 L3 ~4 \; S2 Z8 h/ _4 |6 p

( N9 e( x; R. i0 k" l) U! m6 d" ~

$ v0 [3 j/ I8 T! k5 [/ k# [- S

7 b6 C( Q% r; b& I% _9 i

* w: E+ n& t3 d- s: F) e$ S

5 U R+ v" R( K- H. y x& {

3 l/ f A7 p$ r: x+ v, O7 v' h

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

) G0 X4 |1 j! v4 \

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

1 r3 I1 z2 ]" n* ~

C5 x7 x. D2 o9 O! U# W. E4 [

, O7 `3 v" ]+ |! m' J# s# k, s

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

9 b C- k3 J' `; k& j

Code： 2 C1 ^* p. L( A# ?6 V7 \0 j. z / s2 }+ e& d$ h# E) X

- N. ~) r0 ]) G- y

#include <blitz/array.h> 6 m. K: |. c6 T5 \3 S4 `

) K! S6 v t$ E# {% w

#include <cstdlib>

#include <algorithm> ) A Z6 _$ y- H; X; u

#include <vector> 0 H, c; o5 x) T G9 Y : g: k% W; G4 B# C) q+ S

2 O( [4 e6 z% a# e

using namespace blitz; % a1 n+ E8 J, U! j : @' |" h. S) `! }. a ( q+ d) Q$ D l/ n2 w, j( n7 c

) O% t$ C4 i. @+ F2 ?

) v0 h3 j. n" J6 R& W

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) 4 D/ H' L. [7 z* Q4 g' g* v

7 M4 j2 g. G, M7 \( @

{ ! [* `; Z4 I9 d" G- `! ]: q+ y 6 s T, X0 G, f

" w3 r" I/ O6 S+ A/ m0 T

* V. p ?$ Z# O+ T4 u) V

int n = A.rows(), m = b.cols(); : R9 W/ g' j0 }$ m- w ( T \( f! [9 r$ `, Y

int irow, icol; ( P! p# c4 ]0 w7 ~: G1 Q8 f6 A! x, ] & a3 c0 w1 \7 c# \* u. O8 V

vector<int> indexcol(n), indexrow(n), piv(n);

& V$ _4 T( g' Y0 \

- |) k+ L* t5 j. ]5 Z1 ?3 |4 c

: f6 _: ^) a5 a4 {9 l. U

for (int j=0; j<n; ++j) 0 T3 }3 e4 T6 d: w j# W$ K' _6 \2 O! ~, U

" P& R0 Q# Y0 b/ |1 z, b

piv.at(j) = 0; ' z. A$ |2 f: d9 z+ N" U+ p $ N1 l( P% X* V( O1 R/ k

2 R; T& R% }* g' a

$ n1 c$ _" ]- M* C& A # [% L) T) I+ h0 D& F3 `9 m

//寻找绝对值最大的元素作为主元 * }) l! Q. p" M' N u. ^

for (int i=0; i<n; ++i) { & s# n! J5 P7 q1 w! r' |/ h

' h5 f! b3 X* P! C( G: ~$ ^

double big = 0.0; ( ]% r1 u" \2 e* { 8 J& I* Q! `2 ?8 z1 A& |

- j+ X! d0 {& X- n0 N7 i2 _9 n

, N* G7 F" ]& n* C4 f( M8 n

7 t1 W# u1 b+ o, C

: A. v( C/ ?" W+ i" ^9 Q

for (int j=0; j<n; ++j) & H8 m! I2 Z) q% C6 q0 D 5 Z9 W- x8 E+ E

if (piv.at(j) != 1) 2 K+ d4 C9 X j

: d2 X' a/ ?0 M p: F" Z* R i

' C- p: T' F, |( C5 }. P" Y/ v

for (int k=0; k<n; ++k) { 9 V6 M5 n& `( t

; H5 [2 S/ b$ ~$ D3 h& k

if (piv.at(k) == 0) { ; K" _7 X3 G: A6 t, o! }: Y9 |' P ! }* |' p0 M' A; p3 h

0 o N [6 G' y/ ^# d

if (abs(A(j, k)) >= big) { , K6 E5 C4 r1 A$ {+ ? 9 o- Y5 i1 g1 j0 Q6 w

$ o4 u4 [" I1 N

big = abs(A(j, k)); 9 a$ P! k$ `# Y4 u$ y- B2 x

; Z2 I4 F, O' l& X

$ h8 |5 f( Q* s

irow = j; . ^8 p* N, Q& Y5 \ 6 C3 l' | y! H9 e. z4 T# C" \8 c2 L

) P4 J: x8 L0 e" j3 o

icol = k; ; Q2 z0 M, Z; g. M9 z : G$ c6 J/ N: [3 L( K$ h# P$ @' ^8 v

' s8 d$ r. T/ s" B9 ?* p# w% [

; |: C% O3 J4 W* v5 ]1 {! S$ M

if (irow == icol) break; 5 r5 d* T) m& s4 J S

6 R% c1 r/ L2 m# @* g

} 8 C' \7 M$ U3 r3 H( @

5 \; m7 S4 B" \& q' O, F: W

} C5 e) s2 U: Y% |; C: J' y 6 }3 z( h, M* o5 n) \- ]9 G

' p4 |: }! m; ~& y% s' b

# u7 R; c* `( A! W. _

} ' q. u$ z2 `7 N3 j

1 P3 \; Y8 y! ]. E, y

$ `4 a$ c% p2 e

" |( x0 p* j R' k- H- Z

++piv.at(icol); ( q( c8 z& |; A) f4 k" M 1 a/ n; V- y6 v4 o

, Q" Q- R) ?4 l( @ W; J8 n- g/ w

" \4 `) w. K: l

& O) l) ~2 Y6 t

//进行行交换，把主元放在对角线位置上，列进行假交换， 4 m0 Y( A$ K, d: a5 T' K 7 I, v( F1 C0 ^5 o$ n* {# n

//使用向量indexrow和indexcol记录主元位置， - g+ s* N8 R7 _; y- o3 z4 | ; N4 j% C6 _5 V$ v1 A% e

//这样就可以得到最终次序是正确的解向量。 ; V; U- C4 z1 X& Z" u- k0 ] 0 P3 Q+ H; o9 q- B

' |+ \) l5 |5 l1 m8 Z

if (irow != icol) { ; h4 \' L: ~. E, F3 P- O+ g. _

* ~, E* L2 \0 l/ o m

for (int l=0; l<n; ++l) " m$ E& R1 E* x$ O 2 B2 t+ t" e/ B/ I- _! T" C; }' K

: @0 @& b# e; |* |. }/ i

swap(A(irow, l), A(icol, l)); 6 w+ o/ a) X' }1 W5 k: ?

" l% L' d' Y* K a+ C- P$ n* n

5 ^( z/ y) G% }4 b

8 e: N, m( }8 ~& q2 h/ B

6 k; N# N' v" M! ]: h

for (int l=0; l<m; ++l) - w' M% |1 k; D# V* `" x

4 E% [' v; T9 x: c. t# i: |1 y. J, m

swap(b(irow, l), b(icol, l)); 9 Y! O `9 A7 B! s, _/ D/ [. d

* I, E7 d. T% B6 ]: Y& k

} 0 j9 |9 |0 i, d5 M# \' W9 @5 ~ + [0 f; q- ^( c+ R/ M) d

indexrow.at(i) = irow; " ]' R0 T. U$ H0 g; l1 J

indexcol.at(i) = icol; : r1 a( W0 n* f l( K5 k, r$ |5 G

% x! x! U h. E$ c v

$ P& _; h2 n3 v( S9 _, u K

try {

- c Z q2 f) ~8 X5 D! a4 W

double pivinv = 1.0 / A(icol, icol); + v) t1 g# S2 G* p _

, V2 x, [$ f' A4 x6 `6 S' J

# X7 F! U: d" ~1 a) w8 |

for (int l=0; l<n; ++l) & N# \9 \; x5 x/ U

A(icol, l) *= pivinv; " f. J: _! z0 { 1 m9 X( G/ p) v2 P; ?

for (int l=0; l<m; ++l) P, q; [& ?; E8 [8 K) q

b(icol, l) *= pivinv;

0 v( C( F; K8 S( K% j5 l

9 u4 R0 N; W/ r: v! ~

- ~9 w7 w$ M; r8 A. U: U2 R

//进行行约化 / D* f/ c$ C# {% o. x7 f

' J4 ^6 V m( i# C* x2 {* H# R- W

for (int ll=0; ll<n; ++ll) , ~4 W' p% |7 w

if (ll != icol) { * E% s) B. Y3 s" X* i' v6 V$ r7 o / m* k. z2 c0 h/ w6 k: Z* H( M

. N& K% L; V" i" U p. N

double dum = A(ll, icol); + i8 j& Z: ~9 M0 x b: ~7 I! N2 O9 W1 t u

- a) q# o* l7 g' L/ v

# R" B% S3 H4 c0 q6 C' c

for (int l=0; l<n; ++l) & ?4 n) G; o2 i% M7 F7 J ( k1 v! n* X U) [6 F . @, G4 x- v$ z g J' \

- q5 ^; K3 O# h' n, C# ^

8 v! [* l0 r5 L @8 s

A(ll, l) -= A(icol, l)*dum; 2 B3 I: m# W# T" t) C

g1 n$ [' ]' A& o q. D& Y

for (int l=0; l<m; ++l) % M7 f$ @* v% D. N$ H

# y: H& e. L3 T; Z- m

b(ll, l) -= b(icol, l)*dum; + }( I: V5 O, z

8 w. U& t: l3 {8 Q6 r; L0 U

} # Q8 ` Q7 j2 p5 E5 H# p% ? 0 i) t; h, B$ D% S( e4 S$ a; d

1 i" G; S. d. H6 T0 ]4 S

}

3 J6 u5 [2 X2 _# S* h# ?

catch (...) { * {7 j: G" J' W: y 3 z9 X4 w; }! T9 l# z* Q; P% `

7 ?# ?' D9 |% X" _/ Y! x1 t

2 [2 B4 F- s# P. ^6 y

cerr << "Singular Matrix"; 0 e0 }3 ^ y9 @0 c! B' F5 h

4 F% A2 V$ [! e) G- y

} - ^+ q2 \: S8 A' W; O5 C

' k; k$ i0 [+ g! u

} % M( K2 v! ~, S; d# z, ^

4 l7 v% w# ~) R0 c& \9 \9 S

} N" F, @5 s; @/ q( ~, S

j1 V; h K6 u+ H \. z

* @0 h, A: M& `# \

! d7 H) b" C1 s

int main() " _" a6 |1 n; Q7 n4 i

1 u2 U1 Z% ~! }

{ ! A# S+ s9 j ^7 Y! G+ O, i% r # X9 b% a+ \: E& Y" n

& E; @4 Q2 ~" ^$ g

//测试矩阵 8 P/ G: ] Z0 d! M6 t7 [ 2 j% L' K" E1 P( o

Array<double, 2> A(3,3), b(3,1); 3 l: Z) I5 Y- E4 J8 Y* p. j

9 x5 ~' ?- \* i6 `

A = 10,-19,-2, ) J( L( ~/ e5 C / A- ]% ^" L, n- p9 b

-20, 40, 1, # U1 f$ N* C0 k. D; J0 a % J l3 y& N: x" R( h, ^3 R

1, 4, 5; / X" K: @0 ]* l; t: `! b; n7 P8 j+ Y ! v0 Y# z! ^/ A% A

9 @( O) G! e! K" b0 f# r8 }7 [

, p2 B, l, e) D* \8 \9 E: [, ~2 k n

b = 3, 8 I* |9 r; `/ d7 \, D 1 N7 E/ B( d4 J' j9 ?' d; O1 F! Z 4 Y3 t8 `1 E( B5 y9 Y& {

# _& Z* v/ `- S& r$ D# F8 G8 `

- u; ]8 `) S& @9 ?

4, . ~7 I' L3 C# O% d2 { n# l6 t8 U- T5 W6 s# b3 B* c P

5; 8 ~. w* A+ ~8 i* U

# x: A6 @: A% p8 m: r

$ s! I' U" t7 z. [7 x1 D! i3 g7 E 2 |9 T' J, s4 ^2 T$ j6 M$ t

Gauss_Jordan(A, b); 0 q; z9 `, B2 A u* M 9 b1 L, D+ ^0 C( }

! F! @2 y) l, N6 z' ^

- D" V6 P7 z* J1 |+ J! g( {! v # c- v, `, w# h0 r0 D" `9 k$ l5 e - x1 h' h6 d7 H$ g6 D0 w+ B0 d5 Z

cout << "Solution = " << b <<endl; ; [1 }1 U- m3 G `

% j7 b$ o v3 L

" J& T! p" C( F8 C

} }- F' o( @2 x+ {" }' Z

: P( H, L& ^8 l/ |! c

2 T2 l7 o, J/ F# g) z. P5 {

% B3 a7 F3 X8 e- E! K

: q7 p$ b% d, B" l7 n2 A7 J

Result： ! i6 Q0 S* V" ^ 6 K7 k0 Z+ t% C& `* d0 B; k

% L; g( X/ Z& \) ^2 ~

. _: p- K1 B5 Z v+ E8 a

Solution = 3 x 1 8 `+ u# F/ w% o, f1 A

: Z, D$ ` p) ^# ]

[ 4.41637 : f8 A! c6 P1 R) X( r+ u2 b \8 t+ Z3 x* K2 ^/ Q3 n; N2 i) Q $ S4 d7 J: w1 |& I( z) r' q

2.35231 4 z/ d- }/ m2 ^* Y( R8 i % ]5 ^+ e' P- q! T2 d2 L

3 i( B0 {. V( G2 r* \. \( Y% a0 z N+ n

4 {2 y5 R7 _1 O8 Z& \; E

-1.76512 ] ! u0 k" ]9 v7 z3 a6 K 8 r% {, Y# m9 a( Q# g

: s7 U# f. l! ], }. X

0 o; j, @3 b: B# w3 k

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 " m: N9 n/ _5 K( l% q5 k

/ k6 g9 t( c3 H+ J

9 i- i5 z/ z) s3 ]5 s T

6 `4 }/ t. m, z! \# R4 N* [: n. p" ~; s

注释：[1]主元，又叫主元素，指用作除数的元素 : v; W; z& x: P0 q* L

' O+ y1 N+ x) A$ v0 i( `

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

路过。。。