¿Puedo / debo ejecutar este código en una GPU?

ACTUALIZACIÓN Versión de 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
    }
}

Este solo lo probé para entradas más pequeñas, porque estoy probando mi computadora portátil. Sin embargo, funcionó. Sin embargo, es necesario hacer más testículos.

ACTUALIZACIÓN Versión secuencial

Acabo de hacer esta versión ingenua que realiza su algoritmo para 30,000,000 en menos de 20 segundos (ya cuenta la función para generar datos).

Básicamente, ordena tu matriz de flotadores. Viajará sobre el matriz ordenada, analizando el número de veces que un valor aparece consecutivamente en la matriz y luego coloque este valor en un diccionario junto con el número de veces que aparece.

Puede usar mapa ordenado, en lugar del unordered_map que usé.

Aquí está el código:

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

Si tiene la biblioteca de empuje instalada en su máquina, debe usar esto:

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

En lugar de esto

sort(data, data + size);

Seguro que será más rápido.

Original Post

"Estoy trabajando en una aplicación estadística que tiene una gran matriz que contiene entre 10 y 30 millones de valores de coma flotante".

"Es posible (y tiene sentido) utilizar una GPU para acelerar los cálculos?"

Sí, lo es. Hace un mes puse una simulación Dinámica Molecular completamente en la GPU. Uno de los núcleos, que calcula la fuerza entre pares de partículas, recibe 6 matrices cada una con 500.000 dobles, un total de 3 Millones de dobles (22 MB).

Así que está planeando poner 30 Millones de puntos flotantes esto es aproximadamente 114 MB de memoria global, por lo que esto no es un problema, incluso mi computadora portátil tiene 250 MB.

El número de cálculo puede ser un problema en su caso? Basándome en mi experiencia con la Dinámica Molecular (MD), digo que no. La versión secuencial de MD tarda aproximadamente 25 horas en completarse, mientras que en la GPU tardó 45 minutos. Usted dijo que su aplicación tomó un par de horas, también basado en su ejemplo de código se ve más suave que el Dinámica Molecular.

Aquí está el ejemplo de cálculo de fuerza:

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

Un ejemplo simple de un código en Cuda podría ser la suma de dos matrices 2D:

En c:

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

En 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];
}

En Cuda que básicamente tomó cada iteración y dividir por cada hilo,

1) threadIdx.x + blockIdx.x*blockDim.x;

Cada bloque tiene un Id de 0 a N-1 (N el número máximo de bloques) y cada bloque tiene un número X de hilos con un id de 0 a X-1.

1) Te da el para iteración que cada subproceso calculará en función de su id y el id de bloque en el que se encuentra el subproceso, el blockDim.x es el número de hilos que tiene un bloque.

Así que si tienes 2 bloques cada uno con 10 hilos y un N = 40, el:

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

Mirando su código hice este borrador de lo que podría ser en 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();
}

Tienes que usar atomicAdd porque diferentes hilos de diferentes bloques pueden escribir / leer noOfOccurrences al mismo tiempo, por lo que tienes que inseguro mutuo exclusion.

Este es solo un enfoque que incluso puede dar las iteraciones del bucle externo a los hilos en lugar de los bloques.

Tutoriales

La serie Dr Dobbs Journal CUDA: Supercomputing for the masses de Rob Farmer es excelente y cubre casi todo en sus catorce entregas. También comienza con bastante suavidad y, por lo tanto, es bastante amigable para principiantes.

Y otros:

• Desarrollo Con CUDA-Introducción
• Volumen I: Introducción a la programación CUDA
• Primeros pasos con CUDA
• Lista de Recursos CUDA

Eche un vistazo al último elemento, encontrará muchos enlaces para aprender CUDA.

OpenCL: OpenCL Tutorials / MacResearch