Skip to content

openseespy_solvers.scipy.precond

CPU built-in preconditioners for iterative scipy.sparse.linalg solvers.

Pass these callables through the M= keyword of cg, gmres, or lobpcg. When M is callable, the solver calls it as M(A) after OpenSeesPy supplies the assembled matrix.

Preconditioners

Preconditioner Description
jacobi Diagonal/Jacobi preconditioner
ilu Incomplete LU preconditioner
direct Direct-solver preconditioner

Usage

from openseespy_solvers.scipy import cg
from openseespy_solvers.scipy import precond

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

Function Reference

Preconditioner factories for the scipy backend.

These callables are intended for the M argument of :func:cg, :func:gmres, and :func:lobpcg. Each factory accepts the assembled sparse matrix A from OpenSees and returns a preconditioner object.

jacobi 

jacobi(A: spmatrix) -> spmatrix

Return a Jacobi (diagonal) preconditioner.

Computes M = diag(1 / diag(A)). Zero diagonal entries are left as 1.

Parameters:

Name Type Description Default
A sparse matrix

System matrix assembled by OpenSees.

 required

Returns:

Name Type Description
M spmatrix

Diagonal preconditioner suitable for M= in :func:~openseespy_solvers.scipy.cg.
See Also

ilu direct openseespy_solvers.scipy.cg

Examples:

>>> import scipy.sparse as sp
>>> from openseespy_solvers.scipy import cg, precond
>>> A = sp.diags([2.0, 3.0, 4.0])
>>> M = precond.jacobi(A)
>>> solver = cg(M=precond.jacobi)
Source code in src/openseespy_solvers/scipy/precond.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def jacobi(A: sp.spmatrix) -> sp.spmatrix:
    """Return a Jacobi (diagonal) preconditioner.

    Computes ``M = diag(1 / diag(A))``. Zero diagonal entries are left as ``1``.

    Parameters
    ----------
    A : sparse matrix
        System matrix assembled by OpenSees.

    Returns
    -------
    M : scipy.sparse.spmatrix
        Diagonal preconditioner suitable for ``M=`` in :func:`~openseespy_solvers.scipy.cg`.

    See Also
    --------
    ilu
    direct
    openseespy_solvers.scipy.cg

    Examples
    --------
    >>> import scipy.sparse as sp
    >>> from openseespy_solvers.scipy import cg, precond
    >>> A = sp.diags([2.0, 3.0, 4.0])
    >>> M = precond.jacobi(A)
    >>> solver = cg(M=precond.jacobi)
    """
    diag = A.diagonal()
    inv = np.ones_like(diag, dtype=np.float64)
    nz = diag != 0.0
    inv[nz] = 1.0 / diag[nz]
    return sp.diags(inv)

ilu 

ilu(A: spmatrix, **opts: Any) -> LinearOperator

Return an incomplete LU preconditioner as a LinearOperator.

Parameters:

Name Type Description Default
A sparse matrix

System matrix assembled by OpenSees.

 required
**opts Any

Additional keyword arguments forwarded to :func:scipy.sparse.linalg.spilu (for example drop_tol, fill_factor).

 {}

Returns:

Name Type Description
M LinearOperator

Preconditioner implementing M @ x via forward/back substitution on the ILU factorization.
See Also

scipy.sparse.linalg.spilu jacobi direct

Examples:

>>> from openseespy_solvers.scipy import cg, precond
>>> solver = cg(M=lambda A: precond.ilu(A, drop_tol=1e-4))
Source code in src/openseespy_solvers/scipy/precond.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def ilu(A: sp.spmatrix, **opts: Any) -> spla.LinearOperator:
    """Return an incomplete LU preconditioner as a ``LinearOperator``.

    Parameters
    ----------
    A : sparse matrix
        System matrix assembled by OpenSees.
    **opts
        Additional keyword arguments forwarded to
        :func:`scipy.sparse.linalg.spilu` (for example ``drop_tol``,
        ``fill_factor``).

    Returns
    -------
    M : scipy.sparse.linalg.LinearOperator
        Preconditioner implementing ``M @ x`` via forward/back substitution on
        the ILU factorization.

    See Also
    --------
    scipy.sparse.linalg.spilu
    jacobi
    direct

    Examples
    --------
    >>> from openseespy_solvers.scipy import cg, precond
    >>> solver = cg(M=lambda A: precond.ilu(A, drop_tol=1e-4))
    """
    ilu_factor = spla.spilu(A.tocsc(), **opts)

    def matvec(x: np.ndarray) -> np.ndarray:
        return ilu_factor.solve(x)

    n = A.shape[0]
    return spla.LinearOperator((n, n), matvec=matvec, dtype=A.dtype)

direct 

direct(solver: LinearSolver)

Return a callable M(A) that approximates A^{-1} via a direct solver.

Intended primarily for :func:~openseespy_solvers.scipy.lobpcg as M=. The first application factors A; later applications in the same eigen solve reuse the factorization when the matrix is unchanged.

Parameters:

Name Type Description Default
solver LinearSolver

Direct solver for A x = b, for example :func:~openseespy_solvers.scipy.spsolve or :func:~openseespy_solvers.scipy.umfpack.

 required

Returns:

Name Type Description
preconditioner callable

Function M(A) returning a LinearOperator suitable for M=.
See Also

jacobi ilu openseespy_solvers.scipy.lobpcg openseespy_solvers.hybrid

Examples:

>>> from openseespy_solvers.scipy import lobpcg, precond, spsolve
>>> eigsolver = lobpcg(M=precond.direct(spsolve()), tol=1e-8)
Source code in src/openseespy_solvers/scipy/precond.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
def direct(solver: LinearSolver):
    """Return a callable ``M(A)`` that approximates ``A^{-1}`` via a direct solver.

    Intended primarily for :func:`~openseespy_solvers.scipy.lobpcg` as ``M=``.
    The first application factors ``A``; later applications in the same eigen
    solve reuse the factorization when the matrix is unchanged.

    Parameters
    ----------
    solver : LinearSolver
        Direct solver for ``A x = b``, for example :func:`~openseespy_solvers.scipy.spsolve`
        or :func:`~openseespy_solvers.scipy.umfpack`.

    Returns
    -------
    preconditioner : callable
        Function ``M(A)`` returning a ``LinearOperator`` suitable for ``M=``.

    See Also
    --------
    jacobi
    ilu
    openseespy_solvers.scipy.lobpcg
    openseespy_solvers.hybrid

    Examples
    --------
    >>> from openseespy_solvers.scipy import lobpcg, precond, spsolve
    >>> eigsolver = lobpcg(M=precond.direct(spsolve()), tol=1e-8)
    """
    if not isinstance(solver, LinearSolver):
        raise TypeError("precond.direct() requires a LinearSolver instance.")

    def preconditioner(A: sp.spmatrix) -> spla.LinearOperator:
        needs_refactor = True
        n = A.shape[0]

        def matvec(x: np.ndarray) -> np.ndarray:
            nonlocal needs_refactor
            result = apply_inner_factorization(
                solver,
                A,
                x,
                refactor=needs_refactor,
                on_device=False,
            )
            needs_refactor = False
            return np.asarray(result, dtype=A.dtype).ravel()

        return spla.LinearOperator((n, n), matvec=matvec, dtype=A.dtype)

    return preconditioner