我正在研究GPU算法,它应该做很多模块化计算。特别是对有限域中矩阵的各种运算,从长远来看简化为原始运算,如:(a*b-c*d)modm或(a*b+c)modm,其中a、b、c和d是模m的残差,m是32位素数。
通过实验,我了解到该算法的性能主要受到慢模运算的限制,因为硬件中的GPU不支持整数模(%)和除法运算。
如果有人能告诉我如何使用CUDA实现高效的模块化计算,我将不胜感激?
为了了解这是如何在CUDA上实现的,我使用了以下代码片段:
__global__ void mod_kernel(unsigned *gout, const unsigned *gin) {
unsigned tid = threadIdx.x;
unsigned a = gin[tid], b = gin[tid * 2], m = gin[tid * 3];
typedef unsigned long long u64;
__syncthreads();
unsigned r = (unsigned)(((u64)a * (u64)b) % m);
__syncthreads();
gout[tid] = r;
}
这个代码不应该工作,我只是想看看模块化缩减是如何实现的在CUDA上实现。
当我用cuobjdump--dump sass分解它时(感谢njuffa的建议!),我看到以下内容:
/*0098*/ /*0xffffdc0450ee0000*/ BAR.RED.POPC RZ, RZ;
/*00a0*/ /*0x1c315c4350000000*/ IMUL.U32.U32.HI R5, R3, R7;
/*00a8*/ /*0x1c311c0350000000*/ IMUL.U32.U32 R4, R3, R7;
/*00b0*/ /*0xfc01dde428000000*/ MOV R7, RZ;
/*00b8*/ /*0xe001000750000000*/ CAL 0xf8;
/*00c0*/ /*0x00000007d0000000*/ BPT.DRAIN 0x0;
/*00c8*/ /*0xffffdc0450ee0000*/ BAR.RED.POPC RZ, RZ;
请注意,在对bar.red.popc的两个调用之间,有一个对0xf8过程的调用,该过程实现了一些复杂的算法(大约50条指令甚至更多)。不奇怪的是,mod(%)操作是缓慢的
前段时间,我在GPU上尝试了很多模块运算。在费米GPU上,您可以使用双精度运算来避免昂贵的div和mod运算。例如,模乘可以如下进行:
// fast truncation of double-precision to integers
#define CUMP_D2I_TRUNC (double)(3ll << 51)
// computes r = a + b subop c unsigned using extended precision
#define VADDx(r, a, b, c, subop)
asm volatile("vadd.u32.u32.u32." subop " %0, %1, %2, %3;" :
"=r"(r) : "r"(a) , "r"(b), "r"(c));
// computes a * b mod m; invk = (double)(1<<30) / m
__device__ __forceinline__
unsigned mul_m(unsigned a, unsigned b, volatile unsigned m,
volatile double invk) {
unsigned hi = __umulhi(a*2, b*2); // 3 flops
// 2 double instructions
double rf = __uint2double_rn(hi) * invk + CUMP_D2I_TRUNC;
unsigned r = (unsigned)__double2loint(rf);
r = a * b - r * m; // 2 flops
// can also be replaced by: VADDx(r, r, m, r, "min") // == umin(r, r + m);
if((int)r < 0)
r += m;
return r;
}
然而,这只适用于31位整数模(如果1位对您来说不是关键的)并且您还需要预先计算"invk"。这给出了我能实现的指令的绝对最小值,即:
SHL.W R2, R4, 0x1;
SHL.W R8, R6, 0x1;
IMUL.U32.U32 R4, R4, R6;
IMUL.U32.U32.HI R8, R2, R8;
I2F.F64.U32 R8, R8;
DFMA R2, R2, R8, R10;
IMAD.U32.U32 R4, -R12, R2, R4;
ISETP.GE.AND P0, pt, R4, RZ, pt;
@!P0 IADD R4, R12, R4;
关于算法的描述,你可以看看我的论文:gpu顾问。其他操作,如(xy-zw)mod m也在那里解释。
出于好奇,我比较了所得算法的性能使用模块乘法:
unsigned r = (unsigned)(((u64)a * (u64)b) % m);
针对mulm的优化版本。
具有默认%运算的模块算术:
low_deg: 11; high_deg: 2481; bits: 10227
nmods: 330; n_real_pts: 2482; npts: 2495
res time: 5755.357910 ms; mod_inv time: 0.907008 ms; interp time: 856.015015 ms; CRA time: 44.065857 ms
GPU time elapsed: 6659.405273 ms;
带mul_m:的模运算
low_deg: 11; high_deg: 2481; bits: 10227
nmods: 330; n_real_pts: 2482; npts: 2495
res time: 1100.124756 ms; mod_inv time: 0.192608 ms; interp time: 220.615143 ms; CRA time: 10.376352 ms
GPU time elapsed: 1334.742310 ms;
因此,平均速度大约快5倍。还要注意的是,如果您只是使用带有一组mul_mod操作的内核(如sappy示例)来评估raw算术性能,则可能看不到加速。但在具有控制逻辑、同步屏障等的实际应用中,加速是非常显著的。
高端费米GPU(例如GTX 580)可能会为您提供在运输卡中最好的性能。为了获得最佳性能,您希望所有32位操作数都是"unsigned int"类型,因为在处理有符号除法和模时会有一些额外的开销。
编译器生成了非常有效的带有固定除数的除法和模的代码。我记得它通常是关于费米和开普勒的三到五条机器指令。您可以使用cuobjdump--dump SASS检查生成的SASS(机器代码)。如果只使用几个不同的除数,则可以使用带常除数的模板函数。
您应该看到,在费米和开普勒之间,为具有可变除数的无符号32位运算生成了16条内联SASS指令。该代码受到整数乘法吞吐量的限制,并且对于费米类GPU而言,它与硬件解决方案具有竞争力。由于开普勒级GPU的整数倍吞吐量降低,目前发货的GPU性能有所下降。
[稍后在澄清问题后添加:]
另一方面,无符号64位除法和带可变除数的模被称为费米和开普勒上大约65条指令的子程序。它们看起来接近最优。在Fermi上,这仍然与硬件实现具有相当的竞争力(请注意,在作为内置指令提供64位整数除法的CPU上,64位整数分割并不完全是超快的)。以下是我在NVIDIA论坛上发布的一些代码,用于澄清中描述的任务。它避免了昂贵的除法,但确实假设相当大的操作数批共享相同的除数。它使用双精度算术,这在特斯拉级GPU上尤其快(与消费卡相反)。我只是粗略地测试了一下代码,你可能想在部署它之前更仔细地测试一下
// Let b, p, and A[i] be integers < 2^51
// Let N be a integer on the order of 10000
// for i from 1 to N
// A[i] <-- A[i] * b mod p
/*---- kernel arguments ----*/
unsigned long long *A;
double b, p; /* convert from unsigned long long to double before passing to kernel */
double oop; /* pass precomputed 1.0/p to kernel */
/*---- code inside kernel -----*/
double a, q, h, l, rem;
const double int_cvt_magic = 6755399441055744.0; /* 2^52+2^51 */
a = (double)A[i];
/* approximate quotient and round it to the nearest integer */
q = __fma_rn (a * b, oop, int_cvt_magic);
q = q - int_cvt_magic;
/* back-multiply, representing p*q as a double-double h:l exactly */
h = p * q;
l = __fma_rn (p, q, -h);
/* remainder is double-width product a*b minus double-double h:l */
rem = __fma_rn (a, b, -h);
rem = rem - l;
/* remainder may be negative as quotient rounded; fix if necessary */
if (rem < 0.0) rem += p;
A[i] = (unsigned long long)rem;
有一些技巧可以有效地执行mod运算,但前提是只有m是基数2。
例如,x mod y==x&(y-1),其中y是2^n。执行逐位操作是最快的。
否则,可能是一个查找表?下面是关于高效模实现的讨论链接。你可能需要自己实现它才能最大限度地利用它
mod 的有效计算