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

lckboy

全主元Gauss-Jordan消元法 6 y2 o# p/ R+ M1 w" k2 G

l: @- x( u) Z4 s) }0 j$ d

% |4 W- N$ [+ r" ^% K1 s6 |$ R9 r

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

7 m1 J) i4 X5 D$ ~

h/ r4 l& R j' X( M

; j( z, m: K$ x! X

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

9 ~) L% K) f9 N: Z _# ^, }# e

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。 ! T3 t7 v# ]( W

& b$ W, P p) Y5 D. B# [

8 X7 b0 Q1 k& ^# K. `* X: z0 D

Code： " x3 A2 Y# T9 I" v4 _ % \: k- x( n+ V' P

# J" I5 a: P% m

2 @2 U: A: K' p, l/ |5 I# k" X/ f

, J2 |% ^8 A, g7 U, l

" n8 H) `5 Z {" A/ u) ]

#include <blitz/array.h> - ?: D; Z9 k! m7 [

2 }( A7 f2 w( {/ n$ a- R

#include <cstdlib> ! G# L) k8 ]" |, y2 v: K

0 H# y# G2 V5 @$ n

#include <algorithm>

. y0 ^$ B" d3 j& [8 O7 p; W! T

#include <vector> 4 O9 i: J4 P6 u8 H; T: }

. a6 A' T1 y, x+ ^& v8 G

using namespace blitz; 7 g& g& b: s% `# _0 T

Y1 l% S0 Y) i- o: U. _ m

8 ^1 V2 G/ E# ~

( b, @. M7 |4 R2 X" u

8 }5 s- \" p. o

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) 8 d3 t% x4 b2 i 6 P! O4 Q$ A! T7 z

0 x7 {$ q( A' A2 o* I' M6 H

9 U- C* _3 G: T# I# o4 G! c5 b

{ 6 g5 S/ ^6 m( n! \# W% Z9 B, O8 I 8 H0 s& ?/ d, h; F5 K

int n = A.rows(), m = b.cols(); + P9 w( Q# L3 j. `

C, ]0 B1 [- [# b& ^. V

7 x. [2 j, C2 f; C/ }

int irow, icol;

* ~& j) k$ A: M9 B

9 B9 v @1 z; Q/ M2 v4 G

vector<int> indexcol(n), indexrow(n), piv(n); " [2 U$ Y- e! |. Z$ n# f) k 8 ~6 r; `( P! I# T+ C

! v* ^+ ^3 C( ~7 Q, R

for (int j=0; j<n; ++j) + Q- s* F$ X ^- O: V/ X

piv.at(j) = 0; . l4 Z' p( Z1 P" }. } ' f8 K6 M% N. s* h! ~; S0 z& u! M

( c* }; x( _" d' a! k8 D

7 A( T, `# K0 Z( |) A3 s

//寻找绝对值最大的元素作为主元 * l8 R7 h {, P, a }% R / j$ v- p9 ?/ X

" A, T* N3 F" a2 `" G7 j

for (int i=0; i<n; ++i) { + \1 b( z& X) [+ C9 i3 W' [6 G / R0 v' d/ j! F! L: u, O# @

+ }! y; L6 ]* _: I. M

- V" }. Q, U4 Y. T

double big = 0.0; 6 I5 E6 B) S' l0 P" [ Q( w( r) r9 R* }/ A

4 V0 O2 \4 J Z+ f

9 G. e! a+ V8 w; t+ Y8 l U

0 E' C, U, z8 n6 \" o

for (int j=0; j<n; ++j) 9 i9 d& q/ ^2 ^+ T8 |% u: u ! Z8 j L& q: a: `

; [0 M0 e( n' y- J" Z( ?( J

5 B! @4 P0 U3 ]* _; M4 N

if (piv.at(j) != 1) 1 t) z2 f. D: n, X5 O! m5 M% _ ! d E. k8 R. G* P

for (int k=0; k<n; ++k) { # L: A1 r6 {' M+ G( v1 ~ ' L* k) M5 y6 w+ r

! N+ H* c( Y- O' h4 ~4 a1 ~- k

+ w P: H: y1 Q' f. F- [# E' H

if (piv.at(k) == 0) {

% n4 p. M% s/ o# x* c) A S

6 W5 Z" g! v, }0 s/ `- B% I

if (abs(A(j, k)) >= big) { 1 B- m/ {' g; M7 P ~. I7 {, t' C' p, \8 t; T6 V

E/ X# r% U J/ z) v3 Y

big = abs(A(j, k)); 4 o% ~9 n# y# P; A 8 u* S1 R+ X$ i7 V5 l8 z3 F

irow = j; 1 q! y$ z3 ^; v5 C* x ' i7 N. G5 V, D1 A4 \

8 j. P8 o6 m" U3 _7 S K! c) U

icol = k; : j( L; e' `& K4 x9 C

if (irow == icol) break; , u4 q4 o7 t6 v) F, n' A1 G 4 u: u0 w$ w% O, M3 Y% {/ i

2 p. x9 n/ u# B0 Z" F- j7 L

}

} 2 `% Z! M9 V# Y/ d) R3 m* {4 B . g% l% }1 D% V* f) m- I, M) w# K- b! E

- L5 T9 V& I5 f5 Q+ W

} ) O z K+ A% c- u / t5 h8 ?8 {9 ^8 y 3 @# m6 [% s/ U* t$ d4 p. G7 H

( F9 T0 x* }( B& n- K2 G1 \. G4 l( G

0 G, u% _: z! _8 U; j4 p- b; D

" T' [5 r- n+ ~4 l8 N# w6 l

: J& f3 p# A1 l9 q2 s

++piv.at(icol); / G, h8 k) g9 R0 F3 | ! A- g" r/ O1 ]9 h$ C* A2 i8 M

% U" C7 Z. v6 M5 \0 P. y

* o2 F& J# r& {. j" ^9 L

+ F( Q3 E6 m9 Z' [8 d8 o! J8 g & b6 g, U n& C4 [ x% ^& N" M5 | 5 _. f0 e2 \+ _' {6 z1 b

9 m2 l7 Z& n7 C% O% M* d/ o0 g$ R/ B

//进行行交换，把主元放在对角线位置上，列进行假交换， 7 G, [% i0 O6 F9 q0 n

7 S/ T* x7 W' `; w$ R, ?8 r% I

//使用向量indexrow和indexcol记录主元位置， $ e1 ~/ b* }0 w6 h% \: D4 [ 2 i, F4 g+ p2 `2 r9 z) I

2 E! f \, H F1 T) b$ }

" ?5 W) L4 d6 j, T; N' i. w

//这样就可以得到最终次序是正确的解向量。 3 r$ P* t# Q. a! U. l) @( f 6 t9 j/ N/ a6 m! u4 e2 Y/ ^

! q( ~5 `! ~2 M# O$ v- E

if (irow != icol) { ; P& M! W5 s! r 3 d) F$ l5 v( b2 n

$ u- D7 _9 t% y; r8 D

for (int l=0; l<n; ++l) : j0 J0 G- ^3 u6 K3 C3 q 2 A0 }& P4 E$ U/ m& Z

swap(A(irow, l), A(icol, l)); . F0 c8 C7 u# h. e' e% D, c! S * _# X6 s4 {% ]9 B9 _ O

2 N! j7 C5 a& n; m

3 X/ a* O; s& k; u' p+ Y

( S, L2 |/ L# k- `- M- m( Z, x! _8 i5 O

2 Q- ]3 W3 I1 I& J' r: u* U! F

for (int l=0; l<m; ++l) , U0 m/ S+ u( V2 x* W2 U- a 3 m# Z5 k* o. H" ]- T# l% C' s

7 z3 Y G( B; k- b

; Q1 a+ F$ ~7 ]/ g

swap(b(irow, l), b(icol, l)); 0 }- w1 v2 G: l ) }1 ~/ l$ L" o2 l6 u

; j$ S* L: e$ [' t1 f1 H) k7 ?

} 6 r8 {! N4 F$ [" _& }

- j* k% Q! F' H [* a) a2 y @9 ]

/ s! H+ ~2 J* g( `3 g( M: k

# }- E1 Z+ _4 v, k- k

indexrow.at(i) = irow; 2 @2 t! j0 x3 l+ B" M! @5 X

indexcol.at(i) = icol; * \0 C: C" R, E8 H

0 z( S- h% V3 {0 G

# i3 q$ ~( S0 ?/ h; d& x9 { 5 l" {4 o5 g, V( A! E" m/ {+ D

5 @: I1 x! H7 v- D& Q* k7 o

, ?- s5 V) C' y

try {

4 B) o- h5 j r* M

double pivinv = 1.0 / A(icol, icol);

, ^& F" s/ M, t; s. \8 _

5 W$ ^7 c1 p( t

8 M: K U# a8 y2 V& e( f& s

for (int l=0; l<n; ++l) : `; y% [* u2 s) e! V+ D1 p " \& [6 X: |& k+ ]

6 t' j- m' F/ l4 @5 M

A(icol, l) *= pivinv; 4 E' w8 Y. o$ e6 ?$ [

5 t" T; p3 y$ p4 c' v$ U* L

6 O+ b% U0 B! J0 q3 O( |# R

for (int l=0; l<m; ++l) ) B/ G0 W1 W# A5 n+ H9 z, D

b(icol, l) *= pivinv; ) E0 I0 p0 T! Z( ], W% G. m & H- @, M; y/ \

/ A3 R5 N$ e5 f

: [# o$ l b7 r# p6 N+ J

+ W# ~8 D- k: Z- j$ O+ R: a

5 R- ~- y9 T/ P! h& s/ Y

. E, b' b' p/ I# @

//进行行约化 1 B- @! |/ u9 m3 G9 N+ I; { $ u3 v! ?6 a! P7 Z8 H

* S; e5 ?# n" _; [, O* m8 j; @

/ M9 @. t+ D# y6 C7 N

for (int ll=0; ll<n; ++ll) 3 ]$ _) {& P7 i- g/ F+ _ T H $ a, g' y% `3 D, n3 s$ V* [ u

. t: V# W8 V1 t" |) ^. A

if (ll != icol) { - @, E! G% D7 l1 c2 h5 Y E

* ^! m% m) y8 J

. z( \, L' o. K7 u$ k: }; |# e

double dum = A(ll, icol); " i- s1 U: a& F4 t 8 G- o) H' P1 _% G. I8 T0 @7 l

' L+ W( F; G4 t5 G- L* }

; ^2 k; f* U4 S" K$ b- N

}7 d9 w s8 i3 Y- H! U

- |; q5 b6 f9 a: g

for (int l=0; l<n; ++l) * m- g5 f- E; s 1 s3 H* o* d5 \, s2 j: A) Z # s9 [' H. @! A! V c

4 ]. c4 O0 ?2 F* [

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

4 E) ^$ L: `) }5 d

' K6 t, U7 {( S, O

for (int l=0; l<m; ++l) 5 m$ F$ t. q9 [# H/ U- _: N T

: @+ b& T% U/ E' m' s' K9 g5 F

m4 i2 A- P( u$ |$ ~4 A- ?

b(ll, l) -= b(icol, l)*dum; , P% C6 [: d) S6 U! p# T4 g & u \! Q5 Q0 v% O

+ l$ s' S( j3 l. a. M

+ c, Z6 l; p" j' M7 X7 Q5 h

} 1 z9 ?4 L3 y! [1 F

} * G$ N. L; U) H; Z' D7 {! N5 a3 B/ }; Q ) u/ T' o6 ?+ J8 v; n

% G7 G& R! w( x6 B

catch (...) { / M1 w, z2 h, D3 H* Y& U9 x! m$ K & e3 I) T' L3 I- W

cerr << "Singular Matrix"; 3 J! P/ _/ C6 W+ o" a 1 }" }: f( l0 v, h3 b+ o5 u* J ' C2 X" a/ g5 M" v. i* x( @3 ^

( Z o0 k" d9 A# q# {

} " b. [5 j* G, _$ A A $ F/ s' Q; {# b4 n6 g$ {4 t

- u+ u" b; ^1 C' N/ F6 H; [

} " M& z& L: x5 q. q2 @ 8 E/ U& X: @9 A! F* @) Y# w ( m% h. l. C' X

& k7 w. N( G* r% v/ Z9 i2 \# |. c

} 0 _) c2 z1 Q9 g$ I 1 z! Q* N% @, G z% U% P" T

' I5 h8 k6 a8 z5 X& \4 @) U

' d3 [; E' B4 P

5 @. Y1 p& y* H

int main() / T; v4 [ B Y r

# X& T0 c6 J+ Q8 T4 q

{ 7 F; {! j( E+ N7 O

8 G7 R6 L M4 b5 g: Z! B

//测试矩阵 : f( o" N0 s! m" q 1 z$ s; N( o4 p. F% H) J! D

Array<double, 2> A(3,3), b(3,1); 6 I3 L# Q+ s+ \6 D; H. E0 j

1 [' {7 \$ W# z i8 \

A = 10,-19,-2, 2 X* p7 v( b1 a; v l / K* G" @& o# ], ~" t' o; [7 C

4 M6 i) y' H" m1 g9 i

-20, 40, 1, % m/ l% \# {8 F8 @6 Z ! a# r: A: W0 j1 f% O7 e % @- Q. {8 L% E

1 `. }* r% v4 ~% t* y( G

1, 4, 5; ( T/ P6 U3 q9 v% q$ [. |" f . D- C, S+ e& X+ G. D9 g

, W9 ^! o' x( b8 U/ w5 _: F* w

a$ m3 m" U0 X! e

) |8 }1 Z# Y3 u. v% \4 X5 g

2 E/ i. G% W8 }/ b$ C) Z

b = 3, / S& i7 y# q! P3 k2 |. @) U

& N, B9 y' S. q0 o

4, ( l/ O- m2 Y. s7 }. v& t/ I7 p& i4 R + |7 P6 h' t( `9 A# ]* I$ F 3 Q2 F! _* G! n" d* z

- T" z: a' a& M: e) ?9 F

) C" Y/ z9 ^4 Y) z5 c! a0 d

5; 7 R& p+ ]1 n5 W3 \* I: \9 ? . n1 Y+ A# F( l' {% \# ?5 a

9 C, U$ K; `% U. d5 ~

. L7 q6 \: k0 `* f" ^5 K

Gauss_Jordan(A, b); 2 B% Z, ~" g( j! [ * Y$ H. c* F& m* u. D5 S

, }) V4 T8 N; Y

& B2 |* v/ j; G" S3 ` O: q " D3 t9 M) V T- S7 H5 r

0 v4 U; U! c6 J4 S) v- }, T

0 h' c( d6 D' c4 y1 s1 U8 k) v' ]6 @4 @

cout << "Solution = " << b <<endl; - m& A- \6 \# F: U: f% r , y5 u6 |# o Z6 _

' [- M1 Z6 x) X

; C N/ z+ t9 w2 c# M/ }

} # l: M& }: U# b- L' _% Q

0 p) n* W' T2 i: @. B0 W- y# G

3 ?% O- q: K# S, n( i: ~

8 d# _7 N" V5 T* I! ?$ {& t3 |

Result： + [# Z3 P; z% W( h; E9 F

h4 l1 [: m" I4 A8 @! s5 o( X

; V+ }3 L: ~' i, X1 I% R; O

Solution = 3 x 1 3 H: d' I+ A9 p% y* m6 S# P* | # ?8 R8 F! R! n3 g) B4 v3 x

( v' X2 w1 o5 ?% |* o, z

[ 4.41637 2 o2 C( e" c# y1 x0 d* S

( _+ m2 S* s2 `

4 I2 i. j2 h3 H6 j/ f8 p

2.35231 + W5 Z# a: o) ~( d$ ^* [5 E! t% ~

, p5 h+ t# k: ]2 ]0 t

g0 d5 q3 Q! S: X

-1.76512 ] 3 T' h; n. _3 W M' [ ^. u & b8 s! r0 X6 M# q$ Q+ k " m Y8 d- ~9 l1 q

: `: i' }: ?. |

. @ b/ y& a, P$ _! O

! m9 p$ ~( V& Q2 @1 f# _

& R9 A7 p, F/ z7 e! q( Z

, X/ G# h# t) H+ m1 Q( u3 m$ p

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 . ?4 O) x( i8 a- N3 m" b: S

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

$ X! \) F8 l5 x& @2 n G7 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

路过。。。