Source code for sfepy.linalg.sparse

"""Some sparse matrix utilities missing in scipy."""
from __future__ import absolute_import
import numpy as nm
import scipy.sparse as sp

from sfepy.base.base import assert_
from six.moves import range

[docs] def save_sparse_txt(filename, mtx, fmt='%d %d %f\n'): """Save a CSR/CSC sparse matrix into a text file""" fd = open(filename, 'w') fd.write('%d %d\n' % mtx.shape) fd.write('%d\n' % mtx.size) if mtx.format == 'csr': indptr, indices, data = mtx.indptr, mtx.indices, mtx.data for ir in range(mtx.shape[0]): for ii in range(indptr[ir], indptr[ir+1]): fd.write(fmt % (ir, indices[ii], data[ii])) elif mtx.format == 'csc': indptr, indices, data = mtx.indptr, mtx.indices, mtx.data for ic in range(mtx.shape[0]): for ii in range(indptr[ir], indptr[ir+1]): fd.write(fmt % (indices[ii], ic, data[ii])) else: raise ValueError('matrix format not supported! (%s)' % mtx.format)
[docs] def insert_sparse_to_csr(mtx1, mtx2, irs, ics): """ Insert a sparse matrix `mtx2` into a CSR sparse matrix `mtx1` at rows `irs` and columns `ics`. The submatrix `mtx1[irs,ics]` must already be preallocated and have the same structure as `mtx2`. """ import sfepy.discrete.common.extmods.assemble as asm if isinstance(irs, slice): irs = nm.arange(irs.start, irs.stop, irs.step, dtype=nm.int32) if isinstance(ics, slice): ics = nm.arange(ics.start, ics.stop, ics.step, dtype=nm.int32) n_row, n_col = mtx1.shape assert_((irs.min() >= 0) and (irs.max() < n_row)) assert_((ics.min() >= 0) and (ics.max() < n_col)) aux = mtx2.tocoo() data = nm.ascontiguousarray(aux.data[:,None,None,None]) rows = irs[aux.row[:,None]] cols = ics[aux.col[:,None]] iels = nm.arange(rows.shape[0], dtype=nm.int32) asm.assemble_matrix(mtx1.data, mtx1.indptr, mtx1.indices, data, iels, 1.0, rows, cols)
def _normalize_sizes(sizes): """ Checks whether all the sizes are either slices or not. Transforms slices into their sizes. """ out = [] ns = 0 for size in sizes: if isinstance(size, slice): size = size.stop - size.start ns += 1 else: size = int(size) out.append(size) if ns: if ns != len(sizes): raise ValueError('cannot mix sizes with slices! (%s)' % (sizes,)) is_slice = True else: is_slice = False return out, is_slice
[docs] def compose_sparse(blocks, row_sizes=None, col_sizes=None): """ Compose sparse matrices into a global sparse matrix. Parameters ---------- blocks : sequence of sequences The sequence of sequences of equal lengths - the individual sparse matrix blocks. The integer 0 can be used to mark an all-zero block, if its size can be determined from the other blocks. row_sizes : sequence, optional The required row sizes of the blocks. It can be either a sequence of non-negative integers, or a sequence of slices with non-negative limits. In any case the sizes have to be compatible with the true block sizes. This allows to extend the matrix shape as needed and to specify sizes of all-zero blocks. col_sizes : sequence, optional The required column sizes of the blocks. See `row_sizes`. Returns ------- mtx : coo_matrix The sparse matrix (COO format) composed from the given blocks. Examples -------- Stokes-like problem matrix. >>> import scipy.sparse as sp >>> A = sp.csr_matrix([[1, 0], [0, 1]]) >>> B = sp.coo_matrix([[1, 1]]) >>> K = compose_sparse([[A, B.T], [B, 0]]) >>> print K.todense() [[1 0 1] [0 1 1] [1 1 0]] """ if not len(blocks): raise ValueError('no matrix blocks!') if row_sizes is None: row_sizes = nm.array([-1] * len(blocks)) else: assert_(len(row_sizes) == len(blocks)) if col_sizes is None: col_sizes = nm.array([-1] * len(blocks[0])) else: assert_(len(col_sizes) == len(blocks[0])) rs, is_slice_r = _normalize_sizes(row_sizes) cs, is_slice_c = _normalize_sizes(col_sizes) for ir, row in enumerate(blocks): for ic, mtx in enumerate(row): if isinstance(mtx, int) and (mtx == 0): continue if ic >= len(col_sizes): raise ValueError('invalid row size at (%d, %d)!' % (ir, ic)) if rs[ir] == -1: rs[ir] = mtx.shape[0] elif rs[ir] != mtx.shape[0]: msg = 'incompatible matrix block row size at (%d, %d)!' \ % (ir, ic) raise ValueError(msg) if cs[ic] == -1: cs[ic] = mtx.shape[1] elif cs[ic] != mtx.shape[1]: msg = 'incompatible matrix block column size at (%d, %d)!' \ % (ic, ic) raise ValueError(msg) if nm.any(rs == -1): raise ValueError('incomplete row block sizes! (%s)' % row_sizes) if nm.any(cs == -1): raise ValueError('incomplete column block sizes! (%s)' % col_sizes) if is_slice_r: n_row = row_sizes[-1].stop row_offsets = nm.r_[[ii.start for ii in row_sizes], n_row] else: row_offsets = nm.cumsum(nm.r_[0, rs]) n_row = row_offsets[-1] if is_slice_c: n_col = col_sizes[-1].stop col_offsets = nm.r_[[ii.start for ii in col_sizes], n_col] else: col_offsets = nm.cumsum(nm.r_[0, cs]) n_col = col_offsets[-1] rows = [] cols = [] datas = [] for ir, row in enumerate(blocks): for ic, mtx in enumerate(row): if isinstance(mtx, int) and (mtx == 0): continue aux = sp.coo_matrix(mtx) rows.append(aux.row + row_offsets[ir]) cols.append(aux.col + col_offsets[ic]) datas.append(aux.data) rows = nm.concatenate(rows) cols = nm.concatenate(cols) datas = nm.concatenate(datas) mtx = sp.coo_matrix((datas, (rows, cols)), shape=(n_row, n_col)) return mtx
[docs] def infinity_norm(mtx): """ Infinity norm of a sparse matrix (maximum absolute row sum). Parameters ---------- mtx : spmatrix or array The sparse matrix. Returns ------- norm : float Infinity norm of the matrix. Notes ----- - This serves as an upper bound on spectral radius. - CSR and CSC avoid copying `indices` and `indptr` arrays. - inspired by PyAMG See Also -------- scipy.linalg.norm : dense matrix norms """ ones = nm.ones(mtx.shape[1], dtype=mtx.dtype) if sp.isspmatrix_csr(mtx) or sp.isspmatrix_csc(mtx): # Avoid copying index and ptr arrays. abs_mtx = mtx.__class__((nm.abs(mtx.data), mtx.indices ,mtx.indptr), shape=mtx.shape) norm = (abs_mtx * ones).max() elif sp.isspmatrix(mtx): norm = (abs(mtx) * ones).max() else: norm = nm.dot(nm.abs(mtx), ones).max() return norm