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

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

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

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

Ответ 1

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

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

На мой взгляд, что может быть совершенно неверным, 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 в строке-major-order в качестве результата, blas-процедура должна записать транспонированную матрицу C в столбце major-order (после этого мы не можем изменение) в эту самую память. Однако C^t=(AB)^t=B^t*A^t и B^t an A^t являются исходными матрицами, переинтерпретированными в основном столбце.

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

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

Обратите внимание:

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