Document PythonSparse system and eigen solvers

src/PythonSparse.rst

@@ -0,0 +1,173 @@
+.. include:: sub.txt
+
+.. _systemPythonSparse:
+
+=======================
+ PythonSparse system
+=======================
+
+The **PythonSparse** system delegates linear system solves (:math:`\mathbf{A}\mathbf{x} = \mathbf{b}`) to user-defined Python objects. Sparse matrix buffers (CSR, CSC, or COO format) are exposed to Python via zero-copy memoryviews, enabling efficient integration with SciPy, CuPy, nvmath, or custom solvers without writing C++ wrappers or recompiling OpenSees. See the CuPy and SciPy examples below.
+
+.. function:: system('PythonSparse', config)
+   :noindex:
+
+   Create a linear system that forwards sparse matrix data to a Python solver object. ``config`` is a ``dict`` of solver options; see the table below.
+
+============================   ==================================================================================
+``'solver'``                   Python object with a ``solve(**kwargs)`` method (**required**)
+``'scheme'``                   Sparse storage format: ``'CSR'``, ``'CSC'``, or ``'COO'``. Default: ``'CSR'``.
+``'writable'``                 Buffers the solver may modify: ``'values'``, ``'rhs'``, ``'values,rhs'`` or ``'all'``, or ``'none'``.
+                               Also accepts a list: ``['values', 'rhs']``. Default: ``'none'`` — only the solution buffer ``x`` is writable.
+============================   ==================================================================================
+
+.. note::
+
+   Enabling ``writable`` for ``values`` or ``rhs`` allows in-place updates of the stiffness matrix and right-hand side vector; use only if you know what you are doing.
+
+The solver object must implement a ``solve(**kwargs)`` method. OpenSees passes matrix data as keyword arguments: memoryviews for arrays and scalars for metadata. Use `numpy.frombuffer <https://numpy.org/doc/stable/reference/generated/numpy.frombuffer.html>`_ to interpret memoryviews as NumPy arrays without copying. The keyword structure depends on ``scheme``:
+
+**CSR/CSC format** — compressed storage uses ``index_ptr`` and ``indices``:
+
+============================   =================   ============================================================
+``index_ptr``                  memoryview          Row pointers (CSR) or column pointers (CSC), int32
+``indices``                    memoryview          Column indices (CSR) or row indices (CSC), int32
+``values``                     memoryview          Matrix coefficients (float64)
+``rhs``                        memoryview          Right-hand side vector (float64)
+``x``                          memoryview          Solution buffer (float64, writable). Write the solution *in place*.
+``num_eqn``                    int                 Number of equations
+``nnz``                        int                 Number of nonzeros
+``matrix_status``              str                 ``'UNCHANGED'``, ``'COEFFICIENTS_CHANGED'``, or ``'STRUCTURE_CHANGED'``
+``storage_scheme``             str                 ``'CSR'``, ``'CSC'``, or ``'COO'``
+============================   =================   ============================================================
+
+**COO format** — coordinate storage uses ``row`` and ``col`` instead of ``index_ptr`` / ``indices``:
+
+============================   =================   ============================================================
+``row``                        memoryview          Row indices for each nonzero, int32 (COO only)
+``col``                        memoryview          Column indices for each nonzero, int32 (COO only)
+``values``                     memoryview          Matrix coefficients (float64)
+``rhs``                        memoryview          Right-hand side vector (float64)
+``x``                          memoryview          Solution buffer (float64, writable). Write the solution *in place*.
+``num_eqn``                    int                 Number of equations
+``nnz``                        int                 Number of nonzeros
+``matrix_status``              str                 ``'UNCHANGED'``, ``'COEFFICIENTS_CHANGED'``, or ``'STRUCTURE_CHANGED'``
+``storage_scheme``             str                 ``'COO'``
+============================   =================   ============================================================
+
+The ``solve`` method should return ``0`` on success, or a negative value to indicate failure.
+
+.. warning:: **In-place output required**
+
+   You **must** write the solution directly into the ``x`` buffer. OpenSees uses this buffer; it does **not** use return values or new arrays.
+
+   **Do:** Use in-place assignment, e.g. ``np.frombuffer(x, dtype=np.float64, count=num_eqn)[:] = result``
+
+   **Do not:** Return the solution, assign to a new variable, or write to a copy. The result will be ignored and the analysis will use uninitialized or stale data.
+
+.. _pythonSparseExampleCupy:
+
+Example — CuPy conjugate gradient (GPU)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This example uses `CuPy <https://cupy.dev/>`_ and `cupyx.scipy.sparse.linalg.cg <https://docs.cupy.dev/en/stable/reference/generated/cupyx.scipy.sparse.linalg.cg.html>`_ to solve the linear system on an NVIDIA GPU:
+
+.. code-block:: python
+
+   import numpy as np
+   import cupy as cp
+   import cupyx.scipy.sparse as cupyx_sparse
+   import cupyx.scipy.sparse.linalg as cupyx_splinalg
+   import openseespy.opensees as ops
+
+   class CuPyCGSolver:
+       def __init__(self, rtol=1e-5, atol=1e-12, maxiter=None):
+           self.rtol = rtol
+           self.atol = atol
+           self.maxiter = maxiter
+           self.A = None
+
+       def solve(self, **kwargs):
+           index_ptr = kwargs['index_ptr']
+           indices = kwargs['indices']
+           values = kwargs['values']
+           rhs = kwargs['rhs']
+           x = kwargs['x']
+           num_eqn = kwargs['num_eqn']
+           nnz = kwargs['nnz']
+           matrix_status = kwargs['matrix_status']
+
+           indptr = np.frombuffer(index_ptr, dtype=np.int32, count=num_eqn + 1)
+           idx = np.frombuffer(indices, dtype=np.int32, count=nnz)
+           vals = np.frombuffer(values, dtype=np.float64, count=nnz)
+
+           if matrix_status == 'STRUCTURE_CHANGED' or self.A is None:
+               self.A = cupyx_sparse.csr_matrix(
+                   (cp.asarray(vals), cp.asarray(idx), cp.asarray(indptr)),
+                   shape=(num_eqn, num_eqn)
+               )
+           elif matrix_status == 'COEFFICIENTS_CHANGED':
+               self.A.data[:] = cp.asarray(vals)
+
+           rhs_gpu = cp.asarray(np.frombuffer(rhs, dtype=np.float64, count=num_eqn))
+           x_gpu, info = cupyx_splinalg.cg(
+               self.A, rhs_gpu, tol=self.rtol, atol=self.atol, maxiter=self.maxiter
+           )
+
+           np.frombuffer(x, dtype=np.float64, count=num_eqn)[:] = cp.asnumpy(x_gpu)
+           return -int(info)
+
+   solver = CuPyCGSolver(rtol=1e-8, atol=1e-12, maxiter=1000)
+   ops.system('PythonSparse', {'solver': solver, 'scheme': 'CSR'})
+
+.. _pythonSparseExampleScipy:
+
+Example — SciPy direct solve (CPU)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This example uses `SciPy <https://scipy.org/>`_ `factorized <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.factorized.html>`_ (sparse LU) on CPU:
+
+.. code-block:: python
+
+   import numpy as np
+   import scipy.sparse as sp
+   import scipy.sparse.linalg as sp_linalg
+   import openseespy.opensees as ops
+
+   class SciPyFactorizedSolver:
+       def __init__(self):
+           self.A = None
+           self._solve_func = None
+
+       def solve(self, **kwargs):
+           index_ptr = kwargs['index_ptr']
+           indices = kwargs['indices']
+           values = kwargs['values']
+           rhs = kwargs['rhs']
+           x = kwargs['x']
+           num_eqn = kwargs['num_eqn']
+           nnz = kwargs['nnz']
+           matrix_status = kwargs['matrix_status']
+
+           indptr = np.frombuffer(index_ptr, dtype=np.int32, count=num_eqn + 1)
+           idx = np.frombuffer(indices, dtype=np.int32, count=nnz)
+           vals = np.frombuffer(values, dtype=np.float64, count=nnz)
+
+           if matrix_status != 'UNCHANGED' or self._solve_func is None:
+               self.A = sp.csr_matrix((vals.copy(), idx.copy(), indptr.copy()),
+                                     shape=(num_eqn, num_eqn))
+               self._solve_func = sp_linalg.factorized(self.A)
+
+           rhs_arr = np.frombuffer(rhs, dtype=np.float64, count=num_eqn)
+           x_arr = np.frombuffer(x, dtype=np.float64, count=num_eqn)
+           x_arr[:] = self._solve_func(rhs_arr)
+           return 0
+
+   solver = SciPyFactorizedSolver()
+   ops.system('PythonSparse', {'solver': solver, 'scheme': 'CSR'})
+
+Code developed by: `gaaraujo <https://github.com/gaaraujo>`_
+
+.. seealso::
+
+   * :ref:`eigenPythonSparse` — PythonSparse interface for generalized eigenvalue problems (``eigen('PythonSparse', ...)``)
+   * `EXAMPLES/SolverBenchmark/benchmark_python_sparse.py <https://github.com/OpenSees/OpenSees/blob/master/EXAMPLES/SolverBenchmark/benchmark_python_sparse.py>`_ — benchmark PythonSparse against native solvers

src/eigen.rst

@@ -10,11 +10,147 @@
 
    ================================   ===========================================================================
    numEigenvalues |int|               number of eigenvalues required
-   solver |str|                       optional string detailing type of solver: ``'-genBandArpack'``, ``'-fullGenLapack'``, (optional)
+   solver |str|                       optional string detailing type of solver: ``'-genBandArpack'``, ``'-fullGenLapack'``, ``'-symmBandLapack'``, or ``'PythonSparse'`` (optional)
+   ``config``                         for ``'PythonSparse'`` only: dictionary of solver options; see :ref:`eigenPythonSparse`
    ================================   ===========================================================================
 
+   .. code-block:: python
+
+      eigenvalues = ops.eigen(10)
+      eigenvalues = ops.eigen('-fullGenLapack', 10)
+      eigenvalues = ops.eigen('standard', '-symmBandLapack', 10)
+
 .. note::
 
    #. The eigenvectors are stored at the nodes and can be printed out using a Node Recorder, the nodeEigenvector command, or the Print command.
    #. The default eigensolver is able to solve only for N-1 eigenvalues, where N is the number of inertial DOFs. When running into this limitation the -fullGenLapack solver can be used instead of the default Arpack solver.
    #. The -fullGenLapack option is VERY SLOW for moderate to large models
+
+.. _eigenPythonSparse:
+
+PythonSparse eigen solver
+-------------------------
+
+.. index::
+   single: PythonEigen (see PythonSparse eigen)
+
+The **PythonSparse** eigen interface solves the generalized eigenvalue problem :math:`\mathbf{K}\phi = \lambda \mathbf{M}\phi` using a user-defined Python object. Sparse **K** and **M** buffers are exposed via zero-copy memoryviews so you can call SciPy, CuPy, or a custom eigensolver without recompiling OpenSees.
+
+Call ``ops.eigen('PythonSparse', numEigenvalues, config)``. The ``config`` dictionary accepts:
+
+============================   ==================================================================================
+``'solver'``                   Python object with a ``solve(**kwargs)`` method (**required**)
+``'scheme'``                   ``'CSR'``, ``'CSC'``, or ``'COO'``. Default: ``'CSR'``.
+============================   ==================================================================================
+
+.. note::
+
+   The PythonSparse eigen interface is available only when using the Python interpreter (OpenSeesPy).
+
+The solver object must implement ``solve(**kwargs)``. Keywords depend on ``scheme``:
+
+**CSR/CSC**
+
+============================   =================   ============================================================
+``index_ptr``                  memoryview          Row pointers (CSR) or column pointers (CSC), int32
+``indices``                    memoryview          Column indices (CSR) or row indices (CSC), int32
+``k_values``                   memoryview          Stiffness matrix **K** coefficients (float64)
+``m_values``                   memoryview          Mass matrix **M** coefficients (float64)
+``eigenvalues``                memoryview          Output buffer for eigenvalues (float64, writable)
+``eigenvectors``               memoryview          Output buffer for eigenvectors (float64, writable), row-major by mode
+``num_eqn``                    int                 Number of equations
+``nnz``                        int                 Number of nonzeros (K and M share sparsity pattern)
+``matrix_status``              str                 ``'UNCHANGED'``, ``'COEFFICIENTS_CHANGED'``, or ``'STRUCTURE_CHANGED'``
+``storage_scheme``             str                 ``'CSR'``, ``'CSC'``, or ``'COO'``
+``num_modes``                  int                 Number of modes requested
+``generalized``                bool                ``True``: generalized (K, M); ``False``: standard (K only)
+``find_smallest``              bool                ``True``: smallest eigenvalues; ``False``: largest
+============================   =================   ============================================================
+
+**COO** — uses ``row_indices`` and ``col_indices`` instead of ``index_ptr`` / ``indices`` (same ``k_values``, ``m_values``, outputs, and metadata as above).
+
+The ``solve`` method should return ``None`` on success, or raise an exception on failure.
+
+.. warning:: **In-place output required**
+
+   Write eigenvalues and eigenvectors into the ``eigenvalues`` and ``eigenvectors`` buffers. OpenSees does **not** use return values for the modes.
+
+   **Do:** e.g. ``np.frombuffer(eigenvalues, ...)[:] = eigvals`` and ``np.frombuffer(eigenvectors, ...)[:] = eigvecs.T.flatten()``
+
+   **Do not:** Return arrays only; results must be copied back into the buffers.
+
+
+.. _eigenPythonSparseExample:
+
+Example — SciPy ``eigsh``
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. code-block:: python
+
+   import numpy as np
+   import scipy.sparse as sp
+   import scipy.sparse.linalg as sp_linalg
+   import openseespy.opensees as ops
+
+   class SciPyGeneralizedEigenSolver:
+       def __init__(self, maxiter=None, tol=0.0):
+           self.maxiter = maxiter
+           self.tol = tol
+           self._k_matrix = None
+           self._m_matrix = None
+
+       def solve(self, **kwargs):
+           index_ptr = kwargs['index_ptr']
+           indices = kwargs['indices']
+           k_values = kwargs['k_values']
+           m_values = kwargs['m_values']
+           eigenvalues = kwargs['eigenvalues']
+           eigenvectors = kwargs['eigenvectors']
+           num_eqn = kwargs['num_eqn']
+           nnz = kwargs['nnz']
+           matrix_status = kwargs['matrix_status']
+           num_modes = kwargs['num_modes']
+           find_smallest = kwargs['find_smallest']
+
+           indptr = np.frombuffer(index_ptr, dtype=np.int32, count=num_eqn + 1)
+           idx = np.frombuffer(indices, dtype=np.int32, count=nnz)
+           k_data = np.frombuffer(k_values, dtype=np.float64, count=nnz)
+           m_data = np.frombuffer(m_values, dtype=np.float64, count=nnz)
+
+           if matrix_status == 'STRUCTURE_CHANGED' or self._k_matrix is None:
+               self._k_matrix = sp.csr_matrix((k_data, idx, indptr),
+                                             shape=(num_eqn, num_eqn))
+               self._m_matrix = sp.csr_matrix((m_data, idx, indptr),
+                                             shape=(num_eqn, num_eqn))
+           elif matrix_status == 'COEFFICIENTS_CHANGED':
+               self._k_matrix.data[:] = k_data
+               self._m_matrix.data[:] = m_data
+
+           eigsh_kwargs = {
+               'k': num_modes,
+               'M': self._m_matrix,
+               'which': 'LM',
+           }
+           if find_smallest:
+               eigsh_kwargs['sigma'] = 0.0
+           if self.maxiter is not None:
+               eigsh_kwargs['maxiter'] = self.maxiter
+           if self.tol > 0.0:
+               eigsh_kwargs['tol'] = self.tol
+
+           eigvals, eigvecs = sp_linalg.eigsh(self._k_matrix, **eigsh_kwargs)
+
+           np.frombuffer(eigenvalues, dtype=np.float64,
+                         count=num_modes)[:] = eigvals[:num_modes]
+           np.frombuffer(eigenvectors, dtype=np.float64,
+                         count=num_modes * num_eqn)[:] = eigvecs.T.flatten()
+
+           return None
+
+   solver = SciPyGeneralizedEigenSolver()
+   eigenvalues = ops.eigen('PythonSparse', 10, {'solver': solver, 'scheme': 'CSR'})
+
+.. seealso::
+
+   * :ref:`systemPythonSparse` — PythonSparse linear system (``system('PythonSparse', ...)``)
+   * `EXAMPLES/SolverBenchmark/benchmark_python_sparse_eigen.py <https://github.com/OpenSees/OpenSees/blob/master/EXAMPLES/SolverBenchmark/benchmark_python_sparse_eigen.py>`_

src/system.rst

@@ -25,6 +25,7 @@ The following contain information about available ``systemType``:
 #. :doc:`SparseSYM`
 #. :ref:`PFEM-System`
 #. :doc:`Mumps`
+#. :doc:`PythonSparse`
 
 
 .. toctree::
@@ -39,5 +40,6 @@ The following contain information about available ``systemType``:
    FullGeneral
    SparseSYM
    Mumps
+   PythonSparse