porepy.numerics.linalg.matrix_operations module

module for operations on sparse matrices

block_diag_index(m, n=None)[source]

Get row and column indices for block diagonal matrix

This is intended as the equivalent of the corresponding method in MRST.

Examples: >>> m = np.array([2, 3]) >>> n = np.array([1, 2]) >>> i, j = block_diag_index(m, n) >>> i, j (array([0, 1, 2, 3, 4, 2, 3, 4]), array([0, 0, 1, 1, 1, 2, 2, 2])) >>> a = np.array([1, 3]) >>> i, j = block_diag_index(a) >>> i, j (array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3]), array([0, 1, 1, 1, 2, 2, 2, 3, 3, 3]))

Parameters
Return type

tuple[numpy.ndarray, numpy.ndarray]

block_diag_matrix(vals, sz)[source]

Construct block diagonal matrix based on matrix elements and block sizes.

Parameters

vals: matrix values sz: size of matrix blocks

Returns

sps.csr matrix

Parameters
Return type

spmatrix

copy(A)[source]

Create a new matrix C that is a copy of matrix A This function is equivalent to A.copy(), but does not change the ordering of the A.indices for csc and csr matrices

Parameters:

A (scipy.sparse.spmatrix): A sparce matrix

Return

A (scipy.sparse.spmatrix): A sparce matrix

Parameters

A (spmatrix) –

Return type

spmatrix

csc_matrix_from_blocks(data, block_size, num_blocks)[source]

Create a csc representation of a block diagonal matrix of uniform block size.

The function is equivalent to, but orders of magnitude faster than, the call

sps.block_diag(blocks)

Parameters

  • data (np.array) – Matrix values, sorted column-wise.

  • block_size (int) – The size of all the blocks.

  • num_blocks (int) – Number of blocks to be added.

Returns

csr representation of the block matrix.

Return type

sps.csc_matrix

Raises

ValueError – If the size of the data does not match the blocks size and number of blocks.

Example

>>> data = np.array([1, 2, 3, 4, 5, 6, 7, 8])
>>> block_size, num_blocks = 2, 2
>>> csc_matrix_from_blocks(data, block_size, num_blocks).toarray()
array([[1, 3, 0, 0],
       [2, 4, 0, 0],
       [0, 0, 5, 7],
       [0, 0, 6, 8]])
csr_matrix_from_blocks(data, block_size, num_blocks)[source]

Create a csr representation of a block diagonal matrix of uniform block size.

The function is equivalent to, but orders of magnitude faster than, the call

sps.block_diag(blocks)

Parameters

  • data (np.array) – Matrix values, sorted column-wise.

  • block_size (int) – The size of all the blocks.

  • num_blocks (int) – Number of blocks to be added.

Returns

csr representation of the block matrix.

Return type

sps.csr_matrix

Raises

ValueError – If the size of the data does not match the blocks size and number of blocks.

Example

>>> data = np.array([1, 2, 3, 4, 5, 6, 7, 8])
>>> block_size, num_blocks = 2, 2
>>> csr_matrix_from_blocks(data, block_size, num_blocks).toarray()
array([[1, 2, 0, 0],
       [3, 4, 0, 0],
       [0, 0, 5, 6],
       [0, 0, 7, 8]])
invert_diagonal_blocks(mat, s, method=None)[source]

Invert block diagonal matrix.

Three implementations are available, either pure numpy, or a speedup using numba or cython. If none is specified, the function will try to use numba, then cython. The python option will only be invoked if explicitly asked for; it will be very slow for general problems.

Parameters

mat: sps.csr matrix to be inverted. s: block size. Must be int64 for the numba acceleration to work method: Choice of method. Either numba (default), cython or ‘python’.

Defaults to None, in which case first numba, then cython is tried.

Returns

imat: Inverse matrix

Raises

ImportError: If numba or cython implementation is invoked without numba or

cython being available on the system.

Parameters
Return type

spmatrix

merge_matrices(A, B, lines_to_replace, matrix_format)[source]

Replace rows/coloms of matrix A with rows/cols of matrix B.

If the matrix format is csc, this function is equivalent with

A[:, lines_to_replace] = B

If the matrix format is csr, this funciton is equivalent iwth

A[lines_to_replace, :] = B

Replacement is done in place.

Parameter

A (scipy.sparse.spmatrix): A sparse matrix B (scipy.sparse.spmatrix): A sparse matrix lines_to_replace (ndarray): Lines of A to be replaced by B. matrix_format (str): Should be either ‘csr’ or ‘csc’. Both A and B should adhere

to the respective format.

Return

None

Parameters
Return type

None

optimized_compressed_storage(A)[source]

Choose an optimal storage format (csr or csc) for a sparse matrix.

The format is chosen depending on whether A.shape[0] > A.shape[1] or not.

For very sparse matrices where the number of rows and columns differs significantly (e.g., projection matrices), there can be substantial memory gains by choosing the right storage format, by reducing the number of equal

As an illustration, consider a matrix with shape 1 x N with 1 element: If stored in csc format, this will require an indptr array of size N, while csr format requires only size 2.

Parameters

A (sps.spmatrix) – Matrix to be reformatted.

Returns

The matrix represented in optimal storage format.

Return type

sps.spmatrix

rldecode(A, n)[source]

Decode compressed information in indices.

Example usage: Convert the index pointers in compressed matrix storage (row or column) to a full set of indices.

Acknowledgement: The code is heavily inspired by MRST’s function with the same name, however, requirements on the shape of functions are probably somewhat different.

>>> rldecode(np.array([1, 2, 3]), np.array([2, 3, 1]))
[1, 1, 2, 2, 2, 3]

>>> rldecode(np.array([1, 2]), np.array([1, 3]))
[1, 2, 2, 2]
Parameters

  • A (double, m x k) –

  • 1 (compression should be along dimension) –

  • n (int) – Number of occurences for each element

Returns

The restored array.

Return type

B

See also

rlencode

rlencode(A)[source]

Compress matrix by looking for identical columns.

Example usage: Convert a full set of (row or column) indices of a sparse matrix into compressed storage.

Acknowledgement: The code is heavily inspired by MRST’s function with the same name, however, requirements on the shape of functions are probably somewhat different.

Parameters

A (np.ndarray) – Matrix to be compressed. Should be 2d. Compression will be along the second axis.

Returns

The compressed array, size n x m. np.ndarray: Number of times each row in the first output array should

be repeated to restore the original array.

Return type

np.ndarray

See also

rldecode

slice_indices(A, slice_ind, return_array_ind=False)[source]

Function for slicing sparse matrix along rows or columns. If A is a csc_matrix A will be sliced along columns, while if A is a csr_matrix A will be sliced along the rows.

Parameters

A (scipy.sparse.csc/csr_matrix): A sparse matrix. slice_ind (np.ndarray): Array containing indices to be sliced

Returns

indices (np.ndarray): If A is csc_matrix:

The nonzero row indices or columns slice_ind

If A is csr_matrix:

The nonzero columns indices or rows slice_ind

array_ind (np.ndarray or slice): The indices in the compressed storage format (csc

or csr) corresponding to the slice; so that, if A is csr, A.indices[array_ind] gives the columns of the slice (represented in indices), and the corresponding data can be accessed as A.data[array_ind]. Only returned if return_array_ind is True.

Examples

A = sps.csc_matrix(np.eye(10)) rows = slice_indices(A, np.array([0,2,3]))

Parameters
Return type

Union[ndarray, tuple[numpy.ndarray, Union[numpy.ndarray, slice]]]

slice_mat(A, ind)[source]

Function for slicing sparse matrix along rows or columns. If A is a csc_matrix A will be sliced along columns, while if A is a csr_matrix A will be sliced along the rows.

Parameters

A (scipy.sparse.csc/csr_matrix): A sparse matrix. ind (np.array): Array containing indices to be sliced.

Returns

A_sliced (scipy.sparse.csc/csr_matrix): The sliced matrix

if A is a csc_matrix A_sliced = A[:, ind] if A is a csr_matrix A_slice = A[ind, :]

Examples

A = sps.csc_matrix(np.eye(10)) rows = slice_mat(A, np.array([0,2,3]))

Parameters
Return type

spmatrix

stack_diag(A, B)[source]

Create a new matrix C that contains matrix A and B at the diagonal: C = [[A, 0], [0, B]] This function is equivalent to sps.block_diag((A, B), format=A.format), but does not change the ordering of the A.indices or B.indices

Parameters:

A (scipy.sparse.spmatrix): A sparce matrix B (scipy.sparse.spmatrix): A sparce matrix

Return

None

Parameters
Return type

spmatrix

stack_mat(A, B)[source]

Stack matrix B at the end of matrix A. If A and B are csc matrices this function is equivalent to

A = scipy.sparse.hstack((A, B))

If A and B are csr matrices this function is equivalent to

A = scipy.sparse.vstack((A, B))

Parameters:

A (scipy.sparse.spmatrix): A sparse matrix B (scipy.sparse.spmatrix): A sparse matrix

Return

None

Parameters
zero_columns(A, cols)[source]

Function to zero out columns in matrix A. Note that this function does not change the sparcity structure of the matrix, it only changes the column values to 0.

The matrix is modified in place.

Parameter

A (scipy.sparse.spmatrix): A sparce matrix cols (ndarray): A numpy array of columns that should be zeroed

Return

None

Parameters
Return type

None

zero_rows(A, rows)[source]

Function to zero out rows in matrix A. Note that this function does not change the sparcity structure of the matrix, it only changes the row values to 0.

The matrix is modified in place.

Parameter

A (scipy.sparse.spmatrix): A sparce matrix rows (ndarray): A numpy array of columns that should be zeroed

Return

None

Parameters
Return type

None