低效的核函数

Question

有没有可能加速这个简单的内核函数？我考虑过使用共享内存，但 N 等于 507904，所以它比共享内存数组要多得多。
我的程序每个创建 256 个线程块。

__global__ void compute(COMPLEX_TYPE *a, COMPLEX_TYPE *b,
              FLOAT_TYPE *F, FLOAT_TYPE f, int N) 
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) 
    {
        F[i] = ( a[i].x*a[i].x + a[i].y*a[i].y + b[i].x*b[i].x + b[i].y*b[i].y) / (f);
    }
}

Answer 1

最简单的一般优化是这样的：

__global__ void compute(const COMPLEX_TYPE * __restrict__ a, 
                        const COMPLEX_TYPE * __restrict__ b,
                        FLOAT_TYPE *F, FLOAT_TYPE f, int N) 
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    #pragma unroll 8
    for(; i < N; i += blockDim.x * gridDim.x;)
    {
        COMPLEX_TYPE aval = a[i], bval = b[i]
        FLOAT_TYPE Fval;
        Fval = ( aval.x*aval.x + aval.y*aval.y + bval.x*bval.x + bval.y*bval.y) / (f);
        F[i] = Fval;
    }
}

[免责声明：在浏览器中编写，未经测试，使用风险自负]

这里的想法是只启动与目标 GPU 上并发执行的线程一样多的线程，然后让每个线程执行多个操作而不是一个。这有助于分摊块调度程序和设置代码级别的大量固定开销，并提高整体效率。在大多数架构上，这可能无论如何都会限制内存带宽，因此内存合并和事务优化是您将能够进行的最重要的性能优化。

编辑： 由于这个答案被标记为 CW，我选择在这里添加我的测试，而不是创建我自己的答案。如果有人对此有异议，请将编辑回滚到以前可接受的版本。我没有添加任何新想法，只是测试@talonmies 和@JanLucas

提供的想法

在我的测试用例中，@talonmies 提供的建议（除了 unroll pragma）似乎提高了 ~10% 的性能。 @JanLucas 提出的用浮点乘法代替浮点除法的建议（如果可以接受的话）似乎可以使性能提高一倍。这显然会因 GPU 和其他细节而异。这是我的测试：

$ cat t891.cu
#include <cuComplex.h>
#include <stdio.h>
#include <stdlib.h>

#define DSIZE 507904
#define nTPB 256
#define nBLK 256

#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)

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

typedef cuFloatComplex COMPLEX_TYPE;
typedef float FLOAT_TYPE;

__global__ void compute(COMPLEX_TYPE *a, COMPLEX_TYPE *b,
              FLOAT_TYPE *F, FLOAT_TYPE f, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N)
    {
        F[i] = ( a[i].x*a[i].x + a[i].y*a[i].y + b[i].x*b[i].x + b[i].y*b[i].y) / (f);
    }
}

__global__ void compute_imp(const COMPLEX_TYPE * __restrict__ a,
                        const COMPLEX_TYPE * __restrict__ b,
                        FLOAT_TYPE *F, FLOAT_TYPE f, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
//    #pragma unroll 8
    for(; i < N; i += blockDim.x * gridDim.x)
    {
        COMPLEX_TYPE aval = a[i];
        COMPLEX_TYPE bval = b[i];
        FLOAT_TYPE Fval = ( aval.x*aval.x + aval.y*aval.y + bval.x*bval.x + bval.y*bval.y) / (f);
        F[i] = Fval;
    }
}

__global__ void compute_imp2(const COMPLEX_TYPE * __restrict__ a,
                        const COMPLEX_TYPE * __restrict__ b,
                        FLOAT_TYPE *F, FLOAT_TYPE f, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
//    #pragma unroll 8
    for(; i < N; i += blockDim.x * gridDim.x)
    {
        COMPLEX_TYPE aval = a[i];
        COMPLEX_TYPE bval = b[i];
        FLOAT_TYPE Fval = ( aval.x*aval.x + aval.y*aval.y + bval.x*bval.x + bval.y*bval.y) * (f);
        F[i] = Fval;
    }
}

int main(){

  COMPLEX_TYPE *d_A, *d_B;
  FLOAT_TYPE *d_F, f = 4.0f;

  cudaMalloc(&d_A, DSIZE*sizeof(COMPLEX_TYPE));
  cudaMalloc(&d_B, DSIZE*sizeof(COMPLEX_TYPE));
  cudaMalloc(&d_F, DSIZE*sizeof(FLOAT_TYPE));

  //warm-up
  compute<<<(DSIZE+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_F, f, DSIZE);
  cudaDeviceSynchronize();
  unsigned long long t1 = dtime_usec(0);
  compute<<<(DSIZE+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_F, f, DSIZE);
  cudaDeviceSynchronize();
  t1 = dtime_usec(t1);

  //warm-up
  compute_imp<<<DSIZE/(8*nTPB),nTPB>>>(d_A, d_B, d_F, f, DSIZE);
  cudaDeviceSynchronize();
  unsigned long long t2 = dtime_usec(0);
  compute_imp<<<nBLK,nTPB>>>(d_A, d_B, d_F, f, DSIZE);
  cudaDeviceSynchronize();
  t2 = dtime_usec(t2);

  //warm-up
  compute_imp2<<<(DSIZE+nTPB-1)/nTPB,nTPB>>>(d_A, d_B, d_F, 1/f, DSIZE);
  cudaDeviceSynchronize();
  unsigned long long t3 = dtime_usec(0);
  compute_imp2<<<nBLK,nTPB>>>(d_A, d_B, d_F, 1/f, DSIZE);
  cudaDeviceSynchronize();
  t3 = dtime_usec(t3);
  cudaCheckErrors("some error");
  printf("t1: %fs, t2: %fs, t3: %fs\n", t1/(float)USECPSEC, t2/(float)(USECPSEC), t3/(float)USECPSEC);
}
$ nvcc -O3 -o t891 t891.cu
$ ./t891
t1: 0.000226s, t2: 0.000209s, t3: 0.000110s
$

备注：

unroll pragma 似乎没有帮助（它使运行变慢，对于我尝试的一些测试用例）。在某些情况下，编译器已经会在没有特定提示的情况下展开循环，并且循环展开通常是一种需要调整的优化，也许需要仔细调整。
@talonmies 提议对内核进行修改以创建跨网格循环，这是使特定循环展开行程计数有用需要考虑的因素之一。总体网格维度应至少减少一个等于展开行程计数的因子。但是我找不到 "sweet spot".
我主要在 Quadro5000（Fermi cc2.0 GPU）、CUDA 7.5RC、Fedora20 上进行测试。当然，在不同的 GPU 上，行为会有所不同，尤其是较新的 GPU。
此代码中的 nBLK 参数是另一个 "tunable" 参数，但是当超过 64 左右时我发现它几乎没有变化。最好的情况可能是网格大小与数据相等。

低效的核函数

Inefficient kernel function

optimization

cuda

shared-memory