私はいくつかの研究のために CUDA (GPGPU プログラミング) を使用していますが、新しいハードウェア アーキテクチャが原因で、固有の倍精度のパフォーマンスが単精度のパフォーマンスと比較して (24 倍!) 低下しています。2 つの uint を使用して 1 つの double を表すことにしました。このようにして、パフォーマンスに大きな影響を与えることなく DP 計算を実行できました。
たとえば、このメソッドを使用して double 10.12 を表現したいとします。
uint 実数 = 10、小数 = 12;
したがって; 「real.decimal」は double を視覚的に表します。
このタイプを「シングル」と呼びましょう。
与えられた: シングル a = 10.12, b = 20.24;
2 つの単数を乗算、除算、加算、および減算するための効率的なアルゴリズムは何でしょうか?
シングル c = a * b;
これが機能するためには、任意の DP 計算を実行したり、32 ビットを超えるデータ型を使用したりすることはできないことに注意してください。
@njuffaによって投稿された回答は、ほぼ間違いなく、この欲求を前進させる最も賢い方法です。どれだけ高速になるかはわかりませんが、デバイスで最高のスループットのエンジン（単精度浮動小数点エンジン）を利用するため、競合する代替エンジンの中で最も優れているように思われます。それにもかかわらず、あなたは固定小数点表現を求めていたので、私は先に進んで、この質問をアイデアのリポジトリとして使用して、とにかくこれを投稿すると思いました。
このコードには多くの注意点があります。
だから私はこれをフレームワークとして、娯楽のために提供します、そしておそらくあなたはそれを「真の」固定小数点形式で行うのが例えばSPエンジンを使うのに比べてどれほど遅いかを見ることができます。
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 10000
#define nTPB 512
#define SCALE 1
#define MUL_ERR_LIM (0.0000001 * SCALE * SCALE)
#define ADD_ERR_LIM (0.0000001 * SCALE)
#define cudaCheckErrors(msg) \
do { \
cudaError_t __err = cudaGetLastError(); \
if (__err != cudaSuccess) { \
fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
msg, cudaGetErrorString(__err), \
__FILE__, __LINE__); \
fprintf(stderr, "*** FAILED - ABORTING\n"); \
exit(1); \
} \
} while (0)
// fixed point number representation: 32 bit whole portion and 32 bit decimal
typedef union {
struct {
unsigned d; // decimal portion
unsigned w; // whole number portion
} part;
unsigned long val;
} fxp;
__device__ __host__ inline fxp fxp_add(fxp a, fxp b){
fxp temp;
temp.val = a.val + b.val;
return temp;
}
__device__ __host__ inline fxp fxp_sub(fxp a, fxp b){
fxp temp;
temp.val = a.val - b.val;
return temp;
}
__device__ __host__ inline unsigned fxp_whole(fxp a){
return a.part.w;
}
__device__ __host__ inline unsigned fxp_decimal(fxp a){
return a.part.d;
}
__device__ __host__ inline void fxp_set_whole(fxp *a, unsigned val){
a->part.w = val;
}
__device__ __host__ inline void fxp_set_decimal(fxp *a, unsigned val){
a->part.d = val;
}
__device__ __host__ inline fxp float_to_fxp(float val){
fxp temp;
temp.part.w = (unsigned) truncf(val);
temp.part.d = (unsigned) rintf((val - truncf(val)) * 0x0000000100000000ul);
return temp;
}
__device__ __host__ inline float fxp_to_float(fxp val){
return val.part.w + (val.part.d/(float)0x0000000100000000ul);
}
__device__ __host__ inline fxp double_to_fxp(double val){
fxp temp;
temp.part.w = (unsigned) trunc(val);
temp.part.d = (unsigned) rint((val - trunc(val)) * 0x0000000100000000ul);
return temp;
}
__device__ __host__ inline double fxp_to_double(fxp val){
return val.part.w + (val.part.d/(double)0x0000000100000000ul);
}
__device__ __host__ fxp fxp_mul(fxp a, fxp b){
fxp temp;
unsigned long ltemp = ((unsigned long)a.part.w * (unsigned long)b.part.d) + ((unsigned long)a.part.d * (unsigned long)b.part.w);
unsigned utemp = (unsigned) (ltemp & 0x0FFFFFFFFul);
temp.part.w = (unsigned) (ltemp >> 32);
temp.part.d = (unsigned) (((unsigned long)a.part.d * (unsigned long)b.part.d) >> 32) + utemp;
temp.part.w += (a.part.w * b.part.w) + ((temp.part.d < utemp) ? 1:0);
return temp;
}
__global__ void fxp_test(float *a, float *b, float *s, float *p, double *da, double *db, double *ds, double *dp, unsigned n){
unsigned idx = threadIdx.x + (blockDim.x * blockIdx.x);
if (idx < n){
float la = a[idx];
float lb = b[idx];
double lda = da[idx];
double ldb = db[idx];
s[idx] = fxp_to_float(fxp_add(float_to_fxp(la), float_to_fxp(lb)));
p[idx] = fxp_to_float(fxp_mul(float_to_fxp(la), float_to_fxp(lb)));
ds[idx] = fxp_to_double(fxp_add(double_to_fxp(lda), double_to_fxp(ldb)));
dp[idx] = fxp_to_double(fxp_mul(double_to_fxp(lda), double_to_fxp(ldb)));
}
}
int main(){
fxp a,b,c;
float x,y,z, xp, yp, zp;
double dx, dy, dz, dxp, dyp, dzp;
if (sizeof(unsigned) != 4) {printf("unsigned type size error: %d\n", sizeof(unsigned)); return 1;}
if (sizeof(unsigned long) != 8) {printf("unsigned long type error: %d\n", sizeof(unsigned long)); return 1;}
// test host side
x = 76.705116;
y = 1.891480;
a = float_to_fxp(x);
b = float_to_fxp(y);
// test conversions
xp = fxp_to_float(a);
yp = fxp_to_float(b);
printf("con: a = %f, a should be %f, b = %f, b should be %f\n", xp, x, yp, y);
// test multiply
c = fxp_mul(a, b);
z = x*y;
zp = fxp_to_float(c);
printf("mul: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);
//test add
c = fxp_add(a, b);
z = x+y;
zp = fxp_to_float(c);
printf("add: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);
//test subtract
c = fxp_sub(a, b);
z = x-y;
zp = fxp_to_float(c);
printf("sub: a = %f, b = %f, c = %f, c should be %f\n", x, y, zp, z);
// now test doubles
dx = 6.7877;
dy = 5.2444;
a = double_to_fxp(dx);
b = double_to_fxp(dy);
// test conversions
dxp = fxp_to_double(a);
dyp = fxp_to_double(b);
printf("dbl con: a = %f, a should be %f, b = %f, b should be %f\n", dxp, dx, dyp, dy);
// test multiply
c = fxp_mul(a, b);
dz = dx*dy;
dzp = fxp_to_double(c);
printf("double mul: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);
//test add
c = fxp_add(a, b);
dz = dx+dy;
dzp = fxp_to_double(c);
printf("double add: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);
//test subtract
c = fxp_sub(a, b);
dz = dx-dy;
dzp = fxp_to_double(c);
printf("double sub: a = %f, b = %f, c = %f, c should be %f\n", dx, dy, dzp, dz);
// test device side
float *h_a, *d_a, *h_b, *d_b, *h_s, *d_s, *h_p, *d_p;
double *h_da, *d_da, *h_db, *d_db, *h_ds, *d_ds, *h_dp, *d_dp;
if ((h_a=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_b=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_s=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_p=(float *)malloc(N*sizeof(float))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_da=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_db=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_ds=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
if ((h_dp=(double *)malloc(N*sizeof(double))) == 0) {printf("malloc fail\n"); return 1;}
cudaMalloc((void **)&d_a, N*sizeof(float));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_b, N*sizeof(float));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_s, N*sizeof(float));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_p, N*sizeof(float));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_da, N*sizeof(double));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_db, N*sizeof(double));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_ds, N*sizeof(double));
cudaCheckErrors("cudamalloc fail");
cudaMalloc((void **)&d_dp, N*sizeof(double));
cudaCheckErrors("cudamalloc fail");
for (unsigned i = 0; i < N; i++){
h_a[i] = (float)drand48() * SCALE;
h_b[i] = (float)drand48() * SCALE;
h_da[i] = drand48()*SCALE;
h_db[i] = drand48()*SCALE;
}
cudaMemcpy(d_a, h_a, N*sizeof(float), cudaMemcpyHostToDevice);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(d_b, h_b, N*sizeof(float), cudaMemcpyHostToDevice);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(d_da, h_da, N*sizeof(double), cudaMemcpyHostToDevice);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(d_db, h_db, N*sizeof(double), cudaMemcpyHostToDevice);
cudaCheckErrors("cudamemcpy fail");
fxp_test<<<(N+nTPB-1)/nTPB, nTPB>>>(d_a, d_b, d_s, d_p, d_da, d_db, d_ds, d_dp, N);
cudaCheckErrors("kernel fail");
cudaMemcpy(h_s, d_s, N*sizeof(float), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(h_p, d_p, N*sizeof(float), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(h_ds, d_ds, N*sizeof(double), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudamemcpy fail");
cudaMemcpy(h_dp, d_dp, N*sizeof(double), cudaMemcpyDeviceToHost);
cudaCheckErrors("cudamemcpy fail");
for (unsigned i=0; i<N; i++){
if (fabsf(h_s[i] - (h_a[i] + h_b[i])) > ADD_ERR_LIM) {printf("float add mismatch at %d, a: %f b: %f gpu: %f cpu: %f\n", i, h_a[i], h_b[i], h_s[i], (h_a[i] + h_b[i])); return 1;}
if (fabsf(h_p[i] - (h_a[i] * h_b[i])) > MUL_ERR_LIM) {printf("float mul mismatch at %d, a: %f b: %f gpu: %f cpu: %f\n", i, h_a[i], h_b[i], h_p[i], (h_a[i] * h_b[i])); return 1;}
if (fabs(h_ds[i] - (h_da[i] + h_db[i])) > ADD_ERR_LIM) {printf("double add mismatch at %d, a: %f b: %f gpu: %f cpu: %f\n", i, h_da[i], h_db[i], h_ds[i], (h_da[i] + h_db[i])); return 1;}
if (fabs(h_dp[i] - (h_da[i] * h_db[i])) > MUL_ERR_LIM) {printf("double mul mismatch at %d, a: %f b: %f gpu: %f cpu: %f\n", i, h_da[i], h_db[i], h_dp[i], (h_da[i] * h_db[i])); return 1;}
}
printf("GPU results match!\n");
return 0;
}
余談ですが、分子動力学研究コードAMBERはGPUで高速に実行され、SPとDPの計算の混合や、一部の作業に固定小数点表現を使用するSPFP形式など、速度を達成するためにいくつかの異なる算術モデルを実装します。と他の仕事のためのSP。でも詳細はわかりません。しかし、どうやらそれは良い効果をもたらすことができます。
倍精度を置き換えたい場合は、オペランドが「ヘッド」と「テール」と呼ばれる単精度数のペアであり、正規化要件を満たす「ダブルシングル」表現を使用することを検討できます |しっぽ| <= 0.5 * ulp (|頭|)。CUDA の float2 は、必要な float のペアを保持するのに適した型です。たとえば、重要度の低い x 成分は「尾」を表し、重要度の高い y 成分は「頭」を表します。
これは、倍精度とまったく同じ精度を提供しないことに注意してください。この表現の精度は、49 ビット (単精度の各仮数からの 24 ビットと末尾の符号ビットの 1 ビット) に制限されますが、倍精度の場合は 53 ビットです。演算は IEEE の丸めを提供しません。また、乗算コード内の一時的な量のオーバーフローにより、範囲が多少縮小される場合もあります。非正規量も正確に表現できない場合があります。
GPU のネイティブの倍精度演算を使用するよりもパフォーマンスが大幅に向上するとは思いません。レジスタの負荷が高くなり、各 "double single" 演算に必要な単精度演算の数がかなり多くなるためです。もちろん、両方のバリアントのパフォーマンスを試して測定することもできます。
以下は、「double single」形式での加算の例です。アルゴリズムのソースはコメントに記載されています。引用された研究は、1980 年代後半または 1990 年代初期にこの主題について研究を行った D. Priest の研究に基づいていると私は信じています。「double single」乗算の実際の例については、CUDA に同梱されているファイル math_functions.h の関数 __internal_dsmul() を見ることができます。
CUDA での 64 ビット整数演算に関する簡単なコメント (1 人のコメンターがこれを潜在的な代替手段として指摘したように)。現在の GPU ハードウェアは、ネイティブの 64 ビット整数演算のサポートが非常に限定されており、基本的には浮動小数点型との間の変換、およびロードとストアでの 64 ビット アドレスを使用したインデックス作成を提供しています。それ以外の場合、逆アセンブルされたコード (cuobjdump --dump-sass) を見ると簡単にわかるように、ネイティブの 32 ビット操作を介して 64 ビット整数操作が実装されます。
/* Based on: Andrew Thall, Extended-Precision Floating-Point Numbers for GPU
Computation. Retrieved from http://andrewthall.org/papers/df64_qf128.pdf
on 7/12/2011.
*/
__device__ float2 dsadd (float2 a, float2 b)
{
float2 z;
float t1, t2, t3, t4, t5, e;
t1 = __fadd_rn (a.y, b.y);
t2 = __fadd_rn (t1, -a.y);
t3 = __fadd_rn (__fadd_rn (a.y, t2 - t1), __fadd_rn (b.y, -t2));
t4 = __fadd_rn (a.x, b.x);
t2 = __fadd_rn (t4, -a.x);
t5 = __fadd_rn (__fadd_rn (a.x, t2 - t4), __fadd_rn (b.x, -t2));
t3 = __fadd_rn (t3, t4);
t4 = __fadd_rn (t1, t3);
t3 = __fadd_rn (t1 - t4, t3);
t3 = __fadd_rn (t3, t5);
z.y = e = __fadd_rn (t4, t3);
z.x = __fadd_rn (t4 - e, t3);
return z;
}