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

lckboy

全主元Gauss-Jordan消元法

% q1 k* m* ]2 w }" t$ j

- c6 F7 @" X$ G6 }' X; Y' ]7 o

& M+ m6 G S' S' l

: h' y5 l0 O4 y

3 I" [) e/ D5 Z! |& @8 d. H

% Q0 s, V/ l3 i8 Z; a3 e) c

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

) N7 |" ?9 V; Y9 A/ f/ Z

; W% k9 ~4 z1 R f* m

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

8 ]9 }7 Y/ B2 a# u, M( p8 Z/ H2 y

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。 9 H$ |3 k! E/ [( x0 O, Q, V4 \

U- x. P" ]% [! l9 F* s% c4 g0 {

# ?# n$ V7 f6 s+ e# n0 j) h9 ?

3 ?9 U+ a, y' ]# g

% l) ?# R w3 z! n+ c

Code： 3 I K0 N, U L0 b 4 @/ {( f; q6 [: L/ R8 ?, i

$ d) A8 d: Y2 b3 H* Q

1 l. W" ]" A) C; K/ i! ~' ^

6 Y' V7 y5 L7 s3 S& J

#include <blitz/array.h> 7 B& x! T4 ^$ w% j ' X" J( x& u6 m0 e/ W

/ o7 t& U( u& w: L! c! G! y2 w4 R4 h+ G

" n0 }7 q: _. v

#include <cstdlib> " V# J2 k' p+ y, y u: v( c% b * L2 g* |1 h9 u6 W" a- Z

8 \) @3 U8 Z# Q& o \7 E$ ]7 o

#include <algorithm> 3 v, k6 ~: {6 h [5 h1 W# X) s7 m

4 w1 z6 J/ s Y4 E

* A7 c- @1 b! i1 o5 G8 A

#include <vector>

6 n. I. L2 H& @3 j$ F

using namespace blitz; ; I2 P: C9 I( X# n ^ 2 N, E, K& s6 s8 s

8 \# X/ W x% t+ r7 P

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) # G& K6 [% p* p9 j. } u [

5 b8 _. W+ x# O4 S, u i+ [

{ # [( M1 J( U% i' ?% b) m % s& R9 z" M1 E

- Y: |7 C; l) l4 G& J

int n = A.rows(), m = b.cols(); . \ v5 n1 x& ^) j ; B. }1 L& c" k) k: ^# |6 n6 y 0 c t( R+ k* B/ o

8 \- T9 F* {! x# a1 Z

int irow, icol; 0 C; ]) [ v" {, \% i& Z& f # k4 ]$ d/ c e" J/ Q3 p

@6 e7 K: G& c! ]

vector<int> indexcol(n), indexrow(n), piv(n); 9 k; o" \+ q B$ { 3 d! T- m; V- i$ N8 g t" D2 g

1 l3 Y9 M2 s+ `; I- R5 q6 d

2 f, r: o9 g, z) e- y' q; p

6 l) j( ?' L* _$ h) Y+ U9 F# a/ \

for (int j=0; j<n; ++j) 4 [0 b6 i* U6 x# w 0 |9 E+ k% D& A6 _0 g + D6 T* o) u# E6 S1 u5 Q) ^

2 {. W. K( G# {# l1 \4 t0 W

. G4 \' A S; { f0 F& a

piv.at(j) = 0; / @8 r( Z# O" C: |; C4 {

4 ?1 d& z6 }& m

7 B* h8 P% i& u+ G 8 T* K. S/ C: a+ b

: e8 u. ?# Q5 f: X' z5 \

9 U) `& b* I$ d" V- b9 j- U! V+ G

//寻找绝对值最大的元素作为主元 : R( |' b8 h7 J8 e/ E5 t2 Z$ J9 q

; [4 y( B* F6 }) l/ o+ M. B

for (int i=0; i<n; ++i) {

: I( n! Y4 A) M0 v: J

double big = 0.0; 8 }+ w) Y/ ^" k + S* |5 ]: U; a+ \/ N

# P$ x' G( S) W: |) u$ j8 R, c7 v

( U# {% u5 t& A% @

for (int j=0; j<n; ++j) ' _* T- B" A9 _6 U1 T" J 4 m( b, g1 A: H2 w" ^

+ F* H0 }( o7 V: Y* w

if (piv.at(j) != 1) 9 D$ B" e5 k3 @

5 ]3 @$ L+ Q! j4 Y4 D& e

for (int k=0; k<n; ++k) { 3 F0 [. H4 j: K9 d% F' t ' r, q& Z. o; ^% R: T" h: p, [

if (piv.at(k) == 0) { - ~3 r2 M3 _/ B: X4 J/ N, x: ^0 _* w) e

5 \5 Z. r- s; [0 A& H7 k

/ w/ z, l5 s ^4 X: s1 C

if (abs(A(j, k)) >= big) { - x- Y# d) ?: w) W% y$ y# n , l! @. p+ U( o: c0 m3 H

6 I5 x+ J* L/ c& R* T

+ s+ d9 Y) J0 l

big = abs(A(j, k));

* a. [0 \9 s: E7 ^0 u

9 ^5 O5 F( @) m4 T' a( |4 S; v

irow = j; 8 y' }' u6 ^+ Q* e

' F& O; M7 t: b" K# R7 |

icol = k; ) l% ]4 l9 p) s4 {8 D) u) A + G/ D8 m8 t+ i; U! d, \- H

3 n w0 c6 z& m9 |3 x0 i. z

if (irow == icol) break; ( L- S+ w+ Q; e+ n. H M1 H) l+ D/ J: j+ c( x. Z ! w9 W6 Z5 G: t2 G2 h; w4 L

% q/ G7 y# r8 |. D

4 G# x" m( o1 M- U

} 7 E6 P8 L3 [1 Y$ d- | # Z- |. D% D% p! N6 [6 ^8 V

8 g7 x6 W: A0 x% j" D7 ^

} - O3 }) [0 K' m; \2 g

( q, _8 Q1 ]. j5 Q, L$ a5 h6 {

} 7 J7 f' k( V* X8 W8 w) P+ M

% O& D! K& N6 n3 m Y+ G

- | R% \, X. }/ W. [% S

5 D- V! D: E3 P6 N) J1 _ p

++piv.at(icol); 1 ~ u3 \) K# r/ b/ x( e 2 ^, [1 i8 N: l) L

# }) d- ]! \) v5 {1 ^- n

5 K; y5 G5 |, D4 K: ~) ]1 k. I

//进行行交换，把主元放在对角线位置上，列进行假交换， " \2 W3 R9 \: T 1 B+ E5 J9 k l$ d0 ^2 X* ^" z

//使用向量indexrow和indexcol记录主元位置， * V# `2 K8 p$ v- k8 O 5 }3 v, d' Q2 S$ c

//这样就可以得到最终次序是正确的解向量。

: m# U( D7 b/ T; C1 A

7 y/ m4 N% l/ k+ w0 O

if (irow != icol) { $ s, P; D- N; I: { 8 a6 K) }1 N5 o f

" z# I& q) ?, w; {) x# {2 D

for (int l=0; l<n; ++l) 6 r8 d# x T. n6 u6 u& ]& k9 A F( m ^$ r- i8 p- a! ^ 7 e# d7 J z _. e# g

4 d7 y) g4 E% \0 r9 g

; l8 x( S" L% E8 z6 S

swap(A(irow, l), A(icol, l)); ' W: q: F y5 j" j! C1 ]

6 P9 \; M( S- r0 x5 l# }

$ k- i8 s7 _1 R2 s9 u' o8 T- y

0 p Z9 Z0 Y9 c o8 {4 T% s

5 }0 N& J- ^" U& u5 C3 V1 Q

for (int l=0; l<m; ++l) 7 A2 [5 b: a$ J% ^, l0 V4 u, j9 A $ ], U2 j9 v! W* [9 o) Q

4 U; \8 H% |7 m

swap(b(irow, l), b(icol, l)); % r! d& b, X: M) a) |$ r$ {

$ q5 {9 f0 O8 `7 c. ]$ {: u/ f+ ?" a/ D

& o8 V* Q: u, b g% |6 `- V J

}

% g2 [4 w( L! K! C

indexrow.at(i) = irow; ) {" C* |& O, }1 x l" F7 b ; {! M# s7 p& n$ |7 k

indexcol.at(i) = icol; ( N5 ?3 y$ a$ T7 U9 z2 y! y & U8 L7 M3 w: m8 Q+ A

/ U- k5 n- U. s, B

5 ?7 `8 q9 N) C8 C, M* {

9 s( ?' s$ c8 H8 n6 b& O, h / D) ~+ ]5 R3 W, C- c n2 d 3 w" J) C3 x4 }& T! N- g

9 [# F! V$ s. d6 c2 X

try { / b9 E* n, `# J

double pivinv = 1.0 / A(icol, icol); * I. r" h. _- p) p8 H& A# ~. B

! u7 g5 {8 ]( {% ?' j

for (int l=0; l<n; ++l) 4 `* x7 l/ M% F& j0 G$ [

A(icol, l) *= pivinv; 3 s& {+ L: o) w* ~3 Q! o/ I % ^# r9 R+ a* {7 Q& ~4 G! ^+ t6 T0 u: t# R

# n8 j) Z. B* S

for (int l=0; l<m; ++l) , b, Y. C# m1 [! b ^# o ~% l + q$ x. A# q" U! I

b(icol, l) *= pivinv; 6 r0 _0 L9 k. G8 C5 G 4 ]# C4 T- X" v( i / H# d! d6 x# U8 _ ^' d

0 L' a% X9 P5 ]

( n; h" q6 {5 [( B5 J$ Y

//进行行约化 9 G1 C) u" X- e' }

4 q. e- V, @0 \' w6 Z; u3 h

; u0 `3 _/ f1 w3 A* U1 n6 e

for (int ll=0; ll<n; ++ll)

6 K' \! i" @$ ^5 f3 t( C

if (ll != icol) { 7 N! p# A B2 c

+ X% i& k7 ?9 t0 L

double dum = A(ll, icol); , k, @4 f) s6 |( ` ' i9 n( F1 X' Z7 A

; }2 f; y k+ P0 k: o

/ Z h/ O% I: ?2 \

! k) p0 g5 }* [! `" a. F b

for (int l=0; l<n; ++l) 4 ~ q' V3 N- m 2 B8 i; Q6 y: I, A% r* m* Q- U ! N9 R) k+ k$ S, k/ E, f7 y

r2 @9 ]3 b; B }8 [

A(ll, l) -= A(icol, l)*dum; : [( s/ a9 x* S& t! m" S+ p& w $ @+ `& M) \/ \0 B

1 V, M3 ?$ T5 z. E" g% c

for (int l=0; l<m; ++l) 8 S' p1 K' \+ o" p2 O. U( K 3 N" y# A: V6 Z+ v7 f

6 j* i& K# I% m% u2 ]/ \$ ?

b(ll, l) -= b(icol, l)*dum; " j7 U: c: @" h1 Q6 g t/ p6 C t+ Q0 w! f

} ) m) a# u& C& `( h5 B+ ` 6 B z* t8 ^) G! G G. o+ [( b# j8 M- _. k- _& Z0 H

: K& P! P- b* `" @ z

9 _* n+ F+ W% e( H" e( W

} ! ]3 \3 M5 Z! N* J # D p' n1 [& x# D

catch (...) { ' \, F8 p; ]7 \. y7 e 3 U; c( [% P& b5 a# Y

# T; L5 X; R# b1 v$ h2 X9 {

3 }' y* @% H( d# E

cerr << "Singular Matrix"; ) |6 X# O2 x; ~% v$ Z) A. S6 S1 } 2 ?" y; S' M$ G

5 b6 I2 u. X- |+ c# `! S

} / k, A2 a8 c8 L% `7 j " k4 Q( E% g+ i& Y r! k

% J: K* d2 c D$ C/ K, _ N& o

} 7 C' t" e2 [- }" g 6 l2 g: G# @7 l' l% J2 t

} T6 Z; e5 [5 W + U' n/ _) q, |7 A5 ]; L

/ }" z8 K& ], f

& t" n) L. r( H0 p

0 E, q) T7 d' g* v% W; c# y

int main() I5 q* I6 k6 L' F3 l4 G/ e' Z

/ z( p: n( d0 z% K W& W+ A1 H

5 E, p+ m1 C. ?5 I/ s

{ ' J+ e$ T* u) Y1 `. m/ X& l

( o2 y, k2 y4 Q; s% a

//测试矩阵 0 A0 E# v) X) ^4 M

; ^2 i9 ~3 Z. ?( r4 l+ N% x- o

Array<double, 2> A(3,3), b(3,1); 7 h* j4 u% ~) b4 {) B

' B1 l5 T& I A7 T: V

A = 10,-19,-2, " y1 i8 {2 M8 S ; e3 {2 T- z8 J + k# l0 Y9 {9 h" K( `: E J5 E

/ ?# F6 d! o* U; J6 z

-20, 40, 1, " L# q' l" U$ U& W H- p. J. Z2 [2 P8 p/ k5 F, J

, p4 L8 _. L3 \

; u; t! T' _' F/ s/ ~4 J9 M

1, 4, 5; * L: b) L! W3 U7 E7 [ W

4 T: d9 R0 \5 p% O+ m

. G# q! j# m* M8 R5 I/ P$ R

5 `7 Y4 {' y% R9 E0 U. P- n

b = 3,

4, ! l+ E* M$ P' Y* I$ K1 ~

5 p) F: f1 C, C# W2 j' |( c

7 f' e; C( y/ P

5; % Q- K- @8 F$ j' F

! R! p, ]/ k$ V# Z( B; U9 @3 Y0 C

. M/ N- N# X' X2 S; w

) U& K; ~& }# H! Y

- v- w$ E( a" z

Gauss_Jordan(A, b);

6 d! o! A, n! P5 C7 Q; V+ C

7 |& ~* c1 V3 M4 R, Q+ s' g8 R* \

) l% x# i% @" w- w3 W9 P5 l

6 L' {* a- f* a# \

8 [2 N% k9 ]' j/ b! y: O

cout << "Solution = " << b <<endl;

- }% c2 p1 N' Q" E6 Z5 p

( `* Z# M7 V0 e9 M0 T: @. @. H6 O

} 9 Q( A/ Z! r8 M/ I( m

6 P T5 P3 Q+ T0 R& W2 \* i4 g

& a; e6 Q% T4 E; E6 `% K

8 y0 e1 I7 j0 d: P! Q: O' x1 m

5 W$ b* |2 _8 V; z& X: A* q; D

Result： 8 C! q- j" V: Q + @0 T' l6 o) A/ l# c" l. P

* `) Q8 ^3 o: x8 {' C; q$ ~' t) D

8 k9 j; w) V: X0 {5 N7 ^

9 b+ Q; j6 L& r. ?* K

Solution = 3 x 1 4 C1 H f2 t5 N8 Q- C' S ! H" z1 l( G; w8 \8 w, R

( M- X, u2 W+ W. Q4 Y6 ~; o3 y

[ 4.41637 8 \* l, c: r+ X; T& S: c; s! L: A 1 P* S$ V! ^* P2 r1 s; ~

! g X) Z* o! L

2.35231 1 q8 X, \8 C' K5 h

0 @. R( A$ {& _2 V. A3 c9 m, R$ N

: T, @9 Y( B+ `8 k

-1.76512 ] & \6 ?1 z; C0 x 6 a' `; `& ]5 |2 N) n$ z2 x! y

+ t, `5 s! S1 ^- ?- ? ]5 e

. z$ v: L# r9 H1 K

3 e" h+ H( Q' T5 L* ~& S

/ [) J+ j9 o0 ^ l' {# I. b5 H

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

n2 m2 u( z- Q9 n D9 A2 L

' V( k0 R* I5 Z+ C6 e+ @

1 h4 X) M* F+ b4 j+ _1 G: X- b

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

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

路过。。。