[SciPy-user] Construct sparse matrix from sparse blocks
Neilen Marais
nmarais@sun.ac...
Fri Feb 8 08:15:14 CST 2008
David,
On Thu, 07 Feb 2008 10:45:32 -0500, David Warde-Farley wrote:
> If they're stored as a vector of row indices, a vector of column indices
> and a vector of values (as in the scipy.sparse.coo_matrix ) then
> constructing it should be as straightforward as doing a few array
> concatenations (or copies). This format can then be efficiently
> converted to CSR or CSC with the tocsr() or tocsc() methods, which is
> the format you want it in if you're going to be doing any multiplies,
> etc.
Thanks for the suggestion. I ended up writing the function at the end of
the
message. If anyone else finds it useful I think it may be a good idea to
put
it somewhere in scipy.sparse?
Thanks
Neilen
def merge_sparse_blocks(block_mats, format='coo', dtype=N.float64):
"""
Merge several sparse matrix blocks into a single sparse matrix
Input Params
============
block_mats -- sequence of block matrix offsets and block matrices
such that
block_mats[i] == ((row_offset, col_offset),
block_matrix)
format -- Desired sparse format of output matrix
dtype -- Desired dtype, defaults to N.float64
Output
======
Global matrix containing the input blocks at the desired block
locations. If csr or csc matrices are requested it is ensured that
their
indices are sorted.
Example
=======
The 5x5 matrix A containing a 3x3 upper diagonal block, A_aa and 2x2
lower
diagonal block A_bb:
A = [A_aa 0 ]
[0 A_bb]
A = merge_sparse_blocks( ( ((0,0), A_aa), ((3,3), A_bb)) )
"""
nnz = sum(m.nnz for o,m in block_mats)
data = N.empty(nnz, dtype=dtype)
row = N.empty(nnz, dtype=N.intc)
col = N.empty(nnz, dtype=N.intc)
nnz_o = 0
for (row_o, col_o), bm in block_mats:
bm = bm.tocoo()
data[nnz_o:nnz_o+bm.nnz] = bm.data
row[nnz_o:nnz_o+bm.nnz] = bm.row+row_o
col[nnz_o:nnz_o+bm.nnz] = bm.col+col_o
nnz_o += bm.nnz
merged_mat = sparse.coo_matrix((data, (row, col)), dtype=dtype)
if format != 'coo':
merged_mat = getattr(merged_mat, 'to'+format)()
if format == 'csc' or format == 'csr':
if not merged_mat.has_sorted_indices: merged_mat.sort_indices
()
return merged_mat
>
> David
More information about the SciPy-user
mailing list