Используя два uint для представления двойника, а затем умножить?
Я работаю с CUDA (GPGPU programming) для некоторых исследований, и врожденная производительность двойной точности страдает по сравнению с производительностью одиночной точности (в 24 раза!), благодаря новой аппаратной архитектуре. Я решил попробовать использовать два uint для представления одного двойника. Таким образом, я мог выполнять вычисления DP без существенного снижения производительности.
Например, предположим, что мы хотим представить double 10.12 с помощью этого метода.
Uint real = 10, decimal = 12;
Следовательно; ' реальный.decimal ' визуально представляет нашего двойника.
Назовем этот тип "одиночным".
Дано: одиночный a = 10.12, b = 20.24;
Каков был бы эффективный алгоритм умножения, деления, сложения и вычитания двух одиночных чисел?
Одиночный c = a * b;
Пожалуйста, помните, что для того, чтобы это работало, невозможно выполнить какие-либо вычисления DP или использовать тип данных, превышающий 32 бита.
2 ответа:
Если вы хотите заменить двойную точность, вы можете рассмотреть использование "двойного одинарного" представления, где операндами являются пары чисел с одинарной точностью, обозначаемые как "head" и "tail", и которые удовлетворяют требованию нормализации, что |tail| Обратите внимание, что это не обеспечивает точно такую же точность, как двойная точность. Точность этого представления ограничена 49 битами (24 бита от каждой из мантисс одинарной точности плюс 1 бит знаковый бит хвоста) против 53 для двойной точности. Операции не будут обеспечивать округление IEEE, и диапазон также может быть несколько уменьшен из-за переполнения временных величин внутри кода умножения. Субнормальные величины также могут быть неточными. представимый.
Я не думаю, что производительность будет значительно лучше, чем при использовании собственных операций двойной точности вашего GPU, так как давление регистров будет выше, а количество операций с одной точностью, требуемых для каждой операции "двойной одинарной", довольно существенно. Конечно, вы можете просто попробовать измерить производительность для обоих вариантов.
Ниже приведен пример добавления в формате "double single". Источник алгоритма указан в комментарии. Я полагаю, что цитируемая работа, в свою очередь, основана на работе Д. Приста, который занимался исследованиями в этой области в конце 1980-х или начале 1990-х гг. для работающего примера "двойного одинарного" умножения вы могли бы взглянуть на функцию __internal_dsmul() в файле math_functions. h что корабли с CUDA.
Краткое замечание о 64-битных целочисленных операциях в CUDA (как отметил один комментатор в качестве потенциальной альтернативы). Нынешнее оборудование GPU имеет очень ограниченную поддержку собственных 64-битных процессоров. целочисленные операции, в основном обеспечивающие преобразование из и в типы с плавающей запятой, а также индексирование с использованием 64-разрядных адресов в загрузках и хранилищах. В противном случае 64-битные целочисленные операции реализуются через собственные 32-битные операции, что легко увидеть, посмотрев на дизассемблированный код (cuobjdump --dump-sass).
/* 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; }
Ответ, опубликованный @njuffa, почти наверняка является самым разумным путем для этого желания. Я не знаю, насколько быстро это будет, но поскольку он подключается к самому высокоточному движку на устройстве (движок с плавающей точкой с одной точностью), он, похоже, будет лучшим из любой конкурирующей альтернативы. Тем не менее вы попросили представление с фиксированной точкой, поэтому я подумал, что все равно отправлю это, просто используя этот вопрос как хранилище идей.
Этот код имеет много предостережения:
- не подвергался обширным испытаниям
- вероятно, медленнее, чем большинство альтернатив
- не выполняет арифметику со знаком (вы упомянули uint), поэтому вычитание сомнительно
- не реализовал divide
- нет проверки переполнения / недостаточного потока
- проигнорированный запрос абсолютно не использовать любое 64-битное представление или операцию. Я, конечно, не использовал 64-битные (DP) операции с плавающей запятой, но я использую 64-битные целочисленные операции. Мое обоснование заключается в том ,что (в некоторых случаях) я ломаю по сути 64-битное представление, предложенное в двух 32-битных количествах, делает кучу 32-битных операций, а затем повторная сборка бессмысленна, когда машина (компилятор) будет делать именно это в любом случае, если я запрошу 64-битные (целочисленные) операции.
- не оптимизирован. Вы, вероятно, могли бы сделать лучшую работу с PTX, но это может быть в категории "полировка дерьма".
Поэтому я предлагаю это в качестве основы, для развлечения, и, возможно, чтобы вы могли увидеть, насколько медленнее делать это в "истинный" формат фиксированной точки был бы по сравнению с использованием, например, двигателей 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 быстро работаетна графических процессорах и реализует несколько различных арифметических моделей для достижения своей скорости, включая смесь вычислений SP и DP, а также формат SPFP , который использует представление с фиксированной точкой для одной работы и SP для другой работы. Однако подробностей я не знаю. Но, по-видимому, это может быть сделано с хорошим эффектом.