使用 LAPACK 访问子矩阵

Question

LAPACK 中是否有一个函数可以为我提供特定子矩阵的元素？如果是这样，C++ 中的语法是什么？

或者我需要编码吗？

Answer 1

没有访问子矩阵的功能。但是，由于矩阵数据存储在 LAPACK 例程中的方式，您不需要一个。这样可以节省大量的复制，并且（部分）选择了数据布局是因为这个原因：

回想一下，LAPACK 中的密集（即非带状、三角形、厄米特等）矩阵由四个值定义：

• 指向矩阵左上角的指针
• 矩阵中的行数
• 矩阵中的列数
• 矩阵的“领先维度”；通常这是内存中一行相邻元素之间的距离。

大多数时候，大多数人只使用等于行数的前导维度； 3x3 矩阵通常存储如下：

a[0] a[3] a[6] 
a[1] a[4] a[7]
a[2] a[5] a[8]

假设我们想要一个具有前导维度 lda 的巨大矩阵的 3x3 子矩阵。假设我们特别想要左上角位于 a(15,42) 的 3x3 子矩阵：

         .             .            .
         .             .            .
... a[15+42*lda] a[15+43*lda] a[15+44*lda] ...
... a[16+42*lda] a[16+43*lda] a[16+44*lda] ...
... a[17+42*lda] a[17+43*lda] a[17+44*lda] ...
         .             .            .
         .             .            .

我们可以将这个 3x3 矩阵复制到连续存储中，但是如果我们想将它作为输入（或输出）矩阵传递给 LAPACK 例程，我们不需要这样做；我们只需要适当地定义参数。我们称这个子矩阵为b；然后我们定义：

// pointer to the top-left corner of b:
float *b = &a[15 + 42*lda];
// number of rows in b:
const int nb = 3;
// number of columns in b:
const int mb = 3;
// leading dimension of b:
const int ldb = lda;

唯一可能令人惊讶的是ldb 的值；通过使用“大矩阵”的值lda，我们可以在不复制的情况下对子矩阵进行寻址，并就地对其进行操作。

但是 我撒了谎（有点）。有时您真的无法在适当的位置上操作子矩阵，并且确实需要复制它。我不想谈论这个，因为它很少见，你应该尽可能使用就地操作，但如果不告诉你这是可能的，我会感到很遗憾。套路：

SLACPY(UPLO,M,N,A,LDA,B,LDB)

将左上角为A并以前导维度LDA存储的MxN矩阵复制到左上角为B的MxN矩阵并且具有领先维度LDB。 UPLO参数表示是复制上三角、下三角还是整个矩阵。

在我上面给出的示例中，您可以这样使用它（假设使用 clapack 绑定）：

...
const int m = 3;
const int n = 3;
float b[9];
const int ldb = 3;
slacpy("A",  // anything except "U" or "L" means "copy everything"
       &m,   // number of rows to copy
       &n,   // number of columns to copy
       &a[15 + 42*lda], // pointer to top-left element to copy
       lda,  // leading dimension of a (something huge)
       b,    // pointer to top-left element of destination
       ldb); // leading dimension of b (== m, so storage is dense)
...