ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUD. Учебное пособие казань 2017 2 удк 004. 4, 681. 3 Ббк 32. 81 Печатается по постановлению Редакционноиздательского совета
 Скачать 1.28 Mb.
|
|
Пример 4. Выполнить скалярное произведение двух векторов одинарной точности. Результат должен быть типа double. cublasDotEx(handle, n, dev_x, CUDA_R_32F, 1, dev_y, CUDA_R_32F, 1, &result, CUDA_R_64F, CUDA_R_64F); Приведем пример решения системы линейных уравнений с нижней треугольной матрицей. #include #include #include #include // подключение библиотеки cuBLAS #include // макрос для работы с индексами в стиле FORTRAN 100 #define IDX2C(i,j,ld) (((j)*(ld))+(i)) int main() { const int N = 6; cublasHandle_t handle; float *dev_A, *dev_b; float *x, *A, *b; x = (float *)malloc(N * sizeof(*x)); b = (float *)malloc(N * sizeof(*b)); A = (float *)malloc(N * N * sizeof(*A)); // инициализация матрицы и вектора правой части int ind = 11; for (int j = 0; j < N; j++) { for (int i = 0; i < N; i++) if (i >= j) A[IDX2C(i, j, N)] = (float)ind++; else A[IDX2C(i, j, N)] = 0.0f; b[j] = 1.0f; } // выделяем память на GPU соответствующего размера // для каждой переменной cudaMalloc((void**)&dev_b, N * sizeof(*x)); cudaMalloc((void**)&dev_A, N * N * sizeof(*A)); // инициализируем контекст cuBLAS cublasCreate(&handle); // копируем вектор и матрицу из CPU в GPU cublasSetVector(N, sizeof(*b), b, 1, dev_b, 1); cublasSetMatrix(N, N, sizeof(*A), A, N, dev_A, N); // решаем нижнюю треугольню матрицу 101 cublasStrsv(handle, CUBLAS_FILL_MODE_LOWER, CUBLAS_OP_N, CUBLAS_DIAG_NON_UNIT, N, dev_A, N, dev_b, 1); // копируем результат из GPU в CPU cublasGetVector(N, sizeof(*x), dev_b, 1, x, 1); for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) printf("%3.0f ", A[IDX2C(i, j, N)]); printf(" = %f %4.6f\n", b[i], x[i]); } // освобождаем память в GPU cudaFree(dev_b); cudaFree(dev_A); // уничтожаем контекс cuBLAS cublasDestroy(handle); // освобождаем память в CPU free(x); free(b); free(A); getchar(); return EXIT_SUCCESS; } При компиляции могут возникать следующие ошибки: LNK2019 unresolved external symbol cublasCreate_v2 referenced in function main Для их устранения необходимо подключить библиотеку cublas.lib, добавив ее в Project→Properties→Linker→Input→Additional Dependencies. После этого возможна следующая ошибка для 64-битных систем: LNK1104 cannot open file 'cublas.lib'. Эта проблема устраняется изменением значения Build→Configuration Manager…→Platform на x64. 102 Пакет cuBLAS-Xt cuBLAS-Xt – это набор рутинных подпрограмм, которые ускоряют вызовы в Level 3 BLAS (базовые подпрограммы линейной алгебры) за счет распределения нагрузки на нескольких GPU. Благодаря потоковой структуре, cuBLAS-Xt эффективно автоматически управляет пересылками на шине PCI- Express, что позволяет хранить входные и выходные данные в системной памяти. Это приводит к тому, что операции выполняются с использованием внешней памяти, размер операндов памяти ограничен только размером системной памяти, а не памятью GPU. cuBLAS-XT поддерживается на всех 64- битных платформах. cuBLASXt API заботится о распределении памяти между назначенными графическими процессорами и распределяет нагрузку между ними и, наконец, возвращает результаты обратно на хост. cuBLAS-Xt API поддерживает только вычисление операций вида BLAS3 (матрица-матрица). Начиная с версии CUDA 6.0, бесплатная версия cuBLAS-Xt входит в состав CUDA Toolkit как часть библиотеки cuBLAS. Cublas-Xt API имеет файл заголовка cuBLAS-Xt.h. Примеры cuBLAS-Xt находятся в каталоге: C:\ProgramData\NVIDIA Corporation\CUDA Samples\v8.0\7_CUDALibraries\simpleCUBLASXT Инициализация контекста cuBLAS-XT API Для возможности использования функций библиотеки cuBLAS-Xt требуется проинициализировать контекст. Тип cublasXtHandle_t является типом указателя на cuBLAS-Xt API контекст. cublasXt API контекст инициализируется с помощью cublasXtCreate(), и возвращенный дескриптор должен быть передан на все последующие функции cuBLAS-Xt API. Контекст должен быть уничтожен в конце, используя cublasXtDestroy(). cublasStatus_t status; cublasXtHandle_t handle; status = cublasXtCreate(&handle); 103 ----- status = cublasXtDestroy(handle); После инициализации контекста выбираются графические устройства, между которыми будут распределяться вычисления. int devices[1] = { 0 }; cublasXtDeviceSelect(handle, 1, devices); Функция cublasXtDeviceSelect позволяет пользователю указать количество устройств GPU и соответствующие им идентификаторы, которые будут участвовать в последующих вызовах функций API cuBLAS-Xt Math. cublasXtDeviceSelect(cublasXtHandle_t handle, int nbDevices, int DeviceId[]) Эта функция создает контекст cuBLAS для каждого GPU, представленного в этом списке. Конфигурация устройства не может быть изменена между вызовами функций Math. В этой связи, эта функция должна вызываться только один раз после функции cublasXtCreate. Для того, чтобы иметь возможность запускать несколько конфигураций, требуется несколько контекстов. Распределение нагрузки Для распределения нагрузки между несколькими GPU, API cuBLAS-Xt использует стратегию черепицы (“tiling strategy”) каждая матрица делится на квадратные плитки размерностью, выбранной пользователем. Функция cublasXtSetBlockDim позволяет пользователю установить размер блока, используемого для разбиения матриц для последующих вызовов функций Math. cublasXtSetBlockDim(cublasXtHandle_t handle, int blockDim) 104 Матрицы разделяются в квадратные плитки blockDimxblockDim. Эта функция может быть вызвана в любое время и вступит в силу для следующих вызовов функций Math. Размер блока должен быть выбран таким образом, чтобы оптимизировать работу и убедиться в том, что переводы PCI хорошо покрываются вычислениями. Функция cublasXtGetBlockDim() позволяет пользователю запрашивать размер блока. Функции cuBLAS-Xt API Math Все функции для математических операций имеют вид: cublasStatus_t cublasXt <Т> float 'S' double 'D' cuComplex 'C' cuDoubleComplex 'Z' 1. cublasXt cublasStatus_t cublasXtSgemm(cublasXtHandle_t handle, cublasOperation_t transa, cublasOperation_t transb, size_t m, size_t n, size_t k, const float *alpha, const float *A, int lda, const float *B, int ldb, const float *beta, float *C, int ldc) Эта функция выполняет перемножение матриц. C=αop(A)op(B)+βC, 105 == == == = C, CUBLAS_OP_ transa if , A T, CUBLAS_OP_ transa if , A N, CUBLAS_OP_ transa if A, op(A) H T где α и β – скаляры, A, B и C являются матрицами с размерами m × k , k × n и m × n, соответственно, lda, ldb, ldc – первые размерности двумерных массивов, используются для хранения матриц. 2. cublasXt cublasStatus_t cublasXtChemm(cublasXtHandle_t handle, cublasSideMode_t side, cublasFillMode_t uplo, size_t m, size_t n, const cuComplex *alpha, const cuComplex *A, size_t lda, const cuComplex *B, size_t ldb, const cuComplex *beta, cuComplex *C, size_t ldc) Эта функция выполняет умножение Эрмитовой матрицы на матрицу. == + == + = E_RIGHT, CUBLAS_SID side if C, BA E_LEFT, CUBLAS_SID side if C, AB C β α β α где α и β скаляры, A – Эрмитова матрица, m × n размер матриц B и С соответствующие матрице A, uplo указывает, с какой частью матрицы производить вычисления – нижней треугольной или верхней (т.к нижнюю треугольную матрицу эрмитовой матрицы можно получить из верхней и наоборот). 3. cublasXt cublasStatus_t cublasXtSsymm(cublasXtHandle_t handle, cublasSideMode_t side, cublasFillMode_t uplo, size_t m, size_t n, const float *alpha, const float *A, size_t lda, const float *B, size_t ldb, const float *beta, float *C, size_t ldc) Эта функция выполняет умножение симметричной матрицы на матрицу == + == + = E_RIGHT, CUBLAS_SID side if C, BA E_LEFT, CUBLAS_SID side if C, AB C β α β α 106 где α и β – скаляры, A, B и C являются матрицами с размерами m × k , k × n и m × n, соответственно, lda, ldb, ldc – первые размерности двумерных массивов, используются для хранения матрицы, uplo указывает, с какой частью матрицы производить вычисления – нижней треугольной, либо верхней. 4. cublasXt cublasStatus_t cublasXtStrsm(cublasXtHandle_t handle, cublasSideMode_t side, cublasFillMode_t uplo, cublasOperation_t trans, cublasXtDiagType_t diag, size_t m, size_t n, const float *alpha, const float *A, size_t lda, float *B, size_t ldb) Эта функция решает треугольную линейную систему с несколькими правыми сторонами. == = == = E_RIGHT, CUBLAS_SID side if B, op(A) E_LEFT, CUBLAS_SID side if B, op(A)X α α X где А – треугольная матрица с или без главной диагонали, uplo указывает, с какой частью матрицы производить вычисления – нижней треугольной или верхней. X и В – матрицы размера m×n, и α скаляр. Кроме того, для матрицы А: == == == = C. CUBLAS_OP_ transa if , A T, CUBLAS_OP_ transa if , A N, CUBLAS_OP_ transa if A, op(A) H T Матрица решений X записывается поверх матрицы B на выходе. 5. cublasXt Функция выполняет умножение треугольной матрицы на матрицу. == == = E_RIGHT. CUBLAS_SID side if op(A), B E_LEFT, CUBLAS_SID side if A)B, ( op C α α 107 6. cublasXt Функция выполняет умножение симметрично упакованной матрицы на матрицу. == + == + = E_RIGHT. CUBLAS_SID side if C, BA E_LEFT, CUBLAS_SID side if C, AB C β α β α Гибридные вычисления CPU-GPU В случае очень больших проблем API cuBLAS-Xt предлагает возможность переложить часть вычислений на CPU. Эта функция может быть настроена с процедуры cublasXtSetCpuRoutine() и cublasXtSetCpuRatio(). Эту функцию следует использовать с осторожностью, поскольку он может помешать потокам процессора, работающим с видеокартой. На данный момент только cublasXt Документация библиотеки расположена на сайте http://docs.nvidia.com/cuda/cublas/#using-the-cublasXt-api Приведем листинг кода перемножения матриц: #include #include #include #include /* Подключаем требуемые библиотеки cuda */ #include #include #define N (2) int main(int argc, char **argv) { // статус ошибки возвращается всеми cuBLAS функциями cublasStatus_t status; 108 float *h_A; float *h_B; float *h_C; // Скаляры float alpha = 1.0f; float beta = 0.0f; // Размер матрицы int n2 = N * N; // Указатель на контекст CUBLAS-XT cublasXtHandle_t handle; int devices[1] = { 0 }; // Инициализируем контекст CUBLAS-XT status = cublasXtCreate(&handle); if (status != CUBLAS_STATUS_SUCCESS) { fprintf(stderr, "!!!! CUBLASXT initialization error\n"); return EXIT_FAILURE; } /* Выбираем устройства для запуска функций CUBLAS-XT между которыми будет распределяться нагрузка */ status = cublasXtDeviceSelect(handle, 1, devices); if (status != CUBLAS_STATUS_SUCCESS) { fprintf(stderr, "!!!! CUBLASXT device selection error\n"); return EXIT_FAILURE; } /* Устанавливаем размер блоков (blockDim x blockDim) на которые будут разбиваться матрицы при распределение между устройствами */ status = cublasXtSetBlockDim(handle, 64); if (status != CUBLAS_STATUS_SUCCESS) { fprintf(stderr, "!!!! CUBLASXT set block dimension error\n"); return EXIT_FAILURE; } 109 // Выделяем память для матриц в системной памяти h_A = (float *)malloc(n2 * sizeof(h_A[0])); h_B = (float *)malloc(n2 * sizeof(h_B[0])); h_C = (float *)malloc(n2 * sizeof(h_C[0])); // Заполняем матрицы тестовыми данными srand((unsigned int)time(NULL)); for (int i = 0; i < n2; i++) { h_A[i] = rand() / (float)RAND_MAX; h_B[i] = rand() / (float)RAND_MAX; } // Выполняем операцию перемножения матриц status = cublasXtSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha, h_A, N, h_B, N, &beta, h_C, N); if (status != CUBLAS_STATUS_SUCCESS) { fprintf(stderr, "!!!! kernel execution error.\n"); return EXIT_FAILURE; } // Вывод результата перемножения матриц for (int i = 0; i < N; i++) { for (int j = 0; j < N; j++) printf("%4.4f ", h_B[i*N + j]); printf(" * "); for (int j = 0; j < N; j++) printf("%4.4f ", h_A[i*N + j]); printf(" = "); for (int j = 0; j < N; j++) printf("%4.4f ", h_C[i*N + j]); printf("\n"); } // Очищаем память free(h_A); free(h_B); free(h_C); 110 getchar(); } Контрольные вопросы 1. Как происходит инициализация контекста в библиотеке cuBLAS? 2. Как можно оптимизировать первый пример? Можно ли в этом примере обойтись одним вектором вместо двух? 3. Где должна размещаться (девайс или хост) переменная, в которую возвращается результат вычисления скалярного произведения (cublasDdot())? 4. Перечислите основные этапы работы с cuBLAS-Xt API. 5. Назначение типа cublasXtHandle_t. Возможно ли обойтись без него? 6. Какой вид имеют функции для математических операций? 7. В чем разница библиотек cublas и cuBLAS-Xt? 8. Назначение функции cublasXtSetBlockDim. Лабораторная работа 1. Дополните пример обработкой ошибок. Учитывайте как ошибки при выделении памяти на CPU и GPU (функции CUDA), так и возвращаемый статус функций cuBLAS. 2. Вычислите матрицу С = op(A + B) * (A * B), где A и B квадратные матрицы, используя библиотеку cuBlasXT. Домашнее задание 1. Расширьте пример для решения любых матриц. Для этого использовать LU-разложение из BLAS-like функций. 2. Найдите значения матрицы X. α (A + C) * C * X = B, где α – скаляр, матрицы A и С верхнетреугольные. A=((1,1), (0,1)), C=((1, -1), (0, 1)), B=((12, 10) (6, 8)). 3. Сравните время работы перемножения матриц с использованием библиотек cuBLAS, cuBLAS-Xt и без библиотек. 111 Литература 1. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA. – М.: ДМК Пресс, 2010. – 232 с. 2. Сандерс Дж., Кэндрот Э. Технология CUDA в примерах: Введение в программирование графических процессоров. – М.: ДМК Пресс, 2013. – 232 с. 3. Cheng J., Grossman M., McKercher T. Professional CUDA C Programming. – Indianapolis: John Wiley & Sons, Inc., 2014. – 528 p. Тумаков Дмитрий Николаевич Чикрин Дмитрий Евгеньевич Егорчев Антон Александрович Голоусов Святослав Владимирович ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUDA Учебное пособие Подписано в печать 14.11.2017 Форм. Бум. 60 X 84 1/16. Гарнитура “Таймс”. Печать цифровая. Печ. л. 8,25. Т. 100. Заказ 15. Издательство КФУ, 420008, Казань, ул. пр Нужина, 1/37, +7 (843) 233-73-28, +7 (843) 233-73-59 |