LAPACK zgemm op(A) 尺寸

Question

在来自 netlib 的 this 链接中，它将 M 指定为：

输入时，M 指定矩阵的行数 矩阵 C 的 op( A ) 和 M 必须至少为零。 退出时不变。

因此，如果我想使用 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/宝石

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);

在动词中，如果 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/内存布局

您还必须考虑内存布局。有两种可能：

• 列专业（Fortran 风格），您必须分配一个大小为

  ld*J_size

  的数组，并且 Aᵢⱼ 是

  A[i+ld*j]

  ，其中 0 ≤ i < I_size and 0 ≤ j < J_size

• 行专业（C 风格），您必须分配一个大小为

  I_size*ld

  的数组，并且 Aᵢⱼ 为

  A[j+ld*i]

  ，且 0 ≤ i < I_size and 0 ≤ j < J_size

（其中 ld 是主维）

更新：

• 即使您使用 C++ 编码我建议使用 Fortran 约定（主要列）。 Lapacke 也假装支持行主要模式，但是，在幕后，它只是在调用请求的子例程之前将矩阵“复制”到列主要布局中。所以这个额外的设施只是一个幻觉（关于性能）。更准确地说，这是 LAPACKE_dge_trans() 函数。您可以检查 Lapacke 代码，发现该函数几乎随处使用，就像 Layout=RowMajor 一样（例如，请参阅

  lapacke_dgesv_work()

  代码）。

另请注意，如果您想要通用步幅（I 和 J 方向上的“任意引导尺寸”），您可以使用像
• Blis

  这样的库而不是 Blas。真正的优势是能够创建张量的任意 2D 视图。这个选择取决于你的应用程序，我不知道你是否考虑过张量操作。

B.3/矩阵维度

如果你的矩阵总是小到 3x10 blas/lapack 不是一个好的选择（对于性能）。考虑使用像

Eigen

或 Blaz 这样的库。