使用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]

假设我们想要一个具有领先维度的巨大矩阵的3x3子矩阵lda.假设我们特别想要左上角位于(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)

复制MX N,其左上角是矩阵A,并存储与领先的尺寸LDA为MX N,其左上角是矩阵B,并已导致尺寸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)
...