visualdynamics.core.sparse_solve

sparse_solve

A sparse symmetric matrix factored once and solved many times, by the fastest solver the machine has.

The modal solve of a large model factors (K - sigma M) once and solves with it a few hundred times: every Lanczos step and every polish step is a solve. On macOS the factorization is Apple's, in the Accelerate framework every Mac carries (Sparse Solvers: a supernodal Cholesky, or LDLT when the matrix is not positive definite), reached through ctypes from the layouts in its public header. On Windows and Linux on Intel and AMD processors it is PARDISO, the multithreaded supernodal solver in Intel's oneMKL, through pypardiso (Cholesky, or LDLT with Bunch-Kaufman pivoting). Everywhere else, and wherever those refuse a matrix, it is SuperLU's LU through scipy, told the matrix is symmetric so it orders and pivots for that.

Classes:

Name	Description
SuperLUFactor	SuperLU's factorization of a symmetric matrix, through scipy.
AccelerateFactor	Apple's supernodal factorization of a symmetric matrix: Cholesky,
PardisoFactor	PARDISO's supernodal factorization of a symmetric matrix:

Functions:

Name	Description
probe_backward_error	|A x - b| / |b| for one solve with a fixed right-hand side, the
count_below	How many eigenvalues of K phi = lambda M phi lie below shift,
expected_backend	The backend this machine should factor with: 'accelerate' on a
check	Check that this installation factors with the solver it should.
factor_symmetric	The matrix factored by the fastest solver here that takes it:

Classes

SuperLUFactor

SuperLUFactor(matrix: Any)

SuperLU's factorization of a symmetric matrix, through scipy.

Methods:

Name	Description
solve	X with A X = rhs, for one right-hand side or a block of

Source code in src/visualdynamics/core/sparse_solve.py

def __init__(self, matrix: Any) -> None:
    from scipy.sparse.linalg import splu

    # Told the matrix is symmetric (2026-10-18): the default ordering
    # treats it as general, and on a 135,000-DOF model factored in
    # 14 s with 164 M nonzeros in the factors where this takes 3.8 s
    # and 87 M; every solve is half the work as well
    self._factor = splu(sp.csc_matrix(matrix), permc_spec='MMD_AT_PLUS_A',
                        options={'SymmetricMode': True})
    self.shape: tuple[int, int] = self._factor.shape
    residual = probe_backward_error(matrix, self.solve)
    if not residual < PROBE_LIMIT:
        raise RuntimeError('the matrix is singular to working precision: '
                           f'a probe solve has a backward error of {residual:.1e}')

Methods:

solve

solve(rhs: Any) -> ndarray

X with A X = rhs, for one right-hand side or a block of them, in the shape rhs came in.

Source code in src/visualdynamics/core/sparse_solve.py

def solve(self, rhs: Any) -> np.ndarray:
    """X with A X = `rhs`, for one right-hand side or a block of
    them, in the shape `rhs` came in."""
    return self._factor.solve(np.asarray(rhs, dtype=np.float64))

AccelerateFactor

AccelerateFactor(matrix: Any)

Apple's supernodal factorization of a symmetric matrix: Cholesky, or LDLT when the matrix is not positive definite.

Raises ValueError when Accelerate refuses the matrix, naming its status, so the caller can fall back.

Methods:

Name	Description
solve	X with A X = rhs, for one right-hand side or a block of

Source code in src/visualdynamics/core/sparse_solve.py

def __init__(self, matrix: Any) -> None:
    lib = _library()
    if lib is None:
        raise ValueError('Accelerate is not available here')
    lower = sp.tril(matrix, format='csc')
    lower.sort_indices()
    # kept for the factorization's life: Accelerate reads them in
    # place and the factor may refer back to the structure
    self._starts = np.ascontiguousarray(lower.indptr, dtype=np.int64)
    self._rows = np.ascontiguousarray(lower.indices, dtype=np.int32)
    self._values = np.ascontiguousarray(lower.data, dtype=np.float64)
    n = lower.shape[0]
    self.shape: tuple[int, int] = (n, n)
    self._matrix = _SparseMatrix(
        _MatrixStructure(n, n, self._starts.ctypes.data_as(C.POINTER(C.c_long)),
                         self._rows.ctypes.data_as(C.POINTER(C.c_int)),
                         _LOWER_SYMMETRIC, 1),
        self._values.ctypes.data_as(C.POINTER(C.c_double)))
    self._factor = None
    for kind in (_CHOLESKY, _LDLT):
        del lib['messages'][:]
        factored = lib['factor'](kind, C.byref(self._matrix),
                                 C.byref(lib['symbolic']), C.byref(lib['numeric']))
        if factored.status == 0:
            self._factor = factored
            self.kind: str = 'Cholesky' if kind == _CHOLESKY else 'LDLT'
            # LDLT reports success on an exactly singular matrix and
            # its solves come back NaN (measured 2026-10-18: an empty
            # matrix, a zero row, a rank-one block all 'ok'), or
            # finite garbage on one singular to working precision
            # (2026-10-19, a backward error of 0.25), where Cholesky
            # succeeding proves the matrix definite. One probe solve
            # tells, and a refusal here sends the matrix on to
            # SuperLU, whose own probe the modal solve turns into
            # its message
            if kind == _LDLT and not probe_backward_error(
                    matrix, self.solve) < PROBE_LIMIT:
                lib['destroy'](C.byref(factored))
                self._factor = None
                raise ValueError('Accelerate found the matrix singular')
            return
        lib['destroy'](C.byref(factored))
    why = _STATUS.get(factored.status, f'status {factored.status}')
    said = '; '.join(m.strip() for m in lib['messages'] if m.strip())
    raise ValueError(f'Accelerate refused the matrix: {why}'
                     + (f' ({said})' if said else ''))

Methods:

solve

solve(rhs: Any) -> ndarray

X with A X = rhs, for one right-hand side or a block of them, in the shape rhs came in.

Source code in src/visualdynamics/core/sparse_solve.py

def solve(self, rhs: Any) -> np.ndarray:
    """X with A X = `rhs`, for one right-hand side or a block of
    them, in the shape `rhs` came in."""
    lib = _library()
    given = np.asarray(rhs, dtype=np.float64)
    b = np.asfortranarray(given.reshape(given.shape[0], -1))
    n, count = b.shape
    x = np.zeros((n, count), dtype=np.float64, order='F')
    right = _DenseMatrix(n, count, n, 0, b.ctypes.data_as(C.POINTER(C.c_double)))
    out = _DenseMatrix(n, count, n, 0, x.ctypes.data_as(C.POINTER(C.c_double)))
    needed = (self._factor.solveWorkspaceRequiredStatic
              + count * self._factor.solveWorkspaceRequiredPerRHS)
    workspace = C.create_string_buffer(int(needed) + 64)
    lib['solve'](C.byref(self._factor), C.byref(right), C.byref(out),
                 C.cast(workspace, C.c_void_p))
    return x.reshape(given.shape) if given.ndim == 1 else x

PardisoFactor

PardisoFactor(matrix: Any)

PARDISO's supernodal factorization of a symmetric matrix: Cholesky, or LDLT when the matrix is not positive definite.

Raises ValueError when PARDISO is not available or refuses the matrix, so the caller can fall back.

Methods:

Name	Description
solve	X with A X = rhs, for one right-hand side or a block of

Source code in src/visualdynamics/core/sparse_solve.py

def __init__(self, matrix: Any) -> None:
    pypardiso = _pypardiso()
    if pypardiso is None:
        raise ValueError('PARDISO is not available here')
    from pypardiso.pardiso_wrapper import PyPardisoError

    self._upper = _upper_triangle(matrix)
    n = self._upper.shape[0]
    self.shape: tuple[int, int] = (n, n)
    self._solver = None
    failure: Exception | None = None
    for mtype, kind in ((2, 'Cholesky'), (-2, 'LDLT')):
        solver = _pardiso_solver(mtype)
        try:
            solver.factorize(self._upper)
        except (PyPardisoError, ValueError) as refused:
            solver.free_memory(everything=True)
            failure = refused
            continue
        self._solver, self.kind = solver, kind
        # oneMKL's own memory, which Python does not see, freed when
        # the factor is dropped: without it a factor held its
        # gigabytes until exit. A finalizer rather than __del__,
        # which ran during interpreter shutdown, after the modules
        # pypardiso's release calls on were gone
        self._finalizer = weakref.finalize(self, solver.free_memory, True)
        # LDLT perturbs a pivot within round-off of zero and reports
        # success on a singular matrix, as Accelerate's LDLT does;
        # one probe solve tells, and is cheap beside the factor
        if not probe_backward_error(matrix, self.solve) < PROBE_LIMIT:
            self._release()
            raise ValueError('PARDISO found the matrix singular')
        return
    raise ValueError(f'PARDISO refused the matrix: {failure}')

Methods:

solve

solve(rhs: Any) -> ndarray

X with A X = rhs, for one right-hand side or a block of them, in the shape rhs came in.

Source code in src/visualdynamics/core/sparse_solve.py

def solve(self, rhs: Any) -> np.ndarray:
    """X with A X = `rhs`, for one right-hand side or a block of
    them, in the shape `rhs` came in."""
    return self._solver.solve(self._upper, np.asarray(rhs, dtype=np.float64))

Functions:

probe_backward_error

probe_backward_error(matrix: Any, solve: Any) -> float

|A x - b| / |b| for one solve with a fixed right-hand side, the cosines of the row numbers; infinite when the answer is not.

Source code in src/visualdynamics/core/sparse_solve.py

def probe_backward_error(matrix: Any, solve: Any) -> float:
    """|A x - b| / |b| for one solve with a fixed right-hand side, the
    cosines of the row numbers; infinite when the answer is not."""
    probe = np.cos(np.arange(matrix.shape[0], dtype=np.float64))
    answer = solve(probe)
    if not np.isfinite(answer).all():
        return math.inf
    return float(np.linalg.norm(matrix @ answer - probe) / np.linalg.norm(probe))

count_below

count_below(stiffness: Any, mass: Any, shift: float) -> int

How many eigenvalues of K phi = lambda M phi lie below shift, read off the signs of one factorization of K - shift M.

By Sylvester's law of inertia, a symmetric matrix factored as L D Lᵀ has as many negative entries in D as it has negative eigenvalues, and the negative eigenvalues of K - shift M are the modes below the shift (M positive definite). Where PARDISO is available its LDLT reports the count itself; Bunch-Kaufman's 2 x 2 pivots keep the count. Elsewhere SuperLU in symmetric mode with no pivoting off the diagonal is that factorization, its U's diagonal being D, which is why the threshold is zero here and not the factor the solves use. A pivot within round-off of zero, from a shift on top of an eigenvalue, has a sign that is round-off, so a caller allows the count one either way.

Parameters:

Name	Type	Description	Default
stiffness	scipy.sparse matrix	K and M, symmetric, M positive definite.	required
mass	scipy.sparse matrix	K and M, symmetric, M positive definite.	required
shift	float	The eigenvalue to count below, (rad/s)².	required

Returns:

Type	Description
int	The count, rigid-body modes included when the shift is above zero.

Raises:

Type	Description
RuntimeError	If the factorization fails, had to pivot off the diagonal or perturbed a pivot, so its signs are not the inertia.

References

A standard method, used here to size the search up to a frequency:

1. J. J. Sylvester, A demonstration of the theorem that every homogeneous quadratic polynomial is reducible by real orthogonal substitutions to the form of a sum of positive and negative squares, Philosophical Magazine, 1852.

Source code in src/visualdynamics/core/sparse_solve.py

def count_below(stiffness: Any, mass: Any, shift: float) -> int:
    """How many eigenvalues of K phi = lambda M phi lie below `shift`,
    read off the signs of one factorization of K - shift M.

    By Sylvester's law of inertia, a symmetric matrix factored as
    L D Lᵀ has as many negative entries in D as it has negative
    eigenvalues, and the negative eigenvalues of K - shift M are the
    modes below the shift (M positive definite). Where PARDISO is
    available its LDLT reports the count itself; Bunch-Kaufman's 2 x 2
    pivots keep the count. Elsewhere SuperLU in symmetric mode with no
    pivoting off the diagonal is that factorization, its U's diagonal
    being D, which is why the threshold is zero here and not the
    factor the solves use. A pivot within round-off of zero,
    from a shift on top of an eigenvalue, has a sign that is round-off,
    so a caller allows the count one either way.

    Parameters
    ----------
    stiffness, mass : scipy.sparse matrix
        K and M, symmetric, M positive definite.
    shift : float
        The eigenvalue to count below, (rad/s)².

    Returns
    -------
    int
        The count, rigid-body modes included when the shift is above
        zero.

    Raises
    ------
    RuntimeError
        If the factorization fails, had to pivot off the diagonal or
        perturbed a pivot, so its signs are not the inertia.

    References
    ----------
    A standard method, used here to size the search up to a frequency:

    1. J. J. Sylvester, *A demonstration of the theorem that every
       homogeneous quadratic polynomial is reducible by real orthogonal
       substitutions to the form of a sum of positive and negative
       squares*, Philosophical Magazine, 1852.
    """
    from scipy.sparse.linalg import splu

    shifted = sp.csc_matrix(stiffness - float(shift) * mass)
    # 1.5 s where SuperLU took 29 on a 135,000-DOF frame, 5.2 s to 41 on
    # a 106,000-DOF block of bricks, the same counts (2026-11-02)
    if not _superlu_only() and _pypardiso() is not None:
        return _pardiso_negatives(shifted)
    factor = splu(shifted,
                  permc_spec='MMD_AT_PLUS_A',
                  options={'SymmetricMode': True, 'DiagPivotThresh': 0.0})
    if not np.array_equal(factor.perm_r, factor.perm_c):
        raise RuntimeError('the factorization pivoted off the diagonal')
    return int(np.count_nonzero(factor.U.diagonal() < 0.0))

expected_backend

expected_backend() -> str

The backend this machine should factor with: 'accelerate' on a Mac, 'pardiso' on x86-64 Windows and Linux (where pypardiso is a dependency), 'superlu' anywhere else.

Source code in src/visualdynamics/core/sparse_solve.py

def expected_backend() -> str:
    """The backend this machine should factor with: 'accelerate' on a
    Mac, 'pardiso' on x86-64 Windows and Linux (where `pypardiso` is a
    dependency), 'superlu' anywhere else."""
    import platform

    if sys.platform == 'darwin':
        return 'accelerate'
    if platform.machine() in ('x86_64', 'AMD64'):
        return 'pardiso'
    return 'superlu'

check

check() -> tuple[int, str]

Check that this installation factors with the solver it should.

visualdynamics --check-solver runs this and exits with the returned code. The packaged builds are tested with it, because a bundle that leaves oneMKL out still solves, through SuperLU, many times slower on a large model.

Returns:

Type	Description
tuple of (int, str)	0 and the backend, with the library it loaded, when a small model's factorization took expected_backend and its solve and mode count are right; 1 and what happened otherwise.

Source code in src/visualdynamics/core/sparse_solve.py

def check() -> tuple[int, str]:
    """Check that this installation factors with the solver it should.

    `visualdynamics --check-solver` runs this and exits with the
    returned code. The packaged builds are tested with it, because a
    bundle that leaves oneMKL out still solves, through SuperLU, many
    times slower on a large model.

    Returns
    -------
    tuple of (int, str)
        0 and the backend, with the library it loaded, when a small
        model's factorization took `expected_backend` and its solve
        and mode count are right; 1 and what happened otherwise.
    """
    n = 1000
    matrix = sp.diags([np.full(n - 1, -1.0), np.full(n, 2.5), np.full(n - 1, -1.0)],
                      [-1, 0, 1], format='csc')
    expected = expected_backend()
    try:
        factor = factor_symmetric(matrix)
        rhs = np.cos(np.arange(n, dtype=np.float64))
        residual = float(np.linalg.norm(matrix @ factor.solve(rhs) - rhs)
                         / np.linalg.norm(rhs))
        below = count_below(matrix, sp.identity(n, format='csc'), 1.0)
    except Exception as failure:  # noqa: BLE001 — whatever the packaged app would hit
        return 1, f'the factorization failed: {type(failure).__name__}: {failure}'
    # the eigenvalues are 2.5 - 2 cos(k pi / (n + 1)), k = 1 to n
    eigenvalues = 2.5 - 2.0 * np.cos(np.arange(1, n + 1) * np.pi / (n + 1))
    exact = int(np.count_nonzero(eigenvalues < 1.0))
    where = ''
    if factor.backend == 'pardiso':
        from threadpoolctl import threadpool_info

        layers = {pool.get('threading_layer') for pool in threadpool_info()
                  if pool['internal_api'] == 'mkl'}
        where = f' ({factor._solver.libmkl._name}, threaded by {", ".join(sorted(map(str, layers)))})'
    said = f'{factor.backend}{where}: residual {residual:.1e}, {below} of {n} modes below 1.0'
    if factor.backend != expected:
        why = _pardiso.get('why', '') if expected == 'pardiso' else ''
        return 1, f'expected {expected}, got {said}' + (f'; {why}' if why else '')
    if not (residual < 1e-12 and below == exact):
        return 1, f'wrong answers from {said}; {exact} are below 1.0'
    return 0, said

factor_symmetric

factor_symmetric(matrix: Any) -> AccelerateFactor | PardisoFactor | SuperLUFactor

The matrix factored by the fastest solver here that takes it: Accelerate on a Mac, PARDISO on x86-64 Windows and Linux, else SuperLU. Each answers solve(rhs) and shape.

Parameters:

Name	Type	Description	Default
matrix	sparse matrix	Square and symmetric.	required

Returns:

Type	Description
(AccelerateFactor, PardisoFactor or SuperLUFactor)

Source code in src/visualdynamics/core/sparse_solve.py

def factor_symmetric(matrix: Any) -> AccelerateFactor | PardisoFactor | SuperLUFactor:
    """The matrix factored by the fastest solver here that takes it:
    Accelerate on a Mac, PARDISO on x86-64 Windows and Linux, else
    SuperLU. Each answers `solve(rhs)` and `shape`.

    Parameters
    ----------
    matrix : sparse matrix
        Square and symmetric.

    Returns
    -------
    AccelerateFactor, PardisoFactor or SuperLUFactor
    """
    if not _superlu_only():
        for fast in (AccelerateFactor, PardisoFactor):
            try:
                return fast(matrix)
            except ValueError:
                pass
    return SuperLUFactor(matrix)