cuda - gpu上的模运算

Question

我正在研究应该进行大量模块化计算的 GPU 算法。特别是，从长远来看，对有限域中的矩阵的各种操作
简化为原始操作，例如: (a*b - c*d) mod m 或 (a*b + c) mod m 其中 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 (%) 操作很慢

Answer 1

高端的 Fermi GPU(例如 GTX 580)可能会在这方面为您提供在出货卡中最好的性能。您可能希望所有 32 位操作数都是“无符号整数”类型以获得最佳性能，因为处理有符号除法和模数会产生一些额外的开销。

编译器为具有固定除数的除法和模数生成非常有效的代码我记得在 Fermi 和 Kepler 上通常大约有三到五个机器指令指令。您可以使用 cuobjdump --dump-sass 检查生成的 SASS(机器代码)。如果您只使用几个不同的除数，您也许可以使用带有常量除数的模板函数。

您应该看到在 Fermi 和 Kepler 中为具有变量除数的无符号 32 位操作生成了大约 16 个内联 SASS 指令。代码受到整数乘法吞吐量的限制，对于 Fermi 级 GPU 而言，与硬件解决方案相比具有竞争力。由于整数乘法吞吐量降低，目前出货的 Kepler 级 GPU 的性能有所降低。

[在澄清问题后添加:]

另一方面，无符号 64 位除法和带变量除数的模被称为 Fermi 和 Kepler 上大约 65 条指令的子程序。它们看起来接近最佳状态。在 Fermi 上，这仍然与硬件实现相当有竞争力(请注意，64 位整数除法在作为内置指令提供的 CPU 上并不是非常快)。下面是我在 NVIDIA 论坛上发布的一些代码，用于说明中描述的那种任务。它避免了昂贵的除法，但确实假设相当大批量的操作数共享相同的除数。它使用 double 算术，这在特斯拉级 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;