qscat.linalg¶

Dimension-general sparse linear algebra: Kronecker sums over arbitrary D, a cached sparse LU with a SuperLU and a complex-symmetric MUMPS backend, and the bilinear (non-conjugated) ECS inner product. See docs/physics/mumps-sparse-backend.md.

`kron_sum`	Assemble sum_d I x … x ops[d] x … x I as a CSR matrix.
`c_product`	sum_i a_i b_i -- the bilinear (NOT conjugated) inner product.
`SparseLU`	LU factorization of a square sparse matrix, reusable across solves.
`Ordering`
`ShiftInvertEigs`	The k eigenpairs of a sparse matrix nearest a complex shift.
`default_backend`	Temporarily force the "auto" backend to name within a with block.
`set_default_backend`	Set the default backend SparseLU(backend="auto") resolves to.
`get_default_backend`	The current default backend (see set_default_backend).

Linear-algebra helpers: dimension-general Kronecker sums, cached sparse factorizations, and the exterior-complex-scaling c-product.

Pure linear algebra – nothing here knows about grids, potentials or physics, so it composes with any discretization.

Public API:

kron_sum – sum_d I x … x A_d x … x I for arbitrary D.
c_product – the bilinear (non-conjugated) ECS inner product.
SparseLU – cached sparse LU factorization (factor once, solve many), with fill-in and memory diagnostics and a SuperLU/MUMPS backend switch.
Ordering – the SparseLU/splu fill-reducing ordering literal (“NATURAL”/”MMD_ATA”/”MMD_AT_PLUS_A”/”COLAMD”), re-used by every qscat.core solver’s ordering: kwarg.
ShiftInvertEigs – the k eigenpairs nearest a complex shift (sparse shift-invert Arnoldi on top of SparseLU, reusing its symbolic analysis across shifts).
default_backend / set_default_backend / get_default_backend – process-wide override that SparseLU(backend=”auto”) resolves against, for forcing an internal-SparseLU computation onto one factorization engine (a backend-equivalence check).

See docs/physics/nd-tensor-hamiltonian.md.

class qscat.linalg.ShiftInvertEigs(A, *, k=6, ordering='COLAMD', backend='auto', symmetric=None, ncv=None, tol=0.0, maxiter=None)[source]¶

The k eigenpairs of a sparse matrix nearest a complex shift.

Holds the factorization: the first near call factors A - sigma*I, and every later call refactor`s it, reusing the symbolic analysis (the shifted matrix’s sparsity pattern does not depend on `sigma). On the MUMPS backend that skips the SCOTCH ordering; on scipy it re-runs splu (correct, no reuse) – the same trade SparseLU.refactor documents.

ordering, backend and symmetric are forwarded verbatim to SparseLU. k, ncv, tol and maxiter are ARPACK controls; k may be overridden per call.

Parameters:

A (scipy.sparse.spmatrix) – Square matrix. Promoted to complex CSC internally.
k (int, optional) – Number of eigenpairs to return per call (default 6).
ordering ({"COLAMD", "NATURAL", "MMD_ATA", "MMD_AT_PLUS_A"}, optional) – SuperLU column ordering, forwarded to SparseLU.
backend ({"auto", "scipy", "mumps"}, optional) – Factorization backend, forwarded to SparseLU.
symmetric (bool or None, optional) – Complex-symmetry flag, forwarded to SparseLU (auto-detected if None).
ncv (int or None, optional) – Krylov subspace size. ARPACK’s default is used when None.
tol (float, optional) – ARPACK relative tolerance; 0.0 means machine precision.
maxiter (int or None, optional) – Maximum ARPACK restarts.

Notes

Eigenvectors are returned with numpy’s Euclidean (v^dagger v = 1) normalization, exactly as qscat.dvr.eigen returns them. For exterior-complex-scaling observables re-normalize under the bilinear c_product.

Examples

>>> import numpy as np, scipy.sparse as sp
>>> from qscat.linalg import ShiftInvertEigs
>>> A = sp.diags(np.array([0.0, 1.0, 2.0, 10.0, 11.0], dtype=complex))
>>> vals, vecs = ShiftInvertEigs(A, k=2).near(9.5 + 0.0j)
>>> np.round(np.sort_complex(vals).real, 6)
array([10., 11.])

near(sigma, *, k=None)[source]¶

The k eigenpairs of A nearest sigma, nearest first.

Parameters:

sigma (complex) – The shift: eigenvalues are returned in order of increasing |E - sigma|.
k (int or None, optional) – Number of eigenpairs for this call; the constructor’s k if None.

Returns:

energies (ndarray of complex128, shape (k,)) – Eigenvalues of A, sorted by |E - sigma| ascending.
vectors (ndarray of complex128, shape (n, k)) – vectors[:, i] is the eigenvector of energies[i], Euclidean- normalized.

Raises:

ValueError – If k is not in 1 <= k < n - 1 (an ARPACK constraint).
qscat.exceptions.ConvergenceError – If ARPACK fails to converge.

Return type:

tuple[NDArray[complex128], NDArray[complex128]]

property shape: tuple[int, int]¶: Operator shape (n, n).

property n_factorizations: int¶: How many shifted matrices have been factored (analysis + refactors).

property backend_used: str¶: Which factorization engine ran (“scipy” or “mumps”).

property ordering_used: str¶: The ordering the factorization actually used.

property fill_factor: float¶: Factor nnz relative to the shifted matrix’s nnz.

memory_bytes()[source]¶

Factor memory. NOT cheap – a method, not a property, so the cost is opt-in.

Return type:: int

class qscat.linalg.SparseLU(A, *, ordering='COLAMD', backend='auto', symmetric=None)[source]¶

LU factorization of a square sparse matrix, reusable across solves.

ordering is scipy’s permc_spec: one of “NATURAL”, “MMD_ATA”, “MMD_AT_PLUS_A”, “COLAMD” (the default). For a structurally symmetric pattern – which a Kronecker-sum Hamiltonian has – “MMD_AT_PLUS_A” is often the better choice; measure with fill_factor before assuming.

A real-valued A is silently promoted: the internal CSC conversion always uses dtype=np.complex128, so values are preserved but memory doubles.

fill_factor is cheap at any scale (reads SuperLU’s own nnz count, no array materialization). memory_bytes() is NOT cheap – it is a method, not a property, precisely so that its cost is opt-in rather than hidden behind attribute access – and its cache is permanent for this object’s lifetime. See the module docstring before calling it on a production-size matrix.

backend selects the factorization engine: “scipy” is the SuperLU path above; “mumps” is the complex-symmetric MUMPS factorization (available only where system MUMPS + the qscat[mumps] extra are installed – the Docker image, not a bare Mac – and raising RuntimeError if forced when absent); “auto” (the default) prefers MUMPS when available and falls back to scipy otherwise, so on a MUMPS-less box “auto” and “scipy” are identical in every observable way. backend_used reports which one actually ran. An “auto” call site also consults the context-local override set by set_default_backend / the default_backend context manager (an explicit “scipy”/”mumps” here overrides it) – the seam used to force an entire computation that builds SparseLU internally onto one engine for a backend-equivalence check.

symmetric, if left None, is auto-detected as A == A.T to a SCALED tolerance (an O(nnz) sparse comparison: (abs(A - A.T)).max() <= _SYM_RTOL * abs(A).max(), cheap relative to the factorization itself). The tolerance – not exact equality – is deliberate: the N2 decks this class factors are A = A.T mathematically but symmetric only to round-off (Kronecker-sum float reordering; ~3.6e-17 relative asymmetry), so exact equality would reject them and forfeit the whole point of the MUMPS backend (see _SYM_RTOL). The flag is informational only on the scipy path – SuperLU does not exploit symmetry – but on the MUMPS path it selects the complex-symmetric SYM=2 matrix type (upper triangle only) instead of the general unsymmetric SYM=0 one. Pass an explicit symmetric=True/False to override the auto-detect entirely.

Parameters:

A (sp.spmatrix)
ordering (Ordering)
backend (_Backend)
symmetric (bool | None)

refactor(A_new)[source]¶

Re-factorize A_new reusing this object’s symbolic analysis.

A_new MUST share the original matrix’s sparsity pattern (e.g. a diagonal shift E*I - H across energies). On the MUMPS backend this reuses the analysis (skips re-ordering); on scipy it re-runs splu (correct, no reuse). Keeps the original backend and symmetry decision. Raises ValueError on a shape/pattern mismatch.

Parameters:: A_new (spmatrix)
Return type:: None

property shape: tuple[int, int]¶: Shape of the factored matrix.

property ordering: str¶: The permc_spec column ordering this factorization was built with.

property symmetric: bool¶

Whether A was treated as (complex-)symmetric A == A.T.

Auto-detected from A (to the scaled _SYM_RTOL tolerance) when symmetric=None was passed (the default), or the explicit override. Informational on the scipy path; the MUMPS path uses it to select the complex-symmetric (SYM=2) matrix type instead of the general unsymmetric one.

property backend_used: str¶

“scipy” or “mumps”.

Type:: Which backend actually factorized A

property ordering_used: str¶

The ordering the active backend actually used.

On the scipy path this is scipy’s permc_spec (identical to ordering); on the MUMPS path it is MUMPS’s own chosen ordering (e.g. “scotch”/”metis”/”amd”), read from INFOG(7).

property fill_factor: float¶

(L.nnz + U.nnz) / A.nnz – how much denser the factors are.

Cheap at any scale: on the scipy backend, self._lu.nnz is SuperLU’s own reported L+U nonzero count, read directly off the internal factorization – this NEVER materializes L or U as arrays and costs no extra memory (measured delta < 0.1 MB on an N=6000 matrix with a x300 fill-in). Contrast memory_bytes(), which does materialize them and is priced accordingly – see the module docstring.

memory_bytes()[source]¶

Bytes actually held by the L and U factors (data + index arrays).

NOT CHEAP – a method, not a property, because computing this forces scipy to materialize self._lu.L and self._lu.U as full CSC arrays, which SuperLU then caches for this object’s lifetime (measured: a second call allocates no further memory, and the cache is not released by dropping your own references to the result – only deleting this SparseLU instance does). Read the module docstring’s production-scale estimate (+6 GB on the N2 2-D deck) before calling this on anything but a reduced/test-scale matrix.

Return type:: int

solve(b)[source]¶

Solve A x = b for one (N,) or several (N, k) right-hand sides.

Parameters:: b (NDArray[complex128])
Return type:: NDArray[complex128]

qscat.linalg.c_product(a, b)[source]¶

sum_i a_i b_i – the bilinear (NOT conjugated) inner product.

Shapes are compared BEFORE flattening: c_product(psi_(n0, n1), chi_(n1, n0)) with n0 != n1 raises rather than silently ravelling both down to the same total element count and returning a plausible-looking but physically wrong number (a transposed-axis bug hiding behind a reshape that ravel() alone would never catch, since ravel only cares about total size, not per-axis shape).

Parameters:

a (ArrayLike)
b (ArrayLike)

Return type:

complex

qscat.linalg.default_backend(name)[source]¶

Temporarily force the “auto” backend to name within a with block.

The recommended way to steer internal SparseLU(…) construction – e.g. forcing a whole computation that builds SparseLU internally (projects.n2_2d_cross_section.ve_cross_section_2d) through one specific factorization backend, so two backends can be compared for physics equivalence without threading a backend= kwarg through every call site. Scoped, exception-safe, and context-local, so concurrent threads cannot race each other’s defaults.

Parameters:: name (Literal['auto', 'scipy', 'mumps'])
Return type:: Iterator[None]

qscat.linalg.get_default_backend()[source]¶

The current default backend (see set_default_backend).

Return type:: Literal[‘auto’, ‘scipy’, ‘mumps’]

qscat.linalg.kron_sum(ops)[source]¶

Assemble sum_d I x … x ops[d] x … x I as a CSR matrix.

Each ops[d] must be square. The result is square with dimension prod(n_d). D == 1 returns ops[0] unchanged (as CSR).

The identity blocks used to pad each ops[d] are always built with dtype=complex (np.complex128), regardless of the input operators’ dtype – so a purely real problem (all-real ops, no ECS tail on any axis) still pays complex128 memory for the assembled sum, doubling it relative to a real-only assembly. This mirrors SparseLU’s documented real-to-complex promotion, for the same reason: every real use of this library has at least one ECS-tailed axis, so a real fast path would be untested and this module does not carry one.

Parameters:: ops (Sequence[spmatrix])
Return type:: csr_matrix

qscat.linalg.set_default_backend(name)[source]¶

Set the default backend SparseLU(backend=”auto”) resolves to.

HAZARD: this mutates the CURRENT context for the rest of the process (or thread/task) lifetime and is easy to leave flipped – prefer the default_backend context manager, which restores the previous value on exit (including on exception). Only “auto” call sites consult this; an explicit backend=”scipy”/”mumps” argument always wins.

Parameters:: name (Literal['auto', 'scipy', 'mumps'])
Return type:: None