Enable multiple right-hand sides, fix C-division and n=3-edge case, made tests elaborate

.gitignore

@@ -83,7 +83,7 @@ celerybeat-schedule
 .env
 
 # virtualenv
-.venv
+.venv*
 venv/
 ENV/
 

CHANGELOG.md

@@ -2,6 +2,39 @@
 
 All notable changes to **pentapy** will be documented in this file.
 
+## [1.4.0] - 2024-06
+
+See [#26](https://github.com/GeoStat-Framework/pentapy/pull/26)
+
+### Enhancements
+
+- added support for multiple right-hand sides (currently serial)
+- improved error handling and added debug information to error messages
+
+### Changes
+
+- shotgun refactored and documented the Cython implementation of PTRANS-I and PTRANS-II for single and multiple right-hand sides support
+- fully typed the function ``pentapy.solve``
+- made internal solver alias handling of ``pentapy.solve`` smarter, more robust, and removed all duplicate code
+- gave all solvers a consistent interface
+- made code in ``pentapy.core`` more human-readable and maintainable and added comments
+- fixed typos in documentation
+
+### Bugfixes
+
+- fixed error handling in case of zero-division to trigger dead error handling branch (see [Issue 23](https://github.com/GeoStat-Framework/pentapy/issues/23))
+- fixed edge case error for row/column of 3 (see [Issue 24](https://github.com/GeoStat-Framework/pentapy/issues/24))
+
+### Tests
+
+- transitioned from ``unittest``-based testing to fully ``pytest``-based testing with parametrized and parallelized exhaustive testing (see [Issue 25](https://github.com/GeoStat-Framework/pentapy/issues/25))
+- made actual tests more meaningful by comparing them to LAPACK as reference standard (see [Issue 25](https://github.com/GeoStat-Framework/pentapy/issues/25))
+- included external solver bindings accessible via ``pentapy.solve`` as part of the test suite
+- increased true coverage (not line-hit coverage) close to 100%
+
+### Packaging
+
+- made dependency specification file-based and dynamic
 
 ## [1.3.0] - 2024-04
 
@@ -100,6 +133,7 @@ This is the first release of pentapy, a python toolbox for solving pentadiagonal
 The solver is implemented in cython, which makes it really fast.
 
 
+[1.4.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.3.0...v1.4.0
 [1.3.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.2.0...v1.3.0
 [1.2.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.1.2...v1.2.0
 [1.1.2]: https://github.com/GeoStat-Framework/pentapy/compare/v1.1.1...v1.1.2

examples/README.rst

@@ -88,7 +88,7 @@ If M is a full matrix, you call the following:
 
     X = pp.solve(M, Y)
 
-If M is flattend in row-wise order you have to set the keyword argument ``is_flat=True``:
+If M is flattened in row-wise order you have to set the keyword argument ``is_flat=True``:
 
 .. code-block:: python
 
@@ -99,7 +99,7 @@ If M is flattend in row-wise order you have to set the keyword argument ``is_fla
 
     X = pp.solve(M, Y, is_flat=True)
 
-If you got a col-wise flattend matrix you have to set ``index_row_wise=False``:
+If you got a col-wise flattened matrix you have to set ``index_row_wise=False``:
 
 .. code-block:: python
 

paper/paper.md

@@ -51,8 +51,8 @@ $$
 
 Here, $d_i$ are the diagonal entries and $d_i^{(j)}$ represent the $j$-th minor diagonal.
 
-Recently, @askar presented two algorithms to 
-solve the linear systems of equations for $X$, ``PTRANS-I`` and ``PTRANS-II``, 
+Recently, @askar presented two algorithms to
+solve the linear systems of equations for $X$, ``PTRANS-I`` and ``PTRANS-II``,
 applying first transformation to a triangular matrix and then, respectively, backward and forward substitution.
 ``pentapy`` provides Cython [@cython] implementations of these
 algorithms and a set of tools to convert matrices to row-wise or
@@ -73,7 +73,7 @@ The linear algebra solver of NumPy [@numpy] served as a standard reference, whic
 ``pentapy`` is designed to provide a fast solver for the special case of a
 pentadiagonal linear system. To the best of the author's knowledge,
 this package outperforms the current algorithms for solving pentadiagonal systems in Python.
-The solver can handle different input formats of the coefficient matrix, i.e., a flattend matrix or a
+The solver can handle different input formats of the coefficient matrix, i.e., a flattened matrix or a
 quadratic matrix.
 
 

pyproject.toml

@@ -14,7 +14,7 @@ name = "pentapy"
 authors = [{name = "Sebastian Müller", email = "[email protected]"}]
 readme = "README.md"
 license = {text = "MIT"}
-dynamic = ["version"]
+dynamic = ["version", "dependencies", "optional-dependencies"]
 description = "pentapy: A toolbox for pentadiagonal matrizes."
 classifiers = [
     "Development Status :: 5 - Production/Stable",
@@ -35,32 +35,10 @@ classifiers = [
     "Topic :: Scientific/Engineering",
     "Topic :: Utilities",
 ]
-dependencies = ["numpy>=1.20.0"]
 
-[project.optional-dependencies]
-scipy = ["scipy"]
-umfpack = ["scikit-umfpack"]
-all = [
-    "scipy",
-    "scikit-umfpack",
-]
-doc = [
-    "m2r2>=0.2.8",
-    "scipy>=1.1.0",
-    "matplotlib>=3",
-    "perfplot<0.9",
-    "numpydoc>=1.1",
-    "sphinx>=7",
-    "sphinx-gallery>=0.8",
-    "sphinx-rtd-theme>=2",
-]
-test = ["pytest-cov>=3"]
-check = [
-  "black>=24,<25",
-  "isort[colors]",
-  "pylint",
-  "cython-lint",
-]
+[tool.setuptools.dynamic]
+dependencies = {file = ["requirements/base.txt"]}
+optional-dependencies = {scipy = {file = ["requirements/scipy.txt"]}, umfpack = {file = ["requirements/umfpack.txt"]}, all = {file = ["requirements/all.txt"]}, doc = {file = ["requirements/doc.txt"]}, test = {file = ["requirements/test.txt"]}, check = {file = ["requirements/check.txt"]}}
 
 [project.urls]
 Homepage = "https://github.com/GeoStat-Framework/pentapy"
@@ -103,6 +81,8 @@ max-line-length = 120
         "*examples*",
         "*tests*",
         "*paper*",
+        "src/pentapy/_version.py",
+        "src/pentapy/__init__.py",
     ]
 
     [tool.coverage.report]

pytest.ini

@@ -0,0 +1,2 @@
+[pytest]
+addopts = --doctest-modules

requirements/all.txt

@@ -0,0 +1,2 @@
+scikit-umfpack
+scipy

requirements/base.txt

@@ -0,0 +1 @@
+numpy>=1.20.0

requirements/check.txt

@@ -0,0 +1,4 @@
+black>=24,<25
+isort[colors]
+pylint
+cython-lint

requirements/doc.txt

@@ -0,0 +1,8 @@
+m2r2>=0.2.8
+scipy>=1.1.0
+matplotlib>=3
+perfplot<0.9
+numpydoc>=1.1
+sphinx>=7
+sphinx-gallery>=0.8
+sphinx-rtd-theme>=2

requirements/scipy.txt

@@ -0,0 +1 @@
+scipy

requirements/test.txt

@@ -0,0 +1,4 @@
+pytest>=8
+pytest-cov>=3
+pytest-xdist>=3
+scipy>=1.1.0

requirements/umfpack

@@ -0,0 +1 @@
+scikit-umfpack

setup.py

@@ -1,23 +1,26 @@
-"""pentapy: A toolbox for pentadiagonal matrizes."""
+"""pentapy: A toolbox for pentadiagonal matrices."""
 
 import os
 
+import Cython.Compiler.Options
 import numpy as np
 from Cython.Build import cythonize
 from setuptools import Extension, setup
 
+Cython.Compiler.Options.annotate = True
+
 # cython extensions
 CY_MODULES = [
     Extension(
-        name=f"pentapy.solver",
+        name="pentapy.solver",
         sources=[os.path.join("src", "pentapy", "solver.pyx")],
         include_dirs=[np.get_include()],
         define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")],
     )
 ]
 
 setup(
-    ext_modules=cythonize(CY_MODULES),
+    ext_modules=cythonize(CY_MODULES, nthreads=1, annotate=True),
     package_data={"pentapy": ["*.pxd"]},  # include pxd files
     include_package_data=False,  # ignore other files
     zip_safe=False,

src/pentapy/.gitignore

@@ -0,0 +1,2 @@
+# Cython html files
+*.html

src/pentapy/__init__.py

@@ -44,7 +44,7 @@
 
 try:
     from pentapy._version import __version__
-except ImportError:  # pragma: nocover
+except ImportError:  # pragma: no cover
     # package is not installed
     __version__ = "0.0.0.dev0"
 

src/pentapy/_models.py

@@ -0,0 +1,49 @@
+"""
+Auxiliary models for the pentapy package.
+
+"""
+
+# === Imports ===
+
+from enum import IntEnum
+from typing import Dict
+
+# === Models ===
+
+
+class PentaSolverAliases(IntEnum):
+    """
+    Defines all available solver aliases for pentadiagonal systems, namely
+
+    - ``PTRANS_I``: The PTRANS-I algorithm
+    - ``PTRANS_II``: The PTRANS-II algorithm
+    - ``LAPACK``: Scipy's LAPACK solver :func:`scipy.linalg.solve_banded`
+    - ``SUPER_LU``: Scipy's SuperLU solver :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=False)`
+    - ``UMFPACK``: Scipy's UMFpack solver :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=True)`
+
+    """  # noqa: E501
+
+    PTRANS_I = 1
+    PTRANS_II = 2
+    LAPACK = 3
+    SUPER_LU = 4
+    UMFPACK = 5
+
+
+# === Constants ===
+
+_SOLVER_ALIAS_CONVERSIONS: Dict[str, PentaSolverAliases] = {
+    "1": PentaSolverAliases.PTRANS_I,
+    "ptrans-i": PentaSolverAliases.PTRANS_I,
+    "2": PentaSolverAliases.PTRANS_II,
+    "ptrans-ii": PentaSolverAliases.PTRANS_II,
+    "3": PentaSolverAliases.LAPACK,
+    "lapack": PentaSolverAliases.LAPACK,
+    "solve_banded": PentaSolverAliases.LAPACK,
+    "4": PentaSolverAliases.SUPER_LU,
+    "spsolve": PentaSolverAliases.SUPER_LU,
+    "5": PentaSolverAliases.UMFPACK,
+    "spsolve_umf": PentaSolverAliases.UMFPACK,
+    "umf": PentaSolverAliases.UMFPACK,
+    "umf_pack": PentaSolverAliases.UMFPACK,
+}