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

lckboy

全主元Gauss-Jordan消元法 1 P" l- { ]$ a. c- Q

8 D* g; R Y9 p2 ?" k: `

* c5 ~. r7 p; S. }8 a

/ R3 z s9 b. v) Q* `; \

* W, ^& U8 N2 z; P. b

; p# r; L G/ N& f" e" ?& P

- Q* u0 g: ~) p' K8 n$ U3 k, F5 E' I! {

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

! m! q0 F6 k; a; U6 J ]

0 Z( m* z) q% t, D$ O3 N# `& Z

" q; `% f8 F. e8 F' T' L5 H4 V

7 q9 ~3 e' e/ [% T M' {7 o

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

0 c6 `4 R2 T: |& z9 w

7 @1 M3 N- r3 @% e r

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

Code： ( H* r% h9 }- @$ S ; X X( S, A9 u! Y+ {: u: r : ^ V3 V) X1 U1 U |

! ~: q: ]! i6 k3 t$ O$ O

8 {: P8 z7 `2 J5 p* Z$ `& A3 F& W

#include <blitz/array.h>

0 s! L' p. y! d# f! d7 L3 B

# N. X$ ?9 Q: D; O6 q

#include <cstdlib> ' L; Z! T/ f/ i$ X: |4 B( H 8 e5 K' o' G1 g0 S/ S: v

- K/ F& ?1 n7 p# T7 B

#include <algorithm> 1 {0 @; `7 A) {) j , L" ? L0 {% h4 O

; {$ D0 m& X; @' a

$ R& M- s. m8 x! U! f; W

#include <vector> + L/ t# s2 P1 t* T

- ]6 E: _. K/ @7 A9 o0 F9 R

using namespace blitz;

7 Q% |7 e" ?3 e ^0 E0 X

- ^8 A) {& H2 R: Q7 ~- U

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b)

7 D @, z7 D( p; T

{ 9 W! c+ i/ F: O' |- J , s+ h8 o" n2 _+ {+ k

3 J7 U4 p7 y! g* \3 l0 o

: ~! J$ I3 I) r4 g6 y1 _6 k) r) D

int n = A.rows(), m = b.cols(); 9 M1 l' x% n" l2 `- { , m I. c. ^6 I. q0 D

int irow, icol; , {; ^, W! l2 k6 a

" f3 C' U/ W5 V s

vector<int> indexcol(n), indexrow(n), piv(n); 8 G ^# N( t8 C8 g/ N: W

+ J* x" @* l. v( |& Q+ K) T W/ L

1 f: u4 V! {) A- ] W8 }

% T( M0 w1 ], W& _2 ?

for (int j=0; j<n; ++j) % Q/ ^4 I% Y+ t( H- b" P1 _ x |

- ]( C F# ?( Y* H$ P

1 q9 C8 _2 c w9 M' }

piv.at(j) = 0;

1 {, q+ \0 G! I; }

6 z2 W6 H7 C2 Z: e8 J% Y0 c ( b# e" {5 f. c; k. a

5 S0 j8 e, C: a# n1 g5 N0 c+ u

//寻找绝对值最大的元素作为主元 8 J9 ]- Q; {6 }) O8 R2 f3 M- ]

! i& @7 F! M A( @

4 U7 F2 S3 R9 b R3 ]9 G- z1 v+ U

for (int i=0; i<n; ++i) { 6 B% I$ x o# p4 d( n+ M; e0 u : J: Y% S( n" z& w3 e+ g

double big = 0.0; ' \/ g$ Q4 u" C# Q: c. _) L - h; J; w( B1 _5 [6 p T) L

5 b% B7 L; ^1 F3 p- M( Y1 h- i

3 V2 X8 C4 Y m& E0 ^

1 A* W9 C! p& B

for (int j=0; j<n; ++j) , C) z _* Z$ N% J+ R5 L+ S% F

: e) I1 o' w& S! @ J5 _7 N) x

if (piv.at(j) != 1) 6 `& }; s- s) q; e& H7 U0 e & B/ f7 g3 i( k# d & j" v6 h5 E7 S7 Z; E7 B, D, j0 U

for (int k=0; k<n; ++k) { ?3 r7 Y$ E" b, Y( F) ^ 4 m. D9 B- m; p% M8 t7 }0 b' }% O+ u* [

if (piv.at(k) == 0) { 8 q5 o5 _0 g6 S* ]& E/ a * U/ }0 N- k& f4 J# z( j

8 Y1 \& D# q* u) [+ D

( \# x5 X- s+ y* K q9 [

if (abs(A(j, k)) >= big) { 7 ?$ t7 X7 U% Z p- {, m3 ? 8 g. U: v- ?' a; q ' p% Q3 U& Q4 [2 p4 Z& q4 o' R

$ k8 O/ K6 p. L4 c

big = abs(A(j, k)); 9 E% v. e- C! e; x r% I

* w7 I# m+ ]$ B g( p

irow = j; & J( n `( m: K6 X& _ " ?5 {' d) n( G- C# O& \ 4 q% g9 z! x! F3 _9 v' P# d+ Z: W

icol = k; 8 t4 D, E5 x% j: x. b

" O4 S' ]. R+ j1 c, o# B4 o) G9 ~/ V( v

) }" [! i) E" ^% g7 l3 l9 R

if (irow == icol) break; ! f2 X( G q4 m0 s' ~6 f4 R 5 W4 t5 D, t- O! [

2 B8 C* Q# @ c6 z" ]$ ]7 N

} ; A/ r5 {* t5 U( |+ ?, f. f

* _ o7 F. x% R* A+ C6 r

} & c! U3 }: R* h1 u ( e5 D" K0 i1 _0 m* y

! ?. t+ v, W' G7 u

} 7 |* ^, Q5 e& K) e1 x% z

& Z6 L1 E/ W" ^7 O+ F

6 F' x1 y7 I2 F

++piv.at(icol); : ^$ {8 r: Y3 W2 W) Y 0 r/ R. K( n* P* A, o |; `, k( \6 k! T6 u4 V, _" c- O3 w1 b1 z

6 H0 k) n( q6 h1 r& ~0 w4 F! O

4 c* ]* K" K5 ^2 G) F, L # H; b3 h( ?* i! j. H0 H& T' [2 c3 B

//进行行交换，把主元放在对角线位置上，列进行假交换， " ~0 y9 m( Q9 K, A! [

: K) ]& m8 m5 D

//使用向量indexrow和indexcol记录主元位置， " }' a/ F2 n$ u0 p+ N4 y! \) ~. g: ~! t( w

' }( ~. W, F R: z/ V0 `2 I5 T

//这样就可以得到最终次序是正确的解向量。 ( v' \0 s$ n! X: a, I

if (irow != icol) { & z1 A9 V7 t( L/ \4 k$ a! d

& f$ [9 e" O* b4 g) M

for (int l=0; l<n; ++l) - U- O- I! X( R/ E1 t8 [

. D- D* i5 I# S3 [. {/ P

swap(A(irow, l), A(icol, l)); / Q3 E) {9 Q3 F: _* Y; o9 |4 T 9 y7 V/ c# V% J

. R7 r" `- V( y% V

& f/ N% n+ ?$ V/ l) r

. W0 c2 [, e8 j+ G1 T

4 ]" C$ L5 R5 z- Y# ]" N9 k; N

for (int l=0; l<m; ++l)

swap(b(irow, l), b(icol, l)); & M7 N* ~7 V8 j, `! O! b1 [

} I8 O) Q- L" T, @9 i& n! H S* e j ; O0 W+ C9 k% p/ c d4 I6 [! s 3 i$ c( i3 R, |6 ?& Q5 P2 t

$ w( w/ M. |7 v( K, Y

5 B% H, B4 |2 u; J, f: L/ D

% [& ]$ ^) W$ Y) w) u* N% ?

indexrow.at(i) = irow; d5 e9 b- h9 M1 M

- d+ |4 z8 i8 O" [( ?: j

indexcol.at(i) = icol; % F6 z f- d/ k9 d `" ?8 ?! l( |- z8 o, `9 i

4 `: l; M: t+ F/ G

5 Y( j: |+ G. _7 a' w- t 8 a5 Y7 B0 h( \+ a0 Y

- M0 E1 @9 Y& n4 T

try { : ~* j# n$ q* [) f& ?8 |, A9 E" d- d/ D $ l1 }7 z, k- K1 b# [" Z5 r) T1 P

* v: P2 c! G* \% D. }

double pivinv = 1.0 / A(icol, icol); / e) O% b: d5 @( S9 x3 r' U

, i) o/ f; }( @: o5 I" q `7 _

/ F3 g) z2 h( O8 b1 s+ e# g

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

4 j( l/ E3 g! N; I, b& `

A(icol, l) *= pivinv; ; o$ n/ g0 U3 ?% R4 L * Q% q$ z# s, q) n3 |( ~2 p

: Q6 U0 T0 j0 ]/ B \

for (int l=0; l<m; ++l) 4 ^1 j9 n) x- n0 {% Q' s 4 n0 z/ ]7 p9 S9 p- T. i# r/ }

b(icol, l) *= pivinv; 5 e2 F. U# l% D( ~/ h 8 H& i' G! W$ D

: C5 N. @) o, O7 M* W

+ G# U# I9 e. g: u v3 H

) ^" \2 `$ }# ?: ^

//进行行约化 f7 |% {/ w. `4 E2 L# U 9 t, ~: y+ r0 o7 I

. N( m) }$ Q8 P

3 t @8 a: U6 }5 r! R

for (int ll=0; ll<n; ++ll) % E& B' t" ^7 @6 d, E: y5 E 9 |& A7 ~4 U0 y/ \9 C# `7 D

if (ll != icol) { 9 r7 i2 h/ S& e

1 T- Z) u4 T7 S+ t8 U: b9 V6 M$ n5 I

% C2 {3 a/ l" H r/ J4 T

double dum = A(ll, icol); $ y+ P5 k) K) T

% o& a0 \+ L$ {: n5 p' y3 a b

; s7 i' ^4 z* q; S# P. ]

x/ b2 X1 j- p+ W9 r2 A

+ x% G0 N: L1 F( @! t2 j2 x

for (int l=0; l<n; ++l) 5 ], H7 t$ |' V4 M+ z; @3 N# p) {( r ' [- B( @( i' r- U, ?9 C

. ~; z6 v3 C" g" ^

( v5 u2 P( l7 B, ^7 f: ?# v9 v

A(ll, l) -= A(icol, l)*dum; 0 v2 Y0 R) I% g. |3 v! m * k( U/ i2 D0 {6 u) a

for (int l=0; l<m; ++l)

( ]' f6 G& h7 Q8 v$ Y3 O

b(ll, l) -= b(icol, l)*dum; - F L7 w# Z" d! @) G7 m7 g " V/ s3 u( W: \6 e1 q

} 3 J4 X! y$ O! w8 { & p: W O( Z; |# V, E/ O

; i. x6 R. q6 F e

} # w* F6 l' g8 {) A + `' Z' |* o( {% H( E& X" Y9 M ( e" g. v2 D9 i3 Q

8 q3 q% k8 z* A+ A6 s

catch (...) {

cerr << "Singular Matrix"; 2 Y$ G! X6 x7 p, n

0 g" p8 w4 A1 k Q6 R2 S6 X) I

} / F! D5 W% u+ ^2 H) \; d

} 1 a# [- w9 X; `4 g; G

8 l2 m, |" |" g, V( O1 p+ h0 K/ i( B

}

6 `5 v/ y6 E: @+ \/ q1 c8 {

[ D6 ~* A5 g2 X o, c

% Z0 n4 r3 {8 o i- e4 S

int main() 0 U& T/ |5 a6 v% v+ d 4 \$ F7 f. q) R- u& R4 K

3 N0 U9 s1 b! T$ v/ J9 y: E

{ & T" n& N5 }. K & i8 A' O5 W3 D) n# S+ n$ o

6 X/ t, P3 X) _$ K

//测试矩阵

Array<double, 2> A(3,3), b(3,1); 6 }" z; G" a6 x$ o0 x 6 V r6 ?4 ?# J9 V; R

A = 10,-19,-2, 9 J$ ^8 b! o% O! }# H6 W2 @' X2 G 5 h. e# z8 G: f9 a" S; {

4 e, e- g' N g& j2 i

-20, 40, 1, * _4 e7 S7 s5 R" k# D. p% I) y 8 |; A+ \6 Z c1 J' Y( ~

( {. c+ c `! g, O$ k

1, 4, 5; & H: S% o, D& x1 j* x8 V3 J # x0 t" P+ v" ~) D/ b' E4 f

) g3 g; O2 q1 |5 M9 t

4 @5 {8 ~$ {& H- K2 s6 F

) z* W8 [0 @. p7 y. G) k7 J

b = 3,

4, # P2 j9 T9 u' {3 V ; n+ _% H* y% P4 {3 h5 c5 v i( l

& `. \+ P" _9 w Z& j2 ]( a. A3 Z

1 P% _8 C$ _- f' h

5; 6 I( U6 }6 w% M! [9 o

* F: Z' @* R# m; W

& F+ M6 d; K- z: @+ D

" l+ c+ `* v/ H2 q1 v% [' e( t

% e: l7 E7 ^( H. V! c* W/ k

Gauss_Jordan(A, b); ( f, A) G3 J' B5 y+ R

' I' ?) s; O) V5 b7 c' A' ^

- \, ^0 i& X3 I0 n ) J2 n, u. S9 q( K- G, C

) z. ~7 \9 G4 W/ C7 Z3 \' @

% \) X, N4 l4 E3 M' _

cout << "Solution = " << b <<endl; + Q1 b+ G4 l( J

3 \9 E- ?3 g4 Q4 ~& \0 E5 S9 m/ X

} - p" r0 @4 m# T* V# k" Y: _# p

. p0 G( r% X5 w$ ~: u/ A0 u

3 q3 n+ [& q0 d

4 `5 a) K" F0 q1 V5 [

Result： ; [ `4 j* S$ O' O% s8 J 1 L0 @( U5 Q7 _3 s# i 8 n. x( R! X" E/ }' r

3 ?4 f- \. j* D' a8 f

% {' }6 l+ Q6 H

7 C- f2 @7 G2 L* X/ I

, P, O5 f% f% k, A: O$ n! \, ^! o

Solution = 3 x 1 ; b# w5 r' ?+ n9 H; I+ K8 v* J9 H

; ?! T" q) m3 k

[ 4.41637 & |0 u* n2 Y/ }+ ]9 @5 d! X

2.35231 2 B9 K9 \) G8 m7 M- ^5 L+ d) m" Y ' J' O0 \: F" F$ J* R- M7 t

' A% F1 E# u L1 `

-1.76512 ] : r, u7 K% O1 q& ?- P, J 2 H, a+ U: \& d% @' E5 ]; u) | - x1 U8 C, f0 w% R

% [8 s3 u% b$ H( U: d3 _6 _% I) n

k8 W& V0 U% g; ?& ?

$ p. X& z+ [+ G4 y4 c. l) m

( t8 @ a8 [ {" ?5 o

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 ! _: V" B" W* M- k, F

5 N. s& K/ v/ g ?! L1 Y

% u; E* \" y% u& p! P# \

# ], c" ?, J p8 I! j: S/ l- N

( D. b) A# B& a

% L( i+ W- x: |* ]

, \# Q+ Q* Y* B% [

% N- C% k5 W; J

注释：[1]主元，又叫主元素，指用作除数的元素 ) ~4 s' p( v0 A

3 a! _5 f! N: K- Q$ L& l/ S$ Q

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

路过。。。