C 连续矩阵上的 Fortran gemm 函数

Question

我正在尝试使用 fortran BLAS gemm 函数进行矩阵乘法，请参阅 here。

这个函数的签名是，所有参数的含义都可以在上面找到link。

call sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)

我的问题是，我想使用 C 连续数组而不是 Fortran 连续数组，我已经使用上面的 sgemm 很长一段时间了，仍然很困惑。

请帮我看一些具体的例子。

我所有的输入数组都是 C 连续的。

a = [[0,1],
     [2,3]]
b = [[0,1,2],
     [3,4,5]]
# pre-alloc memory for c
c = [[0,0,0],
     [0,0,0]]

# compute c = a * b, which should be as follows
# c = [[3,4,5],
#      [9,14,19]]

# since sgemm assumes Fortran-contiguous, so I thought it would be
sgemm('T', 'T', 2, 3, 2, 1.0, a, 2, b, 3, 0, c, 2)
      ~~~~~~~   ~~~~~~~         ~~~   ~~~      ~~~
    trans both   m,n,k          lda   ldb      ldc

# HOWEVER, c is not what I expected, 
c = [[3,9,4],
     [14,5,19]]

显然sgemm 以Fortran 连续顺序存储元素，如何解决这个问题？另外我不太明白那些m,n,k,lda,ldb是如何确定的transa/transb='T' or 'N'，希望你能给我一个详细的解释。

注意

我正在使用从 scipy.linalg.cython_blas 导出的这个 gemm 函数，这意味着，除了玩这个 Fortran 排序的东西，我别无选择。

Answer 1

如果你想使用行优先矩阵而不是 Fortran 风格的列优先矩阵，你可以使用 CBLAS API gemm。您可以使用第一个参数选择矩阵存储布局。

https://software.intel.com/en-us/node/520775

或者您仍然可以使用 Fortran API。因为改变矩阵布局等同于矩阵转置。但是，您计算转置 C 的方式有误。

您的代码计算列优先的 C，但您需要行优先的 C。所以需要用FortranAPI计算col-major中的C^T，相当于row-major中的C。

应该是

C^T = B^T * A^T

基本上你需要交换A和B，以及相应的参数。有关这些参数的更多详细信息，您可以查看此答案。

Transpose matrix multiplication in cuBLAS howto

C 连续矩阵上的 Fortran gemm 函数

Fortran gemm function on C-contiguous matrices

c

python

fortran

blas

matrix-multiplication