c++ LAPACK zgemm op(A)尺寸

np8igboo  于 2023-11-19  发布在  其他
关注(0)|答案(1)|浏览(115)

在来自netlib的this链接中,它将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++写代码。

ubby3x7f

ubby3x7f1#

**A/CN.9/2005/10/Add.1/Add.1/Add.2/Add.1/Add.2/Add.2/Add.1/Add.2/Add.1/Add.2/Add.

在给出答案之前,为了更清楚地说明这一点,重要的是要记住这一事实:

  • 在美国,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/内存布局

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

  • column major(Fortran风格),你必须分配一个大小为ld*J_size的数组,A的值为A[i+ld*j],0 ≤ i < I_size,0 ≤ j < J_size
  • row major(C风格)你必须分配一个大小为I_size*ld的数组,A的大小为A[j+ld*i],0 ≤ i < I_size,0 ≤ j < J_size

(其中ld是前导维)

更新内容:

  • 即使您使用C++编写代码,我也建议您使用Fortran约定(column major). Lapacke假装也支持row major模式,但是,在引擎盖下,在调用请求的子程序之前,它只是简单地将矩阵复制到列主布局中。所以这种额外的功能只是一种错觉(关于perfs)。更准确地说,这是LAPACKE_dge_transs()函数。你可以检查Lapacke代码,看看这个函数几乎在Layout=RowMajor的任何地方都可以使用(例如,请参阅lapacke_dgesv_work()代码)。
  • 还要注意的是,如果你想要通用的步幅(在I和J方向上的“任意前导维度”),你可以使用像Blis这样的库来代替Blas。真实的优势是能够创建Tensor的任意二维视图。这个选择取决于你的应用,我不知道你是否有Tensor操作的想法。
    B.3/矩阵维数

如果你的矩阵总是像3x 10那么小,那么lapack不是一个好的选择(为了性能)。考虑使用像EigenBlaz这样的库。

相关问题