使用非均匀节点优化CUDA内核插值

Question

原始问题

我有以下内核，用于在非均匀节点上进行插值，我想要对其进行优化：

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P;
    const double alfa=(2.-1./cc)*pi_double-0.01;
    double phi_cap_s;
    cufftDoubleComplex temp;

    double cc_points=cc*points[i];
    double r_cc_points=rint(cc*points[i]);

    temp = make_cuDoubleComplex(0.,0.);

    if(i<M) {   
        for(int m=0; m<(2*K+1); m++) {
            P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));

            if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
            if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
            if(P==0.) phi_cap_s = alfa/pi_double;        

            PP = modulo((r_cc_points + m -K ),(cc*N)); 
            temp.x = temp.x+phi_cap_s*Uj[PP].x; 
            temp.y = temp.y+phi_cap_s*Uj[PP].y; 
        } 

        result[i] = temp; 
    }
}

K和cc是常数，points包含节点，Uj是要插值的值。modulo是一个函数，基本上像%，但对负值进行了适当扩展。对于某个排列，内核调用需要2.3毫秒。我已经验证了最昂贵的部分是

            if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
            if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
            if(P==0.) phi_cap_s = alfa/pi_double;        

其中大约占总时间的40％，

        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        temp.x = temp.x+phi_cap_s*Uj[PP].x; 
        temp.y = temp.y+phi_cap_s*Uj[PP].y; 

这段代码涉及到IT技术，其中大约60%的计算时间被一个if语句占用。通过可视化分析工具Visual Profiler，我已经验证了前者的性能不会受到if语句存在的影响。请注意，我需要双精度，所以我避免使用__exp()解决方案。我怀疑，对于后者，“随机”内存访问Uj[PP]可能是导致计算百分比如此之高的原因。有什么技巧或建议可以减少计算时间吗？提前感谢。

以下是根据评论和答案进行修改后的版本

根据评论和答案中提供的建议，我最终得到了下面的代码：

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P,tempd;
    const double alfa=(2.-1./cc)*pi_double-0.01;
    cufftDoubleComplex temp = make_cuDoubleComplex(0.,0.);

    double cc_points=cc*points[i];
    double r_cc_points=rint(cc_points);

    cufftDoubleComplex rtemp[(2*K+1)];
    double phi_cap_s[2*K+1];

    if(i<M) {   
     #pragma unroll //unroll the loop
     for(int m=0; m<(2*K+1); m++) {
         PP = modulo(((int)r_cc_points + m -K ),(cc*N)); 
            rtemp[m] = Uj[PP]; //2

         P = (K*K-(cc_points-(r_cc_points+(double)(m-K)))*(cc_points-(r_cc_points+(double)(m-K))));
         if(P<0.) {tempd=rsqrt(-P); phi_cap_s[m] = (1./pi_double)*((sin(alfa/tempd))*tempd);  }
         else if(P>0.) {tempd=rsqrt(P); phi_cap_s[m] = (1./pi_double)*((sinh(alfa/tempd))*tempd); }
         else phi_cap_s[m] = alfa/pi_double;  
     }

     #pragma unroll //unroll the loop
     for(int m=0; m<(2*K+1); m++) {
         temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; 
           temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
     } 

     result[i] = temp; 
     }
 }

特别地： 1）我将全局内存变量Uj移动到大小为2K+1的寄存器rtemp数组中（在我的情况下，K是一个常数，等于6）； 2）我将变量phi_cap_s移动到一个大小为2K+1的寄存器中； 3）我使用if...else语句代替之前使用的三个if语句（条件P<0.和P>0.的出现概率相同）； 4）我为平方根定义了额外的变量； 5）我使用rsqrt而非sqrt（据我所知，CUDA计算sqrt()时采用1/rsqrt()）；

我一次添加每个新功能并验证其与原始版本的改进，但我必须说，它们都没有给我带来任何相关的改进。

执行速度受以下限制： 1）sin / sinh函数的计算（约占时间的40％）；是否有一种方式可以通过某种方式利用内在数学作为“起始猜测”以双精度运算来计算它们？ 2）许多线程最终会访问由于映射索引PP而导致相同的全局内存位置Uj [PP]。避免这种情况的一种可能性是使用共享内存，但这将意味着强烈的线程合作。

我的问题是。我完成了吗？也就是说，有没有改进代码的方法？我通过NVIDIA视觉分析器对代码进行了剖析，以下是结果：

IPC = 1.939 (compute capability 2.1);
Global Memory Load Efficiency = 38.9%;
Global Memory Store Efficiency = 18.8%;
Warp Execution Efficiency = 97%;
Instruction Replay Overhead = 0.7%;

最后，我想指出这个讨论与CUDA中的一维三次样条插值的讨论有关。

使用共享内存的版本

我已经进行了一个使用共享内存的可行性研究。我考虑了 N=64，以便整个Uj适合共享内存。下面是代码（基本上是我的原始版本）。

    __global__ void interpolation_shared(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
 {
         int i = threadIdx.x + blockDim.x * blockIdx.x;

     int PP;
     double P;
     const double alfa=(2.-1./cc)*pi_double-0.01;
     double phi_cap_s;
     cufftDoubleComplex temp;

     double cc_points=cc*points[i];
     double r_cc_points=rint(cc*points[i]);

     temp = make_cuDoubleComplex(0.,0.);

     __shared__ cufftDoubleComplex Uj_shared[128];

     if (threadIdx.x < cc*N) Uj_shared[threadIdx.x]=Uj[threadIdx.x];

     if(i<M) {  
         for(int m=0; m<(2*K+1); m++) {
         P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));

         if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
         if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));  
         if(P==0.) phi_cap_s = alfa/pi_double;        

         PP = modulo((r_cc_points + m -K ),(cc*N)); 
         temp.x = temp.x+phi_cap_s*Uj_shared[PP].x; 
         temp.y = temp.y+phi_cap_s*Uj_shared[PP].y; 
      } 

      result[i] = temp; 
    }
 }

结果再次没有显著提高，尽管这可能取决于输入数组的小尺寸。

详细的 PTXAS 输出：

ptxas : info : Compiling entry function '_Z13interpolationP7double2PdS0_ii' for 'sm_20'
ptxas : info : Function properties for _Z13interpolationP7double2PdS0_ii
  352 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas : info : Used 55 registers, 456 bytes cumulative stack size, 52 bytes cmem[0]

第一层和m=0时，P的值

 0.0124300933082964
 0.0127183892149176
 0.0135847002913749
 0.0161796378170038
 0.0155488126345702
 0.0138890822153499
 0.0121163187739057
 0.0119998374528905
 0.0131600831194518
 0.0109574866163769
 0.00962949548477354
 0.00695850974164358
 0.00446426651940612
 0.00423369284281705
 0.00632921297092537
 0.00655137618976198
 0.00810202954519923
 0.00597974034698723
 0.0076811348379735
 0.00604267951733561
 0.00402922460255439
 0.00111841719893846
 -0.00180949615796777
 -0.00246283218698551
 -0.00183256444286428
 -0.000462696661685413
 0.000725108980390132
 -0.00126793006072035
 0.00152263101649197
 0.0022499598348702
 0.00463681632275836
 0.00359856091027666

模运算

__device__ int modulo(int val, int modulus)
{
   if(val > 0) return val%modulus;
   else
   {
       int P = (-val)%modulus;
       if(P > 0) return modulus -P;
       else return 0;
   }
}

根据答案优化的模数函数

__device__ int modulo(int val, int _mod)
{
    if(val > 0) return val&(_mod-1);
    else
    {
        int P = (-val)&(_mod-1);
        if(P > 0) return _mod -P;
        else return 0;
    }
}

- Vitality

---

当随机访问不可避免时，要做的一件事是同时处理多个值。例如，如果您想将两个东西相加a + b，首先将a1提取到寄存器中，然后是b1，a2，b2，b3，b4...然后按提取顺序将它们相加，a1+b1，a2+b2...在我的应用程序中，我发现每个线程处理4到8个值可以显著提高性能。如果您在这里查看https://docs.google.com/spreadsheet/ccc?key=0Au1A39JI-BwHdHpESVZRODZoM2VoTzFtcG91bW1mVlE#gid=1，我比较块大小和每个线程处理的值。使用处理2个值与4个值相比，我可以获得最多两倍的速度！ - 1-----1

你确定了你的算法是计算受限还是内存受限吗？如果你得到的双精度吞吐量接近芯片的最大值，我认为你已经无能为力了。 - Roger Dahl

占用情况如何？如果您能使用-Xptxas –v编译并向我们展示结果，那就更好了。当进行剖析时，请问您能否向我们展示“指令重放开销”是多少！ - 1-----1

@JackOLantern，“指令重放开销”为0.7％意味着您基本上不会因阻塞的warp而减慢速度，这意味着除非优化指令使用，否则无法真正提高性能。我将在我的帖子中更新一种优化方法。 - 1-----1

@JackOLantern，你的编译命令怎么样了？你是否开启了优化-O3并关闭了调试-G？因为开启调试后，编译器无法对内核进行优化。此外，你真正能做的就是尝试减少指令和寄存器的使用。目前你使用了56个寄存器，这意味着你的占用率可能会受到影响，如果将--maxrregcount设置为26-32，会导致寄存器溢出到本地内存，但可能会加快执行速度。 - 1-----1

Answer 1

//your code above
cufftDoubleComplex rtemp[(2*K+1)] //if it fits into available registers, assumes K is a constant

if(i<M) {   
#pragma unroll //unroll the loop
    for(int m=0; m<(2*K+1); m++) {

        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        rtemp[m] = Uj[PP]; //2
    }
#pragma unroll
    for(nt m=0; m<(2*K+1); m++) {
        P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));
        // 1
        if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
        else if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
        else phi_cap_s = alfa/pi_double;  

        temp.x = temp.x+phi_cap_s*rtemp[m].x; //3
        temp.y = temp.y+phi_cap_s*rtemp[m].y; 
    }

    result[i] = temp; 
}

解释

1. 添加 else if 和 else 语句，因为这些条件是互斥的。如果可能的话，应该按照发生概率对语句进行排序。例如，如果 P＜0 大多数情况下，应该优先评估它。

2. 这将把请求的内存提取到多个寄存器中，您之前所做的可能会导致线程阻塞，因为计算所需的内存未能及时获得。请记住，如果一个 warp 中的线程阻塞，整个 warp 都会被阻塞。如果就绪队列中没有足够的 warp，则程序将阻塞，直到任何一个 warp 就绪。

3. 我们现在将计算向前移动一段时间，以弥补不良内存访问的影响，希望之前所做的计算已经弥补了坏访问模式。

这样做应该有效的原因如下：

从 GMEM 中请求的内存大约需要 400-600 个 tick。如果一个线程尝试对不可用的内存执行操作，它将被阻塞。这意味着如果每个内存请求都不在 L1-L2 中，每个 warp 都必须等待相同或更长的时间才能继续执行。

我怀疑的是temp.x+phi_cap_s*Uj[PP].x正在做这件事。通过将每个内存传输分阶段（步骤2）到寄存器中，并继续进行下一个阶段，您可以隐藏延迟，从而在内存传输时允许您执行其他工作。

当您到达第3步时，内存有望可用，或者您需要等待更短的时间。

如果rtemp无法适应寄存器以实现100％的占用率，则可能需要分批处理。

您还可以尝试将phi_cap_s变成一个数组，并将其放入第一个循环中，如下所示：

#pragma unroll //unroll the loop
    for(int m=0; m<(2*K+1); m++) {
        //stage memory first
        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        rtemp[m] = Uj[PP]; //2

        P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));
        // 1
        if(P>0.)  phi_cap_s[m] = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
        else if(P<0.)  phi_cap_s[m] = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
        else phi_cap_s[m] = alfa/pi_double; 

    }
#pragma unroll
    for(nt m=0; m<(2*K+1); m++) {
        temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; //3
        temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
    }

编辑

表达式

P = (K*K-(cc_points-(r_cc_points+(double)(m-K)))*(cc_points-(r_cc_points+(double)(m-K))));

可以分解为：

const double cc_diff = cc_points-r_cc_points;
double exp = cc_diff - (double)(m-K);
exp *= exp;
P = (K*K-exp);

可能会减少使用的指令数量。

编辑2

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P,tempd;


    cufftDoubleComplex rtemp[(2*K+1)];
    double phi_cap_s[2*K+1];

    if(i<M) {
         const double cc_points=cc*points[i];
         cufftDoubleComplex temp = make_cuDoubleComplex(0.,0.);

         const double alfa=(2.-1./cc)*pi_double-0.01;


         const double r_cc_points=rint(cc_points);
         const double cc_diff = cc_points-r_cc_points;

     #pragma unroll //unroll the loop
         for(int m=0; m<(2*K+1); m++) {
             PP = m-k; //reuse PP
             double exp = cc_diff - (double)(PP); //stage exp to be used later, will explain

             PP = modulo(((int)r_cc_points + PP ),(cc*N)); 
             rtemp[m] = Uj[PP]; //2


             exp *= exp;
             P = (K*K-exp);

             if(P<0.) {tempd=rsqrt(-P); phi_cap_s[m] = (1./pi_double)*((sin(alfa/tempd))*tempd);  }
             else if(P>0.) {tempd=rsqrt(P); phi_cap_s[m] = (1./pi_double)*((sinh(alfa/tempd))*tempd); }
             else phi_cap_s[m] = alfa/pi_double;  
         }

     #pragma unroll //unroll the loop
         for(int m=0; m<(2*K+1); m++) {
             temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; 
             temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
         } 

     result[i] = temp; 
     }
 }

我所做的是将所有计算移至if语句内部，以释放一些资源，包括计算和内存提取。不知道你对第一个if语句if(i<M)有何看法。由于代码中出现了两次m-K，我首先将其放入PP中，在计算exp和PP时使用。

另外，你可以尝试按顺序排列指令，这样如果你设置了一个变量，在下一次使用该变量之前尽可能多地执行指令，因为它需要大约20个时钟周期才能被设置到寄存器中。因此，我将常量cc_diff放在了最上面，但由于这只是一个一次性的指令，可能不会显示任何好处。

模函数

__device__ modulo(int val, int _mod) {
    int p = (val&(_mod-1));// as modulo is always the power of 2
    if(val < 0) {
        return _mod - p;
    } else {
        return p;
    }
}

由于我们总是将_mod作为2的幂次方的整数（cc = 2，N = 64，cc * N = 128），因此我们可以使用此函数代替模运算符。这应该会更快。请检查一下，以确保我的算术正确。这来自Optimizing Cuda - Part II Nvidia第14页。

- 1-----1

---

谢谢你的回答。根据你的回答，我已经修改了我的初始帖子，但是没有显著的改进。请仍然看一下。我有一个问题。根据你的解决方案，在线程可以执行其他操作的同时，Uj正在从全局内存加载到寄存器中，这似乎很合理。但是在使用寄存器值之前，即在最后的for循环之前，我是否还需要__syncthreads()呢？ - Vitality

@JackOLantern 不需要使用__syncthreads()，因为你没有使用共享内存进行线程间通信的操作。也就是说，这些线程之间并不依赖彼此。我会查看更新后的内容。 - 1-----1

如果值尚未加载到寄存器中，则线程和warp将被阻塞。因此，如果值尚未加载到其中，则无需使用__syncthreads，因为它不会进行计算。 - 1-----1

cc_diff已经移出循环。明天会尝试您的新代码并回答其他待定问题。谢谢。 - Vitality

在新的“模数”函数中，您基本上建议使用技巧foo％n等同于（但比）foo＆（n-1）更快，如果n是2的幂。然后我修订了我的“模数”函数（请参见修订后的帖子-您的版本在算术上没有完全等效）。通过您提到的技巧，IPC进一步提高到1.958。谢谢！ - Vitality

Answer 2

一个你可能想要考虑的优化是使用快速数学。使用内置数学函数并使用 -use-fast-math 选项进行编译。 内置数学函数

- Hong Zhou

---

我不使用内置数学函数，因为它们是近似值，而我需要完整的双精度。 - Vitality

Uj[PP]内存的访问是否在块内受限？如果是，那么使用共享内存如何？ - Hong Zhou

我已经修改了原始帖子，并添加了有关共享内存使用的评论。 - Vitality

我已经重新修订了帖子，并添加了使用共享内存进行可行性研究的结果。 - Vitality