matlab - 使用lapack/blas将矩阵的一个子集乘以另一个矩阵

Question

我想使用 dgemm 或任何其他 lapack/blas 函数将矩阵 A 的子集乘以另一个矩阵。我认为由于子矩阵的元素可能不连续，所以如果不将子矩阵复制到另一个空间，我不能直接使用 dgemm。因此，当这个子矩阵本身很大时，它可能非常低效，以至于我认为在 C 中为这个特定问题编写乘法代码可能会更好。由于复制然后使用 lapack/blas 本身，可能根本没有效率。我在 matlab 中使用 lapack/blas 作为 mex 文件。

我的问题是

1- lapack/blas 有没有可以在乘法子矩阵上工作的函数？ 2- 如果不是，是直接编写乘法代码更好还是将子矩阵复制到另一个矩阵并使用 dgemm 更好？

Answer 1

实际上 dgemm 是为子矩阵乘法而设计的。您只需要正确使用每个矩阵的起始指针和参数 LDA、LDB、LDC。

BLAS 的 C 变体是:

void cblas_dgemm (const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE TransA, const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc);

假设你有矩阵:

• A(15x10)
• B(10x20)
• C(15x20)

为Column Major矩阵存储调用dgemm:

cblas_dgemm (CblasColMajor, CblasNoTrans, CblasNoTrans, 15, 20, 10, 1., A, 15, B, 10, 1., C, 15);

假设您需要调用 dgemm 传递子矩阵:

• As(3x2) 从 A 的点 (2,1) 开始
• Bs(2x5) 从 B 的点 (3,5) 开始
• Cs(3x5) 从 C 的点 (4,2) 开始

N、M、K 将变为 3、5、2，但 LDXs 将和上面一样。然后，您必须将正确的指针传递给 dgemm，以便它们指向每个子矩阵的开头。由于您有 C 编号，因此您必须从每个坐标中减去一个。

• As 起点是A + (1+0*15)
• Bs 起点是B + (2+4*10)

• Cs 起点是C + (3+1*15)

  cblas_dgemm (CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 5, 2, 1., A+1, 15, B+42, 10, 1., C+18, 15);
  

N LDA 的想法是说我有一个矩阵 A(LDA,*) 但我将使用上层子矩阵作为(N，*)。在示例中，您不想使用上面的子矩阵，而是使用 A 中的其他子矩阵。在本例中，您创建了一个指向矩阵的新指针 A+1。现在 As 是 A+1 的上子矩阵。

类似地，从 C 调用 Fortran 的原始 dgemm 函数将是

    char NoTrans = `N`;
    int N = 3;
    int M = 5;
    int K = 2;
    int LDA = 15;
    int LDB = 10;
    int LDC = 15;
    double alpha = 1.0;
    double beta = 1.0;
    dgemm (&NoTrans, &NoTrans, N, M, K, alpha, A+1, LDA, B+42, LDB, beta, C+18, LDC);