Accessing submatrices using LAPACK

Question

Is there a function in LAPACK, which will give me the elements of a particular submatrix? If so how what is the syntax in C++?

Or do I need to code it up?

Answer 1

There is no function for accessing a submatrix. However, because of the way matrix data is stored in LAPACK routines, you don't need one. This saves a lot of copying, and the data layout was (partially) chosen for this reason:

Recall that a dense (i.e., not banded, triangular, hermitian, etc) matrix in LAPACK is defined by four values:

• a pointer to the top left corner of the matrix
• the number of rows in the matrix
• the number of columns in the matrix
• the "leading dimension" of the matrix; typically this is the distance in memory between adjacent elements of a row.

Most of the time, most people only ever use a leading dimension that is equal to the number of rows; a 3x3 matrix is typically stored like so:

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

Suppose instead that we wanted a 3x3 submatrix of a huge matrix with leading dimension lda. Suppose we specifically want the 3x3 submatrix whose top-left corner is located at a(15,42):

         .             .            .
         .             .            .
... 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] ...
         .             .            .
         .             .            .

We could copy this 3x3 matrix into contiguous storage, but if we want to pass it as an input (or output) matrix to an LAPACK routine, we don't need to; we only need to define the parameters appropriately. Let's call this submatrix b; we then define:

// 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;

The only thing that might be surprising is the value of ldb; by using the value lda of the "big matrix", we can address the submatrix without copying, and operate on it in-place.

However I lied (sort of). Sometimes you really can't operate on a submatrix in place, and genuinely need to copy it. I didn't want to talk about that, because it's rare, and you should use in-place operations whenever possible, but I would feel bad not telling you that it is possible. The routine:

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

copies the MxN matrix whose top-left corner is A and is stored with leading dimension LDA to the MxN matrix whose top-left corner is B and has leading dimension LDB. The UPLO parameter indicates whether to copy the upper triangle, lower triangle, or the whole matrix.

In the example I gave above, you would use it like this (assuming the clapack bindings):

...
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)
...