使用 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的MxN矩阵，并以前导维度LDA 到左上角为 B 并具有前导维度 LDB 的 MxN 矩阵。 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)
...