ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUD. Учебное пособие казань 2017 2 удк 004. 4, 681. 3 Ббк 32. 81 Печатается по постановлению Редакционноиздательского совета

• Название: Учебное пособие казань 2017 2 удк 004. 4, 681. 3 Ббк 32. 81 Печатается по постановлению Редакционноиздательского совета
• Анкор: ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUD
• Дата: 22.09.2022
• Размер: 1.28 Mb.
• Имя файла: Tumakov___Tekhnologiya_programmirovaniya_CUDA.pdf
• Тип: Учебное пособие; #690313
• страница: 5 из 6

Глава 9. Пакеты для работы с векторами и матрицами
В этой главе рассмотрим два пакета: cuBLAS и cuBLAS-Xt.
Пакет cuBLAS
cuBLAS – реализация интерфейса программирования приложений для создания библиотек, выполняющих основные операции линейной алгебры
BLAS (Basic Linear Algebra Subprograms) для CUDA. Этот пакет позволяет получить доступ к вычислительным ресурсам графических процессоров
NVIDIA. Перед тем, как начать изучение непосредственно библиотеки cuBLAS, стоит ознакомится с интерфейсом BLAS, который реализует данная библиотека.
BLAS (Basic Linear Algebra Subprograms) – стандарт интерфейса программирования приложений (API) для создания библиотек, выполняющих основные операции линейной алгебры, такие как умножение векторов и матриц. Впервые опубликован в 1979 году.
Функциональность BLAS делится на три уровня.
Уровень 1:
Первый уровень содержит векторные операции вида:
y←αx + y,
операции скалярного произведения, взятия нормы вектора и другие операции.
Уровень 2:
Этот уровень содержит операции матрица-вектор вида:
y←αАx + βy,
решение Tx = y для x с треугольной матрицей T и другие операции.
Уровень 3:
Содержит операции матрица-матрица вида:

84
C←αАB + βC,
решение B←αT
-1
B для треугольной матрицы T и другие операции.
Введение в cuBLAS
CuBLAS является самодостаточной библиотекой на уровне API, то есть, прямого взаимодействия с драйвером CUDA не происходит. cuBLAS прикрепляется к одному GPU и автоматически не распараллеливается между несколькими GPU.
Каждая реализуемая функция представлена в нескольких вариантах: для данных одинарной и двойной точности, вещественных и комплексных.
Поскольку пакет BLAS изначально был написан на языке Fortran,
NVIDIA приложила все усилия к тому, чтобы обеспечить максимальную совместимость. Соответственно, массивы в cuBLAS хранятся по столбцам, а не по строкам, как принято в C и C++. При этом сохраняется индексирование с 1.
Чтобы упростить написание программ в рамках тех соглашений, которые приняты в языках C и C++, достаточно ввести макрос для вычисления элемента двумерного массива, хранимого в одномерном виде и в FORTRAN-овской нотации
#define IDX2F(i,j,ld) ((((j)-1)*(ld))+((i)-1))
Здесь ld – это размерности матрицы, и в случае хранения в столбцах, является количеством строк. Для кода изначально написанного на С и С++, можно было бы использовать индексирование с 0, в этом случае макрос выглядит так:
#define IDX2C(i,j,ld) (((j)*(ld))+(i))
В старых версиях библиотеки основные функции cuBLAS (в отличие от вспомогательных функций) не возвращают статус ошибки напрямую (из соображений совместимости с существующими библиотеками BLAS). В этом случае, чтобы помочь в отладке, cuBLAS предоставляет отдельную функцию, которая возвращает последнюю записанную ошибку. Мы будем рассматривать

85 только вторую версию библиотеки. Об особенностях старой (legacy) версии можно прочитать в документации (раздел 1.2. New and Legacy cuBLAS API)
Использование cuBLAS API
Перед использованием функций библиотеки cuBLAS, нужно инициализировать контекст (т.е. содержимое самой библиотеки) и передавать его при каждом вызове функций библиотеки. Контекст хранится в переменной типа cublasHandle_t. Инициализация контекста производиться вызовом функции cublasCreate(), который возвращает значение структуры типа cublasHandle_t. После завершения работы с библиотекой нужно вызвать функцию cublasDestory(), передав переменную, которую вернула функция cublasCreate(), чтобы освободить ресурсы, связанные с контекстом библиотеки cuBLAS.
Такой подход к использованию контекстов позволяет явно контролировать настройки при использовании нескольких нитей в хосте и нескольких графических ускорителей. При этом предполагается, что контекст не будет меняться между вызовами cublasCreate() и cublasDestory(). Для того, чтобы cuBLAS мог использовать разные GPU в одной и той же нити хоста, нужно инициализировать другой контекст, предварительно установив соответствующий девайс вызовом cudaSetDevice().
В дальнейшем, при приведении примеров без полного листинга, будем иметь ввиду, что все необходимые библиотеки подключены, макросы для доступа к элементам матрицы определены и приводимый код находится между вызовами функций создания и уничтожения контекста. Сам контекст хранится в переменной handle. Так же будем считать, что используемые переменные были определены и инициализированы случайными значениями.
Все функции библиотеки cuBLAS возвращают статус ошибки cublasStatus_t. При удачном завершении операции функция возвращает значение CUBLAS_STATUS_SUCCESS. Информацию о других статусах можно получить в официальной документации cuBLAS.

86
Вспомогательные функции cuBLAS
Эти функции не входят в интерфейс BLAS, они реализованы для взаимодействия с графическим ускорителем. Мы здесь приведем лишь краткое описание некоторых функций и их входных параметров. Описание возвращаемых статусов, причины тех или иных ошибок и способы их исправления смотрите в официальной документации. cublasCreate() и cublasDestroy(); cublasStatus_t cublasCreate(cublasHandle_t *handle); cublasStatus_t cublasDestroy(cublasHandle_t handle);
Как было описано ранее, эти функции отвечают за инициализацию контекста. Перед тем, как начать работать с функциями cuBLAS, нужно вызвать cublasCreate(). Так как cublasCreate() выделяет некоторые внутренние ресурсы, нужно в конце работы их освободить при помощи функции cublasDestroy().
Рекомендуется минимизировать использование данных функций в программах (желательно вызывать только один раз для каждого устройства). cublasSetVector(); cublasStatus_t cublasSetVector(int n, int elemSize, const void *x, int incx, void *y, int incy);
Функция копирует n элементов вектора x из памяти хоста в память девайса GPU. Предполагается, что элементы обоих векторов имеют размер elemSize байтов. Интервал между элементами у исходного вектора задается параметром incx, а у вектора получателя – incy.

87
В общем случае, переменная y должна указывать на объект или часть объекта, который был выделен на девайсе до этого. cublasGetVector(); cublasStatus_t cublasGetVector(int n, int elemSize, const void *x, int incx, void *y, int incy);
Функция копирует n элементов вектора x из памяти девайса в память хоста. Предполагается, что элементы обоих векторов имеют размер elemSize байтов. Интервал между элементами у исходного вектора задается параметром incx, а у вектора получателя – incy. cublasSetMatrix(); и cublasGetMatrix(); cublasStatus_t cublasSetMatrix(int rows, int cols, int elemSize, const void *A, int lda, void *B, int ldb); cublasStatus_t cublasGetMatrix(int rows, int cols, int elemSize, const void *A, int lda, void *B, int ldb);
Аналогично функциям cublasSetVector() и cublasGetVector(), cublasSetMatrix() и cublasGetMatrix() производят копирование матриц из хоста на девайс и обратно.
Функции BLAS
Для того, чтобы понять особенности использования библиотеки cuBLAS, опишем некоторые функции базовой подпрограммы линейной алгебры.
Все функции для математических операций имеют следующий вид:

88 cublasStatus_t cublas(cublasHandle_t handle, ...);
Буква T обозначает тип данных
•
S = float;
•
D = double;
•
C = complex float;
•
Z = complex double. и operation – принятое в стандарте BLAS название операции.
Например, функция вычисления нормы вектора с элементами типа float имеет вид: cublasStatus_t cublasSnrm2(cublasHandle_t handle, int n, const float *x, int incx, float *result)
Реализация той же операции для чисел с двойной точностью будет следующей: cublasStatus_t cublasDnrm2(cublasHandle_t handle, int n, const double *x, int incx, double *result)
Описание параметров для функций cublasnrm2().
Параметр Размещение
в памяти
In/out Описание
handle input Контекст библиотеки cuBLAS n input Количество элементов вектора x x device input Вектор типа с n элементами.

89
Параметр Размещение
в памяти
In/out Описание
incx input Расстояние между элементами вектора x. result host или device output Результат вычисления нормы, который равен 0.0 если n,incx<=0.
Также при вычислениях с матрицами, часто к матрицам добавляются некоторые преобразования. То есть вычисления имеют, например, вид y = α
op(A)x + βy, где β и α – скаляры, y и x – векторы, A – матрица, op() - некоторая операция, применяемая к матрице (без преобразования, транспонирование или эрмитово сопряжение). Вид этой операции задается как параметр функции типа cublasOperation_t и может принимать следующие значения:
•
CUBLAS_OP_N – матрица не преобразовывается;
•
CUBLAS_OP_T – к матрице применяется транспонирование;
•
CUBLAS_OP_C – к матрице применяется эрмитово сопряжение.
Также, некоторые функции выполняют лево- или правостороннее умножение матриц, в зависимости от входных параметров. Параметр типа cublasSideMode_t определяет тип умножения и принимает значения:
•
CUBLAS_SIDE_LEFT – матрица с левой стороны;
•
CUBLAS_SIDE_RIGHT – матрица, соответственно, с правой стороны.
cuBLAS функции первого уровня
cublasIamax() функция поиска индекса (наименьшего) максимального элемента; cublasIamin() функция поиска индекса (наименьшего) минимального элемента;

90 cublasasum() функция находит сумму абсолютных значений элементов вектора; cublasaxpy() функция умножает вектор x на скаляр alpha и прибавляет результат к вектору y. Значение вектора y презаписывается; cublascopy() функция поэлементно копирует вектор x в вектор y;
cublasdot() функция производит скалярное умножение векторов x и y;
cublasnrm2() функция вычисляет Евклидову норму вектора; cublasrot() функция применяет поворот Гивенса (поворот вектора на некоторый заданный угол) к векторам x и y;
cublasrotg() функция применяет поворот Гивенса к векторам x и y (с обнулением второго вхождения); cublasrotm() функция применяет модифицированную трансформацию
Гивенса к векторам x и y;
cublasrotmg() функция применяет модифицированную трансформацию
Гивенса к векторам x и y (с обнулением второго вхождения); cublasscal() функция умножает вектор на скаляр и перезаписывает вектор новыми значениями; cublasswap() функция меняет вектора местами.

91
Пример 1. Реализовать следующее вычисление: x = α xx, где x – вектор, α
– скаляр. Вычисления проводятся для чисел с двойной точностью.
// векторы x и y определены на CPU и скопированы
// на GPU (dev_x, dev_y) cublasDcopy(handle, n, dev_x, 1, dev_y, 1); cublasDscal(handle, n, &alpha, dev_y, 1); cublasDdot(handle, n, dev_y, 1, dev_x, 1, &result);
cuBLAS функции второго уровня
cublasgbmv() функция применяет окаймленное умножение матрицы на вектор вида
y = α op(A)x + βy,
где A – ленточная матрица; cublasgemv() функция применяет окаймленное умножение матрицы на вектор вида
y = α op(A)x + βy;
cublasger() функция выполняет операцию вида
A = αxy
T
+ A, где A – матрица размерности n на m;
cublassbmv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – симметричная ленточная матрица; cublasspmv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – симметричная матрица, хранимая в пакетном формате;

92 cublasspr() функция выполняет операцию вида
A = αxx
T
+ A;
cublasspr2() функция выполняет операцию вида
A = α(xy
T
+ yx
T
) + A, где A – симметричная матрица, хранимая в пакетном формате; cublassymv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – симметричная матрица; cublassyr() функция выполняет операцию вида
A = αxx
T
+ A, где A – симметричная матрица; cublassyr2() функция выполняет операцию вида
A = α(xy
T
+ yx
T
) + A, где A – симметричная матрица; cublastbmv() функция выполняет операцию вида
x = op(A)x, где A – треугольная матрица; cublastbsv() функция решает уравнение вида
op(A)x = b, где A – треугольная матрица; cublastpmv() функция выполняет операцию вида
x = op(A)x, где A – треугольная матрица, хранимая в пакетном формате;

93 cublastpsv() функция решает уравнение вида
op(A)x = b, где A – треугольная матрица, хранимая в пакетном формате; cublastrmv() функция выполняет окаймленное умножение матрицы на вектор вида
x = op(A)x ,
где A – треугольная матрица; cublastrsv() функция решает уравнение вида
op(A)x = b, где A – треугольная матрица; cublashemv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – эрмитова матрица; cublashbmv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – эрмитова ленточная матрица; cublashpmv() функция применяет умножение матрицы на вектор вида
y = αAx + βy,
где A – эрмитова матрица, хранимая в пакетном формате; cublasher() функция выполняет операцию вида
A = αxy
H
+ A, где A – эрмитова матрица;

94 cublasher2() функция выполняет операцию вида
A = α(xy
T
+ yx
T
) + A, где A – эрмитова матрица; cublashpr() функция выполняет операцию вида
A = αxy
H
+ A, где A – эрмитова матрица, хранимая в пакетном формате; cublashpr2() функция выполняет операцию вида
A = α(xy
T
+ yx
T
) + A, где A – эрмитова матрица, хранимая в пакетном формате.
Пример 2. Вычислить следующее выражение: y = αAx.
// вектор x и матрица A определены на CPU и скопированы
// на GPU (dev_x, dev_A)
// вектор y инициализирован 0-ми и так же скопирован на
// GPU (dev_y) cublasDgemv()(handle, n, m, &alpha, dev_A, m, dev_x, 1,
&beta, dev_y, 1);
cuBLAS функции третьего уровня
cublasgemm() функция умножения матриц вида
C = α op(A)op(B) + β op(C), где A – m на k, B – k на n, C – m на n;
cublasgemm3m() функция умножения матриц вида
C = α op(A)op(B) + β op(C), где A – m на k, B – k на n, C – m на n. При этом применяется метод Гаусса для ускорения операции;

95 cublasgemmBatched() функция умножения матриц вида
C[i] = α op(A[i])op(B[i]) + β op(C[i]), для массива матриц; cublassymm() функция умножения симметричных матриц вида
C = αAB + βC или C = αBA + βC;
cublassyrk() функция выполняет для симметричных матриц операцию вида
C = α op(A)op(A)
T
+ βC;
cublassyr2k() функция выполняет для симметричных матриц операцию вида
C = α (op(A)op(B)
T
+ op(B)op(A)
T
) + βC;
cublassyrkx() функция выполняет для симметричных матриц операцию вида
C = α op(A)op(B)
T
+ βC;
cublastrmm() функция умножения треугольных матриц вида
C = αop(A)B или C = αBop(A);
cublastrsm() функция решает для треугольных матриц уравнение вида
op(A)X = αB или X op(A)= αB;
cublastrsmBatched() функция решает для массива треугольных матриц уравнение вида
op(A[i])X[i] = αB[i] или X[i] op(A[i])= αB[i];

96 cublashemm() функция умножения для эрмитовых матриц вида
C = αAB + βC или C = αBA + βC;
cublasherk() функция реализует для эрмитовых матриц операцию вида
C = α op(A)op(A)
H
+ βC;
cublasher2k() функция реализует для эрмитовых матриц операцию вида
C = α op(A)op(В)
H
+ α op(В)op(A)
H
+ βC;
cublasherkx() функция выполняет для эрмитовых матриц операцию вида
C = α op(A)op(В)
H
+ βC.
Пример 3. Умножить две матрицы. Результат записать в третью матрицу.
// матрицы A и B определены на CPU и скопированы на GPU
// (dev_A, dev_B)
// матрица C инициализирован 0-ми и скопирован на GPU
// (dev_С)
// alpha и beta равны 1 cublasDdot(handle, CUBLAS_OP_N, CUBLAS_OP_N, n, n, n,
&alpha, dev_A, n, dev_B, n, &beta, dev_A, n);
BLAS-подобное расширение cublasgeam() функция выполняет над матрицами операцию вида
C = α op(A) + β op(B), где A, B, C – m на n;
cublasdgmm() функция выполняет операцию вида
C = A diag(X) или C = diag(X) A ,

97 где A – матрица m на n, X – вектор из n элементов; cublasgetrfBatched() функция выполняет LU-разложение для массива матриц:
P A[i] = LU;
cublasgetrsBatched()функция решает систему линейных уравнений для массива LU-разложенных матриц; cublasgetriBatched() функция вычисляет обратные матрицы для массива матриц; cublasmatinvBatched() функция вычисляет обратные матрицы для массива матриц; cublasgeqrfBatched() функция выполняет для массива матриц QR- разложение; cublasgelsBatched() функция решает для массива матриц следующую задачу:
min||Carr[i] - Aarr[i]*Xarr[i] ||;
cublastpttr() функция выполняет преобразование треугольной матрицы из пакетного формата в «обычный»; cublastrttp() функция выполняет преобразование треугольной матрицы из «обычного» формата в пакетный; cublasgemmEx() функция умножения матриц вида
C = α op(A)op(B) + β op(C),

98 где A – m на k, B – k на n, C – m на n. Входные данные могут иметь более низкую точность, чем выходные; cublasCsyrkEx() функция выполняет для симметричных матриц операцию вида
C = α op(A)op(A)
T
+ βC,
Причем, входные данные могут иметь более низкую точность, но выходные будут типа cuComplex. cublasCsyrk3mEx() функция выполняет модификацию предыдущей операции; cublasCherkEx() функция выполняет для симметричных матриц операцию вида
C = α op(A)op(A)
H
+ βC, входные данные могут иметь более низкую точность , но выходные будут типа cuComplex. cublasCsyrk3mEx() функция выполняет модификацию предыдущей операции; cublasCherkEx() функция выполняет ту же операцию, что и cublasCherk.
Входные данные могут иметь более низкую точность, но выходные будут типа cuComplex. cublasCherk3mEx() функция выполняет модификацию предыдущей операции;

99 cublasNrm2Ex() функция вычисляет Евклидову норму вектора. Входные данные могут иметь более низкую точность, чем выходные (тип входных и выходных параметров передаются дополнительно). cublasAxpyEx() функция является обобщением cublasaxpy, где тип входных и выходных параметров передаются дополнительно. cublasDotEx() функция является обобщением cublasdot, где тип входных и выходных параметров передаются дополнительно. cublasScalEx() функция является обобщением cublasscal, где тип входных и выходных параметров передаются дополнительно.