[CUDA] [C] Суммирование на CUDA

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

ты только учти, что при такой адресации у тебя и проблемы с выравниваением будут и конфиликт банков, поэтому производительность скорее всего будет ниже чем на cpu

печально, если ниже... хотелось бы выше.

но пока хотя бы надо сделать просто чтобы работало %)

Меня этот пример вообще расстроил. Поменял BLOCK_SIZE с 256 на любое другое число (не кратное 256) - и всё(!) результат другой(!) сумма с CPU не сходится! Ну что за подстава?

Не нравится мне такой алгоритм... даже если он быстрый. зато криво работает.

это я к тому, что тот массив, которые мне нужно сложить - его длина ну никак не кратна 256...

Чтобы не было конфликтов, лучше сделать 1 блок, а массив поделить между его потоками. Результаты складывать в массиве разделяемой памяти - каждому потоку свою ячейку, а в конце один поток (например, нулевой) сложит все данные из разделяемой памяти.

Вот, например, такой быдлокод для нормализации вектора:

// нормализация вектора arr
__global__ void normalize_vec(float *arr, int arrsize){
        __shared__ float max[LBLKSZ];
        int idx = threadIdx.x;
        int blksize = (arrsize + blockDim.x - 1) / blockDim.x;
        int b_beg = idx * blksize;
        if(b_beg >= arrsize) return;
        int b_end = b_beg + blksize;
        if(b_end > arrsize) b_end = arrsize;
        int i; float *ptr = &arr[b_beg];
        float mm = *ptr++;
        for(i = b_beg +1 ; i < b_end; i++, ptr++)
                if(mm < *ptr) mm = *ptr;
        max[idx] = mm;
        __syncthreads();
        if(idx == 0){
                mm = max[0];
                for(i = 1; i < LBLKSZ; i++)
                        if(mm < max[i]) mm = max[i];
                max[0] = mm;
        }
        __syncthreads();
        ptr = &arr[b_beg];
        mm = max[0];
        for(i = b_beg ; i < b_end; i++, ptr++)
                *ptr /= mm;
}

http://www.books.ru/shop/books/752123

Примеров из nvidiaSDK + их брошюрок вполне достаточно же.

вот в том-то и дело, что там примеры тоже не очень... %) надо самому много додумывать... чем сейчас и занимаюсь.

спасибо за код... почитаю.

> Чтобы не было конфликтов, лучше сделать 1 блок, а массив поделить между его потоками.

Во-во... точно. Мне приходила в голову такая мысль. Только вот не догнал, почему в примере сделали по-другому... может, от этого производительность страдает?

Короче, попробую. что получитсяю

купить себе что ли... вроде не последний раз на CUDA пишу, наверное ещё много придётся. =)

Я бы лучше купил что по OpenCL на русском. Если бы нашёл. Достать труднее... Видел только на английском, и только за бугром где-то продавалась.

С приведенным выше примером я тоже сначала пробовал делать несколько блоков, но получилось медленнее. Плюс для хранения результатов вычислений каждого блока приходится заводить еще дополнительный массив.

вообще говоря, не очень книга, cuda programming guide для меня показался понятнее и информативнее

хотя я не углублялся в изучение, так что, быть может, я не прав

А как в сравнении GPU/CPU? медленнее или быстрее?

не критично (курсовая работа, думаю, смогу защитить, даже если медленнее будет работать, важно лишь объяснить, почему ;) ), но было бы конечно лучше, если бы алгоритм на GPU был быстрее... «ускорение» всё-таки.

> Чтобы не было конфликтов, лучше сделать 1 блок, а массив поделить между его потоками.

Там ограничение на 65535 нитей в блоке по-моему, да?.. то есть в моём случе (гораздо больше размер массива) мне нужно чтобы одна нить складывала большой кусок? :)

Смотря что за видеокарта и какая версия SDK. У меня, например, ограничение на 512 потоков. Естественно, складывать по кускам.

Книга кстати да не очень. По самой куде и программированию там где то половина. Вторая половина это набор статей на тему решения тех или иных задач с помощью куды.

http://www.nvidia.ru/object/cuda_books.html из этих книг что-нибудь на торентах по OpenCL достать можно ?

Вроде что-то такое получилось...

Кстати, очень много головной боли было из-за того, что применял double, он очень криво к float приводится... то ли с переполнениями, то ли ещё с чем-то...

На всякий случай - если у кого тоже будут какие-то немыслимые глюки вылазить.

Результат не на 100% совпадает с CPU-шным, но уже очень близко... подлатаю потом %) на утро. устал... Должен по идее ведь на 100% совпадать.

И хотелось бы всё-таки чтобы производительность была хоть чуть быстрее ;)

У меня тоже 512 потоков... но это warpSize, он чем-то отличается от blockDim, я не совсем понял, чем...

То есть наверное blockDim на любом железке может быть хоть 65535, но эффекта от этого немного =) только в том, что код одинаковый.

Надо не массив индексов, а массив флагов. Перемножаешь массив флагов на массив значений и затем складываешь.

стоп-стоп-стоп... мысль интересная, но я её не догнал...

массив каких флагов?

в смысле, массив такое же длины, как и исходный, состоящий из 1 и 0? можно подробнее объяснить, как мне это поможет?..

понял, в чём идея... вообще и правда изящное решение... позволяет при суммировании избежать произвольного доступа к памяти, только последовательный.

как я сразу не додумался?

Именно это я и имел в виду :)

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

http://nigma.ru/?s=cuda reduction&t=web

Первая же ссылка выводит на сайт Nvidia с презентацией. Презентация хороша, конечно... да вот только рассматривается там какой-то сферический случай в вакууме - когда исходный массив целиком влазит в разделяемую память. В моём случае этого быть никак не может. И вообще, у меня массив трёхмерный (в частном случае 78x78x78, но мне нельзя на него опираться, размеры могут быть и другими, примерно такого порядка).

Из школьного курса математики можно подсчитать, что если использовать float, который требует 4 байта - то нужно 78*78*78*4=1898208 байт, или ~1854 килобайта, что никак не вяжется в 16 килобайт, или сколько там =)

А разбить один блок по нитям попробую... чтобы каждая нить суммировала свой участок, а потом нулевая нить результаты складывала... напоминает алгоритм, который я писал на MPI когда-то... только у меня подозрение, что CUDA работает по-другому :)

На CPU будет медленнее, потому что в этом случае нужно не тупо сложить весь вектор, а посчитать преобразование фурье, то есть посчитать значение функции в каждой точке и умножить на экспоненту (точнее, на косинус для вещественной части и синус для мнимой части) некоторого числа, которое тоже надо считать «на ходу».

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

Так оно обычно и делается.

посчитать преобразование фурье

В CUDA для этого есть встроенная функция.

Прямо в самом cuda, или в библиотеке cudafft?

cudafft

Так оно обычно и делается.

Так наверное и сделаю. Вариант из доков nvidia мне не понравился... слишком идеализированный.

Ещё вопрос как лучше сделать... Через многомерный массив (трёхмерный) или одномерный? Сейчас сделал «развёртку» трёхмерного массива размера 78^3, соответственно индекс текущего элемента считается как

index = ix + dimx*iy + dimx*dimy*iz

где соответсвенно ix,iy,iz - индексы трёхмерного массива, dimx,dimy - размерность (ну она одинаковая), dimz тут не нужен как таковой.

Все примеры в cuda приводятся с использованием только x-компонента индекса нити (или блока), y и z обычно не используются. Несмотря на то, что иерархия нитей и блоков реально трёхмерная. Мне лучше не париться и использовать одномерную, или я подумал, что удобнее трёхмерную?

Тогда, например, у меня будет gridSize = dim3(1,1,1) и blockSize = dim3(10,10,10), (например, по 10 нитей по каждой координате достаточно? или больше сделать)

А размер массива вот у меня 78x78x78. Мне лучше его увеличивать засчёт нулевых элементов до 80x80x80 (чтобы нити равномерно его считали и не выходили за границы массива, и чтобы код для них писать одинаковый).

Или лучше вообще использовать степени двойки (например, 16x16x16 или 8x8x8)? где-то читал, что так должно быть эффективнее... размер блока выбирать именно такой.

Больше 16x16x16 разделяемой (shared memory) памяти уже не хватит, я подсчитал... её 16 килобайт вроде на всех устройствах...

Я почитаю доки, но у меня подозрение, что применить библиотечную функцию у меня не получится... случай специфический, в формуле я не разобрался

вот под эту не подходит http://upload.wikimedia.org/wikipedia/ru/math/3/b/f/3bfd17b69626915a2de012bd3...

ибо там нет 1/2pi, и в экспоненте не переменная x, а функция, от неё зависящая... и вообще функция двумерная.

может, и можно как-то, но думаю, лучше руками =) просто интеграл сосчитать.

Для больших массивов возможна проблема с размером блока, если использовать трехмерные структуры. Но, если максимальный размер вписывается, почему бы и нет.

Вы формулу приведите, по которой данные надо будет преобразовывать. Может быть и получится через обычное FFT реализовать (это, все-таки, быстрее).

http://dl.dropbox.com/u/17976101/equation.render.png вот такая вот страшная формула.

я не уверен, что это ПФ, мне так научник сказал...

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

под экспонентой скалярное произведение двух векторов. вектора зависят от углов theta и phi.

если надо, могу это зависимость тоже записать в виде формулы... всё это у меня есть (потому и думал о том, чтобы так сосчитать, не используя библиотек, мне почему-то казалось проще)

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

По ссылке ничего нет. А насчет формулы, сдается мне, что это таки Фурье. А мощность из спектра Фурье получить несложно.

только что с дропбоксом что-то было не так - сейчас прошёл по ссылке - всё норм... попробуйте ещё раз.

на всякий случай выложу ещё сюда - http://narod.ru/disk/40690581001/equation.render.png.html

Ан нет, не Фурье... Знакомая формула, не решение Уравнения Шрēдингера, часом?

Внутренний интеграл - дифференциальное сечение рассеяния (электрона, но вроде и необязательно именно его, думаю, если вместо e и m другой заряд/массу подставить), ну а внешний - полное сечение.

Кстати, о птичках. Обнаружил ещё одну интересную библиотеку curand... вообще уже написал код, который использует boost для генерации случайных чисел (думал его и использовать), но раз уж так, не использовать ли мне curand... Наверное, он будет эффективнее. Попробую и то, и то.

Просто мне тяжело представить, как генерацию псевдослучайных чисел можно распараллелить... я читал одну книжку про это (старую, советскую ещё), там каждое следующее число напрямую зависит от предыдущего...

Впрочем, я уверен, что такие алгоритмы устарели, да...