Skip to content

openseespy-solvers

scipy

gaaraujo/openseespy-solvers

openseespy_solvers.scipy¶

CPU sparse linear and eigen solver constructors.

This module follows the naming and keyword conventions of scipy.sparse.linalg where possible. The main difference is that constructors return OpenSeesPy-compatible solver objects; OpenSeesPy supplies A, b, K, and M at solve time.

Solving Linear Problems¶

Direct methods:

Constructor	`scipy` analogue	Description
`spsolve`	`scipy.sparse.linalg.spsolve` / `splu`	Sparse direct solver using SuperLU
`umfpack`	`scikits.umfpack`	Sparse direct solver using UMFPACK

Iterative methods:

Constructor	`scipy` analogue	Description
`cg`	`scipy.sparse.linalg.cg`	Conjugate Gradient
`gmres`	`scipy.sparse.linalg.gmres`	Generalized Minimal Residual

Eigenvalue Problems¶

Constructor	`scipy` analogue	Description
`eigsh`	`scipy.sparse.linalg.eigsh`	ARPACK-based generalized symmetric eigen solve
`lobpcg`	`scipy.sparse.linalg.lobpcg`	Locally Optimal Block Preconditioned Conjugate Gradient

OpenSeesPy Usage¶

from openseespy_solvers.scipy import spsolve

solver = spsolve()
ops.system("PythonSparse", solver.to_openseespy())

For eigen analysis:

from openseespy_solvers.scipy import eigsh

solver = eigsh(tol=1e-8)
lam = ops.eigen("PythonSparse", num_modes, solver.to_openseespy())

Notes¶

spsolve and umfpack are both sparse direct solvers. spsolve uses scipy.sparse.linalg's SuperLU path and is available with the base install. umfpack requires the optional scikit-umfpack package:

python -m pip install "openseespy-solvers[umfpack]"

On Windows, install scikit-umfpack from conda-forge instead.

Parameters documented below follow scipy naming where possible. Extra parameters such as scheme, writable, debug, and dtype control the OpenSeesPy integration.

Function Reference¶

Sparse linear algebra solvers for OpenSeesPy (scipy.sparse.linalg backend).

This module provides solver objects that wrap :mod:scipy.sparse.linalg for OpenSeesPy's PythonSparse linear and eigen commands. Factory signatures match scipy.sparse.linalg except that A and b (or K and M) are supplied by OpenSees at solve time.

Each factory returns a solver object. Register it with OpenSees using :meth:~openseespy_solvers._base.BaseOpenSeesSolver.to_openseespy::

from openseespy_solvers.scipy import cg
solver = cg(rtol=1e-8, M=precond.jacobi)
ops.system("PythonSparse", solver.to_openseespy())

Submodules

precond Preconditioner factories for use with M=.

See Also

scipy.sparse.linalg Underlying sparse linear algebra routines. openseespy_solvers.cupy GPU backend with the same interface.

spsolve ¶

spsolve(*, permc_spec: str = 'COLAMD', scheme: str | None = None, writable: str | list[str] = 'none', debug: bool = False, dtype: Any = float64) -> _SpSolve

Source code in src/openseespy_solvers/scipy/__init__.py

def spsolve(
    *,
    permc_spec: str = "COLAMD",
    scheme: str | None = None,
    writable: str | list[str] = "none",
    debug: bool = False,
    dtype: Any = np.float64,
) -> _SpSolve:
    r"""Configure a sparse direct solver for OpenSees ``PythonSparse``.

    Uses :func:`scipy.sparse.linalg.splu` internally. The LU factorization is
    reused when OpenSees reports ``matrix_status='UNCHANGED'``.

    Parameters
    ----------
    permc_spec : str, optional
        Column permutation applied during LU factorization. Passed to
        :func:`scipy.sparse.linalg.splu`. Default is ``'COLAMD'``.
    """ + _OPENSEES_LINEAR + _LINEAR_RETURNS + _LINEAR_NOTES + """
    See Also
    --------
    scipy.sparse.linalg.spsolve
    scipy.sparse.linalg.splu

    Examples
    --------
    >>> from openseespy_solvers.scipy import spsolve
    >>> solver = spsolve(permc_spec="COLAMD")
    >>> cfg = solver.to_openseespy()
    >>> cfg["scheme"]
    'CSR'
    """
    return _SpSolve(
        permc_spec=permc_spec, scheme=scheme or "CSR", writable=writable, debug=debug, dtype=dtype
    )

umfpack ¶

umfpack(*, scheme: str | None = None, writable: str | list[str] = 'none', debug: bool = False, dtype: Any = float64) -> _Umfpack

Source code in src/openseespy_solvers/scipy/__init__.py

def umfpack(
    *,
    scheme: str | None = None,
    writable: str | list[str] = "none",
    debug: bool = False,
    dtype: Any = np.float64,
) -> _Umfpack:
    r"""Configure a 64-bit UMFPACK direct solver for OpenSees ``PythonSparse``.

    Wraps ``scikits.umfpack.UmfpackContext("dl")`` (the double / 64-bit-integer
    UMFPACK context) to solve ``Ax = b``. The symbolic factorization is reused
    while the sparsity structure is unchanged and is recomputed only when
    OpenSees reports ``matrix_status='STRUCTURE_CHANGED'``; the numeric
    factorization is refreshed on ``'COEFFICIENTS_CHANGED'`` and reused on
    ``'UNCHANGED'``.

    Unlike :func:`spsolve` (SuperLU), this uses the UMFPACK library and 64-bit
    (``int64``) indices, which suits very large systems whose index count
    exceeds the 32-bit range. Requires the optional ``scikit-umfpack`` package
    (``python -m pip install "openseespy-solvers[umfpack]"``); the import is deferred until
    this factory is called.

    Parameters
    ----------
    """ + _OPENSEES_LINEAR + _LINEAR_RETURNS + _LINEAR_NOTES + """
    The OpenSees CSR buffers are assembled into a `scipy.sparse` matrix and converted to
    CSC for UMFPACK, so non-symmetric systems are solved correctly.

    Raises
    ------
    ImportError
        If ``scikit-umfpack`` is not installed (raised when the solver is
        instantiated, not on import).

    See Also
    --------
    spsolve
    scipy.sparse.linalg.spsolve

    Examples
    --------
    >>> from openseespy_solvers.scipy import umfpack  # doctest: +SKIP
    >>> solver = umfpack()  # doctest: +SKIP
    >>> solver.to_openseespy()["scheme"]  # doctest: +SKIP
    'CSR'
    """
    return _Umfpack(scheme=scheme or "CSR", writable=writable, debug=debug, dtype=dtype)

cg ¶

cg(x0: Any = None, *, rtol: float = 1e-05, atol: float = 0.0, maxiter: int | None = None, M: Any = None, callback: Any = None, scheme: str | None = None, writable: str | list[str] = 'none', debug: bool = False, dtype: Any = float64) -> _CG

Source code in src/openseespy_solvers/scipy/__init__.py

def cg(
    x0: Any = None,
    *,
    rtol: float = 1e-5,
    atol: float = 0.0,
    maxiter: int | None = None,
    M: Any = None,
    callback: Any = None,
    scheme: str | None = None,
    writable: str | list[str] = "none",
    debug: bool = False,
    dtype: Any = np.float64,
) -> _CG:
    r"""Configure a Conjugate Gradient solver for OpenSees ``PythonSparse``.

    Uses :func:`scipy.sparse.linalg.cg` internally to solve ``Ax = b`` for
    symmetric, positive-definite ``A``.

    Parameters
    ----------
    x0 : ndarray, optional
        Starting guess for the solution. Forwarded to `scipy.sparse.linalg` when not ``None``.
    rtol, atol : float, optional
        Relative and absolute tolerances for convergence. For convergence,
        ``norm(b - A @ x) <= max(rtol*norm(b), atol)`` must hold. Defaults are
        ``rtol=1e-5`` and ``atol=0.0``.
    maxiter : int, optional
        Maximum number of iterations.
    M : {sparse matrix, LinearOperator, callable}, optional
        Preconditioner for ``A``. If ``M`` is callable, it must accept the
        assembled matrix ``A`` and return a preconditioner; see
        :mod:`openseespy_solvers.scipy.precond`.
    callback : callable, optional
        User-supplied function called after each iteration as ``callback()``.
    """ + _OPENSEES_LINEAR + _LINEAR_RETURNS + _LINEAR_NOTES + """
    See Also
    --------
    scipy.sparse.linalg.cg
    openseespy_solvers.scipy.precond.jacobi
    openseespy_solvers.scipy.precond.ilu

    Examples
    --------
    >>> from openseespy_solvers.scipy import cg, precond
    >>> solver = cg(rtol=1e-8, M=precond.jacobi)
    >>> solver.to_openseespy()["scheme"]
    'CSR'
    """
    return _CG(
        x0,
        rtol=rtol,
        atol=atol,
        maxiter=maxiter,
        M=M,
        callback=callback,
        scheme=scheme or "CSR",
        writable=writable,
        debug=debug,
        dtype=dtype,
    )

gmres ¶

gmres(x0: Any = None, *, rtol: float = 1e-05, atol: float = 0.0, restart: int | None = None, maxiter: int | None = None, M: Any = None, callback: Any = None, scheme: str | None = None, writable: str | list[str] = 'none', debug: bool = False, dtype: Any = float64) -> _GMRES

Source code in src/openseespy_solvers/scipy/__init__.py

def gmres(
    x0: Any = None,
    *,
    rtol: float = 1e-5,
    atol: float = 0.0,
    restart: int | None = None,
    maxiter: int | None = None,
    M: Any = None,
    callback: Any = None,
    scheme: str | None = None,
    writable: str | list[str] = "none",
    debug: bool = False,
    dtype: Any = np.float64,
) -> _GMRES:
    r"""Configure a GMRES solver for OpenSees ``PythonSparse``.

    Uses :func:`scipy.sparse.linalg.gmres` internally to solve ``Ax = b``.

    Parameters
    ----------
    x0 : ndarray, optional
        Starting guess for the solution.
    rtol, atol : float, optional
        Relative and absolute tolerances for convergence. Defaults are
        ``rtol=1e-5`` and ``atol=0.0``.
    restart : int, optional
        Number of iterations between restarts.
    maxiter : int, optional
        Maximum number of iterations.
    M : {sparse matrix, LinearOperator, callable}, optional
        Preconditioner for ``A``.
    callback : callable, optional
        User-supplied function called after each iteration.
    """ + _OPENSEES_LINEAR + _LINEAR_RETURNS + _LINEAR_NOTES + """
    See Also
    --------
    scipy.sparse.linalg.gmres
    cg

    Examples
    --------
    >>> from openseespy_solvers.scipy import gmres
    >>> solver = gmres(rtol=1e-8, restart=30)
    >>> solver.backend
    'scipy'
    """
    return _GMRES(
        x0,
        rtol=rtol,
        atol=atol,
        restart=restart,
        maxiter=maxiter,
        M=M,
        callback=callback,
        scheme=scheme or "CSR",
        writable=writable,
        debug=debug,
        dtype=dtype,
    )

eigsh ¶

eigsh(*, sigma: float | None = None, which: str = 'LM', v0: Any = None, ncv: int | None = None, maxiter: int | None = None, tol: float = 0.0, mode: str = 'normal', scheme: str | None = None, debug: bool = False, dtype: Any = float64) -> _Eigsh

Source code in src/openseespy_solvers/scipy/__init__.py

def eigsh(
    *,
    sigma: float | None = None,
    which: str = "LM",
    v0: Any = None,
    ncv: int | None = None,
    maxiter: int | None = None,
    tol: float = 0.0,
    mode: str = "normal",
    scheme: str | None = None,
    debug: bool = False,
    dtype: Any = np.float64,
) -> _Eigsh:
    r"""Configure a sparse eigenvalue solver for OpenSees ``PythonSparse``.

    Matches OpenSees ``eigsh``: shift-invert with ``sigma=0`` when
    ``find_smallest=True``; plain ARPACK when ``find_smallest=False`` unless
    ``sigma`` is set. ARPACK always uses ``which='LM'`` (OpenSees reference).

    Parameters
    ----------
    sigma : float, optional
        Shift for shift-invert. Default ``0.0`` when ``find_smallest=True``.
    which : {'LM', 'SM', ...}, optional
        Accepted for API compatibility; solves always use ``'LM'`` like OpenSees.
    v0 : ndarray, optional
        Starting vector for ARPACK.
    ncv : int, optional
        Number of Lanczos vectors.
    maxiter : int, optional
        Maximum number of iterations.
    tol : float, optional
        Convergence tolerance. Default is ``0.0``.
    mode : {'normal', 'buckling', 'cayley'}, optional
        Eigenproblem mode. Default is ``'normal'``.
    """ + _OPENSEES_EIGEN + _EIGEN_RETURNS + _EIGEN_NOTES + """
    See Also
    --------
    scipy.sparse.linalg.eigsh
    lobpcg

    Examples
    --------
    >>> from openseespy_solvers.scipy import eigsh
    >>> solver = eigsh(tol=1e-8)
    >>> solver.scheme
    'CSR'
    """
    return _Eigsh(
        sigma=sigma,
        which=which,
        v0=v0,
        ncv=ncv,
        maxiter=maxiter,
        tol=tol,
        mode=mode,
        scheme=scheme or "CSR",
        debug=debug,
        dtype=dtype,
    )

lobpcg ¶

lobpcg(*, X: Any = None, M: Any = None, tol: float | None = None, maxiter: int = 20, largest: bool = False, rng: Any = None, scheme: str | None = None, debug: bool = False, dtype: Any = float64) -> _Lobpcg

Source code in src/openseespy_solvers/scipy/__init__.py

def lobpcg(
    *,
    X: Any = None,
    M: Any = None,
    tol: float | None = None,
    maxiter: int = 20,
    largest: bool = False,
    rng: Any = None,
    scheme: str | None = None,
    debug: bool = False,
    dtype: Any = np.float64,
) -> _Lobpcg:
    r"""Configure a LOBPCG eigenvalue solver for OpenSees ``PythonSparse``.

    Uses :func:`scipy.sparse.linalg.lobpcg` internally to solve the generalized
    eigenproblem ``K x = \lambda M x``.

    Parameters
    ----------
    X : ndarray, optional
        Initial approximation to the eigenvectors. If ``None``, a random matrix
        with ``num_modes`` columns is generated using ``rng``.
    M : {sparse matrix, LinearOperator, callable}, optional
        Optional preconditioner passed to LOBPCG as ``M=``. A callable is
        invoked as ``M(K)``; see :mod:`openseespy_solvers.scipy.precond`.
    tol : float, optional
        Convergence tolerance.
    maxiter : int, optional
        Maximum number of iterations. Default is ``20``.
    largest : bool, optional
        If ``True``, compute largest eigenvalues. Default is ``False``.
    rng : {None, int, `numpy.random.Generator`}, optional
        Random number generator used when ``X`` is ``None``.
    """ + _OPENSEES_EIGEN + _EIGEN_RETURNS + _EIGEN_NOTES + """
    See Also
    --------
    scipy.sparse.linalg.lobpcg
    eigsh

    Examples
    --------
    >>> from openseespy_solvers.scipy import lobpcg
    >>> solver = lobpcg(tol=1e-8, maxiter=40)
    >>> solver.backend
    'scipy'
    """
    return _Lobpcg(
        X=X,
        M=M,
        tol=tol,
        maxiter=maxiter,
        largest=largest,
        rng=rng,
        scheme=scheme or "CSR",
        debug=debug,
        dtype=dtype,
    )