Как работает BLAS sgemm / dgemm?

Question

Я пытаюсь использовать функцию sgemm в BLAS, используя ctypes в python. Пытаясь решить C = A x B , следующий код работает просто отлично:

no_trans = c_char("n")
m = c_int(number_of_rows_of_A)
n = c_int(number_of_columns_of_B)
k = c_int(number_of_columns_of_A)
one = c_float(1.0)
zero = c_float(0.0)

blaslib.sgemm_(byref(no_trans), byref(no_trans), byref(m), byref(n), byref(k),
               byref(one), A, byref(m), B, byref(k), byref(zero), C, byref(m))

Теперь я хотел бы решить это уравнение: C = A 'x A , где A' - транспонирование A , и следующий код выполняется без исключение, но возвращенный результат неверен:

trans = c_char("t")
no_trans = c_char("n")
m = c_int(number_of_rows_of_A)
n = c_int(number_of_columns_of_A)
one = c_float(1.0)
zero = c_float(0.0)

blaslib.sgemm_(byref(trans), byref(no_trans), byref(n), byref(n), byref(m),
               byref(one), A, byref(m), A, byref(m), byref(zero), C, byref(n))

Для теста я вставил матрицу A = [1 2; 3 4] . Правильный результат: C = [10 14; 14 20] , но процедура sgemm выплевывает C = [5 11; 11 25] .

Насколько я понимаю, матрица A не должна быть транспонирована мной, так как алгоритм позаботится об этом. Что не так с моей передачей параметров во втором случае?

Любая помощь, ссылка, статья, совет приветствуется!

Answer 1

Blas обычно использует основные матрицы столбцов (например, Fortran), поэтому A = [1 2; 3 4] означает

    |1 3|   
A = |   |
    |2 4|

, и результат верный (при условии, что ваша библиотека Python делает то же самое).Посмотреть это read-me

Answer 2

Результат, который у вас есть, указывает, что sgemm вычислил A * A ', а не A' * A, как вы хотели.Простое решение состоит в том, чтобы переключить два входа на функцию.