CUDAと浮動小数点の誤差に関する実験メモ
マルチポスト元→http://exth.net/~tgbt/wordpress/2009/08/14/2433/
「CUDAはIEEE754対応だよ！計算精度はCPUと変わらないよ！！！」みたいなことをこれまでに何度か聞いていたような気がするので，実際のところどうなのかを試してみた．これが直接的に論文のネタになるとは思っていないので，公開しておきます．
対象問題
誤差が出る問題として何を使うかちょっと悩んだんだけど，手元の数値計算の教科書（新版 入門数値計算）の最初に掲載されていた「ニュートン・ラフソン法」を実装してみることにした．
もはや授業を受けたときの頃は全然覚えていないんだけど， f(x) = x^3 -2x -5 = 0の実根alphaを求める という問題．「組み立て除法」ってのを使う．
ソースコード
難しいことは考えずに，ガツッと決め撃ちした．
void kernel(TYPE *rs, TYPE *ep)
{
TYPE X3 = 1.0f;
TYPE X2 = 0.0f;
TYPE X1 = -2.0f;
TYPE X0 = -5.0f;
TYPE x3 = 0.0f;
TYPE x2 = 0.0f;
TYPE x1 = 0.0f;
TYPE x0 = 0.0f;
TYPE z0 = INIT;
TYPE z = z0;
TYPE a, b, c, d, e, f, g;
int i;
for(i=0; i<NLOOP; i++){
x3 = X3;
x2 = X2;
x1 = X1;
x0 = X0;
a = x3;
b = a * z;
c = x2 + b;
d = c * z;
e = x1 + d;
f = e * z;
g = x0 + f;
x3 = a;
x2 = c;
x1 = e;
z = g;
a = x3;
b = a * z0;
c = x2 + b;
d = c * z0;
e = x1 + d;
rs[i] = z0 - z/e;
ep[i] = (z/e)*(z/e);
z0 = rs[i];
z = z0;
}
}CPUもGPUも同じ計算カーネル．（GPUのカーネルにはもちろん__global__がつきますが．）
速度はどうでも良いし並列性もどうでも良い．Grid/Blockのパラメタも(1,1)(1,1,1)な設定でOK．
型TYPEはfloatとdoubleとを切り替えてコンパイルできるようにしている．初期値INITも入れ替えられるようにしておいた．
実行環境など
Core i7 + TeslaC1060．CentOS5.x x86_64．cudaのバージョンは2.21．-m 64で64bitのバイナリを吐かせている．doubleを使いたいので，もちろん-arch=sm_13．
実行結果
教科書では初期値2.0でやってたんだけど，収束の具合が見たいので適当に1.0にしてやってみた．ループ回数は特に考えも無しに適当に20回．
はじめにfloatでやった結果から．結果と誤差，それぞれ右下に行くほど繰り返し数が多い際の値．
結果：CPU 7.000000e+00 4.765517e+00 3.348703e+00 2.531600e+00 2.173916e+00 2.097884e+00 2.094558e+00 2.094552e+00 2.094551e+00 2.094552e+00 2.094551e+00 2.094552e+00 2.094551e+00 2.094552e+00 2.094551e+00 2.094552e+00 2.094551e+00 2.094552e+00 2.094551e+00 2.094552e+00
結果：GPU 7.000000e+00 4.765517e+00 3.348703e+00 2.531600e+00 2.173916e+00 2.097884e+00 2.094558e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00 2.094552e+00誤差：CPU 3.600000e+01 4.992913e+00 2.007363e+00 6.676576e-01 1.279377e-01 5.780894e-03 1.106205e-05 3.943923e-11 1.642639e-14 4.562891e-14 1.642639e-14 4.562891e-14 1.642639e-14 4.562891e-14 1.642639e-14 4.562891e-14 1.642639e-14 4.562891e-14 1.642639e-14 4.562891e-14
誤差：GPU 3.600000e+01 4.992913e+00 2.007363e+00 6.676576e-01 1.279377e-01 5.780877e-03 1.106177e-05 3.784591e-11 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
結果を見ると同じに見えるけど，実は途中で少しだけ値が変わってしまっている．
そして誤差を見るとGPUは途中から誤差が無くなってしまっている．なんだか変だ．
とはいえここまではfloat型．本番？はdouble型だよね．
結果：CPU 7.000000e+00 4.765517e+00 3.348703e+00 2.531600e+00 2.173916e+00 2.097884e+00 2.094558e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00
結果：GPU 7.000000e+00 4.765517e+00 3.348703e+00 2.531600e+00 2.173916e+00 2.097884e+00 2.094558e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00 2.094551e+00誤差：CPU 3.600000e+01 4.992913e+00 2.007363e+00 6.676575e-01 1.279377e-01 5.780895e-03 1.106208e-05 3.886632e-11 4.787812e-22 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32 2.532912e-32
誤差：GPU 3.600000e+01 4.992913e+00 2.007363e+00 6.676575e-01 1.279377e-01 5.780895e-03 1.106208e-05 3.886632e-11 4.787810e-22 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33 3.447355e-33
計算結果を見る限りではCPUとGPUの差は見受けられなかった．でも誤差を見ると違いがある．
なんだろうね，これ．
考察
というわけでして，CPUとGPUでは（double型で計算しても）実行結果に差が生じた．
原因はちょっとまだわかっていない．レジスタの違い？
そもそも，こうやって計算結果（誤差計算結果）の違いを見てしまうと，使い方が悪いだけなのか（GPUには問題はないのか）という疑問が生じるところ．実はGPUもちゃんとIEEE754には対応していて，CPUはさらに高精度な計算をしている？それはちょっとおかしいよねえ．
まぁとりあえずこんなところで．
数値計算に詳しい人に聞いてみるか，何も考えずにNVIDIAのForumに投げ込むか．でもForumに投げ込むにはちょっと問題（ソースコード）がシンプルじゃないしなあ．もっとプリミティブな例を作成してみようかしら？
