Gauss Strassen LU LU LU LU 22 5 Gauss LU

Transcription

1 I,

2 Gauss Strassen LU LU LU LU 22 5 Gauss LU Gauss LU LU LDU Gauss LU Gauss Gauss LU LU 34

3 55 LU lu-ver4c () (2) Gauss FORTRAN C () C (2) , LU Gauss LU 60 6 Gauss LU LU Cholesky Cholesky

4 20 (20//6) ( ) LINPACK, LAPACK, MATLAB ( ) ( ) ( ) I ( ) ( ) S I ( ) 3

5 2 (6) (Gauss) (Cramer) (Jordan ) 2x + 3x 2 x 3 = 5 (2) 4x + 4x 2 3x 3 = 3 2x + 3x 2 x 3 = , x 2 x 3 x 3 = x + 3x 2 x 3 = 5, 2x 2 x 3 = 7, 5x 3 = 5 4

6 x 3 = 5 5 = 3, x 2 = 7 + x 3 2 = 2, x = 5 3x 2 + x 3 2 = = 0 (forward elimination) (backward substitution) 3 n + ( ) n ( n ) 5

7 3 3 ( ) A system B ( )

8 a, b R n c = a + b n c i = a i + b i (i =, 2,, n) C : for (i = ; i <= n; i++) c[i] = a[i] + b[i]; a R n λ R c = λa n c i = λa i (i =, 2,, n) for (i = ; i <= n; i++) c[i] = lambda * a[i]; a, b R n c = a b n n n c = a i b i c = a[] * b[]; for (i = 2; i <= n; i++) c += a[i] * b[i]; i= 322 A n b R n c = Ab n 2 n(n ) n c i = a ij b j (i =, 2,, n) for (i = ; i <= n; i++) c[i] = a[i][] * b[]; for (j = 2; j <= n; j++) c[i] += a[i][j] * b[j]; : n i= j= n = n(n ), j=2 ( n + i= ) n = n 2 j=2 7

9 A, B n C = AB n 3 n 2 (n ) n c ij = a ik b kj (i, j =, 2,, n) k= for (i = ; i <= n; i++) for (j = ; j <= n; j++) c[i][j] = a[i][] * b[][j]; for (k = 2; k <= n; k++) c[i][j] += a[i][k] * b[k][j]; : n i= n j= k=2 n = n 2 (n ), n i= ( n + j= ) n = n 3 ( ) A, B M(n; R), c R n ABc (AB)c A(Bc) k=2 323 Gauss Gauss ( ) n n 3 /3 + O(n 2 ) n 3 /6 + O(n 2 ) ( ) C Gauss void gauss(int n, matrix a, vector x) int i, j, k; double q, s; /* */ for (k = ; k <= n-; k++) for (i = k + ; i <= n; i++) q = a[i][k] / a[k][k]; for (j = k + ; j <= n; j++) a[i][j] -= q * a[k][j]; x[i] -= q * x[k]; /* */ x[n] /= a[n][n]; for (k = n - ; k >= ; k--) s = x[k]; for (j = k + ; j <= n; j++) s -= a[k][j] * x[j]; x[k] = s / a[k][k]; 8

10 ( n n n ) + = k= i=k+ j=k+ n = = n n k= j=k+ k= i=k+ n(n )(n + ) n(n ) + = 3 2 n n k= i=k+ n(n )(2n + 5) 6 n n n (n k + ) = (n k)(n k + ) = l(l + ) n(n )(2n ) + 6 n(n )(n + ) 3 k= n(n ) 2 n n = (n k) = l = k= n(n ) 6 l= = n(n ) 6 n(n ) 2 (2(n + ) + 3) = n n = (n k) = l = k= l= n + = + (n ) = n + k= n(n ) 2 n(n + ) 2 + n = n(n + ) 2 n(n ) 2 l= (2n + 3) n(n )(2n + 5) 6 = n ((n )(2n + 5) + 3(n + )) 6 = n ( 2n 2 + 3n 5 + 3n + 3 ) 6 = n(n2 + 3n ) Strassen ( ) n 3 VStrassen n O(n log 2 7 ) (log 2 7 = 28 ) n 98 n 25 ( [3] 3 ) ( : n O(n 3 ) ) 9

11 325 LU Gauss O(n 3 ) 326 n (Ly = b, Ux = y) n 2 3 ( ) 6n trid-luc trilu() 2(n ) trisol() 4n 3 trilu() n 2 = n i=0 trisol() n 2 n 2 + = 2(n ) i=0 i=0 2n 327 LU [3] 5 Ax = b x = A b A A b ( ) 0

12 ( b, b 2,, b m ) x j Ax = b j (Ab b 2 b m ) (Ix x 2 x m ) 2 A LU (A = LU) Ax = b Ly = b Ux = y ( ) 33 (CPU CPU ) ( ) ( ) CORDIC (=coordinate rotation digital computation), STL (=sequential table lookup) ( [0] ) (Taylor ) (FLOPS = floating operations per second) GFLOPS (NEC SX6) 3586 TFLOPS ( ) Fourier (fast Fourier transform, FFT)

13 4 LU (LU ) ( ) (i, j) 0 E ij 42 ( ) E ij E kl = δ jk E il 2

14 42 43 ( ) n (elementary matrix) 3 ) ( (interchange matrix)) i j 0 P n (i, j) def = 0 ( i j ) 2) c 0 Q n (i; c) def = i diag( c ) = c ( (i, i) c ) 3) i j R n (i, j; c) def = I + ce ij ( E ij ) ( (i, j) c ) 44 ( ) ( ) P n (i, j) = E ij + E ji + k i,j E kk = I E ii E jj + E ij + E ji (i j), Q n (i; c) = I + (c )E ii, R n (i, j; c) = I + ce ij (i j) 3

15 45 ( ) A (m, n) () A P m (i, j) A i j (2) A P n (i, j) A i j (3) A Q m (i; c) A i c (4) A Q n (i; c) A i c (5) A R m (i, j; c) A i j c (6) A R n (i, j; c) A i j c 46 ( ) 45 (elementary transformation) (, elementary row operation) (, elementary column operation) 47 ( ) P n (i, j) = P n (i, j) (i j) Q n (i; c) = Q n (i; c ) (c 0) R n (i, j; c) = R n (i, j; c) (i j) 48 A B B A 43 [2] 4

16 49 ( ) ( n ) () i 0 j 0 P (i 0, j 0 )Q(k; c) = Q(l; c)p (i 0, j 0 ), k l = j 0 i 0 (k i 0 k j 0 ) (k = i 0 ) (k = j 0 ) (2) i 0 j 0, k l P (i 0, j 0 )R(k, l; c) = R(k, l ; c)p (i 0, j 0 ), (i) k, l i 0, j 0 = k = k, l = l (ii) k = i 0 l j 0 k = j 0, l = l (iii) k = j 0 l i 0 k = i 0, l = l (iv) k i 0 l = j 0 k = k, l = i 0 (v) k j 0 l = i 0 k = k, l = j 0 (vi) k = i 0 l = j 0 k = l, l = k (vii) k = j 0 l = i 0 k = l, l = k (3) i 0 j 0 R(i 0, j 0 ; c)q(k; d) = Q(k; d)r(i 0, j 0 ; c ), c c = cd c/d (k i 0 k j 0 ) (k = j 0 ) (k = i 0 ) 40 ( ) i 0 j 0 P (i 0, j 0 ) = Q(j 0 ; )R(i 0, j 0 ; )R(j 0, i 0 ; )R(i 0, j 0 ; ) ( [7] ) 5

17 4 ( ) i j 4 ) i j 2) j i 3) i j 4) j ( swap ) 42 ( ) P 0 P (permutation matrix) P P = P T P T (n ) m P = (p ij ) i σ(i) (i =, 2,, m) p ij = (j = σ(i)) 0 (j σ(i)) σ m σ P σ m P σ σ P P = P = P T = P ( ) P P = P P 2 P k P = (P k ) (P k ) (P 2 ) (P ) = (P k ) T (P k ) T (P 2 ) T (P ) T = (P P 2 P k ) T = P T (m, n) A = a a 2 a m 6

18 m P σ A = P σ A A σ : a a σ() a 2 P σ = a σ(2) a m a σ(m) n P τ A = AP τ A τ : ( a a 2 a n )P τ = ( a τ() a τ(2) a τ(n) ) : Octave LU MATLAB Octave lu() A LU P A = LU ( A = P T LU ) yurichan% octave GNU Octave, version 206 (i386-unknown-freebsd34) Copyright (C) 996, 997, 998, 999, 2000 John W Eaton This is free software with ABSOLUTELY NO WARRANTY For details, type warranty octave:> n=0;a=rand(n,n);[l u p]=lu(a);norm(p *l*u-a) ans = 34269e-6 octave:2> p p = octave:3> p*p ans = octave:4> b=rand(n);x=u\(l\(p*b));norm(a*x-b) ans = 756e-5 octave:5> A = U L (P T ) = U L P Ax = b U\(L\(P*b)) 7

19 42 42 A (symmetric partitioning) ( ) A A 2 A = O A 22 A A A 22 A = ( A O A A 2 A 22 A 22 [7] 422 () A = (a ij ), B = (b ij ) n A + B (2) A = (a ij ), B = (b ij ) n AB A, B : a a 2 a n b b 2 b n a b 0 a 22 a 2n 0 b 22 b 2n = 0 a 22 b 22 ) 0 0 a nn 0 0 b nn 0 0 a nn b nn (3) n A = (a ij ) A a ii 0 (i =, 2,, n) A A (4) A, B AB 0 AB = (c ij ) c i,i+ = 0 ( i n ) (5) A n A n = O [7] (3) n A = (a ij ) a ii 0 ( i n) B = (b ij ) AB = I A, B (4) i > j = a ij = b ij =

20 AB = I n a ik b kj = δ ij ( i n, j n) (4) k= i k j a ik b kj = δ ij ( i n, j n) i > j i, j 2 i j i, j i = j a jj b jj = (42) b jj = a jj i < j 0 = j a ik b kj = a ii b ij + k=i (43) b ij = a ii j k=i+ j k=i+ a ik b kj a ik b kj (42), (43) a ij, b ij (* b ij *) for j := to n do begin b jj := / a jj ; for i := j - downto do begin s := 0; for k := i + to j do s := s + a ik * b kj ; b ij := - s / a ii ; end; end; (* b ij *) for i := to n do begin b ii := / a ii ; for j := i - downto do begin s := 0; for k := j + to i do s := s + b ik * a kj ; b ij := - s / a jj ; end; end; 2 i k j k 0 δ ij 0 9

21 (C++ ) MATRIX inv_u(const MATRIX &a) int n = RowDimension(a); MATRIX b = zeros(n,n); for (int j = ; j <= n; j++) b(j,j) = / a(j,j); for (int i = j - ; i >= ; i--) REAL s = 0; for (int k = i + ; k <= j; k++) s = s + a(i,k) * b(k,j); b(i,j) = - s / a(i,i); return b; MATRIX inv_l(const MATRIX &a) int n = RowDimension(a); MATRIX b = zeros(n,n); for (int i = ; i <= n; i++) b(i,i) = / a(i,i); for (int j = i - ; j >= ; j--) REAL s = 0; for (int k = j + ; k <= i; k++) s = s + b(i,k) * a(k,j); b(i,j) = - s / a(j,j); return b; 423 ( ) n 43 0 LU 0 ( ) LU ( ) 43 (Dolittle LU ) n A = (a ij ) A, A 2,, A n = A A L U LU [7] ( ) n n = ( ) A n b A = c T d 20

22 ( A n c T b d ) = ( A n 0 c T d c T A n b ) ( E n A n b 0 T ) A n = L n U n (L n, U n ) ( ) ( ) L n 0 U n U n A n b L = c T (U n ) d c T A, U = n b 0 T L, U A = LU ( ) A = LU = L U L L = U U I L = L, U = U L, U A = LU ( ) ( ) ( L n 0 U n p A n b q T r 0 T = c T a nn ) (44) (45) (46) (47) L n U n = A n, L n p = b, q T U n = c T, q T p + r = a nn (44) L n, U n (45), (46), (47) p = (L n ) b, q T = c T (U n ), r = a nn q T p = a nn c T (U n ) (L n ) b = a nn A n b 432 (Crout LU ) n A = (a ij ) A, A 2,, A n = A A L U LU A T (LDU ) n A = (a ij ) A, A 2,, A n = A A A = LDU L D U 43 A = LU u ij u ij /u ii U D = diag(u, u 22,, u nn ) 2

23 434 (Cholesky ) n A = (a ij ) A, A 2,, A n = A A A = LDL T L D A = LDU LDU A A = A T = (LDU) T = U T DL T A = LDU (LDU ) L T = U, U T = L A = LDL T 44 LU A LU ( ) A ( ) LU [2] m n (m, n) A A = LU (m > n ) L = (l ij ) m l ii = (i =, 2,, m), l ij = 0 (i < j) U = (u ij ) (m, n) u ij = 0 (i > j) Gauss 0 ( 0 ) A LU A rank A = m A A = AP (P ) LU : AP = LU A ( ) P A = LU 22

24 5 Gauss LU ( LU ) 5 Ax = b A x = A b Gauss ( ) Gauss Gauss Gauss LU ( Gauss ) Gauss LU 3 ( ) Gauss LU 52 Gauss Cramer ( ) Gauss Gauss (Gauss- Jordan ) Gauss 23

25 2x + 3x 2 x 3 = 5 4x + 4x 2 3x 3 = 3 2x + 3x 2 x 3 = , x x 2 x 3 = Gauss x + 3x 2 x 3 = 5, 2x 2 x 3 = 7, 5x 3 = 5 x 3 = 5 5 = 3, x 2 = 7 + x 3 2 = 2, x = 5 3x 2 + x 3 2 = (forward elimination) (backward substitution) = 53 LU 53 A =

26 L = , L 2 = , L 3 = L A = , L 2 L A = L := L 3 L 2 L = L := L = U := LU = A, L 3 L 2 L A = , 532 n L = (l ij ) (lower triangular matrix) 0 : l ij = 0 ( i < j n) L = l 0 l 2 l 22 l n l n,n l nn L = 0 l 2 l n l n,n n U = (U n i ) (upper triangular matrix) 0 : u ij = 0 (i > j) 25

27 U = u u 2 u n u 22 u n,n 0 u nn u 2 u n U = u n,n 0 ( ) ( ) 533 LU A L U A = LU, l 0 l 2 l 22 L =, U = u u 2 u n u 22 u n,n l n l n,n l nn (LU ) 0 u nn A LU (LU decomposition) L Crout LU : A = LU, l 2 L = l n l n,n 0, U = (Crout LU ) u u 2 u n u 22 u n,n 0 u nn U Dolittle LU : A = LU, l 0 l 2 l 22 L =, U = l n l n,n l nn (Dolittle LU ) u 2 u n 0 u n,n 26

28 534 LDU A A = LDU, l 2 L = l n l n,n 0, D = d 0 0 d n, U = u 2 u n 0 u n,n L, D, U LDU 53 LDU A = LDU A = L(DU) Crout LU A = (LD)U Dolittle LU Crout LU Dolittle LU LDU Crout LU A = LU, l 2 L = l N l N,N 0, U = u u 2 u N u 22 u N,N 0 u NN D def = ũ kj def = u kj u kk u 0 u 22 0 u NN (k =, 2,, N; j = k +, k + 2,, N),, Ũ = A = LDŨ ũ 2 ũ N 0 ũ N,N 535 A A = LU LU Ax = b ( A = (LU) = U L L, U ) LUx = b Ly = b, Ux = y 27

29 Ly = b l y = b l 2 y +l 22 y 2 = b 2 l 3 y +l 32 y 2 +l 33 y 3 = b 3 l n y +l n2 y 2 +l n3 y 3 +l nn y n = b n y = b /l, y 2 = (b 2 l 2 y ) /l 22, y 3 = (b 3 l 3 y l 32 y 2 ) /l 33, y i = y n = ( ) i b i l ij y j /l ii, j= ( ) n b n l nj y j /l nn n = n(n + )/2 Ux = y j= u x +u,n 2 x n 2 +u,n x n +u n x n = y, u n 2,n 2 x n 2 +u n 2,n x n +u n 2,n x n = y n 2, u n,n x n +u n,n x n = y n, u n,n x n = y n x n = y n /u nn, x n = (y n u n,n x n ) /u n,n, x n 2 = (y n 2 u n 2,n x n u n 2,n x n ) /u n 2,n 2, x i = x = ( ( y i y n j=i+ u ij x j ) /u ii, ) n u j x j /u j= n = n(n + )/2 n(n + )/2 + n(n + )/2 = n(n + ) Gauss n n 3 /3 LU 28

30 A A b n 2 LU (LU ) A LU LU LU A LU ( ) A ( P P A ) LU ( ) Gauss LU ( ) n 3 /3 + O(n 2 ) A n 3 ( ) 532 ( ) A ( ) L, U LU 533 ( ) LU L U n(n + )/2 + n(n )/2 = n 2 A A A LU 54 Gauss LU Gauss LU ( Gauss LU LU Gauss ) Gauss Ax = b (b R n, A M(n; R)) A, b (A b) 2 29

31 54 k =, 2,, n k i (i = k +,, n) (A b) = (A () b () ), (A () b () ) (A (2) b (2) ) (A (k) b (k) ) (A (k+) b (k+) ) (A (n) b (n) ), u u 2 u n A (n) u 22 u 2n def = = U, b (n) = y 0 u nn A (k) = (a (k) ij ), b(k) = (b (k) i ) () a () ij def = a ij, b () def i = b i (i, j =, 2,, n) (2) A (k) = (a (k) ij ), b(k) = (b (k) i ) k =, 2,, n a (k) ij ( i k; j n) (5) (52) a (k+) ij b (k+) i def = def = (3) U = (u ij ), y = (y i ) a (k) ij a(k) ik a (k) kk a (k) kj b (k) i ( i k) b (k) i a(k) ik a (k) kk b (k) k u ij = a (n) ij, (k + i n; j n), (k + i n), y i = b (n) i A (k) = (a (k) ij ) (k =, 2,, n) C for (k = ; k < n; k++) for (i = k + ; i <= n; i++) q = a[i][k] / a[k][k]; for (j = ; j <= n; j++) a[i][j] = a[i][j] - q * a[k][j]; b[i] = b[i] - q * b[k]; 30

32 542 (53) (54) x n = b n, a nn x k = a kk ( b k n j=k+ a kjx j ) C x[n] = b[n] / a[n][n]; for (k = n - ; k >= ; k--) s = b[k]; for (j = k + ; j <= n; j++) s = s - a[k][j] * x[j]; x[k] = s / a[k][k]; (k = n, n 2,, 2, ) 543 Gauss Gauss /* * gauss-verc */ #include <matrixh> void gauss(int n, matrix a, vector x) int i, j, k; double q; /* */ for (k = ; k < n; k++) for (i = k + ; i <= n; i++) q = a[i][k] / a[k][k]; #ifdef ORIGINAL for (j = ; j <= n; j++) #else for (j = k + ; j <= n; j++) #endif a[i][j] -= q * a[k][j]; x[i] -= q * x[k]; /* */ x[n] /= a[n][n]; for (k = n - ; k >= ; k--) for (j = k + ; j <= n; j++) x[k] -= a[k][j] * x[j]; x[k] /= a[k][k]; 3

33 0 j k+ 544 Gauss LU (A (k) b (k) ) (A (k+) b (k+) ) k i (k + i n) L (k) i A (k+) = L (k) n L (k) n L (k) k+ A(k), b (k+) = L (k) n L (k) n L (k) k+ b(k), k i L (k) i = = I q (k) i E ik, q (k) i q (k) i E ik = (i, k) 0 = a(k) ik a (k) kk, L (k) i L (k) def = L (k) n L (k) n L (k) L (k) k : (55) L = k+ A (k+) = L (k) A (k), b (k+) = L (k) b (k) L def = L (n ) L (n 2) L (2) L () A (n) = LA (), b (n) = Lb () U = LA, q () 2 q () 3 q (2) 3 L def = L A = LU y = Lb q n () q n (2) q n (n ) 32

34 L A LU (55) L = L (n ) L (n 2) L (2) L () L (k) = L (k) n L (k) n L (k) k+ L = L = (L () ) (L (2) ) (L (n 2) ) (L (n ) ) (L (k) ) = (L (k) k+ ) (L (k) (L (k) i ) = n ) (L (k) n ) q (k) i (L (k) ) = (L (k) k+ ) (L (k) n ) (L (k) n ) = U L = q () 2 q () 3 q (2) 3 q n () q n (2) q n (n ) q (k), k+ q n (k) q (k) i def = a(k) ik a (k) kk (k =, 2,, n ; i = k +,, n) b Ax = b Gauss A b LU 33

35 545 LU lu-verc /* * lu-verc */ #include <matrixh> #include "lu-verh" void lu(int n, matrix a, matrix L, matrix u) int i, j, k; double q; for (i = ; i <= n; i++) for (j = ; j <= n; j++) u[i][j] = a[i][j]; L[i][j] = 00; L[i][i] = 0; for (k = ; k < n; k++) for (i = k + ; i <= n; i++) q = u[i][k] / u[k][k]; /* j j k OK (k+ ) */ #ifdef ORIGINAL for (j = ; j <= n; j++) #else for (j = k; j <= n; j++) #endif u[i][j] = u[i][j] - q * u[k][j]; L[i][k] = q; /* A LU A x = b */ void solve(int n, matrix L, matrix U, vector b, vector x) /* 434 */ int i, j; double sum; vector y = new_vector(n + ); /* */ for (i = ; i <= n; i++) sum = b[i]; for (j = ; j < i; j++) sum -= L[i][j] * y[j]; y[i] = sum / L[i][i]; /* */ for (i = n; i >= ; i--) sum = y[i]; for (j = i + ; j <= n; j++) sum -= U[i][j] * x[j]; x[i] = sum / U[i][i]; /* */ free_vector(y); A = LU LU 2 L, U 34

36 L (, 0) U ( 0 ) n 2 ( L U ) u u 2 u 3 u n q () 2 u 22 u 23 u 2n q () 3 q (2) 3 u 33 u 3n q n () q n (2) q n (n ) matrix a L, U a (j k+ ) lu-ver2c /* * lu-ver2c */ #include <matrixh> #include "lu-ver2h" void lu(int n, matrix a) int i, j, k; double q; u nn for (k = ; k < n; k++) for (i = k + ; i <= n; i++) q = a[i][k] / a[k][k]; for (j = k + ; j <= n; j++) a[i][j] = a[i][j] - q * a[k][j]; a[i][k] = q; /* A LU A x = b */ void solve(int n, matrix a, vector b) /* 434 */ int i, j; /* */ for (i = ; i <= n; i++) for (j = ; j < i; j++) b[i] -= a[i][j] * b[j]; /* : L_i,i= */ /* */ for (i = n; i >= ; i--) for (j = i + ; j <= n; j++) b[i] -= a[i][j] * b[j]; b[i] /= a[i][i]; 35

37 lu-ver2c (lu-verc) utilc /* * utilc */ #include <stdioh> #include <matrixh> #include "utilh" /* n A */ void print_matrix(int n, matrix A) int i, j; for (i = ; i <= n; i++) for (j = ; j <= n; j++) printf("%0g ", A[i][j]); printf("\n"); /* n x */ void print_vector(int n, vector x) int i; for (i = ; i <= n; i++) printf("%0g ", x[i]); printf("\n"); /* n A, B AB */ void mul_matrix(int n, matrix AB, matrix A, matrix B) /* */ int i, j, k; double sum; for (i = ; i <= n; i++) for (j = ; j <= n; j++) sum = 00; for (k = ; k <= n; k++) sum += A[i][k] * B[k][j]; AB[i][j] = sum; 36

38 test-lu-verc /* * test-lu-verc --- LU * gcc -I/usr/local/include -c utilc * gcc -I/usr/local/include -c lu-verc * ccmg test-lu-verc utilo lu-vero */ #include <stdioh> #include <matrixh> #include "utilh" #include "lu-verh" int main() int n; matrix a, L, u, Lu; vector b, x; n = 3; /* n a, L, u, Lu */ a = new_matrix(n +, n + ); L = new_matrix(n +, n + ); u = new_matrix(n +, n + ); Lu = new_matrix(n +, n + ); /* n b, x */ b = new_vector(n + ); x = new_vector(n + ); /* a ( p27 ) */ a[][] = 2; a[][2] = 3; a[][3] = -; a[2][] = 4; a[2][2] = 4; a[2][3] = -3; a[3][] = -2; a[3][2] = 3; a[3][3] = -; /* LU */ lu(n, a, L, u); /* a, L, u */ printf("a=\n"); print_matrix(n, a); printf("l=\n"); print_matrix(n, L); printf("u=\n"); print_matrix(n, u); /* L u */ mul_matrix(n, Lu, L, u); printf("l U=\n"); print_matrix(n, Lu); /* (p27) */ b[] = 5; b[2] = 3; b[3] = ; solve(n, L, u, b, x); printf("x=\n"); print_vector(n, x); return 0; 37

39 oyabun% gcc -I/usr/local/include -c utilc oyabun% gcc -I/usr/local/include -c lu-verc oyabun% ccmg test-lu-verc utilo lu-vero oyabun% /test-lu-ver A= L= U= L U= x= 2 3 oyabun% LU = A 38

40 (lu-ver2c) test-lu-ver2c /* * test-lu-ver2c --- LU * gcc -I/usr/local/include -c utilc * gcc -I/usr/local/include -c lu-ver2c * ccmg test-lu-ver2c utilo lu-ver2o */ #include <stdioh> #include <matrixh> #include "utilh" #include "lu-ver2h" int main() int n; matrix a; vector b; n = 3; /* n a */ a = new_matrix(n +, n + ); /* n b, x */ b = new_vector(n + ); /* a ( p27 ) */ a[][] = 2; a[][2] = 3; a[][3] = -; a[2][] = 4; a[2][2] = 4; a[2][3] = -3; a[3][] = -2; a[3][2] = 3; a[3][3] = -; /* a */ printf("a=\n"); print_matrix(n, a); /* b */ b[] = 5; b[2] = 3; b[3] = ; printf("b=\n"); print_vector(n, b); /* LU */ lu(n, a); /* a */ printf("lu A=\n"); print_matrix(n, a); /* (p27) */ solve(n, a, b); printf("x=\n"); print_vector(n, b); return 0; 39

41 oyabun% gcc -I/usr/local/include -c utilc oyabun% gcc -I/usr/local/include -c lu-ver2c oyabun% ccmg test-lu-ver2c utilo lu-ver2o oyabun% /test-lui--ver2 A= b= 5 3 LU A= x= 2 3 oyabun% 55 LU 55 lu-ver4c /* * lu-ver4c --- Gauss LU * * A Gauss * * P A = L U, (P, L, U ) * * P, L, U * */ #include <mathh> #include <matrixh> #include "lu-ver4h" void ipvt2perm(int n, ivector perm, ivector ipvt) int i, m, perm_i; for (i = ; i <= n; i++) perm[i] = i; for (i = ; i < n; i++) m = ipvt[i]; if (i!= m) perm_i = perm[i]; perm[i] = perm[m]; perm[m] = perm_i; /* Gauss LU */ void lu(int n, matrix a, ivector ipvt) int i, j, k, m; double q, t, maxpvt; /* */ ipvt[n] = ; /* */ 40

42 for (k = ; k < n; k++) /* a[i][k] (i=k,k+,,n) */ maxpvt = fabs(a[k][k]); m = k; for (i = k + ; i <= n; i++) if (fabs(a[i][k]) > maxpvt) maxpvt = fabs(a[i][k]); m = i; ipvt[k] = m; /* m */ if (k!= m) /* k m */ for (j = ; j <= n; j++) t = a[k][j]; a[k][j] = a[m][j]; a[m][j] = t; /* */ ipvt[n] = - ipvt[n]; /* */ for (i = k + ; i <= n; i++) q = a[i][k] / a[k][k]; for (j = k + ; j <= n; j++) a[i][j] = a[i][j] - q * a[k][j]; a[i][k] = q; /* A LU (lu() ) A x = b */ void solve(int n, matrix a, vector b, ivector ipvt) int i, j; for (i = ; i < n; i++) int m; double t; m = ipvt[i]; if (i!= m) t = b[m]; b[m] = b[i]; b[i] = t; /* */ for (i = ; i <= n; i++) for (j = ; j < i; j++) b[i] -= a[i][j] * b[j]; /* : L_i,i= */ /* */ for (i = n; i >= ; i--) for (j = i + ; j <= n; j++) b[i] -= a[i][j] * b[j]; b[i] /= a[i][i]; /* ipvt[] P A */ void permute_matrix(int n, matrix pa, matrix a, ivector ipvt) ivector perm = new_ivector(n + ); ipvt2perm(n, perm, ipvt); permute_matrix0(n, pa, a, perm); free_ivector(perm); /* perm[] P A */ 4

43 void permute_matrix0(int n, matrix pa, matrix a, ivector perm) int i, j, m; for (i = ; i <= n; i++) m = perm[i]; for (j = ; j <= n; j++) pa[i][j] = a[m][j]; 552 () test-lu-ver4c /* * test-lu-ver4c --- LU * gcc -I/usr/local/include -c utilc * gcc -I/usr/local/include -c lu-ver4c * ccmg test-lu-ver4c utilo lu-ver4o */ #include <stdioh> #include <matrixh> #include "utilh" #include "lu-ver4h" int main() int n; matrix a; vector b; ivector ipvt; n = 3; /* n a */ a = new_matrix(n +, n + ); /* n b, x */ b = new_vector(n + ); /* */ ipvt = new_ivector(n + ); /* a ( p27 ) */ a[][] = 2; a[][2] = 3; a[][3] = -; a[2][] = 4; a[2][2] = 4; a[2][3] = -3; a[3][] = -2; a[3][2] = 3; a[3][3] = -; /* a */ printf("a=\n"); print_matrix(n, a); /* b */ b[] = 5; b[2] = 3; b[3] = ; printf("b=\n"); print_vector(n, b); /* LU */ lu(n, a, ipvt); /* a */ printf("lu A=\n"); print_matrix(n, a); /* (p27) */ solve(n, a, b, ipvt); printf("x=\n"); print_vector(n, b); return 0; 42

44 test-lu-ver4 oyabun% /test-lu-ver4 A= b= 5 3 LU A= x= 2 3 oyabun% 553 (2) test-lu4c /* * test-lu4c --- LU * * * : gcc -o test-lu4 test-lu4o lu-ver4o matrix-baseo * -L/usr/local/lib -lmatrix -lm */ #include <stdioh> /* standard input & output */ #include <mathh> /* */ #include <matrixh> /* matrix */ #include "lu-ver4h" #include "matrix-baseh" int main() /* */ int n; matrix a, LandU, Lu, diff, pa; ivector perm, ipvt; double matrix_norm_2(); printf("n Hilbert LU \n"); printf("n="); scanf("%d", &n); /* ( ) */ a = new_matrix(n +, n + ); LandU = new_matrix(n +, n + ); Lu = new_matrix(n +, n + ); diff = new_matrix(n +, n + ); pa = new_matrix(n +, n + ); perm = new_ivector(n + ); ipvt = new_ivector(n + ); /* */ #ifdef OLD int i, j; for (i = ; i <= n; i++) 43

45 for (j = ; j <= n; j++) printf("a[%d][%d]=", i, j); scanf("%lf", &a[i][j]); #else hilbert(n, a); #endif printf("a =\n"); display_matrix(n, a, "%f "); /* LU */ copy_matrix(n, LandU, a); lu(n, LandU, ipvt); printf("l, U \n"); display_matrix(n, LandU, "%f "); /* L U */ multiply_lumatrix(n, Lu, LandU); printf("l U=\n"); display_matrix(n, Lu, "%f "); /* P A */ permute_matrix(n, pa, a, ipvt); printf("p A=\n"); display_matrix(n, pa, "%f "); /* P A L U */ subtract_matrix(n, diff, pa, Lu); printf("p A - L U=\n"); display_matrix(n, diff, "%g "); printf(" P A - L U =%g\n", matrix_norm_2(n, diff)); return 0; 44

46 test-lu4 oyabun% /test-lu4 n Hilbert LU n=3 A = L, U L U= P A= P A - L U= P A - L U =0 oyabun% 56 Gauss Cramer 2 (Jordan ) 3 Gauss 3 Cramer n + ( ) Cramer n ( n ) Cramer Cramer 2 Cramer ( ) Gauss Gauss 2 ± 45

47 θ 57 3 Gauss LU ( Thomas 3 ) LU LU mathpc00% make test-lu-trid gcc -O -pipe -c test-lu-tridc gcc -O -pipe -c lu-ver2c gcc -O -pipe -c utilc gcc -o test-lu-trid test-lu-trido lu-ver2o utilo -L/usr/local/lib -lmatrix -lm mathpc00% /test-lu-trid n=7 A= LU A= mathpc00% ( [] Fortran ) ( pivoting) ( 0 Smith [8] (Wilkinson ) ( ) ) 3 Smith [8] 46

48 573 FORTRAN 77 [] c c subroutine trid(a,b,c,f,n,id) integer n,id real a(*),b(*),c(*),f(*) integer i (forward elimination) if (id eq 0) then do i=,n- b(i+)=b(i+)-c(i)*a(i+)/b(i) end do endif do i=,n- f(i+)=f(i+)-f(i)*a(i+)/b(i) end do (backward substitution) f(n)=f(n)/b(n) do i=n-,,- f(i)=(f(i)-c(i)*f(i+))/b(i) end do end 574 C () 0 /* trid-luc -- 3 Gauss */ #include "trid-luh" /* 3 ( ) Ax=b * * * n: * al,ad,au: * (al: ie (lower part) * ad: ie (diagonal part) * au: ie (upper part) * * * ad[0] au[0] 0 0 * al[] ad[] au[] 0 0 * 0 al[2] ad[2] au[2] 0 0 * * al[n-2] ad[n-2] au[n-2] * 0 al[n-] ad[n-] * al[i] = A_i,i-, ad[i] = A_i,i, au[i] = A_i,i+, * al[0], au[n-] ) * * b: * ( 0 ie b[0],b[],,b[n-] ) * * al,ad,au: LU * b: * * call LU * * trisol() * 47

49 * Gauss * * */ void trid(int n, double *al, double *ad, double *au, double *b) trilu(n,al,ad,au); trisol(n,al,ad,au,b); /* LU (pivoting ) */ void trilu(int n, double *al, double *ad, double *au) int i, nm = n - ; /* (forward elimination) */ for (i = 0; i < nm; i++) al[i + ] /= ad[i]; ad[i + ] -= au[i] * al[i + ]; /* LU 3 */ void trisol(int n, double *al, double *ad, double *au, double *b) int i, nm = n - ; /* (forward elimination) */ for (i = 0; i < nm; i++) b[i + ] -= b[i] * al[i + ]; /* (backward substitution) */ b[nm] /= ad[nm]; for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + ]) / ad[i]; C trisol() trilu() Fortran C ( LU C ) 575 C (2) /* trid-luc -- 3 Gauss */ #include "trid-luh" /* 3 ( ) Ax=b * * * n: * al,ad,au: * (al: ie (lower part) * ad: ie (diagonal part) * au: ie (upper part) * * * ad[] au[] 0 0 * al[2] ad[2] au[2] 0 0 * 0 al[3] ad[3] au[3] 0 0 * 48

50 * al[n-] ad[n-] au[n-] * 0 al[n] ad[n] * al[i] = A_i,i-, ad[i] = A_i,i, au[i] = A_i,i+, * al[], au[n] ) * * b: * ( ie b[],b[2],,b[n] ) * * al,ad,au: LU * b: * * call LU * * trisol() * * Gauss * * */ void trid(int n, double *al, double *ad, double *au, double *b) trilu(n,al,ad,au); trisol(n,al,ad,au,b); /* LU (pivoting ) */ void trilu(int n, double *al, double *ad, double *au) int i; /* (forward elimination) */ for (i = ; i < n; i++) al[i + ] /= ad[i]; ad[i + ] -= au[i] * al[i + ]; /* LU 3 */ void trisol(int n, double *al, double *ad, double *au, double *b) int i; /* (forward elimination) */ for (i = ; i < n; i++) b[i + ] -= b[i] * al[i + ]; /* (backward substitution) */ b[n] /= ad[n]; for (i = n - ; i >= ; i--) b[i] = (b[i] - au[i] * b[i + ]) / ad[i]; trisol(), trilu(), trid() trisol(), trilu(), trid() ( ) 49

51 void trisol(int n, double *al, double *ad, double *au, double *b) trisol(n, al+, ad+, au+, b+); void trilu(int n, double *al, double *ad, double *au) trilu(n, al+, ad+, au+); void trid(int n, double *al, double *ad, double *au, double *b) trilu(n, al+, ad+, au+, b+); 58 58, N A = (a i,j ) m i j > m = a i,j = 0 a i,j 0 = i j m ( m ) 58 A 3 A 4 ( ) = max i j a ij 0 ( ) ( ) ( ) 0, m N A C 0 A = (a ij ) 0 i,j N m N N 2m + Ã = (ã ik ) 0 i N 0 k 2m 4 50

52 a ij ã ik ã ik a ij 0 i, j N i, j ã i,j i+m ( i j m) a ij = 0 ( i j > m) 0 i N, 0 k 2m i, k a i,i+k m (0 i + k m N ) ã ik = ( ) 583 C 0 A = (a ij ) 0 i,j N m N N m + Ã = (ã ik) 0 i N 0 k m a ij ã ik ã ik a ik 0 i, j N i, j 0 ( i j > m) a ij = ã i,j i (i j i + m) ã j,i j (i m j < i) 0 i N, 0 k m i, k a i,i+k (i + k N ) ã ik = ( ) 584 LU A LU A = LU U U ( (56) ) A det A( : k, : k) 0 (k =, 2,, n) ( A ) A A = LU, L = 0, U = 0 U L U U U = DU, D = 0 0 5, U = 0

53 A = LDU LDU LDU A A = A T = (LDU ) T = (U ) T DL T A LDU U = L T U = DU = DL T l 2 l 3 l n u u 2 u n d u 22 u 2n 0 = d 2 l 32 l n2 l 0 u nn 0 n,n d n 0 U L u ii = d i (i =,, n), u ij = d i lji ( i < j n) (56) l ji = u ij /d i = u ij /u ii ( i < j n) LU /* LU */ void lu(matrix a, int n) int i, j, k; double q; for (i = ; i < n; i++) for (j = i + ; j <= n; j++) q = a[j][i] / a[i][i]; for (k = i + ; k <= n; k++) a[j][k] -= q * a[i][k]; a[j][i] = q; void solve(matrix a, vector b, int n) int i, j; for (i = ; i < n; i++) for (j = i + ; j <= n; j++) b[j] -= a[j][i] * b[i]; b[n] /= a[n][n]; for (i = n - ; i >= ; i--) for (j = i + ; j <= n; j++) b[i] -= a[i][j] * b[j]; b[i] /= a[i][i]; 52

54 L ( ) /* LU ( ) */ void lu2(matrix a, int n) int i, j, k; double q; for (i = ; i < n; i++) for (j = i + ; j <= n; j++) q = a[j][i] / a[i][i]; for (k = i + ; k <= n; k++) a[j][k] -= q * a[i][k]; /* lu() L[j][i] a[j][i] */ /* a[i][j] (i<j) * for (i=; i<=n; i++) for (j=; j<i; j++) a[i][j] = a[j][i] / a[j][j]; */ void solve2(matrix a, vector b, int n) int i, j; for (i = ; i < n; i++) for (j = i + ; j <= n; j++) /* L[j][i] a[i][j] / a[i][i] */ b[j] -= (a[i][j] / a[i][i]) * b[i]; b[n] /= a[n][n]; for (i = n - ; i >= ; i--) for (j = i + ; j <= n; j++) b[i] -= a[i][j] * b[j]; b[i] /= a[i][i]; lu2() a[][] (a[i][j] (i > j) ) void lu3(matrix a, int n) int i, j, k; double q; for (i = ; i < n; i++) for (j = i + ; j <= n; j++) /* a[j][i] a[i][j] */ q = a[i][j] / a[i][i]; /* k i + k j */ for (k = j; k <= n; k++) a[j][k] -= q * a[i][k]; /* solve2() */ 53

55 a[i][j] (i > j) 585 lubandc ( ) /* Gauss * * A m * Ax=b + * * (a_i,j) am, * * am[i][j] = a_i,i+j-m (m i+j < m+n) * = 0 (i+j < m, m+n i+j) * * a_i,j = am[i][j-i+m] (0 j-i+m < 2m+) */ #include <matrixh> /* LU */ void bandlu(int m, int n, matrix am) int i, j, k, m, mm; double aa; m = m - ; for (i = 0; i < m; i++) mm = i + n + ; if (mm > m) mm = m; for (j = i + ; j < mm; j++) /* aa = a[j][i]/a[i][i]; */ aa = am[j][i - j + n] / am[i][n]; for (k = i; k < mm; k++) /* a[j][k] -= aa * a[i][k]; */ am[j][k - j + n] -= aa * am[i][k - i + n]; am[j][i - j + n] = aa; /* LU */ void bandsol(int m, int n, matrix am, vector f) int i, j, jj, m, mm; m = m - ; for (i = ; i < m; i++) jj = i - n; if (jj < 0) jj = 0; for (j = jj; j < i; j++) f[i] -= f[j] * am[i][j - i + n]; f[m] /= am[m][n]; for (i = m - ; i >= 0; i--) mm = i + n + ; 54

56 if (mm > m) mm = m; for (j = i + ; j < mm; j++) f[i] -= f[j] * am[i][j - i + n]; f[i] /= am[i][n]; symbandluc [4] LU ( Fortran ) C /* * symbandluc --- */ #include <matrixh> /* the Gauss elimination method for symmetric band matrices */ /* N m + */ void symbandlu(matrix at, int N, int m) int N, i, j, k, mm; double aa; /* forward elimination; */ N = N - ; for (i = 0; i < N; i++) /* mm = min(i+m,n) */ mm = i + m; if (mm > N) mm = N; for (j = i + ; j < mm; j++) aa = at[i][j-i] / at[i][0]; /* aa = a_i,j / a_i,i */ for (k = j; k < mm; k++) at[j][k-j] -= aa * at[i][k-i]; /* a_j,k -= aa * a_i,k */ void symbandsolve(matrix at, vector f, int N, int m) int N, i, j, k, mm; double aa; /* forward elimination; */ N = N - ; for (i = 0; i < N; i++) mm = i + m; if (mm > N) mm = N; for (j = i + ; j < mm; j++) aa = at[i][j-i] / at[i][0]; /* aa = a_i,j / a_i,i */ f[j] -= aa * f[i]; /* backward substitution */ f[n-] /= at[n-][0]; /* f_n- /= a_n-,n- */ for (i = N - 2; i >= 0; i--) mm = i + m; if (mm > N) mm = N; for (j = i + ; j < mm; j++) f[i] -= at[i][j-i] * f[j]; /* f_i -= at_i,j * f_j */ f[i] /= at[i][0]; /* f_i /= a_i,i */ include symbandluh 55

57 symbandluh /* * symbandluh */ #include <matrixh> void band(matrix at, vector f, int N, int m); void symbandlu(matrix at, int N, int m); void symbandsolve(matrix at, vector f, int N, int m); LAPACK ( ) ( C LAPACK bandluc, symbandluc LAPACK++ ) symbandluc /* * bandtest5c */ #include <stdioh> #include <stdlibh> #include <mathh> #include "matrixh" void lu4(matrix, int, int); void solve4(matrix, vector, int, int); void testmatrix(int, int, matrix, vector, vector); void disp(int, int, matrix, vector, char *, char *); #define FORMAT " %64f " void testmatrix(int n, int m, matrix a, vector x, vector b) int i, j, mm; /* A */ /* A */ for (i = 0; i <= n; i++) for (j = 0; j <= n; j++) a[i][j] = 0; for (i = ; i <= n; i++) mm = i + m; if (mm > n) mm = n; for (j = i + ; j <= mm; j++) a[i][j-i] = rand() / (double) RAND_MAX; a[i][0] = i + ; /* x=(,2,,n)^t b:=a x */ for (i = ; i <= n; i++) x[i] = i; for (i = ; i <= n; i++) b[i] = 0; mm = i - m; if (mm < ) mm = ; for (j = mm; j < i; j++) b[i] += a[j][i-j] * x[j]; 56

58 mm = i + m; if (mm > n) mm = n; for (j = i; j <= mm; j++) b[i] += a[i][j-i] * x[j]; void disp(int n, int m, matrix a, vector x, char *s, char *s2) int i, j; printf(s); for (i = ; i <= n; i++) for (j = 0; j <= m; j++) printf(format, a[i][j]); printf("\n"); printf(s2); for (i = ; i <= n; i++) printf(format, x[i]); printf("\n"); int main() int i, j, k, n, mm, m; matrix a, a0; vector x, b, b0; n = 8; m = 3; a = new_matrix(n +, m + ); b = new_vector(n + ); x = new_vector(n + ); testmatrix(n, m, a, x, b); disp(n, m, a, b, "A=\n", "b=\n"); /* LU */ lu4(a, n, m); solve4(a, b, n, m); disp(n, m, a, b, "LU=\n", "x=\n"); void lu4(matrix a, int n, int m) int i, j, k, mm; double q; for (i = ; i < n; i++) mm = i + m; if (mm > n) mm = n; for (j = i + ; j <= mm; j++) /* q = a_ji / a_ii a_ji a_ij */ q = a[i][j-i] / a[i][0]; /* for (k = i + ; k <= mm; k++) */ for (k = j; k <= mm; k++) a[j][k-j] -= q * a[i][k-i]; void solve4(matrix a, vector b, int n, int m) int i, j, mm; for (i = ; i < n; i++) 57

59 mm = i + m; if (mm > n) mm = n; for (j = i + ; j <= mm; j++) /* b[j] -= L[j][i] * b[i]; */ b[j] -= (a[i][j-i] / a[i][0]) * b[i]; b[n] /= a[n][0]; for (i = n - ; i >= ; i--) mm = i + m; if (mm > n) mm = n; for (j = i + ; j <= mm; j++) b[i] -= a[i][j-i] * b[j]; b[i] /= a[i][0]; oyabun% ccmg bandtest5c oyabun% /bandtest5 A= b= LU= x= oyabun% 586 Poisson Dirichlet N n := N 2 n m N ( ) Gauss n 3 = N 6 5 N = 00 = 0 2 n = N 2 = 0 4, n 2 = N 4 = (0 2 ) 4 = 0 8 = 00M C double 8B 800MB N 6 = (0 2 ) 6 = 0 2 = T = 000G GFLOPS 000 0MFLOPS 5 n 3 /3 n 3 /3 58

60 N = = 64 nm 2 nm = N 3 nm 2 = N 4 8 nmb = 8N 3 B = 8(0 2 ) 3 B = B = 8MB nm 2 = N 4 = (0 2 ) 4 = 0 8 = 00M 0MFLOPS 0 N = = 6 59 Gauss ( ) Gauss LAPACK 59

61 6 Gauss LU ( ) 6 Gauss 6 ( a (k) kk ) 0 A = (a ij ) M(n; R), b R n (6) Ax = b 62 (6) a () x + a () 2 x a () n x n = b () a () 2 x + a () 22 x a () 2n x n = b () 2 a () 3 x + a () 32 x a () 3n x n = b () 3 a () n x + a () n2 x a () nnx n = b () n () 2,, n a () x + a () 2 x a () n x n = b () a (2) 22 x a (2) 2n x n = b (2) 2 a (2) 32 x a (2) 3n x n = b (2) 3 a (2) n2 x a (2) nnx n = b (2) n 60

62 q 2 = a () 2 /a (), q 3 = a () 3 /a (),, q n = a () n /a (), a (2) 22 = a () 22 q 2 a () 2, a (2) 23 = a () 23 q 2 a () 3,, a (2) 2n = a () 2n q 2 a () n, b (2) 2 = b () 2, a (2) 32 = a () 32 q 3 a () 2, a (2) 33 = a () 33 q 3 a () 3,, a (2) 3n = a () 3n q 3 a () n, b (2) 3 = b () 3, a (2) i2 = a() i2 q ia () 2, a (2) i3 = a() i3 q ia () 3,, a (2) in = a() in q ia () n, b (2) i = b () i, a (2) n2 = a () n2 q n a () 2, a (2) n3 = a () n3 q n a () 3,, a (2) nn = a () nn q n a () n, b (2) n = b () n k a () x + a () 2 x a () k + + a () n x n = b () a (2) 22 x a (2) 2k + + a (2) 2n x n = b (2) 2 a (k) kk a k + + a (k) 3n x n = b (k) 3 k k +,, n a (k) nk x k + + a (k) nnx n = b (k) n a () x + a () 2 x a () k x k + a (),k+ x k+ + + a () n x n = b () a (2) 22 x a (2) 2k x k + a (2) 2,k+ x k+ + + a (2) 2n x n = b (2) 2 a (k) kk x k + a (k) k,k+ x k+ + + a (k) 3n x n = b (k) k a (k+) k+,k+ x k+ + + a (k+) 3n x n = b (k+) k+ a (k+) k+,k+ x k+ + + a (k+) 3n x n = b (k+) k+ q ik = a (k) ik /a(k) kk (i = k +, k + 2,, n), a (k+) ij = a (k) ij q ik a (k) kj (i = k +, k + 2,, n; j = k +, k + 2,, n), b (k+) i = b (k) i q ik b (k) k (i = k +, k + 2,, n) n a () x + a () 2 x 2 + a (),n x n + a (),nx n = b (), a (2) 22 x 2 + a (2) 2,n x n + a (2) 2,nx n = b (2) 2, (62) a (n ) n,n x n + a n,nx (n ) n = b (n ) n, a n,nx (n) n = b n (n) 6

63 (62) : x n = b n (n) ( x n = /a nn (n), b (n ) n a (n ) n,nx n ) /a (n ) n,n, x k = x = ( ( b (k) k b () n i=k+ i=2 a (k) ki x i n a () i x i ) ) /a (k) kk, /a () 63 Ax = b A b A () = (A b) = a () a () 2 a () n b () a () 2 a () 22 a () 2n b () 2 a () n a () n2 a () nn b () n a () ij = a ij (i, j =, 2,, n), b () = b i (i =, 2,, n) k k +,, n k A (k) = a () a () 2 a () k a (),k+ a (),n b () a (2) 22 a (2) 2k a (2) 2,k+ a (2) 2,n b (2) 2 a (k) kk a (k) k+,k 0 a (k) n,k a (k) k,k+ a (k) k,n a (k) k+,k+ a (k) k+,n a (k) n,k+ a (k) n,n b (k) k b (k) k+ b (k) n A (k+) = a () a () 2 a () k a (),k+ a (),n b () a (2) 22 a (2) 2k a (2) 2,k+ a (2) 2,n b (2) 2 a (k) kk a (k) k,k+ a (k) k,n a (k+) k+,k+ a (k+) k+,n 0 a (k+) n,k+ a (k+) n,n b (k) k b (k+) k+ b (k+) n 62

64 M k = q k+,k q k+2,k 0 q n,k 0 0 : A (k+) = M k A (k) L def = M n M n 2 M 2 M A () (63) LA () = A (n) = a () a () 2 a () 3 a () n b () a (2) 22 a (2) 23 a (2) 2n b (2) 2 a (3) 33 a (3) 3n b (3) 3 a (n) nn b (n) n LA = a () a () 2 a () 3 a () n a (2) 22 a (2) 23 a (2) 2n a (3) 33 a (3) 3n a (n) nn def = U, L b b 2 b n = b () b (2) 2 b (n) n (M k ) = q k+,k q k+2,k 0 q n,k 0 0 ( ) L def = L = (M ) (M 2 ) (M n ) 63

65 L = q 2 q 3 q 32 0 q 43 q n q n2 q n3 q n,n LA = U A = LU A LU L det A = det L det U = LU ( ) 6 ( ) A GL(n; R) Gauss r a (k) kk (k =, 2,, r) A k δ k ( δ 0 = ) a (k) kk = δ k δ k n k= a (k) kk k k U k : δ k = a () a (2) 22 a (k) kk 62 LU 62 (LU ) n A = (a ij ) LU A 0 : a a k = det δ k def 0 (k =, 2,, n) a k a kk Crout LU Dolittle LU LDU 6 [] [9] ( Rice [6] ) 64

66 63 LU ( [5] ) 63 Cholesky 63 ( Cholesky ) N A = (a ij ) δ k (k =, 2,, N) 0 A = LDL T (L, D ) L = l 2 l 3 l 32 0 l N l N2 l N,N, D = d 0 d 2 0 d N LDL T (i, j) (l i, l i,2 l i,i 0 0) d l j d 2 l j2 d j l j,j d j 0 0 = mini,j k= d k l ik l jk A i j l ii = i i d i + d k l 2 ik (j = i) k= a ij = d k l ik l jk = i k= d k l ik l jk + d i l ji (j = i +, i + 2,, N) k= i =, 2,, N i d i = a ii d k l 2 ik, k= ) l ji = i (a ij d k l ik l jk d i k= (j = i +,, N) d i, l ji i j d i Gauss i d (i) ii 0 65

67 d = a = δ 0, l j = a j /d (j = 2,, N), d 2 = a 22 d l 2 2 = δ 2 /δ 0, l j2 = (a 2j d l 2 l j )/d 2 (j = 3,, N), d 3 = a 33 d l 2 3 d 2 l 2 32 = δ 3 /δ 2 0, l j3 = (a 3j d l 3 l j d 2 l 32 l j2 )/d 3 (j = 4,, N), A Cholesky LU 632 Cholesky A 62 A LDU : A = LDU L = U T A = U T DU D = G = G 2 = D, G T = G a () a (2) 22 a (k) kk > 0 a () a (2) 22 a (N) NN a (N) NN A = U T DU = U T G 2 U = U T G T GU = (GU) T (GU) = S T S S def = LG A = S T S Cholesky Cholesky ( D G ) S S T S = A s i s j + s 2i s 2j + + s ii s ij = a ij ( ) i s ii = aii s 2 kj, s i ij = a ij s ki s kj /s ii k= k= 66

68 64 A LU A A ( ) Gauss ( ) Crout LU A LU P P A LU Gauss ( ) P P A Crout LU LU Crout Dolittle LU 65 ( ) LINPACK, LAPACK ( ) ( ) rnd() ( ) 67

69 [],, 99 [2], 2003, I, II,, (993, 994) [3],, 985 [4], [5] 7, 984 [6] J R Rice Numerical Methods, Software, and Analysis McGraw-Hill, 983 [7], 985 [8] G D Smith Numerical solution of partial differential equations third edition Clarendon Press Oxford, 986 G D,,, (996) [9], 988 [0], [],, 99 68