diff --git a/src/PythonSparse.rst b/src/PythonSparse.rst new file mode 100644 index 0000000..9f56d80 --- /dev/null +++ b/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 `_ 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 `_ and `cupyx.scipy.sparse.linalg.cg `_ 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 `_ `factorized `_ (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 `_ + +.. seealso:: + + * :ref:`eigenPythonSparse` — PythonSparse interface for generalized eigenvalue problems (``eigen('PythonSparse', ...)``) + * `EXAMPLES/SolverBenchmark/benchmark_python_sparse.py `_ — benchmark PythonSparse against native solvers diff --git a/src/eigen.rst b/src/eigen.rst index 250c8ea..26a6939 100644 --- a/src/eigen.rst +++ b/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 `_ \ No newline at end of file diff --git a/src/system.rst b/src/system.rst index cc7061f..dcf4bf6 100644 --- a/src/system.rst +++ b/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