2
$\begingroup$

Let $A$ be an $M \times N$ sparse matrix stored in compressed column format, in a C-like programming environment. I am interested in the best solution to get a sub-matrix of $A$. In MATLAB notation, this would read:

A_sub = A(i:i+m-1, j:j+n-1)

where $A_{sub}$ is again a sparse matrix, but this time being $m \times n$.

$\endgroup$
  • $\begingroup$ I guess solutions should make use of additional data structures such as sparse accumulators to achieve linear time. $\endgroup$ – Tolga Birdal Jan 3 '14 at 10:17
  • $\begingroup$ Linear time in what? Naïve approaches will achieve linear time in the number of nonzero entries. (Not being familiar with sparse accumulators, I'm interested in hearing what advantages they provide.) $\endgroup$ – Geoff Oxberry Jan 3 '14 at 19:02
  • 1
    $\begingroup$ We call this MatGetSubMatrix() in PETSc. The serial algorithm is simple; the parallel implementation is more complicated. $\endgroup$ – Jed Brown Jan 3 '14 at 19:43
4
$\begingroup$

While you mentioned using C, Yousef Saad's package SPARSKIT has a lot of sparse matrix algorithms if you're comfortable with Fortran 77.

The compressed column format consists of three arrays, ja, ia (integers) and a (floats/doubles). For each column index j, ja(j) stores the location in the array ia where the list of all non-zero rows in column j starts.

The first step in extracting a sub-matrix A_sub from A is to count how many non-zero entries there are in each of the n columns; for that purpose, you can make an array count of all zeros and then do something to the effect of:

int k, l;
for (k=0; k<n; k++) {
    for (l=A.ja[k+j]; l<A.ja[k+j+1]; l++) {
        if (A.ia[l]>=i && A.ia[l]<i+m) {
            count[k]++;
        }
    }
}

Next, you'd fill out A_sub.ja as follows:

A_sub.ja[0] = 0;
for (k=0; k<n; k++) {
    A_sub.ja[k+1] = A_sub.ja[k]+count[k];
}

The next step is to fill out the arrays ia and a, which requires a second pass through A.

If you're being space-conscious, instead of making a new array count, you can store the same values starting at the address A_sub.ja[1] with no effect on the output. SPARSKIT does this if you want to see an example.

$\endgroup$
1
$\begingroup$

With the Eigen matrix class library (Eigen), you can efficiently extract a submatrix with a single line:

Eigen::SparseMatrix A(m,n);
Eigen::SparseMatrix B = A.block(startRow, startCol, numRows, numCols);

If A already exists in compressed sparse column (or row) format, you can extract a submatrix with code similar to this:

Eigen::MappedSparseMatrix<double> Amap(A.m, A.n, A.nnz, A.colPtrs, A.rowPtrs, A.vals);
Eigen::SparseMatrix B = Amap.block(startRow, startCol, numRows, numCols);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.