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. Everywhere else, and wherever Accelerate refuses 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,

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,
factor_symmetric	The matrix factored, by Accelerate where it is and takes the

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

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). 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 or had to pivot off the diagonal, 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). 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 or had to pivot off the diagonal,
        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

    factor = splu(sp.csc_matrix(stiffness - float(shift) * mass),
                  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))

factor_symmetric

factor_symmetric(matrix: Any) -> AccelerateFactor | SuperLUFactor

The matrix factored, by Accelerate where it is and takes the matrix, else by SuperLU; both answer solve(rhs) and shape.

Parameters:

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

Returns:

Type	Description
AccelerateFactor or SuperLUFactor

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

def factor_symmetric(matrix: Any) -> AccelerateFactor | SuperLUFactor:
    """The matrix factored, by Accelerate where it is and takes the
    matrix, else by SuperLU; both answer `solve(rhs)` and `shape`.

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

    Returns
    -------
    AccelerateFactor or SuperLUFactor
    """
    if os.environ.get(BACKEND_VARIABLE, '').lower() != 'superlu':
        try:
            return AccelerateFactor(matrix)
        except ValueError:
            pass
    return SuperLUFactor(matrix)