2009年3月13日金曜日

CUDA+Fortran

昨日から今日にかけてホワイトデーのお返しで菓子をあちこちからいっぱいもらったのでなんか嬉しい。もともと甘いもんスキーなので。先生方からもお返しもらってしまった。若干恐縮するなぁ。
研究室の工事があって、夕方から猛烈に接着剤のにおいがしだして、頭痛+アレルギー症状が劇化したので早々に帰ってきた。昔から溶剤なにおいとか防虫剤なにおいとか芳香剤なにおいが嫌いというか嗅ぐと気分悪くなったり頭痛くなる。なわけで今日はあんまり仕事してない。まぁ実は数値計算の新バージョンの計算コードの動作チェック(発狂しそうにめんどくさかった)も完了してあとはひたすら計算機のパワーに任せるだけの段階になったので気分的にのんびりしてたりする。

そんなわけで昨日に引き続きCUDAで遊ぶ。
計算物理のぎょーかいでは速度命だったり昔々の偉人がつくったコードを使いまわしたり、線形代数の神ライブラリLapackを使うことが多いのでいまだにFortranが幅を利かせている。他の分野ではFortranなにそれオイシイの?だろうけどなぁ…。自分もメインはFortran90でありこれまで作ってきたのも大半がそれベース。なのでその一部をCUDAで高速化しようと思うと、FortranからCUDAコードを呼び出す方法をお勉強しとかなあかん。
FortranとCはわりかし互換性がよくて、もともと相互にFortran中でCの関数を呼ぶとかC中でFortranのサブルーチンを呼び出すとかができる。CUDAはCがベースなので、Fortranから呼び出すことも実は比較的容易っぽい。

CUDAに限らず、FortranとCを混在させるときの注意点
・FortranからCの関数を呼ぶときには、呼び出される側のCの関数にexternを付けておく、最後に_を付けておく。
・Fortranはそもそもサブルーチン等の引数が参照渡しなのでC側では引数をポインタとして扱わないといけない。
・Fortranは配列が1から始まる、Cでは0から始まる。Fortranでは列優先なので(a,b)の次は(a+1,b)、Cでは行優先なので(a,b)の次は(a,b+1)。

とりあえず昨日の配列の和のコードをベースに作る
Fortranで配列を準備→CUDA側に引数で渡す、GPUで計算→Fortranサイドで出力チェックてなことをやってみようと。環境はLinux(CentOS5)+ifort

CUDA側(testf.cu)



#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <cuda_runtime.h>


__global__ void wa ( int *array1, int *array2, int len)
{
int i;
for (i=0; i<len ;i++)
array1[i]= array1[i]+array2[i];
return;
}

extern "C" void cudafunc_ (int* arrayH1, int* arrayH2, int* Np)
{
int i;
int N;
int* arrayD1;
int* arrayD2;

size_t array_size;

N=*Np;

printf("inside the C code \n");
printf("input array \n");

printf("input of H1 \n");
for (i=0;i<N;i++)
printf("%d\n",arrayH1[i]);

printf("input of H2\n");
for (i=0;i<N;i++)
printf("%d\n",arrayH2[i]);


array_size = sizeof(int) * N;

cudaMalloc( (void **) &arrayD1, array_size);
cudaMalloc((void **) &arrayD2, array_size);

cudaMemcpy(arrayD1,arrayH1,array_size,cudaMemcpyHostToDevice);
cudaMemcpy(arrayD2,arrayH2,array_size,cudaMemcpyHostToDevice);

wa <<<dim3(1,1),dim3(1,1,1)>>>(arrayD1,arrayD2,N);

cudaMemcpy(arrayH1,arrayD1,array_size,cudaMemcpyDeviceToHost);

printf("output\n");
for (i=0;i<N;i++)
printf("%d\n",arrayH1[i]);


return ;

}




Fortran側(testcuda.f90)



program main
!test for CUDA+Fortran

!the length of array
integer,parameter :: N=8
integer :: arrayH1(1:N)
integer :: arrayH2(1:N)
integer :: i


!setting the input array
do i=1,N
arrayH1(i)=i
arrayH2(i)=i*2
end do

call cudafunc(arrayH1,arrayH2,N)

print*,"output in fortran code"

do i=1,N
print*,arrayH1(i)
end do

end program main


Makefileの例

all: CUDAfortran

# Define Fortran compiler
FC= ifort

CUDAfortran: testcuda.f90 testf.o
$(FC) -o CUDAfortran testcuda.f90 testf.o -L/usr/local/cuda/lib -lcudart

testf.o: testf.cu
nvcc -c testf.cu

clean:
rm testf.o CUDAfortran


うむ、なんとか動いた。これであれやこれやをGPUで計算できるわけだがなんと10秒ルールとやらがあるらしい。GPU上のジョブの時間制限10秒ってマジっすか?10秒で計算終わってもデータ転送に時間かかってるっぽいしなんだか用途の制限多いんじゃ…?まぁ面白そうではあるので引き続きいろいろ調べると思う。

0 コメント: