LAPACK zgemm op(A) 维度

Question

在来自 netlib 的 this link 中，它指定 M 为：

On entry, M specifies the number of rows of the matrix op( A ) and of the matrix C. M must be at least zero. Unchanged on exit.

因此，如果我想使用 3x10 矩阵作为 A，但我想使用它的共轭 zgemm (TRANSA = 'C')，我应该输入什么作为 M？ 3 还是 10？

此外，当我使用其他 LAPACK 例程时，我将 2D 矩阵输入为 1D，如 A[3*3] 而不是 A[3][3] 并且在调用例程时我只是使用 A 作为矩阵，我可以做吗与非方阵一样吗？ A[3*10] 而不是 A[3][10]?

我用 C++ 编写代码。

Answer 1

A/命名convention/clarification

在给出答案之前，为了更清晰，请务必牢记以下事实：

在美国，M 用于行大小，N 用于列大小

而

在其他一些地方，例如欧洲，这是相反的 N 表示行大小，M 表示列大小

评论：

您在 netlib.org 中找到的所有 Blas/Lapack 文档均使用美国公约
我（作为欧洲人）必须承认美国惯例更符合逻辑，例如索引 (i,j) 和 (m,n) 遵循 相同的字母顺序

为了避免这种歧义，我通常使用：

I_size 行大小 和 J_size 列大小

B/答案

B.1/gemm

void cblas_zgemm(CBLAS_LAYOUT layout,
                 CBLAS_TRANSPOSE opA,
                 CBLAS_TRANSPOSE opB,
                 const int M, <-------------- I_Size of op(A) 
                 const int N, <-------------- J_Size of op(B)
                 const int K, <-------------- J_Size of op(A)
                 const void* alpha,
                 const void* A,
                 const int lda,
                 const void* B,
                 const int ldb,
                 const void* beta,
                 void* C,
                 const int ldc);

在动词中 if TRANSA = 'T' 你必须取转置矩阵的维度。

调用 cblas_zgemm 的实现可能如下所示：

const Size_t opA_I_size = (opA == CblasNoTrans) ? A.I_size() : A.J_size();
const Size_t opA_J_size = (opA == CblasNoTrans) ? A.J_size() : A.I_size();

const Size_t opB_I_size = (opB == CblasNoTrans) ? B.I_size() : B.J_size();
const Size_t opB_J_size = (opB == CblasNoTrans) ? B.J_size() : B.I_size();

cblas_zgemm(CblasColMajor,
            opA,
            opB,
            opA_I_size,
            opB_J_size,
            opA_J_size,
            alpha,
            A.data(),
            A.ld(),
            B.data(),
            B.ld(),
            beta,
            C.data(),
            C.ld());

B.2/内存布局

为了 Blas/Lapack 兼容性和更普遍的数字运算...

从不使用 A[I_size][J_size] 但始终使用 A[I_size*J_size]

（原因是：在一种情况下你有一个指针数组，在另一种情况下你有一个连续的内存块，这对于矢量化、缓存友好性等更方便）

更准确地说为

主列（Fortran 风格）你有：A[ld*J_size]
行专业（C风格）你有：A[I_size*ld]

（其中 ld 是主要维度）

更新：

即使您使用 C++ 编写代码我也建议使用 Fortran 约定（主要列）。 Lapacke pretends to support row major mode too however, under the hood, it simply copies your matrix into a column major layout before calling the requested subroutine. So this extra facility is only an illusion (concerning perfs). To be more precise this is the LAPACKE_dge_trans() function. You can check Lapacke code to see that this function is used nearly everywhere as soon as Layout=RowMajor (see the lapacke_dgesv_work() 例如代码）。
另请注意，如果您想要通用步幅（"arbitrary leading dimension" 在 I 和 J 方向），您可以使用像 Blis 这样的库而不是 Blas。真正的优势是能够创建张量的任意二维视图。这个选择取决于你的应用，我不知道你是否考虑过张量操作。

B.3/矩阵维度

如果你的矩阵总是小到 3x10 blas/lapack 不是一个好的选择（为了性能）。考虑使用像 Eigen or Blaz.

这样的库

LAPACK zgemm op(A) 维度

LAPACK zgemm op(A) dimension

c++

matrix

lapack

lapacke