4

CUDA でゼロ次 I0 の修正ベッセル関数を正確に計算する問題を扱っています。

論文によると、長い間、私は合理的なチェビシェフ近似を使用してきました。

JM Blair、「修正ベッセル関数 I_0(x) および I_1(x) の有理チェビシェフ近似」、Math. 、ｖｏｌ．28、n。126、pp. 581-583、1974 年 4 月。

これは、Matlab によって提供された結果と比較して、1e-29 のオーダーの平均誤差を示します。残念ながら、この一見高い精度は、私が取り組んでいる新しいアプリケーションには十分ではありません。

Matlab は、DE Amos によって開発された Fortran ルーチンを使用します。

デ・エイモス、「複雑な引数と非負の順序のベッセル関数のサブルーチン パッケージ」、サンディア国立研究所レポート、SAND85-1018、1985 年 5 月。

Amos、DE、「複雑な引数と非負の順序のベッセル関数のポータブル パッケージ」、Trans。算数。ソフトウェア、1986年。

netlib/amos Web サイトからダウンロードできます。

これらの Fortran ルーチンを C/C++ コードで使用するには、ライブラリ ファイルでコンパイルしてから C/C++ ラッパーを使用する方法があります (たとえば、netlib_wrappingを参照)。これらの Fortran ルーチンからデバイス関数を作成し、CUDA カーネルによって呼び出されるようにする方法があるかどうか疑問に思っています)。

問題の詳細

1 つは Matlab で、もう 1 つは CUDA で記述された 2 つのコードがあります。どちらも次の 3 つのステップで動作します。

1)修正ベッセル関数 I0 によるスケーリングとデータのゼロ パディング。

2) FFT ;

3)補間。

両方を「正確な」結果と比較しています。ステップ 3) の出力として、Matlab は 1e-10 % の相対二乗平均平方根誤差を与え、CUDA は 1e-2% であるため、その理由を調査し始めました。

2 つのコードの最初のステップの二乗平均平方根の差、つまり100*sqrt(sum(abs(U_Matlab_step_1-U_CUDA_step_1).^2))/sqrt(sum(abs(U_Matlab_step_1).^2))は0%(mean(mean(abs(U_Matlab-U_CUDA)))=6e-29代わりに) なので、良いと言えます。残念ながら、ステップ 2 に進むと、エラーが発生し2e-4%ます。最後に、CUDA のステップ 2) に Matlab のステップ 1) の出力を入力すると、ステップ 2) の rms エラー1e-14%は修正されたベッセル関数の。

この議論の興味深い展開について

NVIDIA 開発者ゾーン フォーラムをご覧ください

4

2 に答える 2

4

これは、浮動小数点演算間の精度の違いに起因するのではないかと思います。

チェックアウトすることがいくつかあります

  1. Cuda 5は、計算の形式によりよく一致する可能性のあるいくつかの新しい三角関数を追加します。また、バージョン4以降のCUDA数学ライブラリにはいくつかのベッセル関数があると思いますが、これが正しいかどうか、またはそれらが問題にどの程度関連しているかはわかりません。
  2. テスト用のシリアルCPUバージョンを作成できますか？これにより、精度の問題が、数値の64ビット表現と80ビット表現の使用などの最適化によるものかどうかがわかります。最適化をオフにすると、コンピューターはほとんど80ビット表現を処理します（おそらくmatlabがこれを実行します）が、数学の最適化をオンにすると、コンパイラーは精度の低い64ビット表現を処理する場合があります。これは、x87とSSEの違いに関係しています。
  3. コンピューティング機能のハードウェアが異なれば、精度もわずかに異なります。たとえば、compute 2.0は、より正確で、最適化されたx86により近いFMAを実行します。
  4. Matlabを正しいと見なす物理的な理由はありますか？Matlabがオーバーシュートしているのに、アルゴリズムが結果をアンダーシュートしている可能性があります。この種の状況は、Matlabがグループ化していないときにCUDAが操作をグループ化する場合に発生する可能性があります。
  5. 必要に応じて、Matlabの結果を再作成する必要があります。出力をさまざまな丸めのトリックと照合することで、コードの各ステップを調整することができます。表を参照してください。

丸めテーブル

addition       | x + y        | __dadd_[rn|rz|ru|rd](x, y)
multiplication | x * y        | __dmul_[rn|rz|ru|rd](x, y)
Fused-Mult-Add | fma(x, y, z) | __fma_[rn|rz|ru|rd](x, y, z)
reciprocal     | 1.0 / x      | __drcp_[rn|rz|ru|rd](x)
division       | x / y        | __ddiv_[rn|rz|ru|rd](x, y)
square root    | sqrt(x)      | __dsqrt_[rn|rz|ru|rd](x)

mode | interpretation
rn   | round to nearest, ties to even
rz   | round towards zero
ru   | round towards +∞
rd   | round towards -∞

http://developer.download.nvidia.com/assets/cuda/files/NVIDIA-CUDA-Floating-Point.pdfから

于 2013-01-02T21:30:32.527 に答える
0

あなたの質問に答える入門的なテクニカル トークを見つけました。PDFへのリンクはこちらです。はい、可能ですが、前述のスクリプトを使用してレガシー fortran コードを CUDA C に変換することはできませんでした。おそらく開発者に直接連絡してください。

于 2013-01-02T15:19:03.523 に答える