CUDA search

Всем доброго времени суток!
Почти через месяц защита диссертации. Разбираться во всех дебрях программирования на CUDA уже не осталось времени. Нужна программа, реализующая алгоритм поиска на GPU. Какой именно алгоритм это тоже большой вопрос. Алгоритмов поиска не мало и какой из них лучше распараллелиться на GPU?
Может кто подскажет как реализовать? Буду очень рад любой помощи. Спасибо.

Пока набросал простой линейный поиск. замечаниям и отзывам тоже буду рад.

При создании программы использовались следующие средства и инструменты:
1. Видеокарта из серии nVidia GeForce 8xxx/9xxx или более современная.
2. CUDA Toolkit v.2.3 (скачать можно здесь: www.nvidia.ru/object/cuda_get_ru.html)
3. CUDA SDK v.2.3 (скачать можно там же где Toolkit)
4. Visual Studio 2008
5. CUDA Visual Studio Wizard (скачать можно здесь: sourceforge.net/projects/cudavswizard/)

Текст программы

  1. #include <stdio.h>
  2. #include <time.h>
  3. #define CPU 1
  4. #define GPU 2
  5. // функция выполняемая на центральном процессоре
  6. __host__ void searchCPU ( float * data, int * result )
  7. {
  8.         int len = 16 * 1024 * 1024;    
  9.  
  10.         for ( int j = 0; j < len; j++ )
  11.         {
  12.                 if (data[j] == 2.0f )
  13.                 {
  14.                         result[0] = data[j];
  15.                         result[1] = j;
  16.                 }
  17.         }
  18. }
  19. // функция выполняемая на графическои процессоре
  20. __global__ void searchGPU ( float * data, int * result )
  21. {
  22.    int idx = blockIdx.x * blockDim.x + threadIdx.x;
  23.    //data [idx] = data [idx] + 1.0f;
  24.    if ( data [idx] == 2.0f )
  25.    {
  26.         result[0] = data [idx];
  27.         result[1] = idx;
  28.    }
  29. }
  30.  
  31. int main ( int argc, char *  argv [] )
  32. {
  33.   int n        = 16 * 1024 * 1024;
  34.   int numBytes = n * sizeof ( float );
  35.   int size_int = 2 * sizeof ( int );
  36.         //timer for host
  37.   clock_t start, finish;
  38.   double  duration;
  39.              // allocate host memory
  40.   float * a = new float [n];
  41.   int * host_res = new int[2];
  42.         //init host variables
  43.   for ( int i = 0; i < n-1; i++ )
  44.     a [i] = 0.0f;
  45.     a [n-1] = 2.0f;
  46.   host_res [0] = 0;
  47.   host_res [1] = 0;
  48.         //Выбираем способ расчета
  49.   printf("Select compute mode: 1 - CPU, 2 - GPU\n");
  50.   int mode;
  51.   scanf("%i", &mode);
  52.   if ( mode == CPU )
  53.   {     start = clock(); //start timer
  54.         searchCPU (a, host_res);
  55.         finish = clock(); //stop timer
  56.         duration = (double)(finish - start) / CLOCKS_PER_SEC; //time in seconds
  57.         printf( "time spent executing by the CPU: %.5f seconds\n", duration );
  58.   } else
  59.   if ( mode == GPU )
  60.   {      // allocate device memory
  61.     float * dev = NULL;
  62.         int * res = NULL;
  63.    
  64.     cudaMalloc ( (void**)&dev, numBytes );
  65.     cudaMalloc ( (void**)&res, size_int );
  66.  
  67.                    // set kernel launch configuration
  68.     dim3 threads = dim3(512, 1);
  69.     dim3 blocks  = dim3(n / threads.x, 1);
  70.                   // create cuda event handles
  71.     cudaEvent_t start, stop;
  72.     float gpuTime = 0.0f;
  73.     cudaEventCreate ( &start );
  74.     cudaEventCreate ( &stop );
  75.                  // asynchronously issue work to the GPU (all to stream 0)
  76.     cudaEventRecord ( start, 0 );
  77.         //copy data to device
  78.     cudaMemcpy      ( dev, a, numBytes, cudaMemcpyHostToDevice );
  79.     cudaMemcpy      ( res, host_res, size_int, cudaMemcpyHostToDevice );
  80.                //call device function
  81.     searchGPU<<<blocks, threads>>>(dev,res);
  82. //copy data to host
  83.     cudaMemcpy      ( host_res, res, size_int, cudaMemcpyDeviceToHost );       
  84.     cudaEventRecord ( stop, 0 ); // stop timer for device
  85.     cudaEventSynchronize ( stop );
  86.     cudaEventElapsedTime ( &gpuTime, start, stop ); //gpuTime - time in miliseconds
  87.                // print the gpu times
  88.     printf("time spent executing by the GPU: %.5f seconds\n", gpuTime/1000 );
  89.  // release resources (device)
  90.     cudaEventDestroy ( start );
  91.     cudaEventDestroy ( stop  );
  92.     cudaFree         ( dev   );
  93.   }
  94.   else
  95.   {  printf("not correct input");  return 1;   }
  96.      //print result
  97.      printf ( "value %i, index %i",host_res[0], host_res[1]);
  98. // release resources (host)
  99.     delete a;
  100.     delete host_res;
  101.     return 0;
  102. }

Время выполнения программы на ЦП составило 78 миллисекунд.
При выполнение на графическом процессоре, время выполнения программы составило 47 миллисекунд. Как видим прирост производительности порядка 40%. Не так уж много. Однако время измерения выполнения на графическом процессоре измерялось с учетом копирования исходных данных в память видеокарты и результата обратно в ОЗУ.
Результат выполнения программы на GPU без учета копирования данных составил около 9 миллисекунд. Таким образом, львиная доля времени используется на копирование данных в память GPU и обратно. Но, несмотря на это, все же достигается неплохой прирост производительности.

Программа выполнялась на платформе со следующими характеристиками:

Windows XP Professional (5.1, Build 2600) Service Pack 3
Processor: Intel(R) Pentium(R) 4 CPU 3.00GHz (2 CPUs)
Memory: 1024MB RAM

Display Devices
Name GeForce 8600 GT
Compute Capability 1.1
Clock Rate 1188 MHz
Multiprocessors 4
Warp Size 32
Regs Per Block 8192
Threads Per Block 512
Watchdog Enabled Yes
Threads Dimentions 512 x 512 x 64
Grid Dimentions 65535 x 65535 x 1

Memory InformationTotal Global 511.688 MB
Shared Per Block 16 KB
Pitch 256 KB
Total Constant 64 KB
Texture Alignment 256
GPU Overlap Yes

Performance InformationMemory Copy
Host Pinned to Device 3059.73 MB/s
Host Pageable to Device 1600.54 MB/s
Device to Host Pinned 3062.96 MB/s
Device to Host Pageable 1651.8 MB/s
Device to Device 5295.5 MB/s
GPU Core Performance
Single-precision Float 75641.8 Mflop/s
Double-precision Float Not Supported
32-bit Integer 15184 Miop/s
24-bit Integer 75647.5 Miop/s

Спасибо. Жду ответов и комментариев.

Forums: 

Попробуйте распаралелить

Попробуйте распаралелить алгоритм Бойера-Мура, он и быстрый и некоторые его части хорошо поддадутся параллелизму(его эвристики). Вам необходимо реализовать поиск одной строки в другой при том что у вас по одной копии шаблона и главной строки или их может быть больше одной?

Достаточно будет и по одной

Достаточно будет и по одной строке.
Алгоритмы уже смотрел, у меня возникает трудность имено с реализацией на CUDA.
Еще очень хотелось бы реализовать бинарный поиск. Я так понимаю нужно разбивать исходные данные на части и в каждой из части параллельно осуществлять бинарный поиск. Но опять таки, вопрос, а как это правильно сделать на CUDA?

Спасибо

Бинарный поиск - это же по

Бинарный поиск - это же по сортированым данным? Он сходится за логарифмическое время, если у вас искомое одно (а не много параллельных), то от распараллеливания значимого выигрыша не будет, все уйдет в загрузку данных из случайных мест.

Что касается линейного поиска - то в SDK есть пример scan (хотя я туда не смотрел)

scan из SDK решает совсем

scan из SDK решает совсем другую задачу. В простейшем случае вычисляет по ряду a ряд b по формуле
i
b[i] = SIGMA a[j]
j = 0

Эквивалентный C++ код будет примерно таким

  1. for (int i=0; i<n; ++i)
  2. {
  3.     b[i] = 0;
  4.     for (int j=0; j<i; ++j)    
  5.          b[i]  += a[j];
  6. }

Естествеенно на куда это сложнее, но зато выполняется паралельно, а не производится последовательое накопление результата в b[i].

По поводу паралельного

По поводу паралельного поиска.
Я делаю так:
Пусть имеется массив a в котором ищутся элементы из b.
1. Создаем массив c объединяющий a и b. Затраты O(n).
2. Сортируем c. Затраты O(n log n).
3. Проходим по c сканом и таким образом к каждому b находим ближайший (больший или меньший) a. Затраты O(n).

Т.е. нужно ореинтироваться на обработку на gpu именно МАССОВЫХ запросов. Только так может быть получен выигрышь. А бинарный поиск - это тупик.

Уточнение. По поводу

Уточнение.
По поводу паралельного поиска.
Я делаю так:
Пусть имеется массив a размера n в котором ищутся элементы из b размера m.
1. Создаем массив c объединяющий a и b. Затраты O(n+m).
2. Сортируем c. Затраты O((n+m) log (n+m)).
3. Проходим по c сканом и таким образом к каждому b находим ближайший (больший или меньший) a. Затраты O(n+m).

Т.е. нужно ореинтироваться на обработку на gpu именно МАССОВЫХ запросов. Только так может быть получен выигрышь. А бинарный поиск - это тупик.

Бинарный поиск может и тупик,

Бинарный поиск может и тупик, но на тройке гигабайт данных сотня тысяч лукапов в секунду получается.

А ваш подход имеет смысл, если b известен заранее и в a нету дублей. Так бывает, конечно.

Можно чуть модифициоровать

Можно чуть модифициоровать алгоритм. Допустим a известен и уже отсортирован.
Тогда b можно отдельно накопить и сепаратно отсортировать.
После a и b объеденить в c и после применить к c последнию ступень битонической сортировки (это будет быстрее чем полная сортировка c состоящего из неотсортированных a и b).
Преимущества? b может быть небольшим, значительно меньшим a.

Это сработает примерно со

Это сработает примерно со скоростью O(n), ибо без полного линейного просмотра a не обойдется. Параллельность я не рассмативаю, вы в этом последнем проходе сортировки с неизбежностью упретесь именно в bandwith памяти.

Т.е. выгодно становится если log(n)*m больше чем n (понятно что с изрядными коэффициентами, случайный доступ при бинарном поиске много дороже, чем последовательный доступ при слиянии).

Но мне кажется интересным такой вариант (при сравнительно небольших m): первые фазы бинарного поиска делаются делением на мультипроцессоры: допустим их 30, тогда каждый отвечает только за 1/30-ю входного массива А и за одну итерацию определяет с какими входными диапазонами работает данный конкретный мультипроцессор, после чего отвергает входное данное из b сразу (двумя сравнениями) не выполняя поиск.

>Это сработает примерно со

>Это сработает примерно со скоростью O(n)
В расчете на все m элементов из b. Поэтому b не должен быть слишком маленким. Лично у меня m>>n.

>Но мне кажется интересным такой вариант (при сравнительно небольших m): первые фазы бинарного поиска делаются делением на мультипроцессоры:

Я думал над таким вариантом; это не очень выгодно, дело в том что уже при первой проверке в бинарном поиске отсекается половина элементов.
Нет, все же лучше работать с массовыми запросами.
Если в программе это не так, значит нужно снова продумать архитектуру программы и все же добится перехода от одиночных запросов к массовым, хотя бы и буферизацией.

Ну как

Ну как выгодно-невыгодно.
Если мультипроцессоров 32 (круглое число), то две проверки на больше-меньше (быстрые) отсекают пять лукапов. Не много и не мало, вот столько вот.

Но вы правы в том, что у всех свои задачи, то с чем имею дело я - это обычно маленькие m.

если число элементов 500м =

если число элементов 500м = 2^25 то число уровней сбалансированного бинарного дерева (или в бинарном поиске) будет 25, сокращение 5 из них даст прирост производительности от 32 процессоров по сравнению с одним на 100*5/25=20%.
Это действительно очень мало. И это если еще забыть про проблемы с памятью. Пикапать память в куда не есть гут. Это правда.

Извиняюсь, 500м =

Извиняюсь,
500м = 2^29
100*5/29=17,24%.
Т.е. на само деле еще хуже.

Не, ну 5 бит (лукапов)

Не, ну 5 бит (лукапов) достаются просто на халяву (в предположении о 32 SP, что не так, их может быть или 30 или 48, или меньше, я так просто прикидываю).

Но у нас очень иерархическая память, еще есть регистры и shmem. Предполагая, что половину регистров можно потратить на всякую фигню, у нас 32 килобайт регистров и 32 килобайта shmem на CUDA 1.x (кроме 1.0) - это 16к 4-байтных чисел т.е. записав туда индекс (пары значение-смещение) можно еще 13 лукапов поэкономить. На Fermi всего вдвое больше - значит 14. Вместе с индексом на номере SP - 19.

При этом 500M влезет только на Tesla, реально 250M т.е. 2^28 и ускорились мы раза эдак в три примерно.

Если m относительно небольшие и многократные последовательные чтения из разных блоков - не слишком страшны (понятно что неэкономно, но дофига быстро), то можно сделать индекс на номере блока/треда.

С 500м это я конечно

С 500м это я конечно погорячился, даже и 250м не реально, реально 100м, поскольку элемент меньше 8 байт иметь не будет (4 на ключ, 4 на индекс/указатель). Но это ладно, гораздо интересней поиск.
Бинарный конечно не катит.
А если 32-арный? Т.е. данные представляются как 32-арное дерево, с 32 потомками каждого узла. Тогда каждый процессор (тред) может обрабатывать своего потомка, а результат всего варпа подаем на вход следующей итерации. Таким образом за 2 итерации можно выбрать из 1024 альтернатив, соответственно за 5 итераций можно выбрать из реально максимального 32м массива.
У этой идеи есть одна приятная особенность. Если выделить дополнительную память и в нее перезаписать каждый 32 элемент, то в каждой итерации поиска необходимые для работы данные можно подгружать в шеред память быстро по блокам.

Индексом-указателем может

Индексом-указателем может быть номер в массиве a[], тогда это бесплатно.

Но давайте на проблему с другой стороны посмотрим.
1) Данных много, полный индекс в быстрой памяти не построить. Что бы мы не делали, это будет вариация на тему бинарного поиска, log(n) лукапов на каждый входной элемент, чудес нет.
2) Каждый случайный доступ в a[] - долгий. Т.е. мы можем попрятать latency в то, что на SP будет много блоков, но больших чудес не будет, ограничивает в любом случае лукап в глобальной памяти.
3) Поэкономить эти лукапы можно двумя способами
3a) созданием части индекса в быстрой памяти каким-то образом
3б) группировкой значений b (в смысле последовательности обработки) таким образом, чтобы лукапы в a[] были совсем локальными и случайные чтения заменялись бы линейными (подряд)

3б имеет смысл при относительно больших m, когда полная сортировка слиянием еще невыгодна, но вероятность попасть в близкие куски a[] уже довольно большая (считая что позиционирование в 100 раз медленнее линейного чтения, m должна быть больше чем n/100)

Заметим, что задача очень похожа на типичный базоданновый поиск, только вместо диска у нас - медленная глобальная память (но суть та же - позиционирование медленное, чтение-запись быстрые), а вместо кэша в просто памяти - кэш в быстрой памяти (shared/регистры). Ничто не ново под луною.

>Индексом-указателем может

>Индексом-указателем может быть номер в массиве a[], тогда это бесплатно.

К сожалению это не так просто. При сортировке нужно хранить старое значение этого индекса.

>3б) группировкой значений b (в смысле последовательности обработки) таким образом, чтобы лукапы в a[] были совсем локальными и случайные чтения заменялись бы линейными (подряд)

Это собственно я и предложил, когда писал о вспомогательных массивах, в которых хранятся каждое 32 (32^2, 32^3, 32^4) значение исходного отсортированного.

>Заметим, что задача очень

>Заметим, что задача очень похожа на типичный базоданновый поиск, только вместо диска у нас - медленная глобальная память (но суть та же - позиционирование медленное, чтение-запись быстрые), а вместо кэша в просто памяти - кэш в быстрой памяти (shared/регистры). Ничто не ново под луною.

Как раз сегодня смотрел незабвенного Кнута, 3т. Там была описана работа с внешней памятью (магнитной лентой) при сортировке. Вопрос меня заинтересовал как раз в связи с проблемой глобальной - шеред памяти в куда. Но ничего особенного не вынес. У Кнута это все вариации на тему mergesort. Это мало пригодно для многопроцессорных систем.

Вообще, для линейного поиска

Вообще, для линейного поиска - удивительно что исполнение на GPU оказалось быстрее. Ведь скорость копирования данных на GPU очевидно меньше, чем скорость чтения памяти CPU.

Т.е. для такой простой задачи - если у вас получилась скорость последовательного скана сильно меньше гигабайта в секунду на GPU (16M 4-байтовых элементов, 78msec, но наверное искомый элемент не последний) - это какая-то фигня даже для Pentium4, очень мало.

На такой простой (с такой низкой арифметической интенсивностью) задаче - вы должны получить сильный проигрыш для GPU-решения.

P.S. На правах админа позволил себе поставить тег <code>

Искомы элемент последний

Искомы элемент последний (т.е. наихудший случай для линейного поиска), именно поэтому прирост производительности на GPU, если искомый элемент находиться в начале, тогда GPU сильно проигрывает.

P.S. Cпасибо за тег. я не частый обитатель форумов, подзабыл про то, что можно ставить теги :-)

CPU у вас сильно проигрывает

CPU у вас сильно проигрывает потому что сгенерирован сильно неоптимальный код.

Вы обязаны упираться в скорость линейного чтения из памяти, если это не так, смотрите что вам за код компилятор нагенерировал. Ну или ручками unroll тоже полезно попробовать.

(нет, понятно что CUDA-устройство можно использовать для offload, освободив CPU для чего-то полезного, вроде показа прогресс-бара, но по скорости выигрыша обязано бы не быть)

Искомы элемент последний

Искомы элемент последний (т.е. наихудший случай для линейного поиска), именно поэтому прирост производительности на GPU, если искомый элемент находиться в начале, тогда GPU сильно проигрывает.

P.S. Cпасибо за тег. я не частый обитатель форумов, подзабыл про то, что можно ставить теги :-)

У кого будут какие мнения по

У кого будут какие мнения по поводу следующей статьи?
http://www.usenix.org/event/hotpar09/tech/full_papers/kaldeway/kaldeway_html/