logo标签列表

在CUDA中高效地将无符号值除以2的幂次方并向上取整

optimizationcudaroundinggpgpuinteger-division
3

我刚刚在阅读：

Efficiently dividing unsigned value by a power of two, rounding up

我想知道在CUDA中最快的方法是什么。当然，"快"指的是吞吐量（该问题还涉及到后续调用彼此依赖的情况）。

对于该问题中提到的lg()函数（除数的以2为底的对数），假设我们有：

template <typename T> __device__ int find_first_set(T x);
template <> __device__ int find_first_set<uint32_t>(uint32_t x) { return __ffs(x);   }
template <> __device__ int find_first_set<uint64_t>(uint64_t x) { return __ffsll(x); }

template <typename T> __device__ int lg(T x) { return find_first_set(x) - 1; }

编辑：由于我已经知道PTX中没有find-first-set，而且所有nVIDIA GPU的指令集都没有这个功能，因此让我们用以下内容替换lg()：

template <typename T> __df__ int population_count(T x);
template <> int population_count<uint32_t>(uint32_t x) { return __popc(x);   }
template <> int population_count<uint64_t>(uint64_t x) { return __popcll(x); }

template <typename T>
__device__ int lg_for_power_of_2(T x) { return population_count(x - 1); }

现在我们需要实施

template <typename T> T div_by_power_of_2_rounding_up(T p, T q);

针对 T = uint32_t 和 T = uint64_t，我们需要进行除法运算（p 为被除数，q 为除数）。

注意事项：

  • 与原问题一样，我们不能假设 p <= std::numeric_limits<T>::max() - q 或者 p > 0，因为这将导致各种有趣的情况崩溃 :-)
  • 0 不是 2 的幂，所以我们可以假设 q != 0。
  • 我知道解决方案可能会因为使用 32 位或 64 位而不同；我更关心前者，但也关心后者。
  • 让我们集中在 Maxwell 和 Pascal 芯片上。
是否允许额外的前置条件？如果不会导致溢出，可以编写(被除数+掩码)>>除数的log2。 - harold
@harold：不可以，就像原来的问题一样，这是不允许的，请查看编辑。 - einpoklum
6个回答
3

使用漏斗移位技术，一种可能的32位策略是进行33位移位（基本上）保存加法的进位，因此它可以在移位之前完成，例如：（未经测试）

unsigned sum = dividend + mask;
unsigned result = __funnelshift_r(sum, sum < mask, log_2_of_divisor);

编辑自 @einpoklum:

使用 @RobertCrovella 的程序测试，似乎可以正常工作。SM_61的测试内核PTX如下：

    .reg .pred      %p<2>;
    .reg .b32       %r<12>;


    ld.param.u32    %r5, [_Z4testjj_param_0];
    ld.param.u32    %r6, [_Z4testjj_param_1];
    neg.s32         %r7, %r6;
    and.b32         %r8, %r6, %r7;
    clz.b32         %r9, %r8;
    mov.u32         %r10, 31;
    sub.s32         %r4, %r10, %r9;
    add.s32         %r11, %r6, -1;
    add.s32         %r2, %r11, %r5;
    setp.lt.u32     %p1, %r2, %r11;
    selp.u32        %r3, 1, 0, %p1;
    // inline asm
    shf.r.wrap.b32 %r1, %r2, %r3, %r4;
    // inline asm
    st.global.u32   [r], %r1;
    ret;

而SASS是：

/*0008*/                   MOV R1, c[0x0][0x20];                 /* 0x4c98078000870001 */
/*0010*/                   MOV R0, c[0x0][0x144];                /* 0x4c98078005170000 */
/*0018*/                   IADD R2, RZ, -c[0x0][0x144];          /* 0x4c1100000517ff02 */
                                                                 /* 0x001c4c00fe4007f1 */
/*0028*/                   IADD32I R0, R0, -0x1;                 /* 0x1c0ffffffff70000 */
/*0030*/                   LOP.AND R2, R2, c[0x0][0x144];        /* 0x4c47000005170202 */
/*0038*/                   FLO.U32 R2, R2;                       /* 0x5c30000000270002 */
                                                                 /* 0x003fd800fe2007e6 */
/*0048*/                   IADD R5, R0, c[0x0][0x140];           /* 0x4c10000005070005 */
/*0050*/                   ISETP.LT.U32.AND P0, PT, R5, R0, PT;  /* 0x5b62038000070507 */
/*0058*/                   IADD32I R0, -R2, 0x1f;                /* 0x1d00000001f70200 */
                                                                 /* 0x001fc400fe2007f6 */
/*0068*/                   IADD32I R0, -R0, 0x1f;                /* 0x1d00000001f70000 */
/*0070*/                   SEL R6, RZ, 0x1, !P0;                 /* 0x38a004000017ff06 */
/*0078*/                   MOV32I R2, 0x0;                       /* 0x010000000007f002 */
                                                                 /* 0x0003c400fe4007e4 */
/*0088*/                   MOV32I R3, 0x0;                       /* 0x010000000007f003 */
/*0090*/                   SHF.R.W R0, R5, R0, R6;               /* 0x5cfc030000070500 */
/*0098*/                   STG.E [R2], R0;                       /* 0xeedc200000070200 */
                                                                 /* 0x001f8000ffe007ff */
/*00a8*/                   EXIT;                                 /* 0xe30000000007000f */
/*00b0*/                   BRA 0xb0;                             /* 0xe2400fffff87000f */
/*00b8*/                   NOP;                                  /* 0x50b0000000070f00 */
哦，太棒了！我完全忘记了“漏斗移位”，因为之前我注意到它时没有什么有趣的事情可以做。你能否链接到其他值得注意的用途？ - einpoklum
相同的策略也可以扩展到64位。不过似乎需要针对移位<=32位和> 32位的不同代码路径。 - tera
@tera，你的策略是什么？我没有看到一个好的扩展方式。 - harold
只需使用两个漏斗移位（一个用于高字，一个用于低字）。两个代码路径来自于 __funnelshift_r() 限制移位 <=31 位的事实。 - tera
@einpoklum 我没有任何链接，我所知道的 funnelshift 的唯一用途是旋转和大整数移位（这就是一个例子，因此可以使用它来操纵位流）。 顺便说一下，我就在你对面的 UvA。 - harold
那我们在现实生活中聊聊怎么样？我的个人资料链接到了我的GitHub资料，所以你可以很容易地找到我的姓名和电子邮件... - einpoklum
2
以下是一份针对CPU的高性能答案的改编：原回答。
template <typename T>
__device__ T div_by_power_of_2_rounding_up(T dividend, T divisor)
{
    auto log_2_of_divisor = lg(divisor);
    auto mask = divisor - 1;
    auto correction_for_rounding_up = ((dividend & mask) + mask) >> log_2_of_divisor;

    return (dividend >> log_2_of_divisor) + correction_for_rounding_up;
}

我想知道是否有更好的方法。
SM_61 的 SASS 代码（使用 @RobertCrovella 的测试内核）如下：
code for sm_61
        Function : test(unsigned int, unsigned int)
.headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                           /* 0x001fd400fe2007f6 */
/*0008*/                   MOV R1, c[0x0][0x20];           /* 0x4c98078000870001 */
/*0010*/                   IADD R0, RZ, -c[0x0][0x144];    /* 0x4c1100000517ff00 */
/*0018*/                   MOV R2, c[0x0][0x144];          /* 0x4c98078005170002 */
                                                           /* 0x003fc40007a007f2 */
/*0028*/                   LOP.AND R0, R0, c[0x0][0x144];  /* 0x4c47000005170000 */
/*0030*/                   FLO.U32 R3, R0;                 /* 0x5c30000000070003 */
/*0038*/                   IADD32I R0, R2, -0x1;           /* 0x1c0ffffffff70200 */
                                                           /* 0x001fc400fcc017f5 */
/*0048*/                   IADD32I R3, -R3, 0x1f;          /* 0x1d00000001f70303 */
/*0050*/                   LOP.AND R2, R0, c[0x0][0x140];  /* 0x4c47000005070002 */
/*0058*/                   IADD R2, R0, R2;                /* 0x5c10000000270002 */
                                                           /* 0x001fd000fe2007f1 */
/*0068*/                   IADD32I R0, -R3, 0x1f;          /* 0x1d00000001f70300 */
/*0070*/                   MOV R3, c[0x0][0x140];          /* 0x4c98078005070003 */
/*0078*/                   MOV32I R6, 0x0;                 /* 0x010000000007f006 */
                                                           /* 0x001fc400fc2407f1 */
/*0088*/                   SHR.U32 R4, R2, R0.reuse;       /* 0x5c28000000070204 */
/*0090*/                   SHR.U32 R5, R3, R0;             /* 0x5c28000000070305 */
/*0098*/                   MOV R2, R6;                     /* 0x5c98078000670002 */
                                                           /* 0x0003c400fe4007f4 */
/*00a8*/                   MOV32I R3, 0x0;                 /* 0x010000000007f003 */
/*00b0*/                   IADD R0, R4, R5;                /* 0x5c10000000570400 */
/*00b8*/                   STG.E [R2], R0;                 /* 0xeedc200000070200 */
                                                           /* 0x001f8000ffe007ff */
/*00c8*/                   EXIT;                           /* 0xe30000000007000f */
/*00d0*/                   BRA 0xd0;                       /* 0xe2400fffff87000f */
/*00d8*/                   NOP;                            /* 0x50b0000000070f00 */
                                                           /* 0x001f8000fc0007e0 */
/*00e8*/                   NOP;                            /* 0x50b0000000070f00 */
/*00f0*/                   NOP;                            /* 0x50b0000000070f00 */
/*00f8*/                   NOP;                            /* 0x50b0000000070f00 */

with FLO being the "find leading 1" instruction (thanks @tera). Anyway, those are lots of instructions, even if you ignore the loads from (what looks like) constant memory... the CPU function inspiring this one compiles into just:

    tzcnt   rax, rsi
    lea     rcx, [rdi - 1]
    shrx    rax, rcx, rax
    add     rax, 1
    test    rdi, rdi
    cmove   rax, rdi

（使用clang 3.9.0版本）
1
FLO.U32 -> 查找最高位1。 - tera
2
template <typename T> __device__ T div_by_power_of_2_rounding_up(T p, T q)
{
    return p==0 ? 0 : ((p - 1) >> lg(q)) + 1;
}

如果我数对了的话，比罗伯特之前的答案短一条指令（但请参见他的回应），或者和漏斗移位一样多。虽然有一个分支 - 不确定是否有区别（除非整个warp得到零p输入）：

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

    code for sm_61
        Function : _Z4testjj
    .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                                         /* 0x001fc000fda007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];                                 /* 0x4c98078000870001 */
        /*0010*/                   ISETP.EQ.AND P0, PT, RZ, c[0x0][0x140], PT;           /* 0x4b6503800507ff07 */
        /*0018*/         {         MOV R0, RZ;                                           /* 0x5c9807800ff70000 */
        /*0028*/               @P0 BRA 0x90;        }                                    /* 0x001fc800fec007fd */
                                                                                         /* 0xe24000000600000f */
        /*0030*/                   IADD R0, RZ, -c[0x0][0x144];                          /* 0x4c1100000517ff00 */
        /*0038*/                   LOP.AND R0, R0, c[0x0][0x144];                        /* 0x4c47000005170000 */
                                                                                         /* 0x003fc400ffa00711 */
        /*0048*/                   FLO.U32 R0, R0;                                       /* 0x5c30000000070000 */
        /*0050*/                   MOV R3, c[0x0][0x140];                                /* 0x4c98078005070003 */
        /*0058*/                   IADD32I R2, -R0, 0x1f;                                /* 0x1d00000001f70002 */
                                                                                         /* 0x001fd800fcc007f5 */
        /*0068*/                   IADD32I R0, R3, -0x1;                                 /* 0x1c0ffffffff70300 */
        /*0070*/                   IADD32I R2, -R2, 0x1f;                                /* 0x1d00000001f70202 */
        /*0078*/                   SHR.U32 R0, R0, R2;                                   /* 0x5c28000000270000 */
                                                                                         /* 0x001fc800fe2007f6 */
        /*0088*/                   IADD32I R0, R0, 0x1;                                  /* 0x1c00000000170000 */
        /*0090*/                   MOV32I R2, 0x0;                                       /* 0x010000000007f002 */
        /*0098*/                   MOV32I R3, 0x0;                                       /* 0x010000000007f003 */
                                                                                         /* 0x001ffc00ffe000f1 */
        /*00a8*/                   STG.E [R2], R0;                                       /* 0xeedc200000070200 */
        /*00b0*/                   EXIT;                                                 /* 0xe30000000007000f */
        /*00b8*/                   BRA 0xb8;                                             /* 0xe2400fffff87000f */
        ..........................

我相信通过使用PTX编写它仍然有可能从漏斗移位中削减一到两个指令（早间更新：如罗伯特在此期间已经证明），但我真的需要去睡觉了。
更新2：通过使用哈罗德的漏斗移位并将函数编写为PTX，可以实现上述操作。
_device__ uint32_t div_by_power_of_2_rounding_up(uint32_t p, uint32_t q)
{
  uint32_t ret;
  asm volatile("{\r\t"
               ".reg.u32        shift, mask, lo, hi;\n\t"
               "bfind.u32       shift, %2;\r\t"
               "sub.u32         mask, %2, 1;\r\t"
               "add.cc.u32      lo, %1, mask;\r\t"
               "addc.u32        hi, 0, 0;\r\t"
               "shf.r.wrap.b32  %0, lo, hi, shift;\n\t"
               "}"
                : "=r"(ret) : "r"(p), "r"(q));
  return ret;
}

只需使用与Robert已经用他简单的C代码实现相同的指令次数：

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

    code for sm_61
        Function : _Z4testjj
    .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                            /* 0x001fc000fec007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];                    /* 0x4c98078000870001 */
        /*0010*/                   MOV R0, c[0x0][0x144];                   /* 0x4c98078005170000 */
        /*0018*/         {         IADD32I R2, R0, -0x1;                    /* 0x1c0ffffffff70002 */
        /*0028*/                   FLO.U32 R0, c[0x0][0x144];        }      /* 0x001fc400fec00716 */
                                                                            /* 0x4c30000005170000 */
        /*0030*/                   IADD R5.CC, R2, c[0x0][0x140];           /* 0x4c10800005070205 */
        /*0038*/                   IADD.X R6, RZ, RZ;                       /* 0x5c1008000ff7ff06 */
                                                                            /* 0x003fc800fc8007f1 */
        /*0048*/                   MOV32I R2, 0x0;                          /* 0x010000000007f002 */
        /*0050*/                   MOV32I R3, 0x0;                          /* 0x010000000007f003 */
        /*0058*/                   SHF.R.W R0, R5, R0, R6;                  /* 0x5cfc030000070500 */
                                                                            /* 0x001ffc00ffe000f1 */
        /*0068*/                   STG.E [R2], R0;                          /* 0xeedc200000070200 */
        /*0070*/                   EXIT;                                    /* 0xe30000000007000f */
        /*0078*/                   BRA 0x78;                                /* 0xe2400fffff87000f */
        ..........................
2
借鉴@tera的炫酷答案：
 template <typename T> __device__ T pdqru(T p, T q)
{
    return bool(p) * (((p - 1) >> lg(q)) + 1);
}

11条指令（无分支，无预测），以在R0中获得结果：

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

        code for sm_61
                Function : _Z4testjj
        .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                   /* 0x001fc800fec007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];           /* 0x4c98078000870001 */
        /*0010*/                   IADD R0, RZ, -c[0x0][0x144];    /* 0x4c1100000517ff00 */
        /*0018*/                   LOP.AND R0, R0, c[0x0][0x144];  /* 0x4c47000005170000 */
                                                                   /* 0x003fc400ffa00711 */
        /*0028*/                   FLO.U32 R0, R0;                 /* 0x5c30000000070000 */
        /*0030*/                   MOV R5, c[0x0][0x140];          /* 0x4c98078005070005 */
        /*0038*/                   IADD32I R2, -R0, 0x1f;          /* 0x1d00000001f70002 */
                                                                   /* 0x001fd800fcc007f5 */
        /*0048*/                   IADD32I R0, R5, -0x1;           /* 0x1c0ffffffff70500 */
        /*0050*/                   IADD32I R2, -R2, 0x1f;          /* 0x1d00000001f70202 */
        /*0058*/                   SHR.U32 R0, R0, R2;             /* 0x5c28000000270000 */
                                                                   /* 0x001fd000fe2007f1 */
        /*0068*/                   IADD32I R0, R0, 0x1;            /* 0x1c00000000170000 */
        /*0070*/                   MOV32I R2, 0x0;                 /* 0x010000000007f002 */
        /*0078*/                   MOV32I R3, 0x0;                 /* 0x010000000007f003 */
                                                                   /* 0x001ffc001e2007f2 */
        /*0088*/                   ICMP.NE R0, R0, RZ, R5;         /* 0x5b4b02800ff70000 */
        /*0090*/                   STG.E [R2], R0;                 /* 0xeedc200000070200 */
        /*0098*/                   EXIT;                           /* 0xe30000000007000f */
                                                                   /* 0x001f8000fc0007ff */
        /*00a8*/                   BRA 0xa0;                       /* 0xe2400fffff07000f */
        /*00b0*/                   NOP;                            /* 0x50b0000000070f00 */
        /*00b8*/                   NOP;                            /* 0x50b0000000070f00 */
                ..........................

在研究上述SASS代码之后，似乎很明显这两条指令：
        /*0038*/                   IADD32I R2, -R0, 0x1f;          /* 0x1d00000001f70002 */
                                                                   /* 0x001fd800fcc007f5 */
        ...
        /*0050*/                   IADD32I R2, -R2, 0x1f;          /* 0x1d00000001f70202 */

“不应该真的必要。我没有一个精确的解释，但我的假设是因为SASS指令FLO.U32与__ffs()内置函数并不完全具有相同的语义，编译器显然在使用该内置函数时有一个习惯用法，它包装了正在执行工作的基本FLO指令。在C++源代码级别上如何解决这个问题并不明显，但我能够以一种方式使用bfindPTX指令，进一步减少指令计数，根据我的计算，将答案放入一个寄存器需要7条指令。”
$ cat t107.cu
#include <cstdio>
#include <cstdint>
__device__ unsigned r = 0;


static __device__ __inline__ uint32_t __my_bfind(uint32_t val){
  uint32_t ret;
  asm volatile("bfind.u32 %0, %1;" : "=r"(ret): "r"(val));
  return ret;}

template <typename T> __device__ T pdqru(T p, T q)
{
    return bool(p) * (((p - 1) >> (__my_bfind(q))) + 1);
}

__global__ void test(unsigned p, unsigned q){
#ifdef USE_DISPLAY
  unsigned q2 = 16;
  unsigned z = 0;
  unsigned l = 1U<<31;
  printf("result %u/%u = %u\n", p, q, pdqru(p, q));
  printf("result %u/%u = %u\n", p, q2, pdqru(p, q2));
  printf("result %u/%u = %u\n", p, z, pdqru(p, z));
  printf("result %u/%u = %u\n", z, q, pdqru(z, q));
  printf("result %u/%u = %u\n", l, q, pdqru(l, q));
  printf("result %u/%u = %u\n", q, l, pdqru(q, l));
  printf("result %u/%u = %u\n", l, l, pdqru(l, l));
  printf("result %u/%u = %u\n", q, q, pdqru(q, q));
#else
  r = pdqru(p, q);
#endif
}


int main(){
  unsigned h_r;
  test<<<1,1>>>(32767, 32);
  cudaMemcpyFromSymbol(&h_r, r, sizeof(unsigned));
  printf("result = %u\n", h_r);
}


$ nvcc -arch=sm_61 -o t107 t107.cu -std=c++11
$ cuobjdump -sass t107

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = <unknown>
host = linux
compile_size = 64bit

        code for sm_61

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

        code for sm_61
                Function : _Z4testjj
        .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                        /* 0x001c4400fe0007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];                /* 0x4c98078000870001 */
        /*0010*/         {         MOV32I R3, 0x0;                      /* 0x010000000007f003 */
        /*0018*/                   FLO.U32 R2, c[0x0][0x144];        }  /* 0x4c30000005170002 */
                                                                        /* 0x003fd800fec007f6 */
        /*0028*/                   MOV R5, c[0x0][0x140];               /* 0x4c98078005070005 */
        /*0030*/                   IADD32I R0, R5, -0x1;                /* 0x1c0ffffffff70500 */
        /*0038*/                   SHR.U32 R0, R0, R2;                  /* 0x5c28000000270000 */
                                                                        /* 0x001fc800fca007f1 */
        /*0048*/                   IADD32I R0, R0, 0x1;                 /* 0x1c00000000170000 */
        /*0050*/                   MOV32I R2, 0x0;                      /* 0x010000000007f002 */
        /*0058*/                   ICMP.NE R0, R0, RZ, R5;              /* 0x5b4b02800ff70000 */
                                                                        /* 0x001ffc00ffe000f1 */
        /*0068*/                   STG.E [R2], R0;                      /* 0xeedc200000070200 */
        /*0070*/                   EXIT;                                /* 0xe30000000007000f */
        /*0078*/                   BRA 0x78;                            /* 0xe2400fffff87000f */
                ..........................



Fatbin ptx code:
================
arch = sm_61
code version = [5,0]
producer = cuda
host = linux
compile_size = 64bit
compressed
$ nvcc -arch=sm_61 -o t107 t107.cu -std=c++11 -DUSE_DISPLAY
$ cuda-memcheck ./t107
========= CUDA-MEMCHECK
result 32767/32 = 1024
result 32767/16 = 2048
result 32767/0 = 1
result 0/32 = 0
result 2147483648/32 = 67108864
result 32/2147483648 = 1
result 2147483648/2147483648 = 1
result 32/32 = 1
result = 0
========= ERROR SUMMARY: 0 errors
$

我只展示了32位的例子。
我认为我可以证明，在上述内核SASS中，实际上只有6条指令在“工作”，其余的指令都是内核“开销”和/或将寄存器结果存储到全局内存所需的指令。显然，编译器正是通过这个函数生成这些指令。
        /*0018*/                   FLO.U32 R2, c[0x0][0x144];  // find bit set in q
                                                                        /*  */
        /*0028*/                   MOV R5, c[0x0][0x140];      // load p
        /*0030*/                   IADD32I R0, R5, -0x1;       // subtract 1 from p
        /*0038*/                   SHR.U32 R0, R0, R2;         // shift p right by q bit
                                                                        /*  */
        /*0048*/                   IADD32I R0, R0, 0x1;        // add 1 to result
        /*0050*/                   ...                                  /*  */
        /*0058*/                   ICMP.NE R0, R0, RZ, R5;     // account for p=0 case

然而，这将与我对其他案例的计数方式不一致（它们可能都应该减少1）。
FLO 不应该具有与 __ffs 相同的语义：__ffs 计算 尾部 零位，而不是 前导 零位；你可能想要使用 __clz... 在 PTX 中（对于我的版本），我们实际上有等效于 clz(q & ~p) 的函数，这本身就很奇怪。 - einpoklum
我从未说过FLO应该具有与__ffs相同的语义，而且我说过我的假设是它们是不同的。FLO是查找最高位1的意思，我正在寻找能够执行此操作的PTX指令，似乎是bfind。我假设__ffs()包含额外的不必要指令的原因是由于FLO和__ffs()之间的语义差异。即使这些指令看起来是不必要的，编译器似乎也没有剥离它们，我还没有确定确切的原因。 - Robert Crovella
1
一种可能的简单方法是：
$ cat t105.cu
#include <cstdio>

__device__ unsigned r = 0;

template <typename T>
__device__ T pdqru(T p, T q){

  T p1 = p +  (q-1);
  if (sizeof(T) == 8)
    q = __ffsll(q);
  else
    q = __ffs(q);
  return (p1<p)?((p>>(q-1))+1) :(p1 >> (q-1));
}

__global__ void test(unsigned p, unsigned q){
#ifdef USE_DISPLAY
  unsigned q2 = 16;
  unsigned z = 0;
  unsigned l = 1U<<31;
  printf("result %u/%u = %u\n", p, q, pdqru(p, q));
  printf("result %u/%u = %u\n", p, q2, pdqru(p, q2));
  printf("result %u/%u = %u\n", p, z, pdqru(p, z));
  printf("result %u/%u = %u\n", z, q, pdqru(z, q));
  printf("result %u/%u = %u\n", l, q, pdqru(l, q));
  printf("result %u/%u = %u\n", q, l, pdqru(q, l));
  printf("result %u/%u = %u\n", l, l, pdqru(l, l));
  printf("result %u/%u = %u\n", q, q, pdqru(q, q));
#else
  r = pdqru(p, q);
#endif
}


int main(){
  unsigned h_r;
  test<<<1,1>>>(32767, 32);
  cudaMemcpyFromSymbol(&h_r, r, sizeof(unsigned));
  printf("result = %u\n", h_r);
}


$ nvcc -arch=sm_61 -o t105 t105.cu
$ cuobjdump -sass ./t105

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = <unknown>
host = linux
compile_size = 64bit

        code for sm_61

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

        code for sm_61
                Function : _Z4testjj
        .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                         /* 0x001fc800fec007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];                 /* 0x4c98078000870001 */
        /*0010*/                   IADD R0, RZ, -c[0x0][0x144];          /* 0x4c1100000517ff00 */
        /*0018*/                   LOP.AND R0, R0, c[0x0][0x144];        /* 0x4c47000005170000 */
                                                                         /* 0x005fd401fe20003d */
        /*0028*/                   FLO.U32 R2, R0;                       /* 0x5c30000000070002 */
        /*0030*/                   MOV R0, c[0x0][0x144];                /* 0x4c98078005170000 */
        /*0038*/                   IADD32I R3, -R2, 0x1f;                /* 0x1d00000001f70203 */
                                                                         /* 0x001fd000fc2007f1 */
        /*0048*/                   IADD32I R0, R0, -0x1;                 /* 0x1c0ffffffff70000 */
        /*0050*/                   MOV R2, c[0x0][0x140];                /* 0x4c98078005070002 */
        /*0058*/                   IADD32I R4, -R3, 0x1f;                /* 0x1d00000001f70304 */
                                                                         /* 0x001fd800fe2007f6 */
        /*0068*/                   IADD R5, R0, c[0x0][0x140];           /* 0x4c10000005070005 */
        /*0070*/                   ISETP.LT.U32.AND P0, PT, R5, R0, PT;  /* 0x5b62038000070507 */
        /*0078*/                   SHR.U32 R0, R2, R4;                   /* 0x5c28000000470200 */
                                                                         /* 0x001fd000fc2007f1 */
        /*0088*/                   IADD32I R0, R0, 0x1;                  /* 0x1c00000000170000 */
        /*0090*/                   MOV32I R2, 0x0;                       /* 0x010000000007f002 */
        /*0098*/                   MOV32I R3, 0x0;                       /* 0x010000000007f003 */
                                                                         /* 0x001ffc001e2007f2 */
        /*00a8*/              @!P0 SHR.U32 R0, R5, R4;                   /* 0x5c28000000480500 */
        /*00b0*/                   STG.E [R2], R0;                       /* 0xeedc200000070200 */
        /*00b8*/                   EXIT;                                 /* 0xe30000000007000f */
                                                                         /* 0x001f8000fc0007ff */
        /*00c8*/                   BRA 0xc0;                             /* 0xe2400fffff07000f */
        /*00d0*/                   NOP;                                  /* 0x50b0000000070f00 */
        /*00d8*/                   NOP;                                  /* 0x50b0000000070f00 */
                                                                         /* 0x001f8000fc0007e0 */
        /*00e8*/                   NOP;                                  /* 0x50b0000000070f00 */
        /*00f0*/                   NOP;                                  /* 0x50b0000000070f00 */
        /*00f8*/                   NOP;                                  /* 0x50b0000000070f00 */
                ..........................



Fatbin ptx code:
================
arch = sm_61
code version = [5,0]
producer = cuda
host = linux
compile_size = 64bit
compressed
$ nvcc -arch=sm_61 -o t105 t105.cu -DUSE_DISPLAY
$ cuda-memcheck ./t105
========= CUDA-MEMCHECK
result 32767/32 = 1024
result 32767/16 = 2048
result 32767/0 = 2048
result 0/32 = 0
result 2147483648/32 = 67108864
result 32/2147483648 = 1
result 2147483648/2147483648 = 1
result 32/32 = 1
result = 0
========= ERROR SUMMARY: 0 errors
$

大约有14个SASS指令适用于32位情况，将答案存储在R0中。对于除以零的情况，它会产生虚假结果。

这个答案的等效汇编代码如下：

$ cat t106.cu
#include <cstdio>
#include <cstdint>
__device__ unsigned r = 0;


template <typename T> __device__ int find_first_set(T x);
template <> __device__ int find_first_set<uint32_t>(uint32_t x) { return __ffs(x);   }
template <> __device__ int find_first_set<uint64_t>(uint64_t x) { return __ffsll(x); }

template <typename T>  __device__ T lg(T x) { return find_first_set(x) - 1; }

template <typename T>
__device__ T pdqru(T dividend, T divisor)
{
    auto log_2_of_divisor = lg(divisor);
    auto mask = divisor - 1;
    auto correction_for_rounding_up = ((dividend & mask) + mask) >> log_2_of_divisor;

    return (dividend >> log_2_of_divisor) + correction_for_rounding_up;
}

__global__ void test(unsigned p, unsigned q){
#ifdef USE_DISPLAY
  unsigned q2 = 16;
  unsigned z = 0;
  unsigned l = 1U<<31;
  printf("result %u/%u = %u\n", p, q, pdqru(p, q));
  printf("result %u/%u = %u\n", p, q2, pdqru(p, q2));
  printf("result %u/%u = %u\n", p, z, pdqru(p, z));
  printf("result %u/%u = %u\n", z, q, pdqru(z, q));
  printf("result %u/%u = %u\n", l, q, pdqru(l, q));
  printf("result %u/%u = %u\n", q, l, pdqru(q, l));
  printf("result %u/%u = %u\n", l, l, pdqru(l, l));
  printf("result %u/%u = %u\n", q, q, pdqru(q, q));
#else
  r = pdqru(p, q);
#endif
}


int main(){
  unsigned h_r;
  test<<<1,1>>>(32767, 32);
  cudaMemcpyFromSymbol(&h_r, r, sizeof(unsigned));
  printf("result = %u\n", h_r);
}


$ nvcc -std=c++11  -arch=sm_61 -o t106 t106.cu
$ cuobjdump -sass t106

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = <unknown>
host = linux
compile_size = 64bit

        code for sm_61

Fatbin elf code:
================
arch = sm_61
code version = [1,7]
producer = cuda
host = linux
compile_size = 64bit

        code for sm_61
                Function : _Z4testjj
        .headerflags    @"EF_CUDA_SM61 EF_CUDA_PTX_SM(EF_CUDA_SM61)"
                                                                   /* 0x001fd400fe2007f6 */
        /*0008*/                   MOV R1, c[0x0][0x20];           /* 0x4c98078000870001 */
        /*0010*/                   IADD R0, RZ, -c[0x0][0x144];    /* 0x4c1100000517ff00 */
        /*0018*/                   MOV R2, c[0x0][0x144];          /* 0x4c98078005170002 */
                                                                   /* 0x003fc40007a007f2 */
        /*0028*/                   LOP.AND R0, R0, c[0x0][0x144];  /* 0x4c47000005170000 */
        /*0030*/                   FLO.U32 R3, R0;                 /* 0x5c30000000070003 */
        /*0038*/                   IADD32I R0, R2, -0x1;           /* 0x1c0ffffffff70200 */
                                                                   /* 0x001fc400fcc017f5 */
        /*0048*/                   IADD32I R3, -R3, 0x1f;          /* 0x1d00000001f70303 */
        /*0050*/                   LOP.AND R2, R0, c[0x0][0x140];  /* 0x4c47000005070002 */
        /*0058*/                   IADD R2, R0, R2;                /* 0x5c10000000270002 */
                                                                   /* 0x001fd000fe2007f1 */
        /*0068*/                   IADD32I R0, -R3, 0x1f;          /* 0x1d00000001f70300 */
        /*0070*/                   MOV R3, c[0x0][0x140];          /* 0x4c98078005070003 */
        /*0078*/                   MOV32I R6, 0x0;                 /* 0x010000000007f006 */
                                                                   /* 0x001fc400fc2407f1 */
        /*0088*/                   SHR.U32 R4, R2, R0.reuse;       /* 0x5c28000000070204 */
        /*0090*/                   SHR.U32 R5, R3, R0;             /* 0x5c28000000070305 */
        /*0098*/                   MOV R2, R6;                     /* 0x5c98078000670002 */
                                                                   /* 0x0003c400fe4007f4 */
        /*00a8*/                   MOV32I R3, 0x0;                 /* 0x010000000007f003 */
        /*00b0*/                   IADD R0, R4, R5;                /* 0x5c10000000570400 */
        /*00b8*/                   STG.E [R2], R0;                 /* 0xeedc200000070200 */
                                                                   /* 0x001f8000ffe007ff */
        /*00c8*/                   EXIT;                           /* 0xe30000000007000f */
        /*00d0*/                   BRA 0xd0;                       /* 0xe2400fffff87000f */
        /*00d8*/                   NOP;                            /* 0x50b0000000070f00 */
                                                                   /* 0x001f8000fc0007e0 */
        /*00e8*/                   NOP;                            /* 0x50b0000000070f00 */
        /*00f0*/                   NOP;                            /* 0x50b0000000070f00 */
        /*00f8*/                   NOP;                            /* 0x50b0000000070f00 */
                ..........................



Fatbin ptx code:
================
arch = sm_61
code version = [5,0]
producer = cuda
host = linux
compile_size = 64bit
compressed
$

根据我的计算，这似乎比原来多了1条指令。
这不是一个沿途编辑 - 这是问题的重点（以及原始问题）; 否则，确实，直接的方法相对容易想到，而且似乎没有过于有趣的替代方案。 - einpoklum
我已经修改了它以处理溢出。它添加了一些指令。 - Robert Crovella
SASS指令在32位情况下的计数基本上与您的实现相同。 - Robert Crovella
1

这里提供一种基于计数的替代方案。我只尝试了32位变体，并对其进行了详尽的参考实现测试。由于除数q是2的幂，我们可以轻松通过人口统计操作确定移位计数s。截断除法的余数t可以通过直接从除数q派生的简单掩码m来计算。

// For p in [0,0xffffffff], q = (1 << s) with s in [0,31], compute ceil(p/q)
__device__ uint32_t reference (uint32_t p, uint32_t q)
{
    uint32_t r = p / q;
    if ((q * r) < p) r++;
    return r;
}

// For p in [0,0xffffffff], q = (1 << s) with s in [0,31], compute ceil(p/q)
__device__ uint32_t solution (uint32_t p, uint32_t q)
{
    uint32_t r, s, t, m;
    m = q - 1;
    s = __popc (m);
    r = p >> s;
    t = p & m;
    if (t > 0) r++;
    return r;
}

无论 solution() 是否比之前发布的代码更快，可能取决于特定的GPU架构。使用CUDA 8.0，它编译为以下PTX指令序列：
add.s32         %r3, %r2, -1;
popc.b32        %r4, %r3;
shr.u32         %r5, %r1, %r4;
and.b32         %r6, %r3, %r1;
setp.ne.s32     %p1, %r6, 0;
selp.u32        %r7, 1, 0, %p1;
add.s32         %r8, %r5, %r7;

对于sm_5x，这基本上将其转换为机器代码1:1，除了两个指令SETP和SELP被合并成一个ICMP，因为与0进行比较。
我也没有... 我只有一个5.2，下周可能会用它进行一些性能比较。但你仍然可以发布PTX和/或SASS，使用Robert的内核，这样我们就可以看到它的样子。 - einpoklum
哦，等等，你基本上只是替换了lg()的实现，除此之外，你正在执行更直接的向上舍入修正版本。现在我想想，也许我应该用你基于popcnt的代码替换lg()，这几乎肯定更好。 - einpoklum
@einpoklum，我按照您的建议添加了生成的PTX代码。 - njuffa
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接