[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