Skip to content

openseespy_solvers.cupy.precond

CUDA built-in preconditioners for iterative cupyx.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 on CUDA
ilu Incomplete LU preconditioner
direct Direct-solver preconditioner

Usage

from openseespy_solvers.cupy import cg
from openseespy_solvers.cupy import precond

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

Function Reference

Preconditioner factories for the cupy backend.

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

jacobi 

jacobi(A: Any) -> Any

Return a Jacobi (diagonal) preconditioner on GPU.

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

Parameters:

Name Type Description Default
A spmatrix

System matrix assembled by OpenSees on device.

 required

Returns:

Name Type Description
M LinearOperator

Diagonal preconditioner (y = inv * x on device) for M= in :func:~openseespy_solvers.cupy.cg.
See Also

direct ilu openseespy_solvers.cupy.cg openseespy_solvers.cupy.lobpcg openseespy_solvers.scipy.precond.jacobi

Examples:

>>> from openseespy_solvers.cupy import cg, lobpcg, precond
>>> cg(M=precond.jacobi)
>>> lobpcg(M=precond.jacobi)
Source code in src/openseespy_solvers/cupy/precond.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def jacobi(A: Any) -> Any:
    """Return a Jacobi (diagonal) preconditioner on GPU.

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

    Parameters
    ----------
    A : cupyx.scipy.sparse.spmatrix
        System matrix assembled by OpenSees on device.

    Returns
    -------
    M : cupyx.scipy.sparse.linalg.LinearOperator
        Diagonal preconditioner (``y = inv * x`` on device) for ``M=`` in
        :func:`~openseespy_solvers.cupy.cg`.

    See Also
    --------
    direct
    ilu
    openseespy_solvers.cupy.cg
    openseespy_solvers.cupy.lobpcg
    openseespy_solvers.scipy.precond.jacobi

    Examples
    --------
    >>> from openseespy_solvers.cupy import cg, lobpcg, precond
    >>> cg(M=precond.jacobi)
    >>> lobpcg(M=precond.jacobi)
    """
    cp, _csp, cspla = _import_cupy()
    diag = A.diagonal()
    inv = cp.ones_like(diag)
    mask = diag != 0
    inv[mask] = 1.0 / diag[mask]

    dtype = A.dtype
    n = A.shape[0]

    def matvec(x: Any) -> Any:
        arr = cp.asarray(x, dtype=dtype)
        if arr.ndim == 1:
            return inv * arr
        return inv[:, None] * arr

    def matmat(X: Any) -> Any:
        return inv[:, None] * cp.asarray(X, dtype=dtype)

    return _linear_operator_with_matmat(cspla, n, dtype, matvec, matmat)

ilu 

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

Return an incomplete LU preconditioner as a LinearOperator on GPU.

By default fill_factor=1, which triggers cupy's GPU ILU path with no fill-in and no pivoting. See :func:cupyx.scipy.sparse.linalg.spilu — only fill_factor=1 runs the factorization on GPU; other settings delegate to scipy.sparse.linalg on CPU.

Parameters:

Name Type Description Default
A spmatrix

System matrix assembled by OpenSees on device.

 required
**opts Any

Keyword arguments forwarded to :func:cupyx.scipy.sparse.linalg.spilu. fill_factor defaults to 1 when omitted.

 {}

Returns:

Name Type Description
M LinearOperator

Preconditioner implementing M @ x via solve on the factorization.
See Also

cupyx.scipy.sparse.linalg.spilu jacobi direct openseespy_solvers.scipy.precond.ilu

Notes

With fill_factor=1, factorization and application stay on the GPU but the pattern is restricted to the sparsity structure of A (ILU(0)-like, no fill-in).

Examples:

>>> from openseespy_solvers.cupy import cg, precond
>>> solver = cg(M=precond.ilu)
Source code in src/openseespy_solvers/cupy/precond.py 
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 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
def ilu(A: Any, **opts: Any) -> Any:
    """Return an incomplete LU preconditioner as a ``LinearOperator`` on GPU.

    By default ``fill_factor=1``, which triggers cupy's GPU ILU path with no
    fill-in and no pivoting. See
    :func:`cupyx.scipy.sparse.linalg.spilu` — only ``fill_factor=1`` runs the
    factorization on GPU; other settings delegate to `scipy.sparse.linalg` on CPU.

    Parameters
    ----------
    A : cupyx.scipy.sparse.spmatrix
        System matrix assembled by OpenSees on device.
    **opts
        Keyword arguments forwarded to
        :func:`cupyx.scipy.sparse.linalg.spilu`. ``fill_factor`` defaults to
        ``1`` when omitted.

    Returns
    -------
    M : cupyx.scipy.sparse.linalg.LinearOperator
        Preconditioner implementing ``M @ x`` via ``solve`` on the factorization.

    See Also
    --------
    cupyx.scipy.sparse.linalg.spilu
    jacobi
    direct
    openseespy_solvers.scipy.precond.ilu

    Notes
    -----
    With ``fill_factor=1``, factorization and application stay on the GPU but
    the pattern is restricted to the sparsity structure of ``A`` (ILU(0)-like,
    no fill-in).

    Examples
    --------
    >>> from openseespy_solvers.cupy import cg, precond
    >>> solver = cg(M=precond.ilu)
    """
    cp, _csp, cspla = _import_cupy()
    kwargs = dict(opts)
    kwargs.setdefault("fill_factor", 1)
    ilu_factor = cspla.spilu(A.tocsc(), **kwargs)

    dtype = A.dtype
    n = A.shape[0]

    def matvec(x: Any) -> Any:
        arr = cp.asarray(x, dtype=dtype)
        if arr.ndim == 1:
            return ilu_factor.solve(arr)
        return cp.column_stack([ilu_factor.solve(arr[:, j]) for j in range(arr.shape[1])])

    def matmat(X: Any) -> Any:
        X = cp.asarray(X, dtype=dtype)
        return cp.column_stack([ilu_factor.solve(X[:, j]) for j in range(X.shape[1])])

    return _linear_operator_with_matmat(cspla, n, dtype, matvec, matmat)

direct 

direct(solver: LinearSolver)

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

Intended primarily for :func:~openseespy_solvers.cupy.lobpcg as M= on stiffness K. 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.nvmath.direct_solver or :func:~openseespy_solvers.cupy.spsolve.

 required

Returns:

Name Type Description
preconditioner callable

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

jacobi ilu openseespy_solvers.nvmath.direct_solver openseespy_solvers.cupy.lobpcg openseespy_solvers.hybrid

Examples:

>>> from openseespy_solvers.cupy import lobpcg, precond, spsolve
>>> eigsolver = lobpcg(M=precond.direct(spsolve()), tol=1e-8)
Source code in src/openseespy_solvers/cupy/precond.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def direct(solver: LinearSolver):
    """Return a callable ``M(A)`` that approximates ``A^{-1}`` via a direct solver.

    Intended primarily for :func:`~openseespy_solvers.cupy.lobpcg` as ``M=`` on
    stiffness ``K``. 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.nvmath.direct_solver` or
        :func:`~openseespy_solvers.cupy.spsolve`.

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

    See Also
    --------
    jacobi
    ilu
    openseespy_solvers.nvmath.direct_solver
    openseespy_solvers.cupy.lobpcg
    openseespy_solvers.hybrid

    Examples
    --------
    >>> from openseespy_solvers.cupy 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.")

    cp, _csp, cspla = _import_cupy()

    def preconditioner(A: Any) -> Any:
        needs_refactor = True
        dtype = A.dtype
        n = A.shape[0]

        def _apply_column(rhs: Any) -> Any:
            nonlocal needs_refactor
            result = apply_inner_factorization(
                solver,
                A,
                rhs,
                refactor=needs_refactor,
                on_device=False,
            )
            needs_refactor = False
            return cp.asarray(result, dtype=dtype)

        def matvec(x: Any) -> Any:
            arr = cp.asarray(x, dtype=dtype)
            if arr.ndim == 1:
                return _apply_column(arr)
            return cp.column_stack([_apply_column(arr[:, j]) for j in range(arr.shape[1])])

        def matmat(X: Any) -> Any:
            X = cp.asarray(X, dtype=dtype)
            return cp.column_stack([_apply_column(X[:, j]) for j in range(X.shape[1])])

        return _linear_operator_with_matmat(cspla, n, dtype, matvec, matmat)

    return preconditioner