Неквадратные матрицы C-порядка в cuBLAS (нумба)

Question

Я пытаюсь использовать функции cuBLAS в пакете Anaconda Numba, и у меня возникла проблема. Мне нужно, чтобы входные матрицы были в C-порядке. Вывод может быть в порядке Fortran.

Я могу запустить пример сценария, поставляемого с пакетом, здесь. Скрипт имеет две функции, gemm_v1 а также gemm_v2, В gemm_v1Пользователь должен создать входные матрицы в порядке Фортрана. В gemm_v2, они могут быть переданы в CUDA-реализацию GEMM и перенесены на устройство. Я могу заставить эти примеры работать с квадратными матрицами. Тем не менее, я не могу понять, как получить gemm_v2 работать с неквадратными входными матрицами. Есть ли способ работать с входными матрицами C-порядка, которые не являются квадратными?

Замечания:
В идеале матрицы ввода и вывода должны оставаться на устройстве после вызова GEMM для использования в других вычислениях (это является частью итеративного метода).

Answer 1

Проблема этого примера в том, что он работает только для квадратных матриц. Если матрицы не квадратные, вы не можете рассчитать A^t*B^t из-за несоответствия размеров (при условии, что размеры были правильными для A*B).

У меня под рукой нет работающей установки cuBLAS, так что это своего рода выстрел в темноте, но я был бы очень удивлен, если бы cuBLAS работал иначе, чем обычный BLAS. BLAS ожидает, что матрицы будут в мажорном порядке по столбцам (порядок Fortran), но также могут использоваться для матриц по мажорному порядку строк (C-порядок).

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

Хитрость для вычисления произведения двух C-порядковых матриц с помощью gemm будет работать следующим образом:

Даже если это, вероятно, вам известно, я хотел бы сначала остановиться на главном порядке строк (c-memory-layout) и главном порядке столбцов (fortran-memory-layout), чтобы конкретизировать мой ответ.

Так что, если у нас есть 2x3 (т.е. 2 строки и 3 столбца) матрица Aи сохранить его в некоторой непрерывной памяти, мы получим:

row-major-order(A) = A11, A12, A13, A21, A22, A23
col-major-order(A) = A11, A21, A12, A22, A13, A33

Это означает, что если мы получим непрерывную память, которая представляет матрицу в мажорном порядке строк, и интерпретируем ее как матрицу в мажорном столбце, мы получим совершенно другую матрицу!

Однако, если мы посмотрим на транспонированную матрицу A^t мы можем легко увидеть:

row-major-order(A) = col-major-order(A^t)
col-major-order(A) = row-major-order(A^t)

Это означает, что если мы хотим получить матрицу C в порядке мажорной строки как результат, подпрограмма blas должна записать транспонированную матрицу C в порядке-столбце-мажоре (после всего этого мы не можем изменить) в эту самую память. Тем не мение, C^t=(AB)^t=B^t*A^t а также B^t A^t являются исходными матрицами, переосмысленными в главном порядке столбцов.

Теперь пусть A быть n x k-матрица и B k x m-matrix, вызов функции gemm должен быть следующим:

gemm('N', 'N', m, n, k, 1.0, B, m, A, k, 0.0, C, m)

Пожалуйста, обратите внимание:

1. Нам не нужно транспонировать матрицы A а также B, потому что это обрабатывается путем переосмысления C-порядка как Fortran-порядка.
2. Мы должны поменять местами матрицы A а также B чтобы получить C^t в Фортран-порядке как результат.
3. Полученная матрица C находится в C-порядке (переосмысливая его от Fortran-порядка до C-порядка, мы избавляемся от ^t).