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

(1)

計算情報数学

2008 (

第

13 回

)

C

言語および

LAPACK

(2)

Q & A

Euler 定数の多桁計算が大変積分の計算では精度を出すのが大変 http://phase.hpcc.jp/phase/mppack/long.pdf 参照 KNOPPIX がフリーズ乱数について MT (Mersenne Twister) の改良版の SFMT というのがある. ⇒ 高速で, 零超過状態 (調べてません) からの回復も高速 (だそう) 児玉氏 (神戸高専) について中西先生ではなく細川先生 (中西先生の師匠) の学生だった方です.

(3)

C

言語について

C

はコンパイラ言語

ファイルにプログラムを書いて

,

コンパイルして

から実行

全て関数で記述

main

という関数が必ず最初に呼ばれる

文法は

Asir

類似

(

実際には

Asir

が真似した

)

変数

(

データ

)

には型がある

(

整数

,

浮動小数

,

. . .

)

関数の戻り値

,

変数は全て型宣言が必要

(4)

コンパイル

,

実行

プログラムファイル名は全て拡張子

.c

をつける

例

:

afo.c

main() {

printf("afo\n");

}

コンパイル

,

実行

cc afo.c

./a.out

⇒

afo

が表示される

(5)

プログラム例

授業用ページに掲載 genb.c : M × N 行列データの生成コンパイル : cc -o genb genb.c 実行 : ./genb 1000 1000 a 1000 × 1000 行列 (要素はランダム) を生成して a に書く実行 : ./genb 1000 1 b 1000 × 1 行列 (要素はランダム) を生成して b に書く solveb.c 線形方程式系の LAPACK による求解コンパイル :

(6)

プログラム例

(

つづき

)

checkb.c : 解のチェック

コンパイル :

cc -o checkb checkb.c -llapack -lblas -lg2c

実行 : ./checkb 1000 a x b ||ax − b||₂ を表示 mylu.c 線形方程式系の求解 (非ブロック版; その 1) コンパイル : cc -o mylu mylu.c 実行 : ./mylu 1000 a b x mylu1.c 線形方程式系の求解 (非ブロック版; その 2) コンパイル : cc -o mylu1 mylu1.c 実行 : ./mylu1 1000 a b x

(7)

各プログラムについて

genb m n

が生成するファイル

倍精度浮動小数

(8

バイト

)

が

mn

個並んでいる

solveb

などにおける行列の読み出し

genb

の生成したファイルを

,

列優先で読み込む

実行時間計測

time solveb ...

とすれば表示される

レポートのヒント

プログラムの解読

他のソフトから行列データを生成するプログラ

ムの作成

(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’’); }

(9)

計算したいのは

, LU

分解と

,

それによる求解

LU

分解

,

求解

lapack_linsolve(a,b,n)

内部で

,

dgetrf

(LU

分解

),

dgetrs

(

求解

)

を呼

んでいる

問題

:

行列データ

(

a

),

ベクトルデータ

(

b

)

をどう準

備するか

⇒

ここでは

,

ファイルから読むことにする

(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"

(11)

各部分の説明

変数宣言

double a,b;

int n,nn;

*

がついている変数はポインタを表す

⇒

ここでは

,

配列

(

ベクトル

,

行列

)

を保持すると

思えばよい

a = read matrix(argv[1],nn)

ファイル

argv[1]

に書かれたデータを読んで

,

行列

(

を一次元配列にしたもの

)

を返す

write matrix(argv[3],b,n)

解

b

を

,

argv[3]

で指定されるファイルに書き

出す

.

(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);

これらが呼び出される前に宣言する

⇒

呼び出す関数の引数の型

,

返す値の型をコンパイ

(13)

ファイルから行列を読み込む方法

可読形式ファイルから読む

(Scilab

と同様

)

可読な数字を並べたファイルを作る

⇒

読むのが遅い

,

場合により不正確

バイナリファイルを読む

倍精度浮動小数

: 8

バイトで表現される

バイナリファイル

:

メモリ上の

8 バイトをそのま

まファイルに書く

⇒

読むのが圧倒的に速い

(

変換が不要

),

正確

どうやって用意するかが問題

(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;

}

(15)

バイナリ形式で読み込む例

ファイルは列優先で書かれているとする

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; }

(16)

行列を可読形式で書き出す例

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; }

(17)

行列をバイナリ形式で書き出す例

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; }

(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

に書く

(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; }

(20)

各部分の説明

ipiv = (int )malloc(nsizeof(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

(21)

コンパイル

,

ライブラリのリンク

ライブラリ関数を呼んでいない場合

cc -o mylu mylu.c

ライブラリ関数を呼んでいる場合

LAPACK

を呼んでいるなら

cc -o solve solve.c -llapack -lblas -lg2c

lxxx

は

libxxx.a

の略記

(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

を実行すれば

,

自動コンパイル

,

リンク

ソースが変更を認識してコンパイル

ソースが多数ある場合に特に便利

(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−2ui,j+1+ui+1,j+1)+(u_i−1,j−2ui,j+ui+1,j)

2h2

(24)

Crank-Nicolson

の公式

境界条件

_u

_i,0

_{= u}

_i,n+1

_{= 0}

を考慮して

A =       1 + r −r₂ −r₂ 1 + r . .. . .. . .. ₋r 2 −r₂ 1 + r       C =       1 − r r₂ r 2 1 − r . .. . .. . .. r 2 r 2 1 − r      

により

,

_AU

_j₊₁

_{= CU}

_j

(

_U

_j

₌

t

_(u

_1j

_{, . . . , u}

_mj

₎

)

すなわち

_U

_j₊₁

_{= A}

−1

_CU

j

,

U

0

=

t

(f (h), . . . , h(mh))

特徴

:

誤差は

_O(h

2

_{) + O(k}

2

₎

r

によらず安定

(

振動したりしない

)

⇒ h

,

k

を勝手にとれる

(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

があれば収容可能

)

三重対角行列の

分解を用いた求解

(26)

必要な処理

刻み幅の設定

引数として

,

_m

,

_n

,

_T

を指定

初期値ベクトルをファイルから読む

行列

_A

,

_C

の生成

:

malloc

により動的に確保

U

_j+1

= A

−1

CU

j

の繰り返し

A

−1

C

を計算してはいけない

(

理由は

?)

A

に対し

dgttrf_

を一回だけ実行

,

あとは

dgttrs_

の繰り返し

(

三重対角行列

-

ベクトル積

は自前で用意

)

結果の書き出し

可読形式で出力

_⇒

Scilab

でプロット

(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

:

解を納めたファイル

(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;

(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]; }

(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"); } }