我应该在GPU上运行这个统计应用程序的代码吗?

44

我正在开发一个统计应用程序,其中包含大约 10 到 30 百万个浮点值的数组。

有几种方法在嵌套循环中对数组执行不同但独立的计算,例如:

Dictionary<float, int> noOfNumbers = new Dictionary<float, int>();

for (float x = 0f; x < 100f; x += 0.0001f) {
    int noOfOccurrences = 0;

    foreach (float y in largeFloatingPointArray) {
        if (x == y) {
            noOfOccurrences++;
        }
    }
    noOfNumbers.Add(x, noOfOccurrences);
}

当前应用程序是用C#编写的,运行在Intel CPU上,需要几个小时才能完成。我对GPU编程概念和API没有任何了解,所以我的问题是:

  • 是否有可能(并且是否有意义)利用GPU加速这样的计算?
  • 如果是:是否有人知道任何教程或者有任何示例代码(编程语言无关)?
5个回答

89

更新 GPU 版本

__global__ void hash (float *largeFloatingPointArray,int largeFloatingPointArraySize, int *dictionary, int size, int num_blocks)
{
    int x = (threadIdx.x + blockIdx.x * blockDim.x); // Each thread of each block will
    float y;                                         // compute one (or more) floats
    int noOfOccurrences = 0;
    int a;
    
    while( x < size )            // While there is work to do each thread will:
    {
        dictionary[x] = 0;       // Initialize the position in each it will work
        noOfOccurrences = 0;    

        for(int j = 0 ;j < largeFloatingPointArraySize; j ++) // Search for floats
        {                                                     // that are equal 
                                                             // to it assign float
           y = largeFloatingPointArray[j];  // Take a candidate from the floats array 
           y *= 10000;                      // e.g if y = 0.0001f;
           a = y + 0.5;                     // a = 1 + 0.5 = 1;
           if (a == x) noOfOccurrences++;    
        }                                      
                                                    
        dictionary[x] += noOfOccurrences; // Update in the dictionary 
                                          // the number of times that the float appears 

    x += blockDim.x * gridDim.x;  // Update the position here the thread will work
    }
}

这个我只是测试了一下小的输入,因为我是在我的笔记本电脑上进行测试。尽管如此,它仍然可以工作,但需要进行更多测试。

更新 顺序版本

我刚刚做了这个天真的版本,它可以在不到20秒的时间内执行包含30000000个元素的数组的算法(包括生成数据的函数所花费的时间)。

这个天真的版本首先对浮点数数组进行排序。之后,将遍历排序后的数组,并检查给定的value在数组中出现的次数,然后将这个值与它出现的次数一起放入一个字典中。

您可以使用sorted映射来代替我使用的unordered_map

以下是代码:

#include <stdio.h>
#include <stdlib.h>
#include "cuda.h"
#include <algorithm>
#include <string>
#include <iostream>
#include <tr1/unordered_map>


typedef std::tr1::unordered_map<float, int> Mymap;


void generator(float *data, long int size)
{
    float LO = 0.0;
    float HI = 100.0;
    
    for(long int i = 0; i < size; i++)
        data[i] = LO + (float)rand()/((float)RAND_MAX/(HI-LO));
}

void print_array(float *data, long int size)
{

    for(long int i = 2; i < size; i++)
        printf("%f\n",data[i]);
    
}

std::tr1::unordered_map<float, int> fill_dict(float *data, int size)
{
    float previous = data[0];
    int count = 1;
    std::tr1::unordered_map<float, int> dict;
    
    for(long int i = 1; i < size; i++)
    {
        if(previous == data[i])
            count++;
        else
        {
          dict.insert(Mymap::value_type(previous,count));
          previous = data[i];
          count = 1;         
        }
        
    }
    dict.insert(Mymap::value_type(previous,count)); // add the last member
    return dict;
    
}

void printMAP(std::tr1::unordered_map<float, int> dict)
{
   for(std::tr1::unordered_map<float, int>::iterator i = dict.begin(); i != dict.end(); i++)
  {
     std::cout << "key(string): " << i->first << ", value(int): " << i->second << std::endl;
   }
}


int main(int argc, char** argv)
{
  int size = 1000000; 
  if(argc > 1) size = atoi(argv[1]);
  printf("Size = %d",size);
  
  float data[size];
  using namespace __gnu_cxx;
  
  std::tr1::unordered_map<float, int> dict;
  
  generator(data,size);
  
  sort(data, data + size);
  dict = fill_dict(data,size);
  
  return 0;
}

如果您的计算机安装了thrust库,则应使用以下内容:

#include <thrust/sort.h>
thrust::sort(data, data + size);

改为这样

sort(data, data + size);

肯定会更快。

原始帖子

我正在开发一个统计应用程序,其中包含一个包含10-30百万浮点值的大数组。

是否有可能(并且是否有意义)利用GPU加速这些计算?

是的,可以。一个月前,我在GPU上运行了完全分子动力学模拟。其中一个内核,计算粒子对之间的力量,接收到参数6个每个都有500,000个双精度数组,总共有3百万双精度数(22 MB)。

因此,如果您计划放置30百万个浮点数,总共约为114 MB的全局内存,这将不是问题。

在您的情况下,计算数量是否会成问题? 根据我在分子动力学(MD)方面的经验,我会说不会。 顺序MD版本需要约25小时才能完成,而GPU版本只需要45分钟。 您说您的应用程序花了几个小时,还根据代码示例看起来比MD要“柔和”。

这里是力量计算示例:

__global__ void add(double *fx, double *fy, double *fz,
                    double *x, double *y, double *z,...){
   
     int pos = (threadIdx.x + blockIdx.x * blockDim.x); 
      
     ...
     
     while(pos < particles)
     {
     
      for (i = 0; i < particles; i++)
      {
              if(//inside of the same radius)
                {
                 // calculate force
                } 
       }
     pos += blockDim.x * gridDim.x;  
     }        
  }

在CUDA中,一个简单的代码示例可以是两个2D数组的求和:

使用C语言:

for(int i = 0; i < N; i++)
    c[i] = a[i] + b[i]; 

在CUDA中:

__global__ add(int *c, int *a, int*b, int N)
{
  int pos = (threadIdx.x + blockIdx.x)
  for(; i < N; pos +=blockDim.x)
      c[pos] = a[pos] + b[pos];
}
在CUDA中,您基本上需要将每个for迭代分配给每个线程。
1) threadIdx.x + blockIdx.x*blockDim.x;
每个块都有一个ID,从0到N-1(N是块的最大数量),每个块都有一个“X”线程数,每个线程都有一个ID,从0到X-1。给出了for循环迭代,每个线程将根据其ID和线程所在的块ID计算;blockDim.x是一个块拥有的线程数。因此,如果您有2个块,每个块都有10个线程和N=40,则:
Thread 0 Block 0 will execute pos 0
Thread 1 Block 0 will execute pos 1
...
Thread 9 Block 0 will execute pos 9
Thread 0 Block 1 will execute pos 10
....
Thread 9 Block 1 will execute pos 19
Thread 0 Block 0 will execute pos 20
...
Thread 0 Block 1 will execute pos 30
Thread 9 Block 1 will execute pos 39

看了你目前的代码,我已经为你写出以下使用CUDA的代码草稿:

__global__ hash (float *largeFloatingPointArray, int *dictionary)
    // You can turn the dictionary in one array of int
    // here each position will represent the float
    // Since  x = 0f; x < 100f; x += 0.0001f
    // you can associate each x to different position
    // in the dictionary:

    // pos 0 have the same meaning as 0f;
    // pos 1 means float 0.0001f
    // pos 2 means float 0.0002f ect.
    // Then you use the int of each position 
    // to count how many times that "float" had appeared 


   int x = blockIdx.x;  // Each block will take a different x to work
    float y;
    
while( x < 1000000) // x < 100f (for incremental step of 0.0001f)
{
    int noOfOccurrences = 0;
    float z = converting_int_to_float(x); // This function will convert the x to the
                                          // float like you use (x / 0.0001)

    // each thread of each block
    // will takes the y from the array of largeFloatingPointArray
    
    for(j = threadIdx.x; j < largeFloatingPointArraySize; j += blockDim.x)
    {
        y = largeFloatingPointArray[j];
        if (z == y)
        {
            noOfOccurrences++;
        }
    }
    if(threadIdx.x == 0) // Thread master will update the values
      atomicAdd(&dictionary[x], noOfOccurrences);
    __syncthreads();
}

你需要使用atomicAdd,因为来自不同块的不同线程可能会同时读/写noOfOccurrences,所以你必须确保互斥性

这只是一种方法;你甚至可以将外部循环的迭代分配给线程而不是块。

教程

Rob Farmer撰写的《CUDA:超级计算》系列文章在其十四个部分中涵盖了几乎所有内容,非常适合初学者入门。

还有其他资源:

看最后一个项目,你会发现许多学习CUDA的链接。

OpenCL:OpenCL教程 | MacResearch


11

关于并行处理或GPGPU我不太了解,但是对于这个具体的例子,你可以通过对输入数组进行单一遍历来节省大量时间,而不是重复遍历一百万次。对于大数据集,如果可能的话通常希望进行单一遍历。即使你在进行多个独立的计算,如果它们都是基于相同的数据集,那么在同一次遍历中完成所有计算可能会更快,因为这样会获得更好的引用局部性。但是这可能不值得在代码中增加复杂性。

此外,你真的不想像那样一遍又一遍地给浮点数加上一个小的数,舍入误差会逐渐累积,你将得不到你想要的结果。我已经在下面的示例中添加了一个if语句,以检查输入是否匹配你的迭代模式,但如果你实际上不需要它,请省略。

我不懂C#,但是你的示例的单次遍历实现可能类似于此:

Dictionary<float, int> noOfNumbers = new Dictionary<float, int>();

foreach (float x in largeFloatingPointArray)
{
    if (math.Truncate(x/0.0001f)*0.0001f == x)
    {
        if (noOfNumbers.ContainsKey(x))
            noOfNumbers.Add(x, noOfNumbers[x]+1);
        else
            noOfNumbers.Add(x, 1);
    }
}
希望这有所帮助。

4
您可以通过使用TryGet而不是ContainsKey和noOfNumbers[x]来改进您的代码。使用TryGet可以节省一个字典查找,其时间复杂度为O(1)摊销(即并非始终为O(1)),而且由于字典是一种相对复杂的数据类型,因此这个O(1)操作时间成本较高。总之,+1。 - Eli Algranti
3
感谢你们的帮助,非常感激。你们的建议很快就会被添加到我的应用程序中。但不幸的是,我还有将近100个其他方法,我认为它们不能再优化了。即使我使用代码优化将这些计算加速90%,在快速CPU上完成计算仍可能需要几个小时。 - Mike
3
请将包含有限数据集和您自己的基准测试的完整方法发送给我们。这将使我们能够更好地帮助您。根据我目前在代码中看到的内容,我相信甚至在开始使用GPU之前,我就能将代码的速度提高一倍。 - Martin

9

使用GPU加速此类计算是否可行(并且是否有意义)?

  • 绝对,这种算法通常是大规模数据并行处理的理想选择,而这正是GPU擅长的事情。

如果是:有人知道任何教程或者样例代码吗(编程语言不重要)?

  • 当你想使用GPGPU时,有两种选择: CUDA 或者 OpenCL

    CUDA很成熟,并且有许多工具,但是它以NVIDIA GPU为中心。

    OpenCL是一种标准,可以在NVIDIA和AMD GPU以及CPU上运行。因此,你应该真正偏爱OpenCL。

  • 对于本教程,CodeProject上有一系列由Rob Farber撰写的优秀教程: http://www.codeproject.com/Articles/Rob-Farber#Articles

  • 对于你特定的用例,有许多使用OpenCL构建直方图的示例(注意,许多是图像直方图,但原理相同)。

  • 由于你使用C#,所以可以使用类似OpenCL.NetCloo的绑定。

  • 如果你的数组太大无法存储在GPU内存中,你可以对其进行块分割,并轻松地为每个部分重新运行OpenCL内核。


2
高效直方图算法的另一个资源...http://users.cecs.anu.edu.au/~ramtin/cuda.htm - kineticfocus
2
谢谢你的帮助!非常感激。你对DirectX有什么看法?似乎有一个很好的C# SDK,www.sharpdx.org。 - Mike
2
做了一些额外的研究。OpenCL非常有趣,因为它还支持Xeon Phi和现代Intel CPU集成的GPU,请参见此处http://bit.ly/Ta29ab。 - Mike
1
@Mike:关于DirectCompute,你可以不用考虑它,转而使用更高级别的API,如C++ AMP,但前提是你的计算必须且将永远在Windows上运行。否则,请使用标准API,例如OpenCL,这样如果你需要在Linux集群上运行代码,它将使您的代码具有未来可扩展性。 - Pragmateek

6
我不确定是否使用GPU是一个好的选择,因为需要从内存中检索“largerFloatingPointArray”的值。我的理解是,GPU更适合自包含计算。
我认为将这个单进程应用程序转变为在许多系统上运行的分布式应用程序,并调整算法,应该可以显著加快速度,具体取决于可用的系统数量。
您可以使用经典的“分而治之”方法。我会采取以下一般方法:
使用一个系统对“largeFloatingPointArray”进行预处理,生成哈希表或数据库。这将在单次通过中完成。它将使用浮点值作为键,数组中出现的次数作为值。最坏情况是每个值仅出现一次,但这是不太可能的。如果每次运行应用程序时都会更改largeFloatingPointArray,则在内存中使用哈希表是有意义的。如果它是静态的,则可以将表保存在诸如Berkeley DB之类的键值数据库中。我们称之为“查找”系统。
在另一个系统上,我们称之为“主”系统,创建工作块并将工作项“分散”到N个系统,并在结果可用时“收集”结果。例如,工作项可以简单地表示两个数字,指示系统应处理的范围。当系统完成工作时,它发送回发生次数的数组,并准备好处理另一个工作块。
性能得到改善,因为我们不会一直迭代largeFloatingPointArray。如果查找系统成为瓶颈,则可以在需要时将其复制到尽可能多的系统上。
通过足够数量的并行工作系统,应该可以将处理时间缩短到几分钟。
我正在开发一个针对基于多核的系统(通常称为微型服务器)的C并行编程编译器,该系统将使用多个“系统级芯片”模块构建。 ARM模块供应商包括Calxeda,AMD,AMCC等。英特尔可能也会提供类似的产品。
我有一个编译器版本可用于此类应用程序。基于C函数原型的编译器生成实现跨系统的进程间通信代码(IPC)的C网络代码。可用的IPC机制之一是套接字/ TCP / IP。
如果您需要帮助实施分布式解决方案,我很乐意与您讨论。
添加于2012年11月16日。
我思考了一下算法,我认为这应该可以在单次通过中完成。它是用C编写的,与您所拥有的相比应该非常快。
/*
 * Convert the X range from 0f to 100f in steps of 0.0001f
 * into a range of integers 0 to 1 + (100 * 10000) to use as an
 * index into an array.
 */

#define X_MAX           (1 + (100 * 10000))

/*
 * Number of floats in largeFloatingPointArray needs to be defined
 * below to be whatever your value is.
 */

#define LARGE_ARRAY_MAX (1000)

main()
{
    int j, y, *noOfOccurances;
    float *largeFloatingPointArray;

    /*
     * Allocate memory for largeFloatingPointArray and populate it.
     */

    largeFloatingPointArray = (float *)malloc(LARGE_ARRAY_MAX * sizeof(float));    
    if (largeFloatingPointArray == 0) {
        printf("out of memory\n");
        exit(1);
    }

    /*
     * Allocate memory to hold noOfOccurances. The index/10000 is the
     * the floating point number.  The contents is the count.
     *
     * E.g. noOfOccurances[12345] = 20, means 1.2345f occurs 20 times
     * in largeFloatingPointArray.
     */

    noOfOccurances = (int *)calloc(X_MAX, sizeof(int));
    if (noOfOccurances == 0) {  
        printf("out of memory\n");
        exit(1);
    }

    for (j = 0; j < LARGE_ARRAY_MAX; j++) {
        y = (int)(largeFloatingPointArray[j] * 10000);
        if (y >= 0 && y <= X_MAX) {
            noOfOccurances[y]++;
        }   
    }
}

3
工作可以在第二个时间中分配给机器网络,但在我看来,使用GPU的功率进行廉价(而且通常是巨大的)改进要好得多。至于您的框架,它与MPI相比如何? :) 注:IMHO 表示 “在我看来”;MPI 是一种并行计算库,用于在多台计算机之间共享工作负载。 - Pragmateek
感谢您提供的所有信息和C代码。也许我已经找到了一个解决我的问题的好方法:http://bit.ly/Ta4aSL [PDF]。听起来非常有前途...您认为呢? - Mike
Mike,这是一种有趣的方式,可以利用DirectX而不会被绑定到特定的GPU上。我在考虑任何副作用。当DirectX正在被积极使用时,是否会对其他应用程序渲染图形到显示器产生影响?尝试播放YouTube或Windows Media Player视频,有没有运行您的应用程序,并查看是否注意到正在播放的视频质量有所下降。此外,您是否知道将来是否需要超出工作站的能力范围?由于它都是Windows环境的一部分,我认为值得一试。 - Arun Taylor

6
除了上述帖子中的建议外,适当使用TPL(任务并行库)以在多个核心上并行运行。
上面的示例可以使用Parallel.Foreach和ConcurrentDictionary,但更复杂的映射-归约设置可以将数组拆分为块,每个块生成一个字典,然后将这些字典归约为单个字典,从而获得更好的结果。
我不知道您所有的计算是否正确地映射到GPU功能,但是无论如何,您都必须使用映射-归约算法将计算映射到GPU核心,然后将部分结果归约为单个结果,因此最好在移动到不太熟悉的平台之前在CPU上执行该操作。

3
谢谢您的建议。我已经在较高层次上使用TPL(任务并行库),这意味着我的应用程序并行调用多个方法,看起来效果不错(CPU使用率超过90%)。 - Mike

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接