# Linear Algebra Routines¶

## Hierarchical Direct Solver¶

Warning

All the classes and routines in this module are experimental and the API can change at any point.

### Proxy Point Generation¶

class pytential.linalg.BlockProxyPoints(srcindices: pytential.linalg.utils.BlockIndexRanges, indices: pytential.linalg.utils.BlockIndexRanges, points: numpy.ndarray, centers: numpy.ndarray, radii: numpy.ndarray)[source]
srcindices

A BlockIndexRanges describing which block of points each proxy ball was created from.

indices

A BlockIndexRanges describing which proxies belong to which block.

points

A concatenated list of all the proxy points. Can be sliced into using indices (shape (dim, nproxy_per_block * nblocks)).

centers

A list of all the proxy ball centers (shape (dim, nblocks)).

A list of all the proxy ball radii (shape (nblocks,)).

nblocks
nproxy_per_block
to_numpy(actx: arraycontext.impl.pyopencl.PyOpenCLArrayContext) [source]
class pytential.linalg.ProxyGeneratorBase(places, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf')[source]
nproxy

Number of proxy points in a single proxy ball.

__init__(places, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf')[source]
Parameters
• approx_nproxy – desired number of proxy points. In higher dimensions, it is not always possible to construct a proxy ball with exactly this number of proxy points. The exact number of proxy points can be retrieved with nproxy.

• norm_type – type of the norm used to compute the centers of each block. Supported values are "linf" and "l2".

__call__(actx: arraycontext.impl.pyopencl.PyOpenCLArrayContext, source_dd, indices: pytential.linalg.utils.BlockIndexRanges, **kwargs) [source]

Generate proxy points for each block in indices with nodes in the discretization source_dd.

Parameters

source_dd – a DOFDescriptor for the discretization on which the proxy points are to be generated.

class pytential.linalg.ProxyGenerator(places, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf')[source]

A proxy point generator that only considers the points in the current block when determining the radius of the proxy ball.

Inherits from ProxyGeneratorBase.

class pytential.linalg.QBXProxyGenerator(places, approx_nproxy: Optional[int] = None, radius_factor: Optional[float] = None, norm_type: str = 'linf')[source]

A proxy point generator that also considers the QBX expansion when determining the radius of the proxy ball.

Inherits from ProxyGeneratorBase.

pytential.linalg.partition_by_nodes(actx: arraycontext.impl.pyopencl.PyOpenCLArrayContext, discr: meshmode.discretization.Discretization, *, tree_kind: Optional[str] = 'adaptive-level-restricted', max_particles_in_box: Optional[int] = None) [source]

Generate equally sized ranges of nodes. The partition is created at the lowest level of granularity, i.e. nodes. This results in balanced ranges of points, but will split elements across different ranges.

Parameters
pytential.linalg.gather_block_neighbor_points(actx: arraycontext.impl.pyopencl.PyOpenCLArrayContext, discr: meshmode.discretization.Discretization, pxy: pytential.linalg.proxy.BlockProxyPoints, max_particles_in_box: Optional[int] = None) [source]

Generate a set of neighboring points for each range of points in discr. Neighboring points of a range $$i$$ are defined as all the points inside the proxy ball $$i$$ that do not also belong to the range itself.

### Misc¶

class pytential.linalg.BlockIndexRanges(indices: numpy.ndarray, ranges: numpy.ndarray)[source]

Convenience class for working with subsets (blocks) of an array.

nblocks
indices

An ndarray of not necessarily continuous or increasing integers representing the indices of a global array. The individual blocks are delimited using ranges.

ranges

An ndarray of size (nblocks + 1,) consisting of nondecreasing integers used to index into indices. A block $$i$$ can be retrieved using indices[ranges[i]:ranges[i + 1]].

block_size(i: int) int[source]
block_indices(i: int) [source]
Returns

a view into the indices array for the range corresponding to block i.

block_take(x: numpy.ndarray, i: int) [source]
Returns

a subset of x corresponding to the indices in block i. The returned array is a copy (not a view) of the elements of x.

class pytential.linalg.MatrixBlockIndexRanges(row: pytential.linalg.utils.BlockIndexRanges, col: pytential.linalg.utils.BlockIndexRanges)[source]

Convenience class for working with subsets (blocks) of a matrix.

nblocks
row

A BlockIndexRanges encapsulating row block indices.

col

A BlockIndexRanges encapsulating column block indices.

block_shape(i: int, j: int) Tuple[int, int][source]
Returns

the shape of the block (i, j), where i indexes into the rows and j into the cols.

block_indices(i: int, j: int) Tuple[numpy.ndarray, numpy.ndarray][source]
Returns

a view into the indices that make up the block (i, j).

block_take(x: numpy.ndarray, i: int, j: int) [source]
Returns

a subset of the matrix x corresponding to the indices in the block (i, j). The returned array is a copy of the elements of x.

pytential.linalg.make_block_index_from_array(indices: numpy.ndarray, ranges: Optional[numpy.ndarray] = None) [source]

Wrap a (indices, ranges) tuple into a BlockIndexRanges.

Parameters

ranges – if None, then indices is expected to be an object array of indices, so that the ranges can be reconstructed.

pytential.linalg.make_index_blockwise_product(actx: arraycontext.impl.pyopencl.PyOpenCLArrayContext, idx: pytential.linalg.utils.MatrixBlockIndexRanges) Tuple[Any, Any][source]

Constructs a block by block Cartesian product of all the indices in idx.

The indices in the resulting arrays are laid out in C order. Retrieving two-dimensional data for a block diagonal $$i$$ using the resulting index arrays can be done as follows

offsets = np.cumsum( + [
idx.row.block_size(i) * idx.col.block_size(i)
for i in range(idx.nblocks)
])
istart = offsets[i]
iend = offsets[i + 1]

block_1d = x[rowindices[istart:iend], colindices[istart:iend]]
block_2d = block_1d.reshape(*idx.block_shape(i, i))

assert np.allclose(block_2d, idx.block_take(x, i, i))


The result is equivalent to block_take(), which takes the Cartesian product as well.

Returns

a tuple containing (rowindices, colindices), where the type of the arrays is the base array type of actx.