sorting - ¿Cómo usar Thrust para ordenar las filas de una matriz?

cuda

Question

Tengo una matriz de 5000x500 y quiero ordenar cada fila por separado con cuda. Puedo usar arrayfire pero esto es solo un bucle for sobre Thust :: sort, que no debería ser eficiente.

https://github.com/arrayfire/arrayfire/blob/devel/src/backend/cuda/kernel/sort.hpp

for(dim_type w = 0; w < val.dims[3]; w++) {
            dim_type valW = w * val.strides[3];
            for(dim_type z = 0; z < val.dims[2]; z++) {
                dim_type valWZ = valW + z * val.strides[2];
                for(dim_type y = 0; y < val.dims[1]; y++) {

                    dim_type valOffset = valWZ + y * val.strides[1];

                    if(isAscending) {
                        thrust::sort(val_ptr + valOffset, val_ptr + valOffset + val.dims[0]);
                    } else {
                        thrust::sort(val_ptr + valOffset, val_ptr + valOffset + val.dims[0],
                                     thrust::greater<T>());
                    }
                }
            }
        }

¿Hay alguna manera de fusionar las operaciones de empuje para que las clases se ejecuten en paralelo? De hecho, lo que estoy buscando es una forma genérica de fusionar las iteraciones de bucle.

Answer 1

Puedo pensar en 2 posibilidades, una de las cuales ya ha sido sugerida por @JaredHoberock. No conozco una metodología general para fusionar iteraciones for-loop en empuje, pero el segundo método es el enfoque más general. Supongo que el primer método sería el más rápido de los dos enfoques, en este caso.

1. Utiliza una ordenación vectorizada. Si las regiones que se ordenarán por sus bucles for anidados no se superponen, puede hacer una ordenación vectorizada utilizando 2 operaciones de ordenación estable consecutivas como se describe here .

2. Thrust v1.8 (disponible con CUDA 7 RC, o mediante descarga directa desde el repositorio de empuje github incluye soporte para anidar algoritmos de empuje , al incluir una llamada de algoritmo de empuje dentro de un functor personalizado pasado a otro algoritmo de empuje. Si usa el thrust::for_each operación para seleccionar los thrust::for_each individuales que necesita realizar, puede ejecutar esos ordenamientos con una sola llamada de algoritmo de empuje, incluyendo la operación thrust::sort en el functor que pasa a thrust::for_each .

Aquí hay una comparación completamente trabajada entre 3 métodos:

1. el método original de ordenar en un bucle
2. vectorizado / clasificación por lotes
3. ordenamiento anidado

En cada caso, estamos clasificando los mismos 16000 conjuntos de 1000 ints cada uno.

$ cat t617.cu
#include <thrust/device_vector.h>
#include <thrust/device_ptr.h>
#include <thrust/host_vector.h>
#include <thrust/sort.h>
#include <thrust/execution_policy.h>
#include <thrust/generate.h>
#include <thrust/equal.h>
#include <thrust/sequence.h>
#include <thrust/for_each.h>
#include <iostream>
#include <stdlib.h>

#define NSORTS 16000
#define DSIZE 1000

int my_mod_start = 0;
int my_mod(){
  return (my_mod_start++)/DSIZE;
}

bool validate(thrust::device_vector<int> &d1, thrust::device_vector<int> &d2){
  return thrust::equal(d1.begin(), d1.end(), d2.begin());
}


struct sort_functor
{
  thrust::device_ptr<int> data;
  int dsize;
  __host__ __device__
  void operator()(int start_idx)
  {
    thrust::sort(thrust::device, data+(dsize*start_idx), data+(dsize*(start_idx+1)));
  }
};



#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

int main(){
  cudaDeviceSetLimit(cudaLimitMallocHeapSize, (16*DSIZE*NSORTS));
  thrust::host_vector<int> h_data(DSIZE*NSORTS);
  thrust::generate(h_data.begin(), h_data.end(), rand);
  thrust::device_vector<int> d_data = h_data;

  // first time a loop
  thrust::device_vector<int> d_result1 = d_data;
  thrust::device_ptr<int> r1ptr = thrust::device_pointer_cast<int>(d_result1.data());
  unsigned long long mytime = dtime_usec(0);
  for (int i = 0; i < NSORTS; i++)
    thrust::sort(r1ptr+(i*DSIZE), r1ptr+((i+1)*DSIZE));
  cudaDeviceSynchronize();
  mytime = dtime_usec(mytime);
  std::cout << "loop time: " << mytime/(float)USECPSEC << "s" << std::endl;

  //vectorized sort
  thrust::device_vector<int> d_result2 = d_data;
  thrust::host_vector<int> h_segments(DSIZE*NSORTS);
  thrust::generate(h_segments.begin(), h_segments.end(), my_mod);
  thrust::device_vector<int> d_segments = h_segments;
  mytime = dtime_usec(0);
  thrust::stable_sort_by_key(d_result2.begin(), d_result2.end(), d_segments.begin());
  thrust::stable_sort_by_key(d_segments.begin(), d_segments.end(), d_result2.begin());
  cudaDeviceSynchronize();
  mytime = dtime_usec(mytime);
  std::cout << "vectorized time: " << mytime/(float)USECPSEC << "s" << std::endl;
  if (!validate(d_result1, d_result2)) std::cout << "mismatch 1!" << std::endl;
  //nested sort
  thrust::device_vector<int> d_result3 = d_data;
  sort_functor f = {d_result3.data(), DSIZE};
  thrust::device_vector<int> idxs(NSORTS);
  thrust::sequence(idxs.begin(), idxs.end());
  mytime = dtime_usec(0);
  thrust::for_each(idxs.begin(), idxs.end(), f);
  cudaDeviceSynchronize();
  mytime = dtime_usec(mytime);
  std::cout << "nested time: " << mytime/(float)USECPSEC << "s" << std::endl;
  if (!validate(d_result1, d_result3)) std::cout << "mismatch 2!" << std::endl;
  return 0;
}
$ nvcc -arch=sm_20 -std=c++11 -o t617 t617.cu
$ ./t617
loop time: 8.51577s
vectorized time: 0.068802s
nested time: 0.567959s
$

Notas:

1. Estos resultados variarán significativamente de GPU a GPU.
2. El tiempo / método "anidado" puede variar significativamente en una GPU que pueda soportar paralelismo dinámico, ya que esto afectará la forma en que el empuje ejecuta las funciones de ordenamiento anidado. Para probar con paralelismo dinámico, cambie los -arch=sm_20 compilación de -arch=sm_20 a -arch=sm_35 -rdc=true -lcudadevrt
3. Este código requiere CUDA 7 RC. Usé Fedora 20.
4. El método de clasificación anidada también se asignará desde el lado del dispositivo, por lo tanto, debemos aumentar sustancialmente el montón de asignación de dispositivos usando cudaDeviceSetLimit .
5. Si está utilizando paralelismo dinámico, y dependiendo del tipo de GPU en el que se esté ejecutando, la cantidad de memoria reservada con cudaDeviceSetLimit puede necesitar un aumento tal vez por un factor adicional de 8.