М. В. Ломоносова 2015 Антонюк Валерий Алексеевич Программирование на видеокартах (gpgpu). Спецкурс кафедры мми. М. Физический факультет мгу им. М. В. Ломоносова, 2015. 48 с. Учебное пособие
 Скачать 1.54 Mb.
|
|
STL (например, thrust::sort() и std::sort() ). Существуют в библиотеке также так называемые fancy iterators ("воображаемые" итераторы), ведущие себя в алгоритмах как "настоящие" итераторы: constant_iterator (подобен итератору в бесконечном массиве, заполненном одной величиной) counting_iterator (подобен итератору в бесконечном массиве, заполненном последовательными величинами) transform_iterator (производит итерирование по последовательности, преобразованной с помощью заданной функции) zip_iterator (производит итерирование как бы по массиву структур, хотя принимает массивы отдельных величин, составляющих структуру) permutation_iterator 5 21 Для освоения этих возможностей полезно поэкспериментировать с некоторыми примерами, найденными в Сети. В интересной статье "odeint v2 — Solving ordinary differential equations in C++" ( http://www.codeproject.com/Articles/268589/ ) есть пример решения дифференциальных уравнений в подобном стиле (с использованием thrust::device_vector , thrust::for_each , thrust::get , thrust::make_zip_iterator ). Правда, полного текста программы там вроде бы нет. Пример monte_carlo.cu — это известный способ оценивания константы π методом Монте- Карло. В нём используются формирование равномерно распределённой случайной величины ( thrust::uniform_real_distribution ), редуцирующее преобразование thrust::transform_reduce() , а также упомянутый выше "считающий" итератор thrust::counting_iterator Интересный пример argmin_row_space.cu использования изначально векторной природы операций в Thrust для работы с матрицей (отыскание положения минимальных элементов в её строках) приведён в заметке: http://peterwittek.com/2013/04/argmin-on-the-rows-of-a-matrix-with-thrust/. Там используются thrust::device_vector , thrust::make_zip_iterator() , thrust::counting_iterator() , thrust::make_transform_iterator() Стоит попробовать приводимый в статье "A Brief Test on the Code Efficiency of CUDA and Thrust" ( http://www.codeproject.com/Articles/83757/ ) полезный пример thrustExample.cu (он располагается в прилагаемом к статье архиве с исходным кодом: http://www.codeproject.com/KB/Parallel_Programming/test-on-thrust- efficiency/thrustExample.zip ). Там сравниваются скорости вычисления суммы квадратов вектора просто в CUDA, с Thrust, а также на CPU. Пример вполне рабочий, правда, использует заголовочный файл cutil_inline.h , находящийся в CUDA SDK ( C:\Documents and Settings\All Users\Application Data\NVIDIA Corporation\NVIDIA GPU Computing SDK 4.1 — это обычное место его установки, оно также указано в специальной переменной окружения NVSDKCOMPUTE_ROOT ) в каталоге %NVSDKCOMPUTE_ROOT%\C\common\inc (для Visual Studio указывается как $(NVSDKCOMPUTE_ROOT)\C\common\inc ). Кроме того, пример использует для подсчёта времени средства Windows, поэтому может понадобиться включение заголовочного файла windows.h . Результаты, правда, получаются довольно странные... Библиотека cuBLAS — базовые функции линейной алгебры Эта библиотека содержит основные векторные, векторно-матричные и матричные операции и использует выверенный временем фортрановский код, созданный ранее в рамках библиотеки BLAS. Данная реализация использует оборудование видеокарты. Для того, чтобы ею воспользоваться, приложение должно выделить память для векторов и матриц на GPU, заполнить их, затем вызывать последовательность необходимых функций cuBLAS, после 22 чего результаты надо извлечь с GPU на CPU. Кроме того, предполагается создание (а в конце — разрушение) специального контекста cuBLAS. В библиотеке имеются функции как бы трёх уровней: уровень 1 — работающие только с векторами (скалярное произведение, нормы, AXPY : y[i] = a*x[i]+y[i] ), уровень 2 — работающие с векторами и матрицами (умножение вектора на матрицу общего вида GEMV , решение треугольной системы TRSV ), уровень 3 — работающие только с матрицами (умножение матриц GEMM , ...). Поскольку библиотека основана на старом фортрановском коде, в ней используется хранение матриц по столбцам (а не по строкам, как в языке C), а также индексирование не с нуля, а с единицы. Поддерживаются четыре типа величин в векторах и матрицах: float , double , complex , double complex (условно обозначаются буквами S , D , C , Z ). Именуются функции по схеме: cublas <буква типа данных><тип матриц (2 буквы)><операция>. Например, функция умножения матриц общего вида, содержащих float -величины, называется cublasSgemm (имя составлено из cublas+S+ge+mm ). Возможные типы матриц: ge обычная матрица gb ленточная матрица sy симметричная матрица sp симметричная упакованная матрица sb симметричная ленточная упакованная матрица he эрмитова матрица hp эрмитова упакованная матрица hb эрмитова ленточная матрица tr треугольная матрица tp треугольная упакованная матрица tb треугольная ленточная матрица В примере cublas1Test.cpp мы видим создание контекста ( cublasInit() ), выделение памяти на графической карте ( cublasAlloc() ), копирование информации туда ( cublasSetMatrix() ), вызовы функций cuBLAS ( cublasSscal() в рамках вспомогательной функции modify() ), извлечение результата ( cublasGetMatrix() ), освобождение памяти на карте ( cublasFree() ) и окончательное завершение работы с библиотекой ( cublasShutdown() ). Для доступа к содержимому матриц, хранимых в "фортрановском" виде, уже нельзя применять операции, принятые в языке C для двумерных массивов (т.е., ИмяМатрицы[Индекс][Индекс]), используется индексация одномерного массива, учитывающая расположение элементов по столбцам, с помощью макроопределения: #define IDX2F(i,j,ld) ((((j)-1)*(ld))+((i)-1)) Здесь ld — это число строк в матрице (т.н. leading dimension). Формируется индекс положения элемента в памяти, имитирующий двумерное обращение к "фортрановской" матрице. j — номер столбца (начинается с единицы), i — номер строки (тоже начинается с единицы). 23 Доступ к содержимому матриц, объявленных в C-программе, может быть получен с помощью макроопределения: #define IDX2C(i,j,ld) (((i)*(ld))+(j)) Здесь ld — число столбцов этой матрицы (leading dimension), i , j — номера строки и столбца; оба начинают отсчитываться с нуля. Пример matrix_product_c_view.cu (см. приложение) демонстрирует новую возможность работы с библиотекой (обратите внимание на использование другого заголовочного файла: cublas_v2.h вместо cublas.h ). В этом случае при инициализации cuBLAS создаётся дескриптор (используется вызов cublasCreate(&handle) ), поскольку далее этот дескриптор необходим при вызове функций библиотеки (в данном случае — cublasSgemm() ; ей передаются указатели на память видеокарты (память GPU), полученные с помощью thrust::device_vector и thrust::raw_pointer_cast() ), а результат (вектор F ) прямо из памяти карты выдаётся в стандартный вывод как содержимое матрицы языка C. Завершается всё прекращением работы с контекстом при помощи вызова функции cublasDestroy(handle) Библиотека cuSPARSE — разреженные матрицы Эта библиотека предназначена для работы с разреженными матрицами, где число ненулевых элементов мало. Используемые форматы хранения разреженных матриц: COO, CSR, CSC, ELL, HYB. Последовательность работы с ней — примерно такая же, как и с cuBLAS: создание контекста, выделение ресурсов на GPU, работа с функциями библиотеки, копирование результатов, освобождение ресурсов GPU, освобождение контекста, однако для работы с матрицами нужно также с помощью специальных функций создавать дескрипторы матриц. Поскольку способов хранения у разреженных матриц много, имеет смысл подробнее познакомиться с ними, прежде чем рассматривать примеры работы с библиотекой. 24 Библиотека cuSPARSE Общие сведения Содержит базовые функции работы с разреженными матрицами, которые условно можно разбить на 4 части: 1) операции между разреженным вектором и обычным; 2) операции между разреженной матрицей и обычным вектором; 3) операции между разреженной матрицей и обычной; 4) функции преобразования между различными форматами представления матриц. Остальные действия, необходимые при программировании на GPU (такие как выделение и освобождение памяти, перенос данных в память GPU и обратно), должны выполняться разработчиком. Схема поименования функций библиотеки cuSPARSE для работы с матрицами и векторами аналогична поименованию функций в cuBLAS — с небольшими отличиями: cusparse <ТипД>[<Формат>]<Операция>[<ВыходнойФормат>] где <ТипД> — это буква (может быть S , D , C , Z , или X , что соответствует типу данных float , double , cuComplex , cuDoubleComplex или "обобщённому" типу), <Формат> — dense , coo , csr , csc , или hyb , что соответствует обычным "плотным" матрицам, а также хранящимся в форматах COO (координатный), CSR (компрессированный с разреженными строками), CSC (компрессированный с разреженными столбцами), HYB (гибридный), <Операция> — axpyi , doti , dotci , gthr , gthrz , roti , sctr ; mv , sv ; mm , sm Все функции возвращают тип специальный тип cusparseStatus_t , показывающий, успешно ли завершилось их выполнение. Выполнение функций происходит асинхронно, так что к моменту возврата из них вовращаемый результат может быть ещё не полностью сформирован, поэтому следует использовать функцию cudaDeviceSynchronize() , чтобы гарантированно иметь результат работы функции из библиотеки cuSPARSE. Для этой же цели можно использовать вызов cudaMemcpy() , являющийся блокирующим: он завершается только тогда, когда результаты готовы. Библиотека поддерживает индексирование с нуля и с единицы. Это можно выбрать с помощью величины типа cusparseIndexBase_t , либо передаваемой в виде параметра, либо содержащейся в поле дескриптора матрицы или вектора. Форматы представления данных Векторы в "плотном" формате (Dense Format) хранятся как обычно: все величины (компоненты вектора) последовательно располагаются в памяти. Разреженные векторы (Sparse Format) представлены в виде двух массивов: ненулевых данных и индексов их позиций в обычном, "плотном" формате (вместе с количеством ненулевых элементов). Обычные ("плотные") матрицы хранятся по столбцам и представлены величинами: числом строк, числом столбцов, лидирующим 6 25 размером, который должен быть больше или равен числу строк (если он больше, то это — подматрица), а также указателем на массив с элементами матрицы. Для разреженных матриц придумано много способов хранения, рассмотрим наиболее часто встречающиеся. Форматы представления разреженных матриц MxN Coordinate Format (COO) Разреженная матрица в этом формате хранится по строкам и описывается числом ненулевых элементов ( nnz ) и тремя массивами этой длины: массивом (ненулевых) элементов ( cooValA ), массивом их индексов в строках ( cooRowIndA ) и столбцах ( cooColIndA ). Простой пример — матрица 2x3: 7.0 0.0 5.0 0.0 6.0 0.0 В предположении индексации с нуля матрица хранится в этом формате в таких массивах: cooValA = [7.0 5.0 6.0] cooRowIndA = [ 0 0 1 ] cooColIndA = [ 0 2 1 ] (вместе с количеством ненулевых элементов — 3). Если используется индексация с единицы, то все индексы будут на единицу больше. Предполагается также, что каждая пара индексов встречается только один раз. Compressed Sparse Row Format (CSR) Этот формат отличается от предыдущего тем, что массив индексов в строках сжат, параметр, называемый csrRowPtrA , — это массив из M+1 значения индекса в массивах csrValA и csrColIndA ; эти значения показывают, где каждом из массивов ( csrValA , csrColIndA ) начинаются данные следующей строки. Последнее равно nnz+csrRowPtrA(0) Для нашей матрицы и нулевой индексации: csrValA = [7.0 5.0 6.0] csrRowPtrA = [ 0 2 3 ] csrColIndA = [ 0 2 1 ] Compressed Sparse Column Format (CSC) Этот формат отличается от формата COO двумя вещами: матрица хранится по столбцам, а массив столбцовых индексов сжат аналогично предыдущему формату. Здесь наша матрица (при индексации с нуля) будет храниться так (вместе с nnz = 3 ): cscValA = [7.0 6.0 5.0] cscRowIndA = [ 0 1 0 ] cscColPtrA = [ 0 1 2 3 ] 26 Ellpack-Itpack Format (ELL) Разреженная матрица MxN с максимум K ненулевыми элементами в строке хранится в этом формате в двух "плотных" матрицах размера MxK. Первая содержит величины ненулевых элементов, вторая — соответствующие индексы столбцов, дополняемые при необходимости значениями 0 («невозможное» значение) и -1 («невозможный» индекс) соответственно. Формат предполагает хранение этих матриц по столбцам. Пример хранения для приведённой выше разреженной матрицы (индексы — с нуля): данные 7.0 5.0 хранение : [7.0 6.0 5.0 0.0] 6.0 0.0 Индексы 0 2 хранение : [ 0 1 2 -1] 1 -1 В библиотеке cuSPARSE этот формат непосредственно не используется, но является частью следующего способа хранения регулярной части разреженных матриц. Hybrid Format (HYB) Гибридный формат хранения разреженной матрицы строится из регулярной части, обычно хранимой в ELL-формате, и нерегулярной, обычно хранимой в COO-формате. Обе части используют индексацию с нуля. Существуют также и другие форматы хранения разреженных матриц: Block Compressed Sparse Row Format (BSR), Extended BSR Format (BSRX) и др. Здесь они не рассматриваются. Пример использования библиотеки cuSPARSE Файл cuSparseTest.cu , предлагаемый в качестве примера, взят из описания библиотеки cuSPARSE на сайте NVIDIA. Стоит обратить внимание на процедуру работы с библиотекой. Сначала с помощью вызова функции cusparseCreate(&handle) библиотека инициализируется, причём в результате возвращается дескриптор handle , который потребуется при вызове операций библиотеки. Затем создаётся и заполняется дескриптор разреженной матрицы ( cusparseCreateMatDescr(&descr) ) и указываются её другие свойства (тип: cusparseSetMatType() , начальный индекс: cusparseSetMatIndexBase() ). Делаются необходимые операции с матрицей (в качестве параметра передаётся дескриптор обращения к библиотеке: handle ), после чего результат копируется на компьютер, освобождается дескриптор матрицы ( cusparseDestroyMatDescr() ) и производится завершение работы с библиотекой ( cusparseDestroy(handle) ). Кроме того, в тесте делается и преобразование разреженной матрицы из одного формата (COO) в другой (CSR) с помощью функции cusparseXcoo2csr() 27 Общие замечания о сторонних библиотеках В принципе библиотеки, расширяющие какой-нибудь программный продукт, могут быть трёх видов: в виде заголовочных файлов (содержат только .h или аналогичные им файлы); статические библиотеки ( .lib ) с заголовочными файлами; статические библиотеки ( .lib ) с заголовочными файлами и динамически подгружаемыми компонентами (чаще всего это — .dll - файлы). Соответственно, и подключение этих трёх видов библиотек к разрабатываемой программе осуществляется проще или сложнее. Проще всего работать с библиотеками, содержащими только заголовочные файлы. Они должны подключаться на этапе компиляции, поэтому необходимо указывать в параметрах проекта дополнительные пути поиска этих файлов (C/C++ | General | Additional Include Directories). В случае статических библиотек — помимо путей поиска файлов — нужно также указать, где лежат статические библиотеки, и какие из них надо просматривать для поиска функций, иначе даже если программа будет компилироваться, то собрать её ("слинковать") — не получится. Для этого в Linker | General надо добавить Additional Library Directories, а в Linker | Input | Addtional Dependencies включить необходимые статические библиотеки. Если же помимо статических библиотек существуют ещё и динамические компоненты, то наша программа может успешно откомпилироваться и собраться, но не запустится, поскольку не найдёт динамических библиотек (обычное сообщение в таких случаях при попытке запуска программы: "Не могу найти динамическую библиотеку ..."). Поскольку обеспечить запуск готовой программы, требующей других динамических библиотек, можно разными способами (в зависимости от того, запускается ли она отдельно или в рамках среды редактирования/компиляции/запуска), действия программиста зависят от конкретной ситуации, но всегда следует помнить, что поиск динамических библиотек начинается с текущего каталога и продолжается в тех каталогах, что указаны в переменной окружения PATH Исходя из всего вышесказанного, весьма предусмотрительно начинать сборку программы со сторонними библиотеками с компиляции отдельных файлов, тогда станет ясно, каких заголовочных файлов не хватает или они не найдены. После успешной компиляции отдельных файлов программу можно попытаться собрать, тогда будет понятно, все ли статические библиотеки нашлись (или указаны), в противном случае линкер будет всё время сообщать о каких-нибудь не найденных функциях. Когда сборка завершится успешно, можно пробовать запускать созданную программу. Если у неё нет зависимости от динамических библиотек – либо все они найдены (т.е., находятся в пути поиска), – она запустится. Если же она требует каких-то динамических библиотек — надо решать, что делать, чтобы они при запуске были обнаружены (скажем, добавить каталоги в пути поиска или переместить динамические библиотеки или их копии в нужное место). |