How to efficiently invert only part of an allocated matrix

Question

I have an algorithm that allocates a complex double matrix "A" with a pre-defined size N x N. The elements are initially zero. I also allocated matrix of size N x N to store the inverse, "A_inv". During the algorithm, elements of "A" get filled. At each iteration i, I end up with a submatrix of size i x i. So it looks something like this for the second iteration, where N=4:

| x00 x01 0.0 0.0 |
| x10 x11 0.0 0.0 |
| 0.0 0.0 0.0 0.0 |
| 0.0 0.0 0.0 0.0 |

where the x's indicate some non-zero value. Now I wish to invert the non-zero part of the matrix (a 2x2 matrix in this example). So far I have been doing this in the following way:

1. copy the non-zero elements of "A" to a 2x2 gsl matrix
2. use gsl LU decomposition to invert the 2x2 gsl matrix
3. copy the 2x2 inverted matrix into A_inv

The problem with this approach is that I have to copy over a matrix twice in each iteration. Once to a smaller n x n gsl matrix and once to copy the resulting inverted nxn gsl matrix to A_inv.

I was wondering if anyone knows of a more direct way. Is there a way to use some gsl function to only invert part of a matrix and ignore any zero elements? Say something like this:

A = NxN matrix
A_inv = invert_submatrix(A,n,n)

where n < N. Here invert_submatrix() would only consider the n x n part of A. Furthermore, the original matrix "A" must not be altered by this inversion. Maybe that last demand makes it necessary to copy over the matrix anyway, in which case it would be no more efficient than what I do now. That said, gsl algorithms tend to be a lot more efficient than anything I usually come up with. So thoughts on this are still very welcome.

Answer 1

Sadly as GSL does it's LU decompositions in place I'm not sure that it's possible to do it without copying the submatrix out of A first if you require it not to be modifed. You can, however, use a matrix view to construct the inverse directly from the LU decomposition rather than having to construct it then copy it.

gsl_matrix *invert_submatrix( const gsl_matrix *m, size_t sub_size )
{
    gsl_matrix *inv = gsl_matrix_calloc( m->size1, m->size2 );

    // Create views onto the submatrices we care about 
    gsl_matrix_const_view m_sub_view = 
        gsl_matrix_const_submatrix(m, 0, 0, sub_size, sub_size);

    gsl_matrix_view inv_sub_view = 
        gsl_matrix_submatrix(inv, 0, 0, sub_size, sub_size);

    const gsl_matrix *m_sub = &m_sub_view.matrix;
    gsl_matrix *inv_sub = &inv_sub_view.matrix;

    // Create a matrix for the LU decomposition as GSL does this inplace.
    gsl_permutation *perm = gsl_permutation_alloc(sub_size);
    gsl_matrix *LU = gsl_matrix_alloc(sub_size, sub_size);
    gsl_matrix_memcpy(LU, m_sub);

    int s;
    gsl_linalg_LU_decomp(LU, perm, &s);
    gsl_linalg_LU_invert(LU, perm, inv_sub);

    gsl_matrix_free(LU);
    gsl_permutation_free(perm);

    return inv;
}