BLAS sgemm/dgemm 如何工作？

Question

我正在尝试使用 python 中的 ctypes 来使用 BLAS 中的函数 sgemm 。尝试求解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 库也这样做）。请参阅此自述文件

Answer 2

您得到的结果表明 sgemm 已计算出 A*A' 而不是您想要的 A'*A。简单的解决方案是将两个输入切换到该函数。