(This answer is valid for both MATLAB and Octave, even though I mainly refer to MATLAB)
There are two beasts to slay; but let's first understand the underlying data structure. MATLAB and Octave store sparse matrices in the coordinate (COO) format, i.e. a sparse matrix S is a collection of three arrays of the length equal to the number of nonzero entries of S: row contains the row coordinates of the nonzero entries, col contains the column coordinates of the nonzero entries and finally val contains the values of the nonzero entries. For example, if row(i)=1, col(i)=10 and val(i)=3.14 then we know that the first row - tenth column entry of the sparse matrix S is 3.14 (which also is the i-th nonzero).
We can get this information from MATLAB using the find command:
[i,j,v] = find(S);
Using these arrays, it is trivial to find any banded part of S. We simply look for the entries which are at the distance of n or less from the main diagonal:
idx = find(abs(i-j)<n);
Then the banded part B of S can be constructed:
B = sparse(i(idx),j(idx),v(idx));
So it is just three steps, no conversion to a full matrix is required. These three steps can be replaced by two calls to spdiags: first, use spdiags to extract the diagonal and sub/superdiagonals, then use spdiags together with the extracted data to reconstruct the banded matrix. The spdiags approach works significantly faster in most cases, but the former approach will let us understand the answer to the second part better.
Notice that the n-block diagonal part of a matrix is contained in the n-banded part of the matrix. So a similar strategy may work to efficiently extract the n-block diagonal part of the matrix. The main observation is that depending on the row, there are some column entries we don't want. If it is the first row of the block, we don't want any of the entries to the left of the diagonal entry and if it is the last row of the block, we don't want any of the entries to the right of the diagonal entry. These changes linearly with respect to the row of the block. This should already remind us modular arithmetic, if mod(i(k),n)==0 we are at the last row of a block, so we should not accept any nonzero entry k s.t. j(k)>i(k). We can derive similar conditions for each row, here is a short (somehow cryptic) solution I came up with to summarize these rules:
d = j - i;
idx=find( (d< n-(mod(i-1,n))) .* (d>=-(mod(i-1,n))));
Now, the array idx contains the indices of all the nonzero entries which fall into any diagonal block. We can construct the block diagonal part of S exactly the same way as above:
B = sparse(i(idx),j(idx),v(idx));
I should mention that MATLAB also provides the command blkdiag. This command is more flexible in some sense, because it allows one to create block diagonal matrices with variable block sizes (steps above assume a fixed block size n). However, there is no built-in way to extract diagonal blocks of a sparse matrix. So one would have to extract the block themselves, and then use blkdiag. Even then it is not trivial to get a sparse matrix back; if the input is a cell-structure of dense matrices then MATLAB returns a dense matrix. You should make sure that at least one of them is sparse to guarantee that MATLAB gives you a sparse matrix back. Hence, to my best knowledge, the steps I described above are the best approach to do these operations. Someone might have it implemented in MEX which could run faster or be more memory-efficient, but I couldn't find it on the small section of the internet I looked at.
I would love to hear if anyone has ideas to improve this solution.
