Transcription

1

2

3 iii,,.,,,,.,.,,. =,,.,,.,.,,,,,,, fortran,.,,,,.,,,

4 iv (i),., 10,, (ii),, fortran, fortran, C,.. (i),., (ii),.,, fortran Pasal,, for, if then else. C., fortran C.,.,,.,.,,.,,.,.,,,,.

5 v iii 1 Gauss Gauss Gauss Gauss ( 1) fortran C Gauss (1) fortran C Gauss (2)

6 vi fortran C Gauss fortran C Gauss Gauss fortran C rank fortran C (1) fortran C fortran

7 vii C fortran C fortran C fortran C fortran C fortran C

8 viii fortran C ([D] ) ([B] ) ([D] ) ([D] ) ([B] ) ((i) ) ((ii) ) (i) (ii)

9 ix fortran merge Lagrange serendipity bubble at enter of element quadrati bubble at enter of eah fae and 1 bubble at enter of element

10

11 1 1 Gauss (1) 1. (2) Gauss. 1.1 Gauss 1. [A], {b, { b 1 b 2 b 3 b = 0 0 (1.1) 1. [A (1) ]=[A], { (1) = {.

12 2 1 Gauss step 1. [A (1) ] {( 4) / 5 1 {{{{ A (1) 21 A (1) 11 3 { 1 / 5 1 {{{{ A (1) 31 A (1) 11 4 { 0 / 5 1 {{{{ A (1) 41 A (1) 11 [A (1) ]{b = { (1), [A (2) ]{b = { (2) b 1 b 2 b 3 b = 0 0 (1.2) step 2. [A (2) ] 2 0 3, 4. {( 3 16 ) / {{{{ { 4 A (2) 23 1 {{ A (2) 24 A (2) 22 / 14 5 {{ A (2) 22 2 [A (2) ]{b = { (2), [A (3) ]{b = { (3).

13 1.1. Gauss b 1 b 2 b 3 b 4 = (1.3) step 3. [A (3) ] {( 4 20 ) / {{{{ A (3) 34 A (3) 44 [A (3) ]{b = { (3), [A (4) ]{b = { (4) b 1 b 2 b 3 b = (1.4) (1.1), 1 1. (forward redution) {b b 4, b 3, b 2, b 1. (bakward substitution)

14 4 1 Gauss step 1. step 2. step 3. step 4. { b 2 = b 3 = 1 {{ (4) 2 { 8 7 {{ (4) 3 ( { b 1 = 0 {{ ( 4) {{ (4) 1 A (4) 12 b 4 = 7 6 {{ (4) 4 ( 16 ) 5 {{ A (4) {{ b 2 20 ) {{ 7 A (4) {{ b 3 1 {{ A (4) 13 / 5 = 7 {{ 6 5 A (4) {{ b 4 1 {{ A (4) {{ b 3 / 15 = 12 {{ 7 5 A (4) {{ b 4 0 {{ A (4) 14 / 14 = {{ A (4) 22 7 / 5 = {{{{ b 4 A (4) 11 (1.5) (1.6) (1.7) (1.8) 1.1.3,. (bakward redution).,. (4,0) i = (4) i. (.,,. )

15 1.1. Gauss 5 step 1. step 2. b 4 = 7 6 {{ (4,0) 4 (4,1) 1 = 0 {{ (4,0) 1 (4,1) 2 = 1 {{ (4,0) 2 (4,1) 3 = 8 7 {{ (4,0) 3 / 5 = {{ A (4) ( {{ A (4) 14 {{ A (4) {{ A (4) {{ b {{ b 4 ) 7 5 {{ b 4 (1.9) = 0 (1.10) = 2 5 = 36 7 (1.11) (1.12) b 3 = 36 7 {{ (4,1) 3 (4,2) 1 = 0 {{ (4,1) 1 (4,2) 2 = 2 5 {{ (4,1) 2 / 15 = {{ A (4) 33 1 {{ A (4) 13 ( 16 5 {{ A (4) {{ b 3 ) 12 5 {{ b 3 = 12 5 = (1.13) (1.14) (1.15) step 3. b 2 = {{ (4,2) 2 (4,3) 1 = 12 5 {{ (4,2) 1 / 14 = {{ A (4) 22 ( 4) {{ A (4) {{ b 2 (1.16) = 8 (1.17)

16 6 1 Gauss step 4. / b 1 = 8 5 = 8 5 {{{{ (4,3) 1 A (4) 11 (1.18),, 6, 6, 4. step,,. loop index,. skyline Gauss, skyline Gauss,, , 1, i, j (i j) C, 0 [P (i, j : C)].

17 j [P (i, j : C)] = (1.19). 1. i C [P (i, j : C)] [A], [A] i [A] j C. 1.1,,, [A], [A] [S], [P n ]. [P n ]...[P 2 ][P 1 ][A] =[S] (1.20), (1.1) (1.20) {{{{{{ step3 step2 step = (1.21).

18 8 1 Gauss step 1. step 2. step (1.22) (1.23) (1.24) [P (i, j : C)] j, i = j +1 n [P (i +1,i: C i+1 )][P (i +2,i: C i+2 )]...[P (u, i : C n )] = i i C i C i C n 1 0 (1.25)

19 1.2. 9, 1, i i +1 C i+1,c i+2,...c n 0.. [P (i, j : C)] [P (i, j : C)], C i+1,c i+2,...c n C i C n 1 0,, [L i ] = (1.26) C i C n 1 0 [L i ]= i i L i+1,i L n,i 1 0 (1.27), L i+j,i = A(i) i+j,i A (i) i,i (1.28)

20 10 1 Gauss [L i ] [L 1 i ] [L 1 i ]. [L 1 i ]= i i L i+1,i L n,i 1 0 (1.29) [A] n, [L 1 i ],, [A (n) ]. [A (2) ]=[L 1 1 ][A(1) ] [S] =[A (n) ],.. [A (n) ]=[L 1 n 1 ][A(n 1) ] (1.30) [L 1 n 1 ]...[L 1 2 ][L 1 1 ][A] =[S] (1.31) (1.1) (1.29), (1.31) [S] [A (1) ]= [L 1 1 ]= 1 0 L 2,1 1 L 3,1 0 1 L 4, L 2,1 = A(1) 2,1 = 4 A (1) 5 1,1 L 3,1 = A(1) 3,1 = 1 5 A (1) 1,1 L 4,1 = A(1) 4,1 = 0 5 =0 A (1) 1,1 (1.32) (1.33)

21 [L 1 1 ][A(1) ]= [L 1 2 ]= [L 1 2 ][A(2) ]= [L 1 3 ]= [L 1 3 ][A(3) ]= =[A (2) ] (1.34) L 3,2 = A(2) 3,2 = 8 A (2) 7 2,2 0 L 3,2 1 L 4,2 = A(2) 4,2 = 5 (1.35) 0 L 4,2 0 1 A (2) 14 2, =[A (3) ] (1.36) L 4,3 = A(3) 4,3 = 4 (1.37) A (3) 3 3,3 0 0 L 4, =[A (4) ]=[S] (1.38) (1.31) [L n 1 ], [L n 2 ],...,[L 2 ], [L 1 ]. [L]. [A] =[L 1 ][L 2 ]...[L n 1 ][S] (1.39) [L] =[L 1 ][L 2 ]...[L n 1 ] (1.40)

22 12 1 Gauss [L],. 1 0 L 2,1 1. L 3,1 L 3,2.. [L] = 1... L i+1,i... 1 L n,1 L n,2... L n,i... L n,n 1 1 (1.41) [L] (1.39) [A] =[L][S] (1.42) [L] [S].. [S],,,, [S] [D] [S] =[D][U]

23 [D], [U],. S 11 S S 1n S [S] = S 0 nn S 11 S 0 11 /S 11 S 12 /S S 1n /S 11 S = 22 S 22 /S S 2n /S S 0 nn S 0 nn /S nn S 11 1 S 0 12 /S S 1n /S 11 S = 22 1 S 2n /S S 0 nn 1 0 (1.43) S 11 0 S [D] = S 0 nn (1.44) 1 S 12 /S S 1n /S 11 1 S [U] = 2n /S (1.45) [U] [L] 1. [A]. [A] =[L][D][U] (1.46), loop index, [A] [U] =[L T ]. [A]

24 14 1 Gauss. 2 [A] =[L][D][L T ]=[U T ][D][U] (1.47) [A] =[L][D][U], 1 [A]{b = {,. [L][D][U]{b = { [D][U]{b =[L 1 ]{ [U]{b =[D 1 ][L 1 ]{ {b =[U 1 ][D 1 ][L 1 ]{ (1.48) [L 1 ]{. {x =[L 1 ]{, {x 1. [L]{x = { (1.49) 1, (forward substitution).. x 1 = 1 (1.50) x 2 = 2 L 21 x 1 (1.51) x 3 = 3 L 31 x 1 L 32 x 2 (1.52)., skyline. 2 [L][D][L T ],, [U T ][D][U].

25 [D 1 ][L 1 ]{. {x, [D 1 ][L 1 ]{ =[D 1 ]{x {y =[D 1 ]{x {y 1. [D]{y = {x (1.53) 1.. y i = x i /D ii (i =1 n) (1.54) [U 1 ][D 1 ][L 1 ]{. {y {b =[U 1 ]{y, {b 1. [U]{b = {y (1.55) 1., (1.48).,,. 1.1 [A (4) ] [S], { (4) {x =[L 1 ]{. (1.1),. 1 0 [L] = [D] = (1.56) (1.57)

26 16 1 Gauss [U] =[L T ]= [L 1 ]=[L 1 n 1 ] [L 1 2 ][L 1 1 ]= (1.58) (1.59) [L 1 ]{ = = 8 = { (4) (1.60) Gauss Gauss.., , ,,, ( ) 3 15 ( ) , 5,, 5.,.

27 1.3. Gauss Gauss , skyline.. 1.2, L [L] = 21 1 (1.61) L 31 L 32 1 L 41 L 42 L 43 1 D 11 0 D [D] = 22 (1.62) D 33 D U 12 U 13 U 14 1 U [U] = 23 U 24 (1.63) 1 U [L], [D], [U] [A] =[L][D][U]. [A] = D 11 D 11 U 12 D 11 U 13 D 11 U 14 L 21 D 11 L 21 D 11 U 12 + D 22 L 21 D 11 U 13 + D 22 U 23 L 21 D 11 U 14 + D 22 U 24 L 31 D 11 L 31 D 11 U 12 + L 32 D 22 L 31 D 11 U 13 + L 32 D 22 U 23 + D 33 L 31 D 11 U 14 + L 32 D 22 U 24 + D 33 U 34 L 41 D 11 L 41 D 11 U 12 + L 42 D 22 L 41 D 11 U 13 + L 42 D 22 U 23 + L 43 D 33 L 41 D 11 U 14 + L 42 D 22 U 24 + L 43 D 33 U 34 + D 44 (1.64)

28 18 1 Gauss. A 11 = D 11 (1.65) A 1i = D 11 U 1i i =2 n (1.66) A i1 = L i1 D 11 i =2 n (1.67) i 2 i 1 A ij = D ii U ij + L ik D kk U kj (i<j: ) (1.68) k=1 j 1 A ij = L ij D jj + L ik D kk U kj (i>j: ) (1.69) k=1 i 1 A ii = D ii + L ik D kk U ki (i = j) (1.70) k= [L], [D], [U]. j =1 D 11 = A 11 (1.71) j =2 U 12 = A 12 /D 11 (1.72) L 21 = A 21 /D 11 D 22 = A 22 L 21 D 11 U 12 (1.73)

29 1.3. Gauss 19 j =3 U 13 = A 13 /D 11 (1.74) L 31 = A 31 /D 11 U 23 =(A 23 L 21 D 11 U 13 )/D 22 (1.75) L 32 =(A 32 L 31 D 11 U 12 )/D 22 D 33 = A 33 L 31 D 11 U 13 L 32 D 22 U 23 (1.76).. for j =1 D 11 = A 11 for j =2 U 12 = A 12 /D 11 L 21 = A 21 /D 11 D 22 = A 22 L 21 D 11 U 12 for j =3 n U 1j = A 1j /D 11 L j1 = A j1 /D 11 for i =2 j 1 for k =1 i 1 temp U = temp U + L ik D kk U kj temp L = temp L + L jk D kk U ki end for

30 20 1 Gauss U ij =(A ij temp U)/D ii L ji =(A ji temp L)/D ii end for for k =1 j 1 temp = temp + L jk D kk U kj end for D jj = A jj temp end for.,., 1,. ( ) for i =1 x 1 = 1 for i =2 n for j =1 i 1 temp = temp + L ij x j end for x i = i temp end for for i =1 n end for y i = x i /D ii

31 1.3. Gauss 21 ( ) for j = n b n = y n for j = n 2 step = 1 for i =1 j 1 y i = y i U ij b j end for b j 1 = y j 1 end for 1.1,,. skyline. 1.1:

32 22 1 Gauss Gauss ( 1). (1). [A] n n n..... n n n... n (2) {b = {1, 2,...,n T [A]{b = { {. (3) [A]{b = { 1 {b, {1, 2,...,n T. (4) , pu time.. n.. A a 2 n,n b b 1 n 1 n L al 2 n,n D ad 2 n,n U au 2 n,n x x 1 n y y 1 n

33 1.3. Gauss 23,. 2, fortran fortran. init gauss gauss ver1. Makefile FC = g77 F_OPT = -O2 OBJS = main.o mis.o gauss_ver1.o TARGET = gauss_ver1 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f gauss elimination ver. 1 by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 1000) parameter(i_hek = 1) dimension a(n,n),b(n),(n),al(n,n),au(n,n),ad(n,n),x(n),y(n) all init(a,b,,n) if (i_hek.ge. 2) all hek_matrix(a,b,,n) all gauss_ver1(a,b,,al,au,ad,x,y,n)

34 24 1 Gauss if (i_hek.ge. 1) all hek_solution(b,n) stop end mis.f ####################################################################### subroutine init(a,b,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),b(*),(*) do i = 1,n do j = 1,i a(i,j) = dfloat(i) a(j,i) = dfloat(i) do i = 1,n b(i) = dfloat(i) do i = 1,n temp = 0.D0 do j = 1,n temp = temp + a(i,j) * b(j) (i) = temp return end ####################################################################### subroutine hek_matrix(a,b,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),b(*),(*) write(*,*) matrix A do i = 1,n

35 1.3. Gauss 25 write(*,1000) (a(i,j),j=1,n) write(*,*) vetor b write(*,1000) (b(i),i=1,n) write(*,*) vetor write(*,1000) ((i),i=1,n) 1000 format(8e10.3) return end ####################################################################### subroutine hek_solution(b,n) ####################################################################### impliit real*8(a-h,o-z) dimension b(*) write(*,*) solution vetor b write(*,1000) (b(i),i=1,n) 1000 format(8e10.3) return end gauss ver1.f ####################################################################### subroutine gauss_ver1(a,b,,al,au,ad,x,y,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),b(*),(*),al(n,*),au(n,*),ad(n,*),x(*),y(*) ad(1,1) = a(1,1) au(1,2) = a(1,2) / ad(1,1) al(2,1) = a(2,1) / ad(1,1) ad(2,2) = a(2,2) - al(2,1) * ad(1,1) * au(1,2) do j = 3,n au(1,j) = a(1,j) / ad(1,1) al(j,1) = a(j,1) / ad(1,1)

36 26 1 Gauss do i = 2, j-1 temp_u = 0.d0 temp_l = 0.d0 do k = 1, i-1 temp_u = temp_u 1 + al(i,k) * ad(k,k) * au(k,j) temp_l = temp_l 1 + al(j,k) * ad(k,k) * au(k,i) au(i,j) = (a(i,j) - temp_u) / ad(i,i) al(j,i) = (a(j,i) - temp_l) / ad(i,i) temp = 0.d0 do k = 1,j-1 temp = temp + al(j,k) * ad(k,k) * au(k,j) ad(j,j) = a(j,j) - temp x(1) = (1) do i = 2,n temp = 0.d0 do j = 1,i-1 temp = temp + al(i,j) * x(j) x(i) = (i) - temp do i = 1,n y(i) = x(i) / ad(i,i) b(n) = y(n) do j = n,2,-1 do i = 1,j-1 y(i) = y(i) - au(i,j) * b(j) b(j-1) = y(j-1) return end

37 1.3. Gauss C C, 1, C 0., n x i x, x 1 = x[0], x 2 = x[1],, x n = x[n-1],.,. (1), 0. (2) 0. (1) C, 0,.,.,, 1, 1,,. (2), fortran,,.,.,, C,.,,,,,,., x i (x 1 =1.0,x 2 = 2.0,x 3 =3.0,,x n = n) 1 x.

38 28 1 Gauss double *x; x = (double *)mallo(sizeof(double)*n); for (i = 0; i < n; ++i){ x[i] = (double)(i+1);, x i = i = x[i-1]... mallo x, x, *(x), x[0], x., xx, x 1. *(xx+1) *(x)., xx[1], x[0]. xx 1 n.. xx x, x i = i = x[i-1] = xx[i]. double *x, *xx; x = (double *)mallo(sizeof(double)*n); xx = x - 1; for (i = 1; i <= n; ++i){ xx[i] = (double)(i); 2, 1..

39 1.3. Gauss [A] = , fortran,, 1 a[0] a[n-1], 2 a[0 + n] a[n-1 + n]., i j i 1+(j 1) n, a[i-1+(j-1)*n]. double *a; a = (double *)mallo(sizeof(double)*n*n); for (j = 0; j < n; ++j){ for (i = 0; i < n; ++i){ a[i + j * n] = (double)((i+1)*10 + (j+1));, aa, a 1+n. *(aa+1+1*n) *(a),, aa[i+j*n], a[(i-1)+(j-1)*n].. double *a, *aa; a = (double *)mallo(sizeof(double)*n*n); aa = a - n - 1; for (j = 1; j <= n; ++j){ for (i = 1; i <= n; ++i){ aa[i + j * n] = (double)(i*10 + j);

40 30 1 Gauss,. init gauss gauss ver1. Makefile CC = g C_OPT = -O2 OBJS = main.o mis.o gauss_ver1.o TARGET = gauss_ver1 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main. /* gauss elimination ver. 1 by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 1000 #define I_CHECK 1 void init(double *a, double *b, double *, int n); void hek_matrix(double *a, double *b, double *, int n); void gauss_ver1(double *a, double *b, double *, double *al, double *au, double *ad, double *x, double *y, int n); void hek_solution(double *b, int n); main() { double *a,*b,*,*al,*au,*ad,*x,*y; int n;

41 1.3. Gauss 31 n = N; a = (double *)mallo(sizeof(double)*n*n); al = (double *)mallo(sizeof(double)*n*n); au = (double *)mallo(sizeof(double)*n*n); ad = (double *)mallo(sizeof(double)*n*n); b = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); x = (double *)mallo(sizeof(double)*n); y = (double *)mallo(sizeof(double)*n); init(a,b,,n); if (I_CHECK >= 2) hek_matrix(a,b,,n); gauss_ver1(a,b,,al,au,ad,x,y,n); if (I_CHECK >= 1) hek_solution(b,n); mis. void init(double *arg_a, double *arg_b, double *arg_, int n) { double *a, *b, *; int i,j, ij, ji; double temp; a = arg_a-1-n; b = arg_b-1; = arg_-1; for (i = 1; i <= n; ++i){ for (j = 1; j <= n; ++j){ a[i+j*n] = (double)(i); a[j+i*n] = (double)(i); for (i = 1; i <= n; ++i){ b[i] = (double)(i); for (i = 1; i <= n; ++i){ temp = 0.0;

42 32 1 Gauss for (j = 1; j <= n; ++j){ temp = temp + a[i+j*n] * b[j]; [i] = temp; void hek_matrix(double *arg_a, double *arg_b, double *arg_, int n) { double *a, *b, *; int i,j; a = arg_a-1-n; b = arg_b-1; = arg_-1; printf("matrix A\n"); for (i = 1; i <= n; ++i){ for (j = 1; j <= n; ++j) printf("% 10.3E",a[i+j*n]); printf("\n"); printf("vetor b\n"); for (i = 1; i <= n; ++i) printf("% 10.3E",b[i]); printf("\n"); printf("vetor \n"); for (i = 1; i <= n; ++i) printf("% 10.3E",[i]); printf("\n"); void hek_solution(double *arg_b, int n) { double *b; int i; b = arg_b-1; printf("solution vetor b\n"); for (i = 1; i <= n; ++i) printf("% 10.3E",b[i]); printf("\n"); gauss ver1. void gauss_ver1(double *arg_a, double *arg_b, double *arg_, double *arg_al, double *arg_au, double *arg_ad, double *arg_x, double *arg_y, int n) {

43 1.3. Gauss 33 double *a, *b, *, *al, *au, *ad, *x, *y; int i, j, k; double temp_u, temp_l, temp; a = arg_a -1-n; al = arg_al-1-n; au = arg_au-1-n; ad = arg_ad-1-n; b = arg_b-1; = arg_-1; x = arg_x-1; y = arg_y-1; ad[1+1*n] = a[1+1*n]; au[1+2*n] = a[1+2*n] / ad[1+1*n]; al[2+1*n] = a[2+2*n] / ad[1+1*n]; ad[2+2*n] = a[2+2*n] - al[2+1*n] * ad[1+1*n] * au[1+2*n]; for (j = 3; j <= n; ++j){ au[1+j*n] = a[1+j*n] / ad[1+1*n]; al[j+1*n] = a[j+1*n] / ad[1+1*n]; for (i = 2; i <= j-1; ++i){ temp_u = 0.0; temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_u = temp_u + al[i+k*n] * ad[k+k*n] * au[k+j*n]; temp_l = temp_l + al[j+k*n] * ad[k+k*n] * au[k+i*n]; au[i+j*n] = (a[i+j*n] - temp_u) / ad[i+i*n]; al[j+i*n] = (a[j+i*n] - temp_l) / ad[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + al[j+k*n] * ad[k+k*n] * au[k+j*n]; ad[j+j*n] = a[j+j*n] - temp; x[1] = [1]; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + al[i+j*n] * x[j];

44 34 1 Gauss x[i] = [i] - temp; for (i = 1; i <= n; ++i){ y[i] = x[i] / ad[i+i*n]; b[n] = y[n]; for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ y[i] = y[i] - au[i+j*n] * b[j]; b[j-1] = y[j-1]; 1.4 Gauss (1),,.., , 1.

45 1.4. Gauss (1) 35 x(1) = (1) (1.77) temp = 0.D0 (1.78) temp = al(2, 1)x(1) (1.79) x(2) = (2) temp (1.80) temp = 0.D0 (1.81) temp = al(3, 1)x(1) + al(3, 2)x(2) (1.82) x(3) = (3) temp (1.83). x(i), x(1),...,x(i 1) (i), (1),...,(i 1)., x(i) (i) x. y,. 2.1 ( { ) Gauss, pu time OS fortran fortran. gauss ver1 gauss ver21. Makefile

46 36 1 Gauss FC = g77 F_OPT = -O2 OBJS = main.o mis.o gauss_ver21.o TARGET = gauss_ver21 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f gauss elimination ver. 2.1 by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 1000) parameter(i_hek = 1) dimension a(n,n),b(n),(n),al(n,n),au(n,n),ad(n,n) all init(a,b,,n) if (i_hek.ge. 2) all hek_matrix(a,b,,n) all gauss_ver21(a,,al,au,ad,n) if (i_hek.ge. 1) all hek_solution(,n) stop end mis.f gauss ver1 gauss ver21.f ####################################################################### subroutine gauss_ver21(a,,al,au,ad,n) #######################################################################

47 1.4. Gauss (1) 37 impliit real*8(a-h,o-z) dimension a(n,*),(*),al(n,*),au(n,*),ad(n,*) ad(1,1) = a(1,1) au(1,2) = a(1,2) / ad(1,1) al(2,1) = a(2,1) / ad(1,1) ad(2,2) = a(2,2) - al(2,1) * ad(1,1) * au(1,2) do j = 3,n au(1,j) = a(1,j) / ad(1,1) al(j,1) = a(j,1) / ad(1,1) do i = 2, j-1 temp_u = 0.d0 temp_l = 0.d0 do k = 1, i-1 temp_u = temp_u 1 + al(i,k) * ad(k,k) * au(k,j) temp_l = temp_l 1 + al(j,k) * ad(k,k) * au(k,i) au(i,j) = (a(i,j) - temp_u) / ad(i,i) al(j,i) = (a(j,i) - temp_l) / ad(i,i) temp = 0.d0 do k = 1,j-1 temp = temp + al(j,k) * ad(k,k) * au(k,j) ad(j,j) = a(j,j) - temp do i = 2,n temp = 0.d0 do j = 1,i-1 temp = temp + al(i,j) * (j) (i) = (i) - temp do i = 1,n (i) = (i) / ad(i,i)

48 38 1 Gauss do j = n,2,-1 do i = 1,j-1 (i) = (i) - au(i,j) * (j) return end C C. gauss ver1 gauss ver21. Makefile CC = g C_OPT = -O2 OBJS = main.o mis.o gauss_ver21.o TARGET = gauss_ver21 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main. /* gauss elimination ver. 2.1 by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 1000 #define I_CHECK 1 void init(double *a, double *b, double *, int n); void hek_matrix(double *a, double *b, double *, int n);

49 1.4. Gauss (1) 39 void gauss_ver21(double *a, double *, double *al, double *au, double *ad, int n); void hek_solution(double *b, int n); main() { double *a,*b,*,*al,*au,*ad; int n; n = N; a = (double *)mallo(sizeof(double)*n*n); al = (double *)mallo(sizeof(double)*n*n); au = (double *)mallo(sizeof(double)*n*n); ad = (double *)mallo(sizeof(double)*n*n); b = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); init(a,b,,n); if (I_CHECK >= 2) hek_matrix(a,b,,n); gauss_ver21(a,,al,au,ad,n); if (I_CHECK >= 1) hek_solution(b,n); mis. gauss ver1 gauss ver21. void gauss_ver21(double *arg_a, double *arg_, double *arg_al, double *arg_au, double *arg_ad, int n) { double *a, *, *al, *au, *ad; int i, j, k; double temp_u, temp_l, temp; a = arg_a -1-n; al = arg_al-1-n; au = arg_au-1-n; ad = arg_ad-1-n; = arg_-1;

50 40 1 Gauss ad[1+1*n] = a[1+1*n]; au[1+2*n] = a[1+2*n] / ad[1+1*n]; al[2+1*n] = a[2+2*n] / ad[1+1*n]; ad[2+2*n] = a[2+2*n] - al[2+1*n] * ad[1+1*n] * au[1+2*n]; for (j = 3; j <= n; ++j){ au[1+j*n] = a[1+j*n] / ad[1+1*n]; al[j+1*n] = a[j+1*n] / ad[1+1*n]; for (i = 2; i <= j-1; ++i){ temp_u = 0.0; temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_u = temp_u + al[i+k*n] * ad[k+k*n] * au[k+j*n]; temp_l = temp_l + al[j+k*n] * ad[k+k*n] * au[k+i*n]; au[i+j*n] = (a[i+j*n] - temp_u) / ad[i+i*n]; al[j+i*n] = (a[j+i*n] - temp_l) / ad[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + al[j+k*n] * ad[k+k*n] * au[k+j*n]; ad[j+j*n] = a[j+j*n] - temp; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + al[i+j*n] * [j]; [i] = [i] - temp; for (i = 1; i <= n; ++i){ [i] = [i] / ad[i+i*n]; for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ [i] = [i] - au[i+j*n] * [j];

51 1.5. Gauss (2) Gauss (2) Gauss,. 1.2(a) A 11,. j D kk j U kj k =1 k =2. j i L ik k =1 k = i 1 U ij k = i 1 (a) (b) 1.2: U ij (i j) A ij, 1.2(b) U kj (k =1 i 1),L ik (k =1 i 1),D kk (k =1 i 1), [A].,. U kj,l ik,d kk, [U], [L]. [A] [D], [U], [L]

52 42 1 Gauss [U], [D], [L] pu time OS fortran fortran. gauss ver21 gauss ver22. Makefile FC = g77 F_OPT = -O2 OBJS = main.o mis.o gauss_ver22.o TARGET = gauss_ver22 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f gauss elimination ver. 2.2 by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 1000) parameter(i_hek = 1) dimension a(n,n),b(n),(n)

53 1.5. Gauss (2) 43 all init(a,b,,n) if (i_hek.ge. 2) all hek_matrix(a,b,,n) all gauss_ver22(a,,n) if (i_hek.ge. 1) all hek_solution(,n) stop end mis.f gauss ver1 gauss ver22.f ####################################################################### subroutine gauss_ver22(a,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*) a(1,1) = a(1,1) a(1,2) = a(1,2) / a(1,1) a(2,1) = a(2,1) / a(1,1) a(2,2) = a(2,2) - a(2,1) * a(1,1) * a(1,2) do j = 3,n a(1,j) = a(1,j) / a(1,1) a(j,1) = a(j,1) / a(1,1) do i = 2, j-1 temp_u = 0.d0 temp_l = 0.d0 do k = 1, i-1 temp_u = temp_u 1 + a(i,k) * a(k,k) * a(k,j) temp_l = temp_l 1 + a(j,k) * a(k,k) * a(k,i) a(i,j) = (a(i,j) - temp_u) / a(i,i) a(j,i) = (a(j,i) - temp_l) / a(i,i)

54 44 1 Gauss temp = 0.d0 do k = 1,j-1 temp = temp + a(j,k) * a(k,k) * a(k,j) a(j,j) = a(j,j) - temp do i = 2,n temp = 0.d0 do j = 1,i-1 temp = temp + a(i,j) * (j) (i) = (i) - temp do i = 1,n (i) = (i) / a(i,i) do j = n,2,-1 do i = 1,j-1 (i) = (i) - a(i,j) * (j) return end C C. gauss ver21 gauss ver22. Makefile CC = g C_OPT = -O2 OBJS = main.o mis.o gauss_ver22.o

55 1.5. Gauss (2) 45 TARGET = gauss_ver22 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main. /* gauss elimination ver. 2.2 by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 1000 #define I_CHECK 1 void init(double *a, double *b, double *, int n); void hek_matrix(double *a, double *b, double *, int n); void gauss_ver22(double *a, double *, int n); void hek_solution(double *b, int n); main() { double *a,*b,*; int n; n = N; a b = (double *)mallo(sizeof(double)*n*n); = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); init(a,b,,n); if (I_CHECK >= 2) hek_matrix(a,b,,n); gauss_ver22(a,,n); if (I_CHECK >= 1) hek_solution(b,n); mis. gauss ver1

56 46 1 Gauss gauss ver22. void gauss_ver22(double *arg_a, double *arg_, int n) { double *a, *; int i, j, k; double temp_u, temp_l, temp; a = arg_a -1-n; = arg_-1; a[1+2*n] = a[1+2*n] / a[1+1*n]; a[2+1*n] = a[2+2*n] / a[1+1*n]; a[2+2*n] = a[2+2*n] - a[2+1*n] * a[1+1*n] * a[1+2*n]; for (j = 3; j <= n; ++j){ a[1+j*n] = a[1+j*n] / a[1+1*n]; a[j+1*n] = a[j+1*n] / a[1+1*n]; for (i = 2; i <= j-1; ++i){ temp_u = 0.0; temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_u = temp_u + a[i+k*n] * a[k+k*n] * a[k+j*n]; temp_l = temp_l + a[j+k*n] * a[k+k*n] * a[k+i*n]; a[i+j*n] = (a[i+j*n] - temp_u) / a[i+i*n]; a[j+i*n] = (a[j+i*n] - temp_l) / a[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + a[j+k*n] * a[k+k*n] * a[k+j*n]; a[j+j*n] = a[j+j*n] - temp; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + a[i+j*n] * [j]; [i] = [i] - temp; for (i = 1; i <= n; ++i){ [i] = [i] / a[i+i*n];

57 1.6. Gauss 47 for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ [i] = [i] - a[i+j*n] * [j]; 1.6 Gauss ,,. i =1 U ij. U 1j = A 1j /D 11 U 2j =(A 2j L 21 D 11 U 1j )/D 22 U 3j =(A 3j L 31 D 11 U 1j L 32 D 22 U 2j )/D 33. U ij =(A ij L i1 D 11 U 1j L i2 D 22 U 2j... L ii 1 D i 1i 1 U i 1j )/D ii U i+1j =(A i+1j L i+11 D 11 U 1j L i+12 D 22 U 2j... L i+1i D ii U ij )/D i+1i+1 (1.84), i =2, U ij U kj (k =1 i 1) D kk, U kj., (D kk U kj ), (L jk D kk ), U ij,l ji,. for j =1 D 11 = A 11

58 48 1 Gauss j k =1. i i +1 D kk L ik L i+1k U kj U ij U i+1j k = i 1 1.3: for j =2 U 12 = A 12 /D 11 L 21 = A 21 /D 11 D 22 = A 22 L 21 D 11 U 12 for j =3 n for i =2 j 1 for k =1 i 1 temp U = temp U + L ik (D kk U kj ) temp L = temp L +(L jk D kk )U ki end for

59 1.6. Gauss 49 (D ii U ij )=A ij temp U (L ji D ii )=A ji temp L end for for i =1 j 1 U ij =(D ii U ij )/D ii L ji =(L ji D ii )/D ii end for for k =1 j 1 temp = temp + L jk D kk U kj end for D jj = A jj temp end for pu time (D kk U kj ), (L jk D kk ) U kj, L jk, A. OS fortran fortran. gauss ver22 gauss ver3. Makefile FC = g77 F_OPT = -O2

60 50 1 Gauss OBJS = main.o mis.o gauss_ver3.o TARGET = gauss_ver3 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f gauss elimination ver. 3 by WATANABE Hiroshi (2000 May) impliit real*8(a-h,o-z) parameter(n = 1000) parameter(i_hek = 1) dimension a(n,n),b(n),(n) all init(a,b,,n) if (i_hek.ge. 2) all hek_matrix(a,b,,n) all gauss_ver3(a,,n) if (i_hek.ge. 1) all hek_solution(,n) stop end mis.f gauss ver1 gauss ver3.f ####################################################################### subroutine gauss_ver3(a,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*)

61 1.6. Gauss 51 a(1,2) = a(1,2) / a(1,1) a(2,1) = a(2,1) / a(1,1) a(2,2) = a(2,2) - a(2,1) * a(1,1) * a(1,2) do j = 3,n do i = 2, j-1 temp_u = 0.d0 temp_l = 0.d0 do k = 1, i-1 temp_u = temp_u + a(i,k) * a(k,j) temp_l = temp_l + a(j,k) * a(k,i) a(i,j) = a(i,j) - temp_u a(j,i) = a(j,i) - temp_l do i = 1,j-1 a(i,j) = a(i,j) / a(i,i) a(j,i) = a(j,i) / a(i,i) temp = 0.d0 do k = 1,j-1 temp = temp + a(j,k) * a(k,k) * a(k,j) a(j,j) = a(j,j) - temp do i = 2,n do j = 1,i-1 (i) = (i) - a(i,j) * (j) do i = 1,n (i) = (i) / a(i,i) do j = n,2,-1 do i = 1,j-1 (i) = (i) - a(i,j) * (j) return end

62 52 1 Gauss C C. gauss ver22 gauss ver3. Makefile CC = g C_OPT = -O2 OBJS = main.o mis.o gauss_ver3.o TARGET = gauss_ver3 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main. /* gauss elimination ver. 3 by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 1000 #define I_CHECK 1 void init(double *a, double *b, double *, int n); void hek_matrix(double *a, double *b, double *, int n); void gauss_ver3(double *a, double *, int n); void hek_solution(double *b, int n); main() { double *a,*b,*; int n;

63 1.6. Gauss 53 n = N; a b = (double *)mallo(sizeof(double)*n*n); = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); init(a,b,,n); if (I_CHECK >= 2) hek_matrix(a,b,,n); gauss_ver3(a,,n); if (I_CHECK >= 1) hek_solution(b,n); mis. gauss ver1 gauss ver3. void gauss_ver3(double *arg_a, double *arg_, int n) { double *a, *; int i, j, k; double temp_u, temp_l, temp; a = arg_a -1-n; = arg_-1; a[1+2*n] = a[1+2*n] / a[1+1*n]; a[2+1*n] = a[2+2*n] / a[1+1*n]; a[2+2*n] = a[2+2*n] - a[2+1*n] * a[1+1*n] * a[1+2*n]; for (j = 3; j <= n; ++j){ for (i = 2; i <= j-1; ++i){ temp_u = 0.0; temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_u = temp_u + a[i+k*n] * a[k+j*n]; temp_l = temp_l + a[j+k*n] * a[k+i*n]; a[i+j*n] = a[i+j*n] - temp_u; a[j+i*n] = a[j+i*n] - temp_l;

64 54 1 Gauss for (i = 1; i <= j-1; ++i){ a[i+j*n] = a[i+j*n] / a[i+i*n]; a[j+i*n] = a[j+i*n] / a[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + a[j+k*n] * a[k+k*n] * a[k+j*n]; a[j+j*n] = a[j+j*n] - temp; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + a[i+j*n] * [j]; [i] = [i] - temp; for (i = 1; i <= n; ++i){ [i] = [i] / a[i+i*n]; for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ [i] = [i] - a[i+j*n] * [j]; 1.7 Gauss Gauss [A] [A] =[L][D][L T ]=[U T ][D][U]. 3,. [U].

65 1.7. Gauss Gauss 55 for j =1 D 11 = A 11 for j =2 U 12 = A 12 /D 11 D 22 = A 22 U 12 D 11 U 12 for j =3 n for i =1 j 1 for k =1 i 1 temp = temp + U ki (D kk U kj ) end for D ii U ij = A ij temp end for for i =1 j 1 U ij = D ii U ij /D ii end for for i =1 j 1 temp = temp + U ij D ii U ij end for D jj = A jj temp end for 4, Gauss. [U] [L], pu time. 3.

66 56 1 Gauss fortran fortran. gauss_ver4l [L], gauss_ver4u [U]., pu time..,, CPU., CPU,,..,. fortran 2 1. a(i, j), a n (j 1) n + i. a(i +1,j) a(i, j +1) a(i, j) a(i +1,j) a(i, j +1) n.. Makefile FC = g77 F_OPT = -O2 OBJS = main.o mis.o gauss_ver4l.o gauss_ver4u.o TARGET = gauss_ver4 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $<

67 1.7. Gauss Gauss 57 j a(1,j) a(1,j). a(i, j) a(i, j +1) a(i +1,j). (j 1) n n n a(1,j) a(1,j). a(i, j) a(i +1,j).. a(i, j +1). 1.4: main.f gauss elimination ver. 4l, 4u by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 1000) parameter(i_hek = 1) dimension a(n,n),b(n),(n) all init(a,b,,n) if (i_hek.ge. 2) all hek_matrix(a,b,,n) all gauss_ver4l(a,,n) all gauss_ver4u(a,,n)

68 58 1 Gauss if (i_hek.ge. 1) all hek_solution(,n) stop end mis.f gauss ver1 gauss ver4l.f ####################################################################### subroutine gauss_ver4l(a,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*) a(2,1) = a(2,1) / a(1,1) a(2,2) = a(2,2) - a(2,1) * a(1,1) * a(2,1) do j = 3,n do i = 2, j-1 temp_l = 0.d0 do k = 1, i-1 temp_l = temp_l + a(j,k) * a(i,k) a(j,i) = a(j,i) - temp_l do i = 1,j-1 a(j,i) = a(j,i) / a(i,i) temp = 0.d0 do k = 1,j-1 temp = temp + a(j,k) * a(k,k) * a(j,k) a(j,j) = a(j,j) - temp do i = 2,n do j = 1,i-1 (i) = (i) - a(i,j) * (j)

69 1.7. Gauss Gauss 59 do i = 1,n (i) = (i) / a(i,i) do j = n,2,-1 do i = 1,j-1 (i) = (i) - a(j,i) * (j) return end gauss ver4u.f ####################################################################### subroutine gauss_ver4u(a,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*) a(1,2) = a(1,2) / a(1,1) a(2,2) = a(2,2) - a(1,2) * a(1,1) * a(1,2) do j = 3,n do i = 2, j-1 temp_u = 0.d0 do k = 1, i-1 temp_u = temp_u + a(k,i) * a(k,j) a(i,j) = a(i,j) - temp_u do i = 1,j-1 a(i,j) = a(i,j) / a(i,i) temp = 0.d0 do k = 1,j-1 temp = temp + a(k,j) * a(k,k) * a(k,j) a(j,j) = a(j,j) - temp

70 60 1 Gauss do i = 2,n do j = 1,i-1 (i) = (i) - a(j,i) * (j) do i = 1,n (i) = (i) / a(i,i) do j = n,2,-1 do i = 1,j-1 (i) = (i) - a(i,j) * (j) return end C fortran. gauss_ver4l [L], gauss_ver4u [U]. Makefile CC = g C_OPT = -O2 OBJS = main.o mis.o gauss_ver4l.o gauss_ver4u.o TARGET = gauss_ver4 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main.

71 1.7. Gauss Gauss 61 /* gauss elimination ver. 4l, 4u by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 1000 #define I_CHECK 1 void init(double *a, double *b, double *, int n); void hek_matrix(double *a, double *b, double *, int n); void gauss_ver4l(double *a, double *, int n); void gauss_ver4u(double *a, double *, int n); void hek_solution(double *b, int n); main() { double *a,*b,*; int n; n = N; a b = (double *)mallo(sizeof(double)*n*n); = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); init(a,b,,n); if (I_CHECK >= 2) hek_matrix(a,b,,n); /* gauss_ver4l(a,,n); */ gauss_ver4u(a,,n); if (I_CHECK >= 1) hek_solution(b,n); mis. gauss ver1 gauss ver4l. void gauss_ver4l(double *arg_a, double *arg_, int n) { double *a, *; int i, j, k; double temp_u, temp_l, temp;

72 62 1 Gauss a = arg_a -1-n; = arg_-1; a[2+1*n] = a[2+2*n] / a[1+1*n]; a[2+2*n] = a[2+2*n] - a[2+1*n] * a[1+1*n] * a[2+1*n]; for (j = 3; j <= n; ++j){ for (i = 2; i <= j-1; ++i){ temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_l = temp_l + a[j+k*n] * a[i+k*n]; a[j+i*n] = a[j+i*n] - temp_l; for (i = 1; i <= j-1; ++i){ a[j+i*n] = a[j+i*n] / a[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + a[j+k*n] * a[k+k*n] * a[j+k*n]; a[j+j*n] = a[j+j*n] - temp; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + a[i+j*n] * [j]; [i] = [i] - temp; for (i = 1; i <= n; ++i){ [i] = [i] / a[i+i*n]; for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ [i] = [i] - a[j+i*n] * [j];

73 1.7. Gauss Gauss 63 gauss ver4u. void gauss_ver4u(double *arg_a, double *arg_, int n) { double *a, *; int i, j, k; double temp_u, temp_l, temp; a = arg_a-1-n; = arg_-1; a[1+2*n] = a[1+2*n] / a[1+1*n]; a[2+2*n] = a[2+2*n] - a[1+2*n] * a[1+1*n] * a[1+2*n]; for (j = 3; j <= n; ++j){ for (i = 2; i <= j-1; ++i){ temp_u = 0.0; temp_l = 0.0; for (k = 1; k <= i-1; ++k){ temp_u = temp_u + a[k+i*n] * a[k+j*n]; a[i+j*n] = a[i+j*n] - temp_u; for (i = 1; i <= j-1; ++i){ a[i+j*n] = a[i+j*n] / a[i+i*n]; temp = 0.0; for (k = 1; k <= j-1; ++k){ temp = temp + a[k+j*n] * a[k+k*n] * a[k+j*n]; a[j+j*n] = a[j+j*n] - temp; for (i = 2; i <= n; ++i){ temp = 0.0; for (j = 1; j <= i-1; ++j){ temp = temp + a[j+i*n] * [j]; [i] = [i] - temp; for (i = 1; i <= n; ++i){ [i] = [i] / a[i+i*n];

74 64 1 Gauss for (j = n; j >= 2; --j){ for (i = 1; i <= j-1; ++i){ [i] = [i] - a[i+j*n] * [j]; ,., x 1 n +1 ( m i (i =1 n + 1)) n ( k i (i =1 n)). x i, f i (i =1 n +1),. m 1 ü 1 = f 1 + k 1 (u 2 u 1 ) m 2 ü 2 = f 2 + k 1 (u 1 u 2 ) + k 2 (u 3 u 2 ) m 3 ü 3 = f 3 + k 2 (u 2 u 3 ) + k 3 (u 4 u 3 ). m n ü n = f n + k n 1 (u n 1 u n ) + k n (u n+1 u n ) m n+1 u n+1 = f n+1 + k n (u n u n+1 ) u i (i =1 n +1) m i x 1..

75 x 3 x n n+1 x 1 1.5: - k 1 k 1 k 1 k 1 + k 2 k 2 k 2 k 2 + k 3 k 3... k n 1 k n 1 + k n k n k n k n+1 u 1 u 2 u 3. u n u n+1 = f 1 m 1 ü 1 f 2 m 2 ü 2 f 3 m 3 ü 3. f n m n ü n f n+1 m n+1 ü n+1 (1.85), ü 1 =ü 2 = =ü n =ü n+1 =0, u 1 = u 2 = = u n = u n+1 =0.,,., rank n. ( ), u 1 =0,,., v 1 = u 2 u 1,v 2 = u 3 u 2, v n = u n+1 u n.,,,.,, u 1 =0.,,, u j = a,. u j, f j.,, 0.

76 66 1 Gauss. A 11 A A 1n b 1 C 1 A 21. b 2 C 2 = (1.86).... A n A nn b n C n, C j b j.. A 11 + A 12 b A 1j b j + + A 1n b n = C 1 A 21 + A 22 b A 2j b j + + A 2n b n = C 2.. (1.87) A j1 + A j2 b A jj b j + + A jn b n = C j.. A n1 + A n2 b A nj b j + + A nn b n = C n j C j, C j b i (i =1 j 1,j+1 n), C j = A j1 b 1 + A j2 b Ajjb j + + A jn b n (1.88)., b i (i =1 j 1,j+1 n). j n 1. A 11 + A 12 b A 1j b j + + A 1n b n = C 1 A 21 + A 22 b A 2j b j + + A 2n b n = C 2.. A j A j 1 2 b A j 1 j b j + + A j 1 n b n = C j 1 A j A j+1 2 b A j+1 j b j + + A j 1 n b n = C j+1.. A n1 + A n2 b A nj b j + + A nn b n = C n (1.89)

77 A 1j b j,a 2j b j,,a nj b j. A 11 + A 12 b A 1 j 1 b j 1 + A 1 j+1 b j A 1n b n = C 1 A 1j b j A 21 + A 22 b A 2 j 1 b j 1 + A 2 j+1 b j A 2n b n = C 2 A 2j b j.. A j A j 1 2 b A j 1 j 1 b j 1 + A j 1 j+1 b j A j 1 n b n = C j 1 A j 1 j b j A j A j+1 2 b A j+1 j 1 b j 1 + A j+1 j+1 b j A j+1 n b n = C j+1 A j+1 j b j. A n1 + A n2 b A nj 1 b j 1 + A nj+1 b j A nn b n = C n A nj b j (1.90),., j j, j,, j b j j.. A A 1 j 1 A 1 j+1... A 1n A A 2 j 1 A 2 j+1... A 2n. A j A j 1 j 1 A j+1 j+1... A j 1 n A j A j+1 j+1 A j+1 j+1... A j+1 n. A n1... A nj 1 A nj+1... A nn b 1 b 2. = b j 1 b j+1. b n C 1 A 1j C 2 A 2j. b j C j 1 C j+1 A j+1 j. C n. A j 1 j. A nj (1.91) b j 0 2., b j. b j(1),b j(2),,b j(m). j(1),j(2),,j(m). b j(1),b j(2),,b j(m)..

78 68 1 Gauss rank, n =2. [P (i, j; C)], [R(i; C)]. [P (i, j; C)] = i j C (1.92) [R(i; C)] = i i C (1.93) (1.94). k 1 k 1 0 [K] = k 1 k 1 + k 2 k 2 0 k 2 k 2 (1.95)

79 k 1 k 1 0 k 1 k 1 0 [P (2, 1; 1)][K] = k 1 k 1 + k 2 k 2 = 0 k 2 k 2 (1.96) k 2 k 2 0 k 2 k k 1 k 1 0 k 1 k 1 0 [P (3, 2; 1)][P (2, 1; 1)][K] = k 2 k 2 = 0 k 2 k 2 (1.97) k 2 k k 1 k k [P (3, 2; 1)][P (2, 1; 1)][K][P (1, 2; 1)] = 0 k 2 k = 0 k 2 k (1.98) k k [P (3, 2; 1)][P (2, 1; 1)][K][P (1, 2; 1)][P (2, 3; 1)] = 0 k 2 k = 0 k (1.99) [R(2; 1 )][R(1; 1 )][P (3, 2; 1)][P (2, 1; 1)][K][P (1, 2; 1)][P (2, 3; 1)] k 2 k k k = 1 0 k k 2 0 = (1.100) [K] rank 2. (n +1) (n +1) [K], [K] rank n.

80 70 1 Gauss [R(n; 1 k n )] [R(1; 1 k 1 )][P (n +1,n; 1)] [P (i +1,i; 1)] [P (2, 1; 1)][K][P (1, 2; 1)] [P (j, j + 1; 1)] [P (n, n + 1; 1)] n = (1.101) {{ n skyline 1 b j(1),,b j(m) b j(1),, b j(m) j(1),,j(m) b j(1) = b j(1),,b j(m) = b j(m) j(1),,j(m) j(1),,j(m) 0 1 j(1),,j(m) b j(1),, b j(m) 6,m=3,j(1) = 1,j(2) = 3,j(3) = 6

81 A 11 A 12 A 13 A 14 A 15 A 16 A 21 A 22 A 23 A 24 A 25 A 26 A 31 A 32 A 33 A 34 A 35 A 36 A 41 A 42 A 43 A 44 A 45 A 46 A 51 A 52 A 53 A 54 A 55 A 56 A 61 A 62 A 63 A 64 A 65 A 66 b 1 b 2 b 3 b 4 b 5 b 6 = (1.102) b 1 b A 22 0 A 24 A 25 0 b 2 2 A 21 A 23 A b 3 b = b 1 b 3 b 6 0 A 42 0 A 44 A 45 0 b 4 4 A 41 A 43 A 46 0 A 52 0 A 54 A 55 0 b 5 5 A 51 A 53 A b b 6 (1.103)

82 72 1 Gauss. for i =1 m for k =1 n end for end for k = k b j(i) A kj(i) for i =1 m end for j(i) = b j(i) for i =1 m for k =1 n A kj(i) =0 A j(i)k =0 end for A j(i)j(i) =1 j(i) =0 end for

83 m j(m) b j(m) n A nn b n n n b given i b given(n b given) v b given(n b given) n a(n,n) b(n) (n) Gauss 4 k 1 = k 2 = = k m = k 1 n =4 k =1 u 1 =0 f 5 =1 2 n =4 k =1 u 1 = u 5 =0 f 3 =1 3 n =4 k =1 u 1 = u 5 =0 u 3 =1 m 1 m 2 m 3 m 4 m 5 1.6: fortran fortran Makefile

84 74 1 Gauss FC = g77 F_OPT = -O2 OBJS = main.o stiff.o ex11-13.o gauss_ver4u.o bound1.o mis2.o TARGET = bound1 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f bound operation ver.1 with gauss elimination ver. 4u by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 5) parameter(i_hek = 3) parameter(i_type = 3) dimension a(n,n),b(n),(n) dimension i_b_given(n),v_b_given(n) nk = 4 ak = 1.D0 all stiff(a,n,nk,ak) if (i_hek.ge. 3) all hek_stiff(a,n) if (i_type.eq. 1) all ex11(,n_b_given,i_b_given,v_b_given,n) if (i_type.eq. 2) all ex12(,n_b_given,i_b_given,v_b_given,n) if (i_type.eq. 3) all ex13(,n_b_given,i_b_given,v_b_given,n) all bound1(a,,n,n_b_given,i_b_given,v_b_given) if (i_hek.ge. 2) all hek_matrix2(a,,n) all gauss_ver4u(a,,n) if (i_hek.ge. 1) all hek_solution(,n) stop end

85 stiff.f ####################################################################### subroutine stiff(a,n,nk,ak) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*) do i = 1,n do j = 1,n a(j,i) = 0.D0 do i = 1,nk a(i,i ) = a(i,i ) + ak a(i,i+1) = a(i,i+1) - ak a(i+1,i ) = a(i+1,i ) - ak a(i+1,i+1) = a(i+1,i+1) + ak return end ex11-13.f ####################################################################### subroutine ex11(,n_b_given,i_b_given,v_b_given,n) ####################################################################### impliit real*8(a-h,o-z) dimension (*),i_b_given(*),v_b_given(*) do i = 1,n (i) = 0.D0 (n) = 1.D0 n_b_given = 1 i_b_given(1) = 1 v_b_given(1) = 0.D0 return end

86 76 1 Gauss ####################################################################### subroutine ex12(,n_b_given,i_b_given,v_b_given,n) ####################################################################### impliit real*8(a-h,o-z) dimension (*),i_b_given(*),v_b_given(*) do i = 1,n (i) = 0.D0 (3) = 1.D0 n_b_given = 2 i_b_given(1) = 1 i_b_given(2) = 5 v_b_given(1) = 0.D0 v_b_given(2) = 0.D0 return end ####################################################################### subroutine ex13(,n_b_given,i_b_given,v_b_given,n) ####################################################################### impliit real*8(a-h,o-z) dimension (*),i_b_given(*),v_b_given(*) do i = 1,n (i) = 0.D0 n_b_given = 3 i_b_given(1) = 1 i_b_given(2) = 3 i_b_given(3) = 5 v_b_given(1) = 0.D0 v_b_given(2) = 1.D0 v_b_given(3) = 0.D0 return end gauss ver4u.f gauss ver4

87 bound1.f ####################################################################### subroutine bound1(a,,n,n_b_given,i_b_given,v_b_given) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*),i_b_given(*),v_b_given(*) do i = 1,n_b_given do k = 1,n (k) = (k) - v_b_given(i) * a(k,i_b_given(i)) do i = 1,n_b_given (i_b_given(i)) = v_b_given(i) do i = 1,n_b_given do k = 1,n a(k,i_b_given(i)) = 0.D0 a(i_b_given(i),k) = 0.D0 a(i_b_given(i),i_b_given(i)) = 1.D0 return end mis2.f ####################################################################### subroutine hek_stiff(a,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*) write(*,*) matrix A do i = 1,n write(*,1000) (a(i,j),j=1,n) 1000 format(8e10.3)

88 78 1 Gauss return end ####################################################################### subroutine hek_matrix2(a,,n) ####################################################################### impliit real*8(a-h,o-z) dimension a(n,*),(*) write(*,*) matrix A do i = 1,n write(*,1000) (a(i,j),j=1,n) write(*,*) vetor write(*,1000) ((i),i=1,n) 1000 format(8e10.3) return end ####################################################################### subroutine hek_solution(b,n) ####################################################################### impliit real*8(a-h,o-z) dimension b(*) write(*,*) solution vetor b write(*,1000) (b(i),i=1,n) 1000 format(8e10.3) return end C C Makefile CC = g

89 C_OPT = -O2 OBJS = main.o stiff.o ex11-13.o gauss_ver4u.o bound1.o mis2.o TARGET = bound1 $(TARGET):$(OBJS) $(CC) $(C_OPT) -o $(TARGET) $(OBJS)..o: $(CC) - $(C_OPT) $< main. /* bound operation ver.1 with gauss elimination ver. 4u by WATANABE Hiroshi (2000 July) */ #inlude <stdio.h> #define N 5 #define I_CHECK 3 #define I_TYPE 3 void stiff(double *a, int n, int nk, double ak); void hek_stiff(double *a, int n); void ex11(double *, int *n_b_given, int *i_b_given, double *v_b_given, int n); void ex12(double *, int *n_b_given, int *i_b_given, double *v_b_given, int n); void ex13(double *, int *n_b_given, int *i_b_given, double *v_b_given, int n); void bound1(double *a, double *, int n, int n_b_given, int *i_b_given, double *v_b_given); void hek_matrix2(double *a, double *, int n); void gauss_ver4u(double *a, double *, int n); void hek_solution(double *b, int n); main() { double *a,*b,*,*v_b_given; int *i_b_given; double ak; int n,nk,n_b_given; n = N;

90 80 1 Gauss nk = 4; ak = 1.0; a = (double *)mallo(sizeof(double)*n*n); b = (double *)mallo(sizeof(double)*n); = (double *)mallo(sizeof(double)*n); v_b_given = (double *)mallo(sizeof(double)*n); i_b_given = (int *)mallo(sizeof(int)*n); stiff(a,n,nk,ak); if (I_CHECK >= 3) hek_stiff(a,n); if (I_TYPE == 1) ex11(,&n_b_given,i_b_given,v_b_given,n); if (I_TYPE == 2) ex12(,&n_b_given,i_b_given,v_b_given,n); if (I_TYPE == 3) ex13(,&n_b_given,i_b_given,v_b_given,n); bound1(a,,n,n_b_given,i_b_given,v_b_given); if (I_CHECK >= 2) hek_matrix2(a,,n); gauss_ver4u(a,,n); if (I_CHECK >= 1) hek_solution(,n); stiff. void stiff(double *arg_a, int n, int nk, double ak) { double *a; int i,j; a = arg_a -1-n; for (i = 1; i <= n; ++i){ for (j = 1; j <= n; ++j){ a[i+j*n] = 0.0; for (i = 1; i <= nk ; ++i){ a[ i + i *n] = a[ i + i *n] + ak; a[ i +(i+1)*n] = a[ i +(i+1)*n] - ak; a[(i+1)+ i *n] = a[(i+1)+ i *n] - ak; a[(i+1)+(i+1)*n] = a[(i+1)+(i+1)*n] + ak;

91 ex void ex11(double *arg_, int *n_b_given, int *arg_i_b_given, double *arg_v_b_given,int n) { double *, *v_b_given; int *i_b_given; int i; = arg_ -1; v_b_given = arg_v_b_given -1; i_b_given = arg_i_b_given -1; for (i = 1; i <= n; ++i){ [i] = 0.0; [n] = 1.0; *n_b_given = 1; i_b_given[1] = 1; v_b_given[1] = 0.0; void ex12(double *arg_, int *n_b_given, int *arg_i_b_given, double *arg_v_b_given,int n) { double *, *v_b_given; int *i_b_given; int i; = arg_ -1; v_b_given = arg_v_b_given -1; i_b_given = arg_i_b_given -1; for (i = 1; i <= n; ++i){ [i] = 0.0; [3] = 1.0; *n_b_given = 2; i_b_given[1] = 1; i_b_given[2] = 5; v_b_given[1] = 0.0; v_b_given[2] = 0.0;

92 82 1 Gauss void ex13(double *arg_, int *n_b_given, int *arg_i_b_given, double *arg_v_b_given,int n) { double *, *v_b_given; int *i_b_given; int i; = arg_ -1; v_b_given = arg_v_b_given -1; i_b_given = arg_i_b_given -1; for (i = 1; i <= n; ++i){ [i] = 0.0; *n_b_given = 3; i_b_given[1] = 1; i_b_given[2] = 3; i_b_given[3] = 5; v_b_given[1] = 0.0; v_b_given[2] = 1.0; v_b_given[3] = 0.0; gauss ver4u. gauss ver4 bound1. void bound1(double *arg_a, double *arg_, int n, int n_b_given, int *arg_i_b_given, double *arg_v_b_given) { double *a,*,*v_b_given; int *i_b_given,i,k; a = arg_a -1-n; = arg_ -1; v_b_given = arg_v_b_given -1; i_b_given = arg_i_b_given -1; for (i = 1; i <= n_b_given; ++i){ for (k = 1; k <= n; ++k){ [k] = [k] - v_b_given[i] * a[k+i_b_given[i]*n];

93 for (i = 1; i <= n_b_given; ++i){ [i_b_given[i]] = v_b_given[i]; for (i = 1; i <= n_b_given; ++i){ for (k = 1; k <= n; ++k){ a[k+i_b_given[i]*n] = 0.0; a[i_b_given[i]+k*n] = 0.0; a[i_b_given[i]+i_b_given[i]*n]= 1.0; mis2. void hek_stiff(double *arg_a, int n) { double *a; int i,j; a = arg_a -1-n; printf("matrix A\n"); for (i = 1; i <= n; ++i){ for (j = 1; j <= n; ++j) printf("% 10.3E",a[i+j*n]); printf("\n"); void hek_matrix2(double *arg_a, double *arg_, int n) { double *a,*; int i,j; a = arg_a -1-n; = arg_ -1; printf("matrix A\n"); for (i = 1; i <= n; ++i){ for (j = 1; j <= n; ++j) printf("% 10.3E",a[i+j*n]); printf("\n"); printf("vetor \n"); for (i = 1; i <= n; ++i) printf("% 10.3E",[i]); printf("\n");

94 84 1 Gauss void hek_solution(double *arg_b, int n) { double *b; int i; b = arg_b-1; printf("solution vetor b\n"); for (i = 1; i <= n; ++i) printf("% 10.3E",b[i]); printf("\n"); (1) A 11 A 12 A 13 A 14 A 15 A 16 A 21 A 22 A 23 A 24 A 25 A 26 A 31 A 32 A 33 A 34 A 35 A 36 A 41 A 42 A 43 A 44 A 45 A 46 A 51 A 52 A 53 A 54 A 55 A 56 A 61 A 62 A 63 A 64 A 65 A 66 b 1 b 2 b 3 b 4 b 5 b 6 = (1.104) b 1 b A 22 0 A 24 A 25 0 b 2 2 A 21 A 23 A b 3 b = b 1 b 3 b 6 0 A 42 0 A 44 A 45 0 b 4 4 A 41 A 43 A 46 0 A 52 0 A 54 A 55 0 b 5 5 A 51 A 53 A b b 6 (1.105) b 1 = b 6 =0 2,4

95 m z j z (m) b jz(m) n b nonzero i b nonzero(n b nonzero) v b nonzero(n b nonzero) bound fortran fortran Makefile FC = g77 F_OPT = -O2 OBJS = main.o stiff.o ex21-23.o gauss_ver4u.o bound2.o mis2.o TARGET = bound2 $(TARGET):$(OBJS) $(FC) $(F_OPT) -o $(TARGET) $(OBJS).f.o: $(FC) - $(F_OPT) $< main.f bound operation ver.2 with gauss elimination ver.4 by WATANABE Hiroshi (2000 June) impliit real*8(a-h,o-z) parameter(n = 5) parameter(i_hek = 3) parameter(i_type = 3) dimension a(n,n),b(n),(n) dimension i_b_given(n),v_b_given(n)