低效的内核功能

最简单的常规优化如下所示：

__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提出的建议（除展开实用性外）似乎使性能提高了约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
$

笔记：

1. 展开实用程序似乎无济于事（对于我尝试的一些测试用例，它会使运行速度变慢）。在某些情况下，编译器将已经展开了没有特定提示的循环，并且循环展开通常是需要进行调整（可能需要仔细调整）的优化。
2. @talonmies建议对内核进行的修改以创建网格跨越循环，这是为了使特定的循环展开行程计数有用而需要考虑的因素之一。整个网格的尺寸至少应减少等于展开行程数的一个因数。但是我找不到一个“最佳地点”。
3. 我主要在Quadro5000（Fermi cc2.0 GPU），CUDA 7.5RC，Fedora20上进行过测试。当然，在不同的GPU（尤其是较新的GPU）上，行为将有所不同。
4. nBLK这段代码中的参数是另一个“可调”参数，但是在大约64左右时，此参数几乎没有变化。最好的情况是使网格的大小等于数据的大小。