3

在netlib 的这个链接中，它将 M 指定为：

输入时，M 指定矩阵 op(A) 和矩阵 C 的行数。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++ 编写代码。

4

1 回答 1

6

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'你必须取转置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() 函数。您可以检查 Lapacke 代码，看看这个函数几乎无处不在（例如，Layout=RowMajor参见lapacke_dgesv_work()代码）。

  • 另请注意，如果您想要通用步幅（I 和 J 方向上的“任意前导维度”），您可以使用Blis之类的库而不是Blas。真正的优势是能够创建张量的任意 2D 视图。这个选择取决于您的应用程序，我不知道您是否考虑过张量操​​作。

B.3/ 矩阵尺寸

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

于 2017-09-22T08:36:41.237 回答