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

lckboy

全主元Gauss-Jordan消元法 0 L% Z2 B" |7 W7 |

p+ B0 L7 o! w# n5 E

8 i+ F. l8 \, D3 E+ b9 J

0 T# E) r8 L" A& N2 U' u7 P: A

; \/ K2 s/ A! d9 p* [# H% U

( U7 }7 C$ M" n% G9 U$ j2 T# ~

Gauss-Jordan消元法是经典的线性方程组A·X=b求解方法，该方法的优点是稳定，而全主元[1]法可能更加稳定，同时，它也存在着弱点——求解速度可能较慢。 - L K! x2 ]8 Q6 z0 V2 j; c! R" q8 ?

+ G1 b2 P! z" O! V X/ } s$ N

: v, `7 ]$ F) }: i6 V

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

2 W5 R' e% h |0 C, ?

1 M6 b8 ]; D6 z; S h* C1 |5 k

下面，我使用Blitz++的Array库，编写全主元Gauss-Jordan消元法。 E1 Y9 S9 l8 B' S

; _5 E+ Q" f, g

8 y( |6 F a1 f, ?+ n$ c

Code： R9 O% U* K9 i2 |* J

2 H+ D# B( K& u7 N+ W

1 q# ^+ [1 X* G; v

4 A9 Q% D0 J9 V& g, I+ |& J

2 c8 F8 k2 m% T' q8 {; ^! R' y* \' B

- h$ m$ F8 c$ I# R4 x. i

#include <blitz/array.h> $ X! R/ ^6 ?8 I5 k o

#include <cstdlib> ) ?$ x; V7 r5 f6 x1 u1 k # `9 V, P M% u* e& l! F3 U ' B( e- R4 G* F z

#include <algorithm> . z$ t6 ]9 T! M; T- ~ 0 K) S7 K, o( B1 z% t3 m1 B

9 w0 ~9 k* c$ R0 @& h

# Q5 G9 w" p3 m! I5 Q' e

#include <vector> ; b0 D! `( K) ~# R, W4 O. i$ _

2 @+ A) T5 {/ ^2 u7 f* y: Q# M

using namespace blitz; 4 \ v8 [3 ?, y7 W9 A9 d2 d

* |: q' y; t5 h( Q

6 k& u! e3 {/ a2 B, y! N" q! J

" l8 G3 p) f) B: |' ]' I! r

4 P* h; F% q# y, Q( _% m

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) 6 r) C! r) l9 p9 w% J7 n/ ]! H* ] 6 g7 }' \9 F o. D

& A! X. H2 d+ p. o5 {

{ % Q! B, \. N! A

int n = A.rows(), m = b.cols(); , v% z$ g0 N) M

. _+ K$ z! e/ e& G

3 b! X& ~' K5 w, r1 r/ V$ M8 ?/ \

int irow, icol;

( @( g5 l7 H0 d4 N" h. p9 d r' F

4 |6 [. Z4 Y; i* b0 x( V7 I

vector<int> indexcol(n), indexrow(n), piv(n); & ~# F) a; M" \3 _8 t . C1 r, g2 c4 j! Y

. @' I% B3 s& z4 ~! F ?: H

; D+ U1 I8 U# J) Y j

9 |! `( K- Q0 Y. x

6 J4 o. ?& W! u# t0 ?

for (int j=0; j<n; ++j) 7 x, U1 G2 Q& D& z# J ( v$ H9 g1 z& d8 p8 o

piv.at(j) = 0; 9 M6 b) p ?$ [" \8 y- M0 N* h1 U 8 j- ^4 k0 L$ h' K' C$ _

! I% N; K$ m, U2 m

5 Z, ^! L( D: x ( J. ^# o/ c. P) N/ e0 G. ]$ {

//寻找绝对值最大的元素作为主元

* v: N: g" t- F. u' @5 ]/ c

S! ^- p; `# }

for (int i=0; i<n; ++i) { 5 t' a3 F) H" M' g 2 G5 q5 s6 E0 a3 T2 ?0 h+ F

# v' z5 Y5 O% R* z5 M# a4 A' `9 I' ]

' ^* E. m% ?+ s1 b5 S

double big = 0.0; ' E7 `5 X$ t ^# M 3 h3 V Y( c% v: P. A S" c4 M

4 @% p- q+ W) ^! F% K7 j" w, M; a$ f

# R& {9 F! b, \3 g) G9 U- T: m6 O

) ?) O1 n7 v+ H

% x, I! A2 S) h6 B

for (int j=0; j<n; ++j) & o- m6 a$ _# I1 l! w ) }+ r% A& y! X7 g" m6 u

if (piv.at(j) != 1) * w3 |0 C; [$ D$ i! g {! z

$ M* ^' k) h+ U7 ^

for (int k=0; k<n; ++k) { " e+ V* s! s$ B8 z I$ e/ F$ | & m9 A" I% u/ b6 Z& S/ L" ^8 A ; [+ B/ Z, S# t" m% w

if (piv.at(k) == 0) { $ ?9 U6 ^1 z9 ?# Y' m4 I* | - q$ C, E+ `* V

7 p" ^2 J6 V/ Q

: o E2 ~% W3 M( A9 Y3 c7 f1 G

if (abs(A(j, k)) >= big) { 2 z4 @* \# ]: z* G3 a

big = abs(A(j, k)); 3 G6 H0 K: ^; V3 e4 x# d

% M! u3 e% R3 h9 v, w

. o% u5 ~) F# ^1 |: S3 y, f6 ~

irow = j; 9 T) @. W+ O4 k

8 {# z. @2 M4 y6 a0 q" G

icol = k; ( a7 Z5 a$ x6 A

4 v! I& ~3 L6 H7 x Z

if (irow == icol) break;

: |3 h- l3 D, V

} + l8 {+ j) s4 S4 V& E $ }* I4 F0 l* k$ q5 f' j

- X' f6 n( m- T8 ^1 S5 m

" c. t) d0 `4 f1 h

} 0 X- Z4 M5 r# I4 _- M( x [- x

/ V' i- {7 }! z% t0 |. {1 b

, j' N2 |$ h- q

} % m$ S# z1 A; L

2 J1 g8 ~& V- y' ?

9 W2 k( w$ y1 s# B( q/ ?

2 ^4 ~& ^% A/ h% B

++piv.at(icol); 6 v* _3 Q' u2 p0 V1 |5 d \" B

! A+ w- B2 y& M9 ` 6 L9 i& N, H+ Z. y& P9 p

//进行行交换，把主元放在对角线位置上，列进行假交换， : N- Z, s. \: n$ O3 ?- t% E 5 v5 c% B" } x6 H# M- `8 W

: o& l4 Z, [. O2 I

//使用向量indexrow和indexcol记录主元位置， 3 }8 b. T' x9 f4 F2 a' U' {1 L 7 f6 Y3 u6 z8 @ s1 T

//这样就可以得到最终次序是正确的解向量。 : t* s! N/ r9 |: T$ \ 9 q7 w. Z/ I& B

( C7 s4 W- ^: h

( } H8 a3 e# b4 D

if (irow != icol) { ; n- T d3 f/ F) T8 U 6 v" @1 r$ n" G1 G

5 B" U7 w# K7 Q, o( o

& x- V b1 b5 n" ]6 P

for (int l=0; l<n; ++l) 7 C4 v* Q+ W, z1 a2 |* L+ u / { t- m8 a5 m

5 R- A# Y- `" d( W1 y |1 o

swap(A(irow, l), A(icol, l)); 5 X# N. V x( M7 D8 o

" i; I+ H" a& `& O

! x- m* A$ i# P0 J/ Y5 `+ _

for (int l=0; l<m; ++l) 9 o- K! X6 H* U K, |0 | 7 S) ~% m, G/ R T7 g! i$ L8 [

# h0 O5 E( V5 w2 ~% z B

" K; J# z, ~8 @. W- S4 l

swap(b(irow, l), b(icol, l)); 0 w& @4 I* N2 h3 K% n! ^2 L

} ! f# o. a2 E5 }6 ~+ g: n( L 7 G+ N: V7 G9 } . L6 s, R: ^8 m1 O

" g9 R8 \* U7 [, |: |3 x. t+ [( y" r

0 o/ F: n3 t/ n

$ S& l# [ V1 m, ]. \! w# e& t) T

indexrow.at(i) = irow; % @3 ?/ v- g: C ]0 s0 P# s# B5 D 1 J; H$ r" O1 m V

- Z |4 h2 Q9 t- R; z2 U7 l

indexcol.at(i) = icol; : q, d% V; ^; ]1 Y f; r! i7 N5 V8 [! l) R' q6 T

5 c. i3 t$ J, O# _. h

) J. v Z! j5 a F5 I$ ~+ t2 u

# ?! d4 p* ^! Z5 G6 K1 M" ]

3 u0 i$ Z( l4 w* K$ X- N

try { ( ~4 x; l1 Q. `/ k" t $ [8 \0 k0 T- ^1 e8 m4 F' \

+ g1 w+ @2 E0 P+ }9 p( `2 j' b- E

( q6 t. j! \* j+ f5 C7 i- v1 M3 b0 Y

double pivinv = 1.0 / A(icol, icol); 8 O2 R# d7 A9 q T - D. s, u* g" Q: B " y G& V4 c% \$ U6 n8 d4 l' u7 u

1 I3 Z9 @# G8 B/ W1 P

; O d) f) q. K& ? I4 a! e8 `( H

8 q6 L! Q# P" ~9 k8 c8 F

for (int l=0; l<n; ++l) * e/ y% I$ ]+ H/ n& D' K1 N3 W3 a0 j; M

A(icol, l) *= pivinv;

5 @' E( B! E: O. B- G

for (int l=0; l<m; ++l) - k" M, `" X9 r7 Y F0 d, g& ~0 ]

- t6 t* K: K1 ]; {% W- ^1 H4 ]

b(icol, l) *= pivinv; * S" G3 {6 U3 q+ ?! b7 x7 K: T. ] * l5 g: `0 L8 d) z# R3 R

: }8 ~" D0 T( d$ _0 h

/ ~& v* |6 ?/ K: R' f# d

7 m& I7 n/ X. }# t* e: Q+ e. Q

/ u5 n, x" i( l' X" D

//进行行约化 ( A- Y1 H, j9 c# x/ G# _ } 3 u8 j" b: k1 |# n4 M : f4 }& B# }0 m$ x4 P7 g: o

( C) r3 r( x2 j$ k# X

( x0 R( I$ Q# E$ K' Q5 A; D# u; u

for (int ll=0; ll<n; ++ll) 1 G+ }' B& w& X5 u 4 a. w' q, g) J. P/ X

' P5 Q# Y c8 S2 Q! t: i6 g5 a) R

if (ll != icol) { % O9 b7 D \# C: u/ y& G

' o% W0 I( {! G. J7 Q2 { G# D

double dum = A(ll, icol); + u, p2 U. K/ j& |8 k& i # h5 @; I* I6 |5 j* L% V! _ 0 z5 ]$ U" n; \4 a, K

1 `; ~" Q5 d( v

. U' q; k6 F0 A) l: _5 O M: o

% i: P% D- _( I

for (int l=0; l<n; ++l) J) \: M ^, I 9 x0 R7 C1 @7 F5 q

3 M* W/ ? R4 X; }" g; Y3 O6 G7 R

A(ll, l) -= A(icol, l)*dum; ; J; e' v% f2 L+ W, y$ K

for (int l=0; l<m; ++l) ! K% Z' V" l; U: q; N4 R

b(ll, l) -= b(icol, l)*dum; ; Y0 ]! T Z, |3 ?+ j# f& |. L : d( } |, |# K" C3 h# ~3 T

2 m8 j) [* u( i0 |3 j7 _: \ b* M

} 1 A( Z8 F6 ~0 S5 g0 k

# N6 ~! l$ j3 z6 I- s% y

( G0 N$ q* c {

} 6 s9 U: _' H* ]- H! S3 h5 J

4 L/ ~) U' J6 @. `$ \

catch (...) {

4 S& ?# t0 s+ t9 U

cerr << "Singular Matrix"; ( n- r. t5 r1 ^" a) C/ D) F * |. D" g# I0 {% }4 c; ? - \" @( C$ c- ?% w8 k, h

} # j; V3 a6 [& B; d# J' G0 Y 8 { O' y: j( E* f

} + m g8 O% Q( O: Y0 U' @ ; S' `+ q, G3 ?; a9 |

} 1 P% U5 S# Z' S3 h9 p

) |- e h7 D$ D5 I7 x; j

' t/ i, _3 p+ N3 O, Q: m' q

int main() 2 t% Z( @) j* C' ]! l" e 9 e0 x3 F O5 m; g1 @" v/ |

{ 8 l* D) M0 P0 H: n . h" s+ k3 Q: A8 m1 Y9 ? 1 W7 w) n. \5 d+ c& Z

+ Y" W/ V+ B0 z. o o. L

//测试矩阵 ! u0 R& y" d* ]1 r$ i' l' Y

) C8 d V$ p( Y+ x- h6 S7 w

& Z3 U* @1 m0 h; S' \8 a' Z

Array<double, 2> A(3,3), b(3,1); $ ]' n) f: h9 p+ t6 W2 V6 p # J! ?/ x" X: K. J1 R1 n3 @( @

A = 10,-19,-2, " k$ T. g1 H; S [8 u. \6 H" m, p ' d( V/ Q; l5 I8 `; e* v6 K

; _ N" F3 P5 f2 Q n1 Z( N* [

-20, 40, 1, s1 J% g8 N/ I

! @" \$ A& _& }; @" z6 I

1, 4, 5; 3 s0 k; i; _/ i; R! w* F& z ( x7 b: x+ U# } 4 L8 H1 A( q% r Q) n* t

" }9 C/ ^+ D' v; @. t- r4 b9 Y

: h5 |1 M5 A0 d. G m

+ L. ~. g. b$ |

b = 3,

6 G! q7 ?" U9 b0 ?( J+ ^+ O: h

* `4 ~ m- D/ O

4, ( h) x6 K8 S$ q1 f* S9 q 0 b% F" ?( L% j. Z) Q7 R( A7 Q

" g, H+ j; D$ S1 Y3 I

5; 3 R# e! G2 I2 i) q( S9 j

& L+ T# I" _, S1 y0 ~2 T1 d3 [

8 {- J9 z6 c0 P- A" q, L

Gauss_Jordan(A, b); # U4 D1 ` s. D( A1 M6 a5 @

8 m$ o% A& w/ G

) i' n0 W# @* l" a8 r

9 g/ T% A$ G2 Z1 ^8 A

cout << "Solution = " << b <<endl; ! K5 K }& @( x' ?& x' _ 4 n' v: n& a) I0 Q( c/ H3 o 0 B) j! B0 L1 V B4 N. U, J# u- ?

6 r, p1 Y6 C2 S, [/ r+ z$ @2 T

} ! W. T) s. F1 R" ` , S3 `( C J" ? ! E- V3 w: Q1 @4 O) i. a7 K

$ Q4 W) \' R; C* W, m3 J

+ T4 }3 T+ D- l% i

; G% J4 d+ V3 l5 M' b( h

t7 c ^ |' O) E, t

' Q8 h2 Y# O* a( e. n5 p/ s

Result：

8 k% i" b4 ]% M2 x, H1 e

$ D, R1 A5 v) m9 L6 I0 V5 L0 C

Solution = 3 x 1 . q! L+ ` F% b- h* y8 X3 d* z& W) t

7 h7 Y: f3 H: |- A# L

[ 4.41637 ! `& I2 P. h: B7 U6 l, \# E6 n ^

2.35231 ! O9 c# d4 L* y$ R0 `$ q# d" N : V; Z- U& ]- d/ \% W

: {, u- d+ d% j0 u( w0 l

$ j5 [" \7 E$ V

-1.76512 ] $ x" u; k6 B" ^' M" |. ^6 ? - A) @! V |9 w0 {

3 q2 W$ |, ^, x

1 R, c$ o8 W! Z* n1 }0 Y/ f* f, i

- d# G0 [' S; E4 A, X

6 H$ V) E8 R5 `. k4 u- Q

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 9 e# Y) ?2 J: |( p& S. T

6 P( N" B5 F$ J. u/ J# m( i

# j+ t) [1 k$ @1 a4 O

6 \' @. m) ]9 S8 l* K7 N$ c! `

8 O4 ^0 G2 d- O+ N# P

注释：[1]主元，又叫主元素，指用作除数的元素 V: @& q0 x. O1 g

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

路过。。。