LAPACK zgemm op(A)维度

Question

在netlib的这个链接中,它将M指定为:

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

因此,如果我想使用3x10矩阵作为A,但我想使用它的zgmm共轭(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);

在动词中,如果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这样的库.真正的优势是能够创建Tensors的任意2D视图.这个选择取决于你的应用程序,我不知道你是否考虑到张量操作.

B.3 /矩阵尺寸

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