程序员经验分享-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 内核插值 的相关文章

  • __device__ __constant__ 常量

    有什么区别吗 在 CUDA 程序中定义设备常量的最佳方法是什么 在 C 主机 设备程序中 如果我想将常量定义在设备常量内存中 我可以这样做 device constant float a 5 constant float a 5 问题 1

  • IE10中的图像插值

    这是我的用例 我有一个采用响应式设计的网页 该页面垂直分成两半 我想在右侧显示图像 呈现为 PNG 或 JPG 的 PDF 页面 调整窗口大小后 图像的大小应立即更改 我以为我已经解决了这个问题 我将服务器上的图像渲染得足够大 以适应最大可

  • 内联 PTX 汇编代码强大吗？

    我看到一些代码示例 人们在 C 代码中使用内联 PTX 汇编代码 CUDA工具包中的文档提到PTX很强大 为什么会这样呢 如果我们在 C 代码中使用这样的代码 我们会得到什么好处 内联 PTX 使您可以访问未通过 CUDA 内在函数公开的指

  • CUDA 中指令重放的其他原因

    这是我从 nvprof CUDA 5 5 获得的输出 Invocations Metric Name Metric Description Min Max Avg Device Tesla K40c 0 Kernel MyKernel do

  • 如何运行和理解CUDA Visual Profiler？

    我已经设置了 CUDA 5 0 并且我的 CUDA 项目运行良好 但我不知道如何使用 Visual Profiler 分析我的 CUDA 项目 如何运行它 我还需要安装更多吗 又该如何做呢 我的电脑使用Window 7 64位 CUDA 5

  • 使用 data.table 对分组数据进行插值

    这是我最初发布的问题的延续http r 789695 n4 nabble com subset Between data table list and single data table object tp4673202 html http

  • cuda-gdb 错误消息

    我尝试使用 cuda gdb 调试我的 CUDA 应用程序 但遇到了一些奇怪的错误 我设置了选项 g G O0构建我的应用程序 我可以在没有 cuda gdb 的情况下运行我的程序 但没有得到正确的结果 因此我决定使用 cuda gdb 但

  • 大型跨平台软件项目的技巧/资源

    我将开始一个大型软件项目 涉及跨平台 GUI 和大量的数字运算 我计划用 C 和 CUDA 编写大部分应用程序后端 并用 Qt4 编写 GUI 我计划使用 Make 作为我的构建系统 这将是一个只有两名开发人员的项目 一旦我相对深入地了解它

  • cuda 文件组织的有效方式：.cpp .h .cu .cuh .curnel 文件

    cuda最容易理解 最高效的代码组织是什么 经过一番调查后 我发现 cuda 函数声明应位于 cuh 文件中 实现位于 cu 文件中 内核函数实现位于 curnel 文件中 其他 C 内容通常在 cpp 和 h 文件中 最近我发布了一个问题

  • CUDA 5.0错误LNK2001：cuda方法无法解析的外部符号

    我的链接器有错误 1 gt ManifestResourceCompile 1 gt All outputs are up to date 1 gt kernel cu obj error LNK2001 unresolved extern

  • Python/Scipy 2D 插值（非均匀数据）

    这是我上一篇文章的后续问题 Python Scipy 插值 地图坐标 https stackoverflow com questions 5124126 python scipy interpolation map coordinates

  • 如何选择图像插值方法？ （Emgu/OpenCV）

    Emgu OpenCV的 net包装器 提供的图像调整大小功能可以使用四种插值方法中的任意一种 http www emgu com wiki files 1 4 0 0 html 596dd03d 301e d3c6 4c53 c42855

  • 用于计算邻居列表的最佳 GPU 算法

    给定 3D 中数千个点的集合 我需要获取落在某个截止值 以欧几里得距离而言 内的每个粒子的邻居列表 并且如果可能的话 从最近到最远排序 在 CUDA 或 OpenCL 语言中 哪种 GPU 算法最快 我所知道的最快的 GPU MD 代码之一

  • OpenCV 2.4.3rc 和 CUDA 4.2：“OpenCV 错误：没有 GPU 支持”

    我在这张专辑中上传了几张截图 https i stack imgur com TELST jpg https i stack imgur com TELST jpg 我正在尝试在 Visual Studio 2008 中的 OpenCV 中

  • C 中的 CUDA：如何使用 cudaMemcpyAsync 修复错误 11

    我目前正在尝试使用 CUDA 运行一个简单的多 GPU 程序 它的基本作用是将一个包含一些虚拟数据的大型数组复制到 GPU GPU 进行一些数学计算 然后将结果数组复制回来 我在 VS2017 的输出中没有收到任何错误 但我设置的一些错误消

  • CUDA 添加矩阵的行

    我试图将 4800x9600 矩阵的行加在一起 得到一个 1x9600 的矩阵 我所做的是将 4800x9600 分成 9 600 个矩阵 每个矩阵长度为 4800 然后我对 4800 个元素进行缩减 问题是 这真的很慢 有人有什么建议吗

  • 布尔实现的atomicCAS

    我想弄清楚是否存在错误答案 https stackoverflow com a 57444538 11248508 现已删除 关于Cuda like的实现atomicCAS for bool是 答案中的代码 重新格式化 static inl

  • 完全禁用 NVCC 优化

    我正在尝试测量 GPU 上的峰值单精度触发器 为此我正在修改 PTX 文件以在寄存器上执行连续的 MAD 指令 不幸的是 编译器正在删除所有代码 因为它实际上没有做任何有用的事情 因为我没有执行任何数据的加载 存储 是否有编译器标志或编译指

  • CUDA cutil.h 在哪里？

    有谁知道包含 cutil h 的 SDK 工具包在哪里 我尝试了 CUDA toolkits3 2 和 toolkits5 0 我知道这个版本已经不支持 cutil h 我还注意到一些提到的如何在 Linux 中包含 cutil h htt

  • goto 指令对 CUDA 代码中扭曲内发散的影响

    对于CUDA中简单的warp内线程发散 我所知道的是SM选择一个重新收敛点 PC地址 并在两个 多个路径中执行指令 同时禁用未采用该路径的线程的执行效果 例如 在下面的代码中 if threadIdx x lt 16 A do someth

随机推荐

热门标签