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

lckboy

全主元Gauss-Jordan消元法

; K R1 F/ E- V2 v( [9 }

( C- E( A! x3 O5 q, C

% x! i. c+ {& j" \+ o# P

: |* W1 R4 C8 g1 _/ f- E

" E0 N" ] F! X, q5 _9 V. s: n& \

9 ~- c$ y& u7 N/ o0 U4 K' {$ n! R

) ?+ P" v4 S% F, T! \( e

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

& L& R9 M8 \3 h G& O/ {( H

5 T$ ^4 l5 ^2 {4 N' u

) G( _4 q2 J; y) W/ s9 e

* b# R' v+ M# I8 G; j& \

8 k/ \+ w5 X' B) f% Z

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

+ J( G# G( b) ]& Q4 H9 L# g

8 y# A4 p+ O" h8 \

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

. Q' o1 B9 k0 m

* i D* |3 P/ u* C

" q! T o% O5 @7 a$ r( R

Code： 4 p0 a# ?4 M6 C! W9 M - C& [& L8 |5 I! S7 a0 C1 {

+ ~+ S8 p9 f l5 p( f, N

" H4 |' S8 j" a; D1 k0 r

& m3 U) K- M6 J) ]; W

#include <blitz/array.h>

0 h' R0 j9 b4 J G1 \

#include <cstdlib> 2 g$ D i- z4 Y D. G9 ^! y( j$ m) @

: I" `) B9 {- F& }( K

#include <algorithm> & b1 K: Z E9 f ! z) K' U/ B) I5 ~6 g" k4 L

) X1 e: ?: f3 f9 P" Q) z

#include <vector>

|' `" i+ J7 b% ^) z

9 I& ^& ]' i) T- h! E

using namespace blitz; 9 |# v. P1 v& b/ B3 ?$ W5 T

' V( r3 I& x2 U( A

2 k# Q ] l% _, D2 _4 y

( V: {2 c; O6 m4 D; J8 {) q

void Gauss_Jordan (Array<double, 2>& A, Array<double, 2>& b) + d, ~4 `8 G) v6 U, S; Q 7 ~2 f/ s# F6 U1 M6 _

! _- Z$ K) A8 x t- }7 n% U# ^$ `! d

{ 6 Y! l5 g( J! w, R# s7 Y5 u/ t' t $ _" W3 I4 G x: i( V& g

int n = A.rows(), m = b.cols(); 5 R2 `" l& ], \- ?8 k

5 \% b7 v% ?* l2 L, @5 A

) L7 g2 O# U5 l- W: B5 W i6 s% T

int irow, icol; " j" g+ L: V" @0 T8 s! w7 _, `) T; t 6 R! K" x# k6 ?1 r

3 o6 L$ s' o% Q/ b

7 G; }2 {* H! v( I3 I3 ~2 A

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

0 @. K; D8 {' l* x# [" }7 O5 ^; o) Y

! o% s7 n' q% C8 n, C

* }6 N( i- a9 f+ Y7 @) N

6 w1 j! e3 e. K% A; Y

for (int j=0; j<n; ++j) / N* Q4 C) s/ y n6 T" [5 y, a! X

piv.at(j) = 0; `, H6 A" L" n+ q1 O, B- [. N+ D3 A : u" @% d( @2 k. ^# U

1 F. |3 e- D' z. ~

8 ~3 w3 ^' l- b+ z, G! y* g/ h & P( L. m1 t3 k/ K3 f/ y8 t

6 i# J2 p. m. W l% ^ |6 I3 ~

9 E; d* j* h1 P5 u8 ]; G

//寻找绝对值最大的元素作为主元 7 y" r, V' B1 x9 {7 q+ g: B$ C! c$ }& _

for (int i=0; i<n; ++i) { 4 n5 s1 X" Y$ K; F0 ^9 T* q ; c# ?/ L' A+ q0 B! J

/ \3 K, }9 j! Y; C+ a

double big = 0.0;

( I3 H0 E9 o9 O; m8 } E

0 T W! q( M% l/ V9 t7 r* C7 k

for (int j=0; j<n; ++j) : ] l/ q0 k" h2 l2 d& ^

, e! i" e" |% M- \5 i/ m1 K

if (piv.at(j) != 1) / J; ^1 Y- }$ Z. }2 n + X7 K& h& }! }( K: }9 W

for (int k=0; k<n; ++k) { 5 ?# F q- U2 ]6 G / K( k( Y/ C+ g& G( L

5 }/ q, m" F5 L/ U+ i

if (piv.at(k) == 0) { X1 y! s) E, m8 z) y- r$ |4 b & O% F/ C& B# s* k9 r, M& }

- y* @- N( q$ V0 H

if (abs(A(j, k)) >= big) { : ^& R4 [( P3 ~5 g8 N5 R1 ^

. g, J+ C3 f/ @2 z' ~: @

big = abs(A(j, k)); 1 M8 P8 g/ u9 [1 ` $ I8 x2 g7 |1 d4 Q% p

irow = j; . D& S! X7 z( j& Y! b- R" t+ T

icol = k; + M( ]( c& I. [" D8 p ; ?9 y: s8 y# ~+ m

6 l& O# R& u/ `5 `1 I0 A, {

if (irow == icol) break; $ @' `+ L, {" j( q, k7 e& J 6 ]& l! K7 Q e6 V3 L1 `) H3 F

. S {* |6 g. Z

7 m) T ^# v8 t& C! ~+ I

} 0 ~+ T% f7 k" G) d1 _! z: l- M I& d( c2 j1 c3 N2 P# j0 c$ @ ( \; p x r0 [( d2 o

} 6 {# ?1 e" V: b$ I; _ p& Z3 |" d) R! _/ G 5 t. [/ _; k" ^

} " Q [. Q4 C& S/ Q0 ~. I) ?

* y9 d1 Z/ U; b2 ?: n s

- ^" ~3 r1 k4 W' |# M7 D) y2 |

++piv.at(icol);

# @% m$ o! _8 I9 b% N& F. V& O/ ]

2 T5 v3 n3 P |$ m4 y

//进行行交换，把主元放在对角线位置上，列进行假交换， 8 v! p n6 D; C

2 I1 D7 o5 Z9 c

//使用向量indexrow和indexcol记录主元位置， Y3 {/ ]5 Q# r! u4 S; R

2 @) r+ q1 K8 w6 ^7 c

//这样就可以得到最终次序是正确的解向量。 0 @8 ~3 Y2 X( F0 S9 S) y; ? ' }) z9 U7 k4 s% I( \% m: J' a

7 ^' D% V. S: s4 j. N8 \

/ ?/ ~, u9 z( w6 b! h

if (irow != icol) { 2 d" K U$ _. w' q# h . y: u5 i0 |( S6 @4 h! Y 1 O7 P+ n4 d0 S) ^' Z q+ u

for (int l=0; l<n; ++l) 7 S; T3 b3 [* P6 ~8 o4 i s8 ] 1 I$ H" W1 b6 n& Z

5 G( m; i. @4 |7 G" {! W( q

swap(A(irow, l), A(icol, l)); 8 {- H; i1 d4 {, E5 w 2 S) }5 _7 O/ t, t& h 3 ~, |' b5 `6 J, m9 @! m2 G

, S' @. e# X* }7 b3 g4 y5 j. |

( N0 w3 F, E7 t2 m8 e

1 e: Y6 `" [3 h, `8 @

: D' ~' f' x7 H3 T5 h- f: S$ b

% J" D8 @- c% E& ]% J

for (int l=0; l<m; ++l) 0 m" v2 _' Z8 m( m

6 f; }; h: q5 l6 L4 s# Q. B" k9 [

. X6 u6 }, ~& k" }1 s8 T

swap(b(irow, l), b(icol, l));

} ! k, u( B t3 K. c/ V7 E% h5 h& a0 a

" L$ L/ G& f, w$ T; x2 y8 Z

indexrow.at(i) = irow; , d8 P7 M8 a. ^6 F/ c' i

5 M* l$ `" l3 i1 j9 \

indexcol.at(i) = icol; 6 e) R6 L$ t0 y+ n) a1 e ' F/ t8 m& H3 x+ T : b+ ]; B! A2 C7 s) ]: Y* X

7 I( J; [ R' w$ [' H! X' S p

/ v: i+ i9 s( b9 d$ E* ]9 b

9 U5 `* o& p6 R; y

try {

% j2 H2 ? `; W2 G4 j \

double pivinv = 1.0 / A(icol, icol); 8 k; P) f, \4 Y8 F* I' \) b , s4 S- L# i2 k6 \. B- `$ j; t0 u

1 k" O: l/ E& w

2 g; b; c4 s L! }; }! V& ^4 j! k

( d4 j3 v8 S3 L; O

for (int l=0; l<n; ++l) 5 |) }& T2 Y' r/ j6 L ; O6 N ?0 D; H

9 U, K2 j/ X7 K& q* R" r

A(icol, l) *= pivinv; 8 {( G I5 ~6 ?9 v" O5 L

for (int l=0; l<m; ++l) A h9 Z% e" L. O8 x5 O ) _ [0 R" T3 T, Q7 L( y

b(icol, l) *= pivinv; ; [- e+ K3 M4 N, v8 b4 ]4 {, g6 ^) L

$ S( b! y# s; l }

0 e2 O+ j- o4 w8 g0 f H8 U

//进行行约化 & t% N' |4 U4 [* {: R2 ^

1 Z. V1 |- @( W/ w

for (int ll=0; ll<n; ++ll) D" n+ l1 e4 m& D

. B9 r& ^6 l7 n* f5 Z. w/ J

if (ll != icol) { & W: t$ n6 i C0 p% x$ M) M

/ K5 W# t& u: `$ ]5 w

, C, O+ a8 ]- R1 S

double dum = A(ll, icol); 4 d% }3 Y1 X( T# }6 S4 f

% O. _) N J& x

9 m' `5 ^# `1 o! q" V5 q+ \

$ ^* }! |0 Z3 g3 m$ @* O5 f

for (int l=0; l<n; ++l) * G2 J, p6 Y/ T9 ~ ! f8 f! M$ k A0 W# t! D- A7 L+ z

. }3 o$ O5 i! B5 T4 O7 b0 n( K$ }

- w; x `" ^. D, I5 u4 m

A(ll, l) -= A(icol, l)*dum; ! f; B" u! z( [/ K9 ?

9 L2 C2 Z' j# z

for (int l=0; l<m; ++l) 6 E. H. j7 t9 x7 E( J& y& p% Q {

5 D" B6 T9 A" h+ u0 t/ w! d$ b9 u

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

8 p( W' u# L2 ]* u( L

} 4 z- k& }: P8 r# i, C( {, Y6 h 2 z3 x; l+ |) s$ I: _, H5 |; @. C # C6 c: @ ~1 b* {6 J* d

5 |8 n, P; S* W. @; h

}

catch (...) { * J; H- T' K* k% W' \8 |& B _ . ~% Q8 P( B S# V; U) A9 B; u0 p/ ]$ K/ G

/ a. S# `' N# O8 K4 I H1 O0 ^5 z8 i

cerr << "Singular Matrix";

} $ ~" |5 w2 U, n( x1 A

} 6 P8 t/ ^! U+ l' g: q# ^ . F, l. |* M" d- _ , r$ E# Z @& g

3 J5 }( J! @2 d8 f2 {7 i

} # f! y6 k2 n( G2 C, ? $ o8 H |- l+ A) G [+ s8 l & Z Z& W$ ]2 M; h1 K

- W0 Q5 l! ]0 _7 G4 Q

. \. |- o& z% g# v

8 \# M4 h3 a3 P4 o& u1 H2 D

/ E, ~3 S7 ?& ~; L0 J: ^

int main() , q) }- a3 |0 _2 p3 P u. P9 \ H8 b - d" z/ r/ j1 [5 z; @' P

5 y( z( f) V# L

) O0 E6 x8 v$ F6 h! Z/ }

{ . P' ~0 u |0 Y }

6 z ]+ H2 M0 i# ]8 b

//测试矩阵 7 y, _% w$ f8 x. d8 m R0 c) \" J$ {2 V8 O3 E' @

! u5 |8 c6 K0 ]3 T8 X9 w

: L+ |1 J$ F: Z) z j

Array<double, 2> A(3,3), b(3,1); 0 Q' c5 j5 w) C1 Y! C3 j2 y& }7 C& Y ( L! j# W7 S D& ]% x

# E+ ~# C2 M, b) m; }7 \

" I' ?# B" m0 s( q h- {; G

A = 10,-19,-2, 3 Z% n- o' R' S& h7 ~: @ . x, U8 l# Q, ^# P$ I8 d " M8 U. U2 f4 j& t/ L0 C4 D

4 M' H# u# o1 Z

-20, 40, 1,

+ h8 z" K2 F" K; c; J( [4 l

1, 4, 5; $ ~. ~7 P* F4 V3 x0 V+ `

+ k$ O$ e; P" L) l

b = 3, 5 ?: `# K; o* S7 p; d* k& z7 K" W

" T, y- L8 V2 }8 G l6 B3 g

2 T. t+ V6 p# P5 g c

4, $ c7 t9 O# ]% ]8 c+ S

5 p) V( }" X4 k+ d4 n

5; # v L1 s0 \! t3 n* _, C$ ^2 c 7 ]) j0 ]8 U6 v9 c

, h) a0 l) S4 S% V; u

& `/ s& k- Y9 b( } - q& v1 G2 ~9 I: v

9 y& ~) N, m3 Z& X6 A k

Gauss_Jordan(A, b); ' t" P! {. j$ y + O3 |4 D+ ]* _- a 6 G: Z2 x# w7 t O+ \' F4 s

1 ^, [9 d8 e, N2 K0 G3 k

1 n( z4 z: k0 s2 s% E, V; D

/ W* q6 _2 ^; D4 d# `

cout << "Solution = " << b <<endl; . r) S" W- T1 C0 q; M) E# M$ i 4 Q# \/ Q4 p' x) y

, z& T4 T+ s0 n0 T

} 0 M. d+ ~+ q# K. Z3 N* B

6 h& v; J: S- c9 L. P! j

; a' i* ~7 e* j% J5 s6 ?5 j# g7 G

Result： + X8 m7 I5 D! j9 a+ n / w# n) V0 |4 _ k+ R6 v' Z& v

& g% `( F% n f( s

% p7 a: I3 t8 Z

Solution = 3 x 1 $ ^+ E' L8 q) e4 e$ X" t" Z) t7 S % V) N# Q3 i; L

9 I! F8 I$ h5 M

[ 4.41637 % L2 \, z* v1 P9 \: z , e& y2 X1 [, c5 W; T

" ~) E) p% P) A4 j$ y' T0 R

2.35231 $ U. K2 e: R* _- J) u- L

-1.76512 ] 1 h" W- I/ F* u5 I 6 e7 M+ A& s/ _4 { ! B! j3 e1 T L- a. h1 x( i

, d& \+ m l, y% N1 a1 K5 y( s

0 k% E: G( c& M: {4 _: S9 r

6 i( Z7 [5 y3 u9 w

9 R( g" n2 Z4 d; a6 [9 o+ S

$ h7 ~! K# D3 W* |) Q4 ]( F4 U

& ?# i, V* D" ~) a9 G- ?2 A

从代码的过程可以看出，矩阵A的逆在A中逐步构造，最终矩阵A演变成单位矩阵，解向量X也逐步替代右端项向量，且使用同一存储空间。 9 ~$ H; ~9 F- g% d3 s7 O6 Z

; T9 z6 n1 E8 u6 p1 ?+ r. q

) P. S/ }! Z; A9 j

3 K1 }9 y9 z7 P7 d

$ ~, G) y, a# _

; r R3 R7 M( v2 c

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

* ^# Q t. r2 h8 A' s: @ A B r, Q) V% Z" l' ~2 P( N

, B/ O7 F- {9 `) }6 |6 D

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

路过。。。