C 连续矩阵上的 Fortran gemm 函数

Question

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

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

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

我的问题是，我想使用 C-contiguous 数组而不是 Fortran-contiguous 数组，并且我已经使用上述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-contiguous 顺序存储元素，如何解决这个问题？另外我不太明白那些m,n,k,lda,ldb是怎么确定的transa/transb='T' or 'N'，希望你能详细解释一下。

注意

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

Answer 1

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

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

---

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

您的代码以列专业计算 C，但您需要以行专业计算 C。所以你需要通过 Fortran API 计算 col-major 中的 C^T，相当于 row-major 中的 C。

应该是

C^T = B^T * A^T

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

Transpose matrix multiplication in cuBLAS howto