程序员经验分享-hwhale

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

2024-03-22

原问题

我有以下内核使用非均匀节点执行插值，我想对其进行优化：

__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.3ms。我已经证实最昂贵的零件是

            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;

大约占 60%。通过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移动到大小为2*K+1的寄存器rtemp数组（在我的例子中K是等于6的常数）； 2）我将变量phi_cap_s移动到2*K+1大小的寄存器； 3）我使用了if ... else语句代替了之前使用的三个if（条件P0.具有相同的出现概率）； 3）我为平方根定义了额外的变量； 4）我使用rsqrt而不是sqrt（据我所知，sqrt()由CUDA计算为1/rsqrt()）；

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

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

我的问题是。我完成了吗？即，有什么办法可以改进代码吗？我使用 NVIDIA Visual Profiler 对代码进行了分析，结果如下：

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：CUDA 中的一维三次样条插值 https://stackoverflow.com/questions/10859322/cuda-1-dimensional-cubic-spline-interpolation-in-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;
    }
}


//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

  2. 这会将请求的内存获取到多个寄存器，您之前所做的操作肯定会由于没有及时可用的内存用于计算而导致该线程阻塞。请记住，如果一个线程在扭曲中阻塞，则整个扭曲都会被阻塞。如果就绪队列中没有足够的扭曲，程序将阻塞，直到任何扭曲准备就绪。

  3. 现在，我们将计算时间进一步向前推进，以补偿错误的内存访问，希望之前完成的计算能够补偿错误的访问模式。

这应该起作用的原因如下：

GMEM 中的内存请求大约 >~400-600 个时钟周期。如果线程尝试对当时不可用的内存执行操作，它将阻塞。这意味着，如果每个内存请求不在 L1-L2 中，则每个 warp 都必须等待该时间或更长时间，直到它可以继续。

我怀疑的是temp.x+phi_cap_s*Uj[PP].x就是这么做的。通过将每个内存传输暂存（步骤 2）到寄存器，然后继续进行下一个暂存，您将允许您在内存传输时执行其他工作，从而隐藏延迟。

当您到达步骤 3 时，内存有望可用，否则您必须等待更少的时间。

If 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; 
    }

Edit

表达

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

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

Edit 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). As m-K代码中出现了两次，我先放入PP计算时使用exp and 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），我们可以使用这个函数来代替 mod 运算符。这应该“快得多”。检查一下，这样我的算术就正确了。它来自第 14 页。

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

CUDA

interpolation

使用非均匀节点优化 CUDA 内核插值 的相关文章

  • 构建 Erlang 服务器场（用于业余爱好项目）最便宜的方法是什么？ [关闭]

    Closed 这个问题是无关 help closed questions 目前不接受答案 假设我们有一个 本质上并行 的问题需要用 Erlang 软件来解决 我们有很多并行进程 每个进程都执行顺序代码 不是数字运算 并且我们向它们投入的 C

  • 指定 NVCC 用于编译主机代码的编译器

    运行 nvcc 时 它始终使用 Visual C 编译器 cl exe 我怎样才能让它使用GCC编译器 设置CC环境变量到gcc没有修复它 我在可执行文件帮助输出中也找不到任何选项 在 Windows 上 NVCC 仅支持 Visual C

  • 如何找到具有负值的索引并将该值替换为最接近的正值索引值？

    我知道如何从矩阵中查找具有负值的索引 matrix matrix lt 0 something should be done 但不知道如何用最接近的恰好为正的索引值替换它们的值 这里最近的索引应该与观察到的索引位于同一行 如果该行中没有具有

  • 如何在 gitlab-ci docker 执行器中使用 cuda

    我们正在使用 gitlab 持续集成来构建和测试我们的项目 最近 其中一个项目添加了 CUDA 的要求以启用 GPU 加速 我不想改变我们的管道 docker 和 gitlab ci 对我们来说运行良好 所以我想以某种方式让 docker

  • 在 R 中有效地插值网格中的值

    我有一个按位置排列的海洋深度数据网格 并且正在尝试为选定的 GPS 点插入深度值 我们一直在使用 RSAGA pick from points 它对于小数据集效果很好 require RSAGA depthdata lt cbind dat

  • Ubuntu 11.10/12.04 上的 CUDA“无兼容设备”错误

    一段时间以来 我一直在尝试在我的笔记本电脑上设置 Ubuntu 环境来进行 CUDA 编程 我目前双启动 Windows 8 和 Ubuntu 12 04 并想在 Ubuntu 上安装 CUDA 5 该笔记本电脑配有 GeForce GT

  • CUDA程序导致nvidia驱动程序崩溃

    当我超过大约 500 次试验和 256 个完整块时 我的 monte carlo pi 计算 CUDA 程序导致我的 nvidia 驱动程序崩溃 这似乎发生在 monteCarlo 内核函数中 任何帮助都会受到赞赏 include

  • 为什么numba cuda调用几次后运行速度变慢？

    我正在尝试如何在 numba 中使用 cuda 然而我却遇到了与我预想不同的事情 这是我的代码 from numba import cuda cuda jit def matmul A B C Perform square matrix m

  • R 中的线性插值

    我有一个真实数据的数据集 例如如下所示 Dataset 1 with known data known lt data frame x c 0 6 y c 0 10 20 23 41 39 61 plot known x known y t

  • CUDA Thrust 和 sort_by_key

    我正在寻找 CUDA 上的排序算法 它可以对元素数组 A 双精度 进行排序 并返回该数组 A 的键 B 数组 我知道sort by keyThrust 库中的函数 但我希望元素数组 A 保持不变 我能做些什么 我的代码是 void sort

  • 如何对重新采样的音频数据进行双三次（或其他非线性）插值？

    我正在编写一些以不同速度播放 WAV 文件的代码 以便波形要么更慢 音调更低 要么更快 音调更高 我目前正在使用简单的线性插值 如下所示 int newlength int Math Round rawdata Length lengthM

  • 如何在 Visual Studio 2010 中设置 CUDA 编译器标志？

    经过坚持不懈的得到error identifier atomicAdd is undefined 我找到了编译的解决方案 arch sm 20旗帜 但是如何在 VS 2010 中传递这个编译器标志呢 我已经尝试过如下Project gt P

  • CUDA 估计 2D 网格数据的每块线程数和块数

    首先我要说的是 我已经仔细阅读了所有类似的问题 确定每个块的线程和每个网格的块 https stackoverflow com questions 4391162 cuda determining threads per block blo

  • 如何在 CUDA 中执行多个矩阵乘法？

    我有一个方阵数组int M 10 以便M i 定位第一个元素i th 矩阵 我想将所有矩阵相乘M i 通过另一个矩阵N 这样我就收到了方阵数组int P 10 作为输出 我看到有不同的可能性 分配不同元素的计算M i 到不同的线程 例如 我

  • 设置最大 CUDA 资源

    我想知道是否可以设置 CUDA 应用程序的最大 GPU 资源 例如 如果我有一个 4GB GPU 但希望给定的应用程序只能访问 2GB 如果它尝试分配更多 就会失败 理想情况下 这可以在进程级别或 CUDA 上下文级别上设置 不 目前没有允

  • Python 3D 插值加速

    我有以下用于插入 3D 体积数据的代码 Y X Z np shape volume xs np arange 0 X ys np arange 0 Y zs np arange 0 Z points list zip np ravel re

  • Yocto for Nvidia Jetson 由于 GCC 7 而失败 - 无法计算目标文件的后缀

    我正在尝试将 Yocto 与 meta tegra 一起使用 https github com madisongh meta tegra https github com madisongh meta tegra 为 Nvidia Jets

  • cuda中有模板化的数学函数吗？ [复制]

    这个问题在这里已经有答案了 我一直在寻找 cuda 中的模板化数学函数 但似乎找不到 在普通的 C 中 如果我调用std sqrt它是模板化的 并且将根据参数是浮点数还是双精度数执行不同的版本 我想要这样的 CUDA 设备代码 我的内核将真

  • 有没有一种有效的方法来优化我的序列化代码？

    这个问题缺乏细节 因此 我决定创建另一个问题而不是编辑这个问题 新问题在这里 我可以并行化我的代码吗 还是不值得 https stackoverflow com questions 17937438 can i parallelize my

  • CUDA、NPP 滤波器

    CUDA NPP 库支持使用 nppiFilter 8u C1R 命令过滤图像 但不断出现错误 我可以毫无问题地启动并运行 boxFilterNPP 示例代码 eStatusNPP nppiFilterBox 8u C1R oDeviceS

随机推荐

热门标签