vetschn / bsparse Goto Github PK

To achieve better parallelization of our applications, it would be nice to include contiguous tensors in bsparse:
The first two dimensions allow any sparsity format like csr, coo, dia. The sparsity in the higher dimensions is constant with respect to the first dimensions which can be understood as multiple sparse matrices with the same sparsity pattern. This would result in that the datastructure of csr/coo becomes 1*n dimensional and once the corresponding index pointer and indices. Similarly BSR would hold tensors instead of single 2D blocks. The accessing of the third dimension could be handelt in the following way A[i, j, k, ...] where i,j are the indices in the 2D sparse space and k in the first higher dimension.

• extending data structure of fundamental sparse formats, correct handling of slicing and operations
• extending BDIA/BCSR/BCOO

The individual blocks of the FEM matrices are ~99% sparse. After a quick and rough evaluation of the speed of solve and inv, as well as sparse --> dense conversions, the idea should be:

1. Store FEM matrix blocks as scipy.sparse.csr_array
2. Convert block to dense every time an inversion or solve is required (sp.linalg is very slow)

Right now there are a few things in here that require some refactoring. There is a bunch of duplicated code that could probably partly be moved to the base classes. Also, the tests are not done in the best way possible.

• Reduce code duplication
• Abstract things that can be abstracted
• Refactor the unit test suite (code duplication ...)

There are a few mistakes/inconsistencies in the docstrings of various classes and methods. It would be nice to have some pictures demonstrating some of the algorithms etc. It would be worthwhile to look into also including some examples (and check with doctest!)

• Check and update documentation of all documented objects.
• Think of ways to better explain what's going on in the matrices.
• Add examples and doctest.

In bdia:~543 np.find_common_type() depreciated from np 1.25
See np.result_type()

Currently, the input sparray is converted to LIL format as a first step, to allow slicing. The resulting data arrays are thus also in LIL format.

This is unexpected behavior and should be adjusted to convert the subarray slices back to the original sparse format for all matrices.

• Adjust behavior for all BSparse types
• Test

It would be convenient to have a method that converts all saved matrix blocks to dense matrices.

Add transpose (.T) and Hermitian conjugate (.H) to bsparse.

import numpy as np
from quasi import vbsr, vbdia


OFFSET = 2
BLOCKSIZES = [3, 2, 1, 2, 3]
BLOCKSIZES = [size + OFFSET for size in BLOCKSIZES]

vbdia.diag([np.random.rand(size, size) for size in BLOCKSIZES], overlap=OFFSET).show()

This gives unexpected behaviour. Problem is with the .from_spmatrix constructor.

Update

Since transferring this issue here, almost nothing has remained unchanged. The overlap functionality needs to be reimplemented.

• Implement overlap functionality.

The block sparse matrices we deal with usually exhibit some sort of symmetry across the main diagonal, i.e. they are symmetric/Hermitian
$$\mathbf{A}^{T} = \mathbf{A}$$ $$\mathbf{A}^{H} = \mathbf{A}$$

or skew-symmetric/skew-Hermitian
$$\mathbf{A}^T = -\mathbf{A}$$ $$\mathbf{A}^{H} = -\mathbf{A}$$

By implementing a proper data structure, we can save a bunch of memory.

After bit of evaluation it seems like C++ bindings using nanobind with a CMake build system are a sensible/straightforward option. The idea should be to realize the most demanding and repetitive computations (e. g. matrix multiplications and index computations) in the C++ back-end.

• Set up the build system
• Implement C++ routines
  • matmul
  • conversion routines
  • validation routines
• Link to the generated library in the Python front-end

The bsparse interface should support writing matrices to disk in an efficient binary format.

The bsparse.diag constructor is pretty confusing. You should have the option to specify the row_sizes and col_sizes or a uniform block_shape or something.

• Think about how to best integrate this
• Implement it
• Tests it