[SciPy-user] Construct sparse matrix from sparse blocks

Neilen Marais nmarais@sun.ac...
Fri Feb 8 08:15:14 CST 2008


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 
message. If anyone else finds it useful I think it may be a good idea to 
it somewhere in scipy.sparse?


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), 
    format -- Desired sparse format of output matrix
    dtype -- Desired dtype, defaults to N.float64


    Global matrix containing the input blocks at the desired block
    locations. If csr or csc matrices are requested it is ensured that 
    indices are sorted.


    The 5x5 matrix A containing a 3x3 upper diagonal block, A_aa and 2x2 
    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