Linear Algebra Routines¶

In the linear algebra parts of pytential, the following naming scheme is used:

block refers to a piece of a vector operator, e.g. the \(S_{xx}\) component of the Stokeslet.
cluster refers to a piece of a block as used by the recursive proxy-based skeletonization of the direct solver algorithms. Clusters are represented by a TargetAndSourceClusterList.

GMRES¶

pytential.linalg.gmres.gmres(op: Callable[[ArrayOrContainerT], ArrayOrContainerT], rhs: ArrayOrContainerT, restart: int | None = None, tol: float | None = None, x0: ArrayOrContainerT | None = None, inner_product: Callable[[ArrayOrContainerT, ArrayOrContainerT], float] | None = None, maxiter: int | None = None, hard_failure: bool | None = None, no_progress_factor: float | None = None, stall_iterations: int | None = None, callback: Callable[[ArrayOrContainerT], None] | None = None, progress: bool = False, require_monotonicity: bool = True) → GMRESResult[source]¶

Solve a linear system \(Ax = b\) using GMRES with restarts.

Parameters:

op – a callable to evaluate \(A(x)\).
rhs – the right hand side \(b\).
restart – the maximum number of iteration after which GMRES algorithm needs to be restarted
tol – the required decrease in residual norm (relative to the rhs).
x0 – an initial guess for the iteration (a zero array is used by default).
inner_product – a callable with an interface compatible with numpy.vdot() that returns a host scalar.
maxiter – the maximum number of iterations permitted.
hard_failure – if True, raise GMRESError in case of failure.
stall_iterations – number of iterations with residual decrease below no_progress_factor indicates stall. Set to 0 to disable stall detection.

Returns:

a GMRESResult.

class pytential.linalg.gmres.GMRESResult(solution: ArrayContainer, residual_norms: Sequence[float], iteration_count: int, success: bool, state: str)[source]¶

solution¶

residual_norms¶

iteration_count¶

success¶: A bool indicating whether the iteration succeeded.

state¶: A description of the outcome.

exception pytential.linalg.gmres.GMRESError[source]¶

class pytential.linalg.gmres.ResidualPrinter(inner_product=<function structured_vdot>)[source]¶

Hierarchical Direct Solver¶

Note

High-level API for direct solvers is in progress.

Low-level Functionality¶

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.proxy.ProxyPointSource(lpot_source: QBXLayerPotentialSource, proxies: Array)[source]¶

get_expansion_for_qbx_direct_eval(base_kernel: Kernel, target_kernels: Sequence[Kernel]) → ExpansionBase[source]¶: Wrapper around pytential.qbx.QBXLayerPotentialSource.get_expansion_for_qbx_direct_eval to allow this class to be used by the matrix builders.

class pytential.linalg.proxy.ProxyPointTarget(lpot_source: QBXLayerPotentialSource, proxies: Array)[source]¶

class pytential.linalg.proxy.ProxyClusterGeometryData(places: GeometryCollection, dofdesc: DOFDescriptor, srcindex: IndexList, pxyindex: IndexList, points: Array, centers: Array, radii: Array, cluster_radii: Array)[source]¶

The proxy information generated by ProxyGeneratorBase.

places: GeometryCollection¶

dofdesc: DOFDescriptor¶: A descriptor for the geometry used to compute the proxy points.

nclusters¶: Number of clusters.

srcindex: IndexList¶: A IndexList describing which cluster of points each proxy ball was created from.

pxyindex: IndexList¶: A IndexList describing which proxies belong to which cluster.

points: Array¶: A concatenated array of all the proxy points. Can be sliced into using pxyindex (shape (dim, nproxies)).

centers: Array¶: An array of all the proxy ball centers (shape (dim, nclusters)).

radii: Array¶: An array of all the proxy ball radii (shape (nclusters,)).

cluster_radii: Array¶: An array of all the cluster ball radii (shape (nclusters,)). The proxy radii are essentially just the cluster radii multiplied by the radius_factor.

__init__(places: GeometryCollection, dofdesc: DOFDescriptor, srcindex: IndexList, pxyindex: IndexList, points: Array, centers: Array, radii: Array, cluster_radii: Array) → None¶

to_numpy(actx: PyOpenCLArrayContext) → ProxyClusterGeometryData[source]¶

as_sources() → ProxyPointSource[source]¶

as_targets() → ProxyPointTarget[source]¶

class pytential.linalg.proxy.ProxyGeneratorBase(places: GeometryCollection, approx_nproxy: int | None = None, radius_factor: float | None = None, norm_type: str = 'linf', _generate_ref_proxies: Callable[[int], NDArray[floating]] | None = None)[source]¶

places: GeometryCollection¶

radius_factor: float¶: Factor multiplying the cluster radius. This is currently fixed for all clusters.

norm_type: str¶: The type of the norm used to compute the bounding box of the cluster. The “radius” is also computed in terms of this norm. Can be one of "l2" or "linf".

ref_points: NDArray[floating]¶: Reference proxy points on the unit sphere. The reference points are translated and scaled for each cluster.

nproxy¶: Number of proxy points in a single proxy ball.

__init__(places: GeometryCollection, approx_nproxy: int | None = None, radius_factor: float | None = None, norm_type: str = 'linf', _generate_ref_proxies: Callable[[int], NDArray[floating]] | None = None) → None[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.
radius_factor – Factor multiplying the cluster radius (i.e radius of the bounding box) to get the proxy ball radius.
norm_type – type of the norm used to compute the centers of each cluster. Supported values are "linf" and "l2".

__call__(actx: PyOpenCLArrayContext, source_dd: DOFDescriptorLike | None, dof_index: IndexList, **kwargs: ArrayContainer) → ProxyClusterGeometryData[source]¶

Generate proxy points for each cluster in dof_index_set 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.proxy.ProxyGenerator(places: GeometryCollection, approx_nproxy: int | None = None, radius_factor: float | None = None, norm_type: str = 'linf', _generate_ref_proxies: Callable[[int], NDArray[floating]] | None = None)[source]¶

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

Inherits from ProxyGeneratorBase.

class pytential.linalg.proxy.QBXProxyGenerator(places: GeometryCollection, approx_nproxy: int | None = None, radius_factor: float | None = None, norm_type: str = 'linf', _generate_ref_proxies: Callable[[int], NDArray[floating]] | None = None)[source]¶

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

Inherits from ProxyGeneratorBase.

pytential.linalg.proxy.partition_by_nodes(actx: PyOpenCLArrayContext, places: GeometryCollection, *, dofdesc: DOFDescriptorLike | None = None, tree_kind: TreeKind | None = 'adaptive-level-restricted', max_particles_in_box: int | None = None) → IndexList[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:

dofdesc – a DOFDescriptor for the geometry in places which should be partitioned.
tree_kind – if not None, it is passed to boxtree.TreeBuilder.
max_particles_in_box – value used to control the number of points in each partition (and thus the number of partitions). See the documentation in boxtree.TreeBuilder.

pytential.linalg.proxy.gather_cluster_neighbor_points(actx: PyOpenCLArrayContext, pxy: ProxyClusterGeometryData, tgtindex: IndexList | None = None, *, max_particles_in_box: int | None = None) → IndexList[source]¶

Generate a set of neighboring points for each cluster of points in pxy.

Neighboring points of a cluster \(i\) are defined as all the points from tgtindex that are inside the proxy ball \(i\) but outside the cluster itself. For example, given a cluster with radius \(r_s\) and proxy radius \(r_p > r_s\), then we gather all points such that \(r_s < \|\mathbf{x}\| <= r_p\).

Skeletonization¶

class pytential.linalg.skeletonization.SkeletonizationWrangler(exprs: ObjectArray1D[ArithmeticExpression], input_exprs: tuple[var, ...], domains: tuple[DOFDescriptor, ...], context: dict[str, Any], neighbor_cluster_builder: type[ClusterMatrixBuilderBase], weighted_targets: bool, target_proxy_exprs: ObjectArray1D[ArithmeticExpression], proxy_target_cluster_builder: type[ClusterMatrixBuilderBase], weighted_sources: bool, source_proxy_exprs: ObjectArray1D[ArithmeticExpression], proxy_source_cluster_builder: type[ClusterMatrixBuilderBase])[source]¶

nrows¶: Number of output exprs in the operator.

ncols¶: Number of input_exprs in the operator.

exprs¶: An ndarray of shape (nrows,) of expressions (layer potentials) that correspond to the output blocks of the matrix. These expressions are tagged for nearfield neighbor evaluation.

source_proxy_exprs¶

target_proxy_exprs¶: Like exprs, but stripped down for farfield proxy evaluation.

input_exprs¶: A tuple of size (ncols,) of densities that correspond to the input blocks of the matrix.

domains¶: A tuple of the same length as input_exprs defining the domain of each input.

context¶: A dict with additional parameters required to evaluate the expressions.

The following attributes and methods are internal and used for skeletonization.

weighted_sources¶: A flag which if True adds a weight to the source to proxy evaluation. This can only be meaningfully set to False when skeletonizing direct P2P interactions.

weighted_targets¶: A flag which if True adds a weight to the proxy to target evaluation. This can only be meaningfully set to False when skeletonizing direct P2P interactions.

proxy_source_cluster_builder¶

proxy_target_cluster_builder¶: A callable that is used to evaluate farfield proxy interactions. This should follow the calling convention of the constructor to pytential.symbolic.matrix.P2PClusterMatrixBuilder.

neighbor_cluster_builder¶: A callable that is used to evaluate nearfield neighbour interactions. This should follow the calling convention of the constructor to pytential.symbolic.matrix.QBXClusterMatrixBuilder.

evaluate_source_proxy_interaction(actx: PyOpenCLArrayContext, places: GeometryCollection, pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, ibrow: int, ibcol: int) → tuple[NDArray[Any], TargetAndSourceClusterList][source]¶

evaluate_target_proxy_interaction(actx: PyOpenCLArrayContext, places: GeometryCollection, pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, ibrow: int, ibcol: int) → tuple[NDArray[Any], TargetAndSourceClusterList][source]¶

evaluate_source_neighbor_interaction(actx: PyOpenCLArrayContext, places: GeometryCollection, pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, ibrow: int, ibcol: int) → tuple[NDArray[Any], TargetAndSourceClusterList][source]¶

evaluate_target_neighbor_interaction(actx: PyOpenCLArrayContext, places: GeometryCollection, pxy: ProxyClusterGeometryData, nbrindex: IndexList, *, ibrow: int, ibcol: int) → tuple[NDArray[Any], TargetAndSourceClusterList][source]¶

class pytential.linalg.skeletonization.make_skeletonization_wrangler(places: GeometryCollection, exprs: ArithmeticExpression | Sequence[ArithmeticExpression], input_exprs: var | Sequence[var], *, domains: Sequence[Hashable] | None = None, context: dict[str, Any] | None = None, auto_where: Hashable | tuple[Hashable, Hashable] | None = None, _weighted_proxy: bool | tuple[bool, bool] | None = None, _proxy_source_cluster_builder: type[ClusterMatrixBuilderBase] | None = None, _proxy_target_cluster_builder: type[ClusterMatrixBuilderBase] | None = None, _neighbor_cluster_builder: type[ClusterMatrixBuilderBase] | None = None)[source]¶

class pytential.linalg.skeletonization.SkeletonizationResult(L: ObjectArray1D[NDArray[Any]], R: ObjectArray1D[NDArray[Any]], tgt_src_index: TargetAndSourceClusterList, skel_tgt_src_index: TargetAndSourceClusterList, _src_eval_result: _ProxyNeighborEvaluationResult | None = None, _tgt_eval_result: _ProxyNeighborEvaluationResult | None = None)[source]¶

Result of a skeletonization procedure.

A matrix \(A\) can be reconstructed using:

\[A \approx L S R\]

where \(S = A_{I, J}\) for a subset \(I\) and \(J\) of the rows and columns of \(A\), respectively. This applies to each cluster in tgt_src_index. In particular, for a cluster pair \((i, j)\), we can reconstruct the matrix entries as follows

Aij = tgt_src_index.cluster_take(A, i, j)
Sij = skel_tgt_src_index.cluster_take(A, i, j)
Aij_approx = L[i] @ Sij @ R[j]

nclusters¶: Number of clusters that have been skeletonized.

L¶: An object ndarray of size (nclusters,).

R¶: An object ndarray of size (nclusters,).

tgt_src_index¶: A TargetAndSourceClusterList representing the indices in the original matrix \(A\) that have been skeletonized.

skel_tgt_src_index¶: A TargetAndSourceClusterList representing a subset of tgt_src_index, i.e. the skeleton of each cluster of \(A\). These indices can be used to reconstruct the \(S\) matrix.

pytential.linalg.skeletonization.skeletonize_by_proxy(actx: PyOpenCLArrayContext, places: GeometryCollection, tgt_src_index: TargetAndSourceClusterList, exprs: var | Sequence[var], input_exprs: var | Sequence[var], *, domains: Sequence[Hashable] | None = None, context: dict[str, Any] | None = None, approx_nproxy: int | None = None, proxy_radius_factor: float | None = None, id_eps: float | None = None, id_rank: int | None = None, rng: Generator | None = None, max_particles_in_box: int | None = None) → ObjectArray2D[NDArray[Any]][source]¶

Evaluate and skeletonize a symbolic expression using proxy-based methods.

Parameters:

tgt_src_index – a TargetAndSourceClusterList indicating which indices participate in the skeletonization.
exprs – see make_skeletonization_wrangler().
input_exprs – see make_skeletonization_wrangler().
domains – see make_skeletonization_wrangler().
context – see make_skeletonization_wrangler().
approx_nproxy – see ProxyGenerator.
proxy_radius_factor – see ProxyGenerator.
id_eps – a floating point value used as a tolerance in skeletonizing each block in tgt_src_index.
id_rank – an alternative to id_eps, which fixes the rank of each skeletonization.
max_particles_in_box – passed to boxtree.TreeBuilder as necessary.

Returns:

an ndarray of SkeletonizationResult of shape (len(exprs), len(input_exprs)).

Internal Functionality and Utilities¶

Warning

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

Misc¶

class pytential.linalg.utils.InexactT¶: alias of TypeVar(‘InexactT’, bound=inexact)

class pytential.linalg.utils.IndexList(indices: NDArray[integer], starts: NDArray[integer])[source]¶

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

nclusters¶: Number of clusters in the index list.

nindices¶: Total number of indices in the index list.

indices: NDArray[integer]¶: An ndarray of not necessarily continuous or increasing integers representing the indices of a global array. The individual cluster slices are delimited using starts.

starts: NDArray[integer]¶: An ndarray of size (nclusters + 1,) consisting of nondecreasing integers used to index into indices. A cluster \(i\) can be retrieved using indices[starts[i]:starts[i + 1]].

cluster_size(i: int | integer) → int[source]¶

cluster_indices(i: int | integer) → NDArray[integer][source]¶

Returns:: a view into the indices array for the range corresponding to cluster i.

cluster_take(x: NDArray[InexactT], i: int | integer) → NDArray[InexactT][source]¶

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

class pytential.linalg.utils.TargetAndSourceClusterList(targets: IndexList, sources: IndexList)[source]¶

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

nclusters¶: Number of clusters in the index list.

shape¶: Shape of the Cartesian product of the targets and sources.

targets: IndexList¶: An IndexList encapsulating target cluster indices.

sources: IndexList¶: An IndexList encapsulating source cluster indices.

cluster_shape(i: int | integer, j: int | integer) → tuple[int, int][source]¶

Returns:: the shape of the cluster (i, j), where i indexes into the targets and j into the sources.

cluster_indices(i: int | integer, j: int | integer) → tuple[NDArray[integer], NDArray[integer]][source]¶

Returns:: a view into the indices that make up the cluster (i, j).

cluster_take(x: NDArray[InexactT], i: int | integer, j: int | integer) → NDArray[InexactT][source]¶

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

flat_cluster_take(x: NDArray[InexactT], i: int) → NDArray[InexactT][source]¶

Returns:: a subset of an array x corresponding to the indices in the cluster i. Unlike cluster_take(), this method indexes into flattened cluster arrays (see also make_index_cluster_cartesian_product()).

pytential.linalg.utils.make_index_list(indices: NDArray[integer] | ObjectArray1D[NDArray[integer]], starts: NDArray[integer] | None = None) → IndexList[source]¶

Wrap a (indices, starts) tuple into an IndexList.

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

pytential.linalg.utils.make_index_cluster_cartesian_product(actx: PyOpenCLArrayContext, mindex: TargetAndSourceClusterList) → tuple[Array, Array][source]¶

Constructs a cluster by cluster Cartesian product of all the indices in mindex.

The indices in the resulting arrays are laid out in C order. Retrieving two-dimensional data for a cluster \(i\) using the resulting index arrays can be done as follows

offsets = np.cumsum([0] + [
    mindex.targets.cluster_size(i) * mindex.sources.cluster_size(i)
    for i in range(mindex.nclusters)
    ])
istart = offsets[i]
iend = offsets[i + 1]

cluster_1d = x[tgtindices[istart:iend], srcindices[istart:iend]]
cluster_2d = cluster_1d.reshape(*mindex.cluster_shape(i, i))

assert np.allclose(cluster_2d, mindex.cluster_take(x, i, i))

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

Returns:: a tuple containing (tgtindices, srcindices), where the type of the arrays is the base array type of actx.

class pytential.linalg.direct_solver_symbolic.KernelTransformationRemover[source]¶

A mapper that removes the transformations from the kernel of all IntGs in the expression.

This includes source and target derivatives and other such transformations. Any unnecessary kernel arguments are also removed from kernel_arguments.

This mapper is meant to be used in the directs solver for proxy interaction, where it is not possible to evaluate source or target directional derivatives.

class pytential.linalg.direct_solver_symbolic.IntGTermCollector[source]¶

A mapper that removes all non-IntG terms from the expression and all their non-constant factors.

In particular, an expression of the type

\[\sum_{i = 0}^N f_i(\mathbf{x}, \sigma) + \sum_{i = 0}^M c_i g_i(\mathbf{x}) \mathrm{IntG}_i(\mathbf{x})\]

is reduced to

\[\sum_{i = 0}^M c_i \mathrm{IntG}_i(\mathbf{x}).\]

The intended used of this transformation is to allow the evaluation of the proxy interactions in the direct solver for a given expression meant for self-evaluation.

class pytential.linalg.direct_solver_symbolic.DOFDescriptorReplacer(source: DOFDescriptorLike, target: DOFDescriptorLike)[source]¶

Mapper that replaces all the DOFDescriptors in the expression with the given ones.

This mapper is meant to allow for evaluation of proxy interactions in the direct solver when the given expression is already partially (or fully) tagged.

__init__(source: DOFDescriptorLike, target: DOFDescriptorLike) → None[source]¶

Parameters:

source – a descriptor for all expressions to be evaluated on the source geometry.
target – a descriptor for all expressions to be evaluate on the target geometry.