c++ - LAPACK zgemm op(A) 维度

Question

在this来自 netlib 的链接将 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/命名约定/说明

在给出答案之前以及为了更清楚地了解这一点很重要:

• 在美国，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' 您必须采用转置 A 矩阵的维度。

调用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假装也支持行主要模式，但是，在幕后，它只是在调用请求的子例程之前将您的矩阵复制到列主要布局中。所以这个额外的设施只是一种幻觉(关于性能)。更准确地说，这是 LAPACKE_dge_trans() function .您可以检查 Lapacke 代码，看看这个函数在 Layout=RowMajor 后几乎无处不在(例如，参见 lapacke_dgesv_work() 代码)。

• 另请注意，如果您想要通用步幅(I 和 J 方向上的“任意前导尺寸”)，您可以使用类似 Blis 的库而不是 Blas。真正的优势是能够创建张量的任意二维 View 。这个选择取决于你的应用，我不知道你是否考虑过张量操​​作。

B.3/矩阵维度

如果你的矩阵总是小到 3x10 blas/lapack 不是一个好的选择(为了性能)。考虑使用类似 Eigen 的库或 Blaz .