2008 ( 13 ) C LAPACK 2008 ( 13 )C LAPACK p. 1

Size: px

Start display at page:

Download "2008 ( 13 ) C LAPACK 2008 ( 13 )C LAPACK p. 1"

あかりとべ
5 years ago
Views:

1 2008 ( 13 ) C LAPACK LAPACK p. 1

2 Q & A Euler KNOPPIX MT (Mersenne Twister) SFMT., ( ) ( ) ( ) ( ). LAPACK p. 2

3 C C, main Asir ( Asir ) ( ) (,,...), LAPACK p. 3

4 ,.c : afo.c } main() { printf("afo\n");, cc afo.c./a.out afo LAPACK p. 4

5 genb.c : M N : cc -o genb genb.c :./genb a ( ) a :./genb b ( ) b solveb.c LAPACK : cc -o solveb solveb.c -llapack -lblas -lg2c :./solveb 1000 a b x ax = b (b ), x LAPACK p. 5

6 ( ) checkb.c : : cc -o checkb checkb.c -llapack -lblas -lg2c :./checkb 1000 a x b ax b 2 mylu.c ( ; 1) : cc -o mylu mylu.c :./mylu 1000 a b x mylu1.c ( ; 2) : cc -o mylu1 mylu1.c :./mylu a b x LAPACK p. 6

7 genb m n (8 ) mn solveb genb, time solveb... LAPACK p. 7

8 main ( ) #include <stdio.h> #include <stdlib.h> #include <f2c.h> /* */ main(int argc,char **argv) { double *a,*b; } int n,nn; n = atoi(argv[1]); nn = n*n; a = read_matrix(argv[2],nn); b = read_matrix(argv[3],n); ret = lapack_linsolve(a,b,n); if ( ret == 0 ) { write_matrix(argv[4],b,n); } else fprintf(stderr, error\n ); LAPACK p. 8

9 , LU, LU, lapack_linsolve(a,b,n), dgetrf (LU ), dgetrs ( ) : (a), (b), (read_matrix). LAPACK p. 9

10 #include <stdio.h> #include <stdlib.h> #include <f2c.h> ( f2c.h LAPACK ) main main(int argc,char **argv) argc : ( ) argv : a.out 123 abc argc=3 argv[0]="a.out" argv[1]="123" argv[2]="abc" LAPACK p. 10

11 double *a,*b; int n,nn; *, (, ) a = read matrix(argv[1],nn) argv[1], ( ) write matrix(argv[3],b,n) b, argv[3]. LAPACK p. 11

12 int write_matrix(char *fname,double *a, int len); double *read_matrix(char *fname,int len); int lapack_linsolve(double *m_array, double *b_array,int n);, LAPACK p. 12

13 (Scilab ), : 8 : 8 ( ), LAPACK p. 13

14 #include <stdio.h> #include <stdlib.h> #include <f2c.h> double *read_matrix(char *fname,int len){ double *a; FILE *fp; int i; } if (!(fp = fopen(fname,"r")) ) return 0; a = (double *) malloc(len*sizeof(double)); for ( i = 0; i < len; i++ ) fscanf(fp,"%lf",&a[i]); return a; LAPACK p. 14

15 #include <stdio.h> #include <stdlib.h> #include <f2c.h> double *read_matrix(char *fname,int len){ double *a; FILE *fp; } if (!(fp = fopen(fname,"rb")) ) return 0; a = (double *) malloc(len*sizeof(double)); fread(a,len,sizeof(double),fp); return a; LAPACK p. 15

16 #include <stdio.h> #include <stdlib.h> #include <f2c.h> int write_matrix(double *a,char *fname, int len){ FILE *fp; int i; } if (!(fp = fopen(fname,"wb")) ) return 0; for ( i = 0; i < len; i++ ) fprintf(fp,"%lf",a[i]); fclose(fp); return 1; LAPACK p. 16

17 #include <stdio.h> #include <stdlib.h> #include <f2c.h> int write_matrix(double *a,char *fname, int len){ FILE *fp; } if (!(fp = fopen(fname,"wb")) ) return 0; fwrite(a,len,sizeof(double),fp); fclose(fp); return 1; LAPACK p. 17

18 fp=fopen(fname,"r"), fp=fopen(fname,"w"), open fclose(fp) fscanf(fp,"%lf",&a[i]),, a[i] fread(a,len,sizeof(double),fp) len, a LAPACK p. 18

19 lapack linsolve int lapack_linsolve(double *a,double *b, int n) { int i,j; int nrhs,lda,ldb,info; int *ipiv; char trans; ipiv = (int *)malloc(n*sizeof(int)); trans = N ; nrhs = 1; lda = n; ldb = n; dgetrf_(&n,&n,a,&lda,ipiv,&info); if ( info ) return -1; dgetrs_(&trans,&n,&nrhs,a,&lda,ipiv, b,&ldb,&info); if ( info ) return -1; return 0; } LAPACK p. 19

20 ipiv = (int *)malloc(n*sizeof(int)) ipiv dgetrf (&n,&n,a,&lda,ipiv,&info) LU. a. dgetrs (&trans,&n,&nrhs,a,&lda,...) LU ( piv) ax=b. b. lda, ldb : n nrhs ( ) : 1 LAPACK p. 20

21 , cc -o mylu mylu.c LAPACK cc -o solve solve.c -llapack -lblas -lg2c lxxx libxxx.a, Makefile LAPACK p. 21

22 Makefile all: sol.c.o: $(CC) -c -g $< clean: $(RM) -rf *.o sol sol: sol.o $(CC) -o sol sol.o -llapack -lblas -lg2c make, LAPACK p. 22

23 : 1 u = u(x,t) ((x,t) [0, 1] [0,T]), u t = u xx, u(x, 0) = f(x), u(0,t) = 0 h = 1/(m + 1), k = T/(n + 1), r = k/h 2 u ij = u(ih,jk) (i = 0,...,n + 1, j = 0,...,m + 1) u t (ih,jk) u i,j+1 ui,j k, u xx (ih,jk) (u i 1,j+1 2u i,j+1 +u i+1,j+1 )+(u i 1,j 2u i,j +u i+1,j ) 2h 2 ( j, j + 1 u xx ) LAPACK p. 23

24 Crank-Nicolson u i,0 = u i,n+1 = r r 2 A = r r r C = r r 2 1 r r 2 r. 2 1 r r 2 r 2 1 r, AU j+1 = CU j (U j = t (u 1j,...,u mj )) U j+1 = A 1 CU j, U 0 = t (f(h),...,h(mh)) : O(h 2 ) + O(k 2 ) r ( ) h, k LAPACK p. 24

25 LAPACK LU dgttrf_(int *n,double *l,double *d,double *u,double *u2, int *ipiv,int *info) l : ( n-1), d : ( n), u : ( n-1), u2 : ( n-2; ) l, d, u l, d, u, u2 (, u2 ) LU dgttrs_(int *trans,...) LAPACK p. 25

26 , m, n, T A, C : malloc U j+1 = A 1 CU j A 1 C (?) A dgttrf_, dgttrs_ ( - ) Scilab LAPACK p. 26

27 heat1.c cc -o heat1 heat1.o -llapack -latlas -lg2c : heat1 n m T f sol n :, m : T : t [0,T] f : (genb ) sol : : heat.c LAPACK p. 27

28 , : A,,, 2 malloc,. (C ) al=(double *)malloc((nx-1)*sizeof(double)); ad=(double *)malloc((nx)*sizeof(double)); au=(double *)malloc((nx-1)*sizeof(double)); au2=(double *)malloc((nx-2)*sizeof(double)); for ( i = 0; i < nx-1; i++ ) { al[i] = au[i] = -r/2; ad[i] = 1+r; } ad[nx-1] = 1+r; LAPACK p. 28

29 - BLAS void mul_td(double *l,double *d,double *u, double *x,double *y,int n) { int i; y[0] = d[0]*x[0]+u[0]*x[1]; for ( i = 1; i < n-1; i++ ) { y[i] = l[i-1]*x[i-1]+d[i]*x[i] +u[i]*x[i+1]; } y[n-1] = l[n-2]*x[n-2]+d[n-1]*x[n-1]; } LAPACK p. 29

30 U j+1 = A 1 CU j, u ipiv = (int *)malloc(nx*sizeof(int)); dgttrf_(&nx,al,ad,au,au2,ipiv,&info); trans = N ; nrhs = 1; ldb = nx; out = fopen(argv[5],"wb"); for ( i = 0; i < nt; i++ ) { mul_td(cl,cd,cu,u,w,nx); dgttrs_(&trans,&nx,&nrhs,al,ad,au,au2, ipiv,w,&ldb,&info); memcpy(u,w,nx*sizeof(double)); if (!(i%10) ) { /* 10 1 */ for ( j = 0; j < nx; j++ ) fprintf(out,"%lf ",w[j]); fprintf(out,"\n"); } } LAPACK p. 30

31 Scilab plot, Scilab U j Scilab plot(x) x, x plot, clf : sol, 1000 t=read("sol",-1,1000) ; plot(x); (A ) LAPACK p. 31

32 Scilab f(x) Fourier, S 1 ( ), etc. LAPACK 8/12 ( ) LAPACK p. 32

file"a" file"b" fp = fopen("a", "r"); while(fgets(line, BUFSIZ, fp)) {... fclose(fp); fp = fopen("b", "r"); while(fgets(line, BUFSIZ, fp)) {... fclose

filea fileb fp = fopen(a, r); while(fgets(line, BUFSIZ, fp)) {... fclose(fp); fp = fopen(b, r); while(fgets(line, BUFSIZ, fp)) {... fclose I117 9 2 School of Information Science, Japan Advanced Institute of Science and Technology file"a" file"b" fp = fopen("a", "r"); while(fgets(line, BUFSIZ, fp)) {... fclose(fp); fp = fopen("b", "r"); while(fgets(line,