Add mumps solver

make use of python-mumps package

.github/workflows/python-package-conda.yml

@@ -35,7 +35,7 @@ jobs:
         conda info
         conda list
         conda config --show
-        conda install --quiet --yes -c conda-forge pip numpy scipy pytest pytest-cov pydiso
+        conda install --quiet --yes -c conda-forge pip numpy scipy pytest pytest-cov pydiso python-mumps
 
     - name: Install Our Package
       run: |

pymatsolver/__init__.py

@@ -34,7 +34,7 @@
   :toctree: generated/
 
   SolverCG
-  BicgJacobi
+  BiCGJacobi
 
 Direct Solvers
 ==============
@@ -45,15 +45,8 @@
   Solver
   SolverLU
   Pardiso
+  Mumps
 """
-from pymatsolver.solvers import Diagonal, Forward, Backward
-from pymatsolver.wrappers import WrapDirect
-from pymatsolver.wrappers import WrapIterative
-from pymatsolver.wrappers import Solver
-from pymatsolver.wrappers import SolverLU
-from pymatsolver.wrappers import SolverCG
-from pymatsolver.wrappers import SolverBiCG
-from pymatsolver.iterative import BicgJacobi
 
 SolverHelp = {}
 AvailableSolvers = {
@@ -66,8 +59,24 @@
     "Mumps": False
 }
 
+# Simple solvers
+from .solvers import Diagonal, Forward, Backward
+from .wrappers import WrapDirect
+from .wrappers import WrapIterative
+
+# Scipy Iterative solvers
+from .iterative import SolverCG
+from .iterative import SolverBiCG
+from .iterative import BiCGJacobi
+
+# Scipy direct solvers
+from .direct import Solver
+from .direct import SolverLU
+
+BicgJacobi = BiCGJacobi  # backwards compatibility
+
 try:
-    from pymatsolver.direct import Pardiso
+    from .direct import Pardiso
     AvailableSolvers['Pardiso'] = True
     PardisoSolver = Pardiso  # backwards compatibility
 except ImportError:
@@ -77,27 +86,17 @@
 to be installed through conda.
 """
 
-# try:
-#     from pymatsolver.mumps import Mumps
-#     AvailableSolvers['Mumps'] = True
-#     MumpsSolver = Mumps  # backwards compatibility
-# except Exception:
-#     SolverHelp['Mumps'] = """Mumps is not working.
-#
-# Ensure that you have Mumps installed, and know where the path to it is.
-# Try something like:
-#     cd pymatsolver/mumps
-#     make
-#
-# When that doesn't work, you may need to edit the make file, to modify the
-# path to the mumps installation, or any other compiler flags.
-# If you find a good way of doing it, please share.
-#
-# brew install mumps --with-scotch5 --without-mpi
-# mpicc --showme
-# """
+try:
+    from .direct import Mumps
+    AvailableSolvers['Mumps'] = True
+except ImportError:
+    SolverHelp['Mumps'] = """Mumps is not working.
+
+Ensure that you have python-mumps installed, which may also require Python
+to be installed through conda.
+"""
 
 __version__ = '0.2.0'
-__author__ = 'Rowan Cockett'
+__author__ = 'SimPEG Team'
 __license__ = 'MIT'
-__copyright__ = 'Copyright 2017 Rowan Cockett'
+__copyright__ = '2013 - 2024, SimPEG Team, https://simpeg.xyz'

pymatsolver/direct/__init__.py

@@ -0,0 +1,18 @@
+from ..wrappers import WrapDirect
+from scipy.sparse.linalg import spsolve, splu
+
+Solver = WrapDirect(spsolve, factorize=False, name="Solver")
+SolverLU = WrapDirect(splu, factorize=True, name="SolverLU")
+
+__all__ = ["Solver", "SolverLU"]
+try:
+    from .pardiso import Pardiso
+    __all__ += ["Pardiso"]
+except ImportError:
+    pass
+
+try:
+    from .mumps import Mumps
+    __all__ += ["Mumps"]
+except ImportError:
+    pass

pymatsolver/direct/mumps.py

@@ -0,0 +1,79 @@
+from pymatsolver.solvers import Base
+from mumps import Context
+
+class Mumps(Base):
+    """
+    Mumps solver
+    """
+
+    _factored = False
+    _transposed = False
+    ordering = ''
+
+    def __init__(self, A, context=None, **kwargs):
+        self.set_kwargs(**kwargs)
+        if context is None:
+            self.solver = Context()
+            self._analyzed = False
+            self._set_A(A)
+        else:
+            self.solver = context
+        self.A = A
+
+    def _set_A(self, A):
+        self.solver.set_matrix(
+            A,
+            symmetric=self.is_symmetric,
+            # positive_definite=self.is_positive_definite  # doesn't (yet) support setting positive definiteness
+        )
+
+    @property
+    def ordering(self):
+        return getattr(self, '_ordering', "metis")
+
+    @ordering.setter
+    def ordering(self, value):
+        self._ordering = value
+
+    @property
+    def _factored(self):
+        return self.solver.factored
+
+    @property
+    def transpose(self):
+        trans_obect = Mumps(
+            self.A,
+            self.solver,
+            is_symmetric=self.is_symmetric,
+            is_positive_definite=self.is_positive_definite,
+            ordering=self.ordering,
+        )
+        trans_obect._transposed = not self._transposed
+
+        return trans_obect
+
+    T = transpose
+
+    def factor(self, A=None):
+        reuse_analysis = False
+        if A is not None:
+            self._set_A(A)
+            self.A = A
+            # if it was previously factored then re-use the analysis.
+            reuse_analysis = self._factored
+        if not self._factored:
+            pivot_tol = 0.0 if self.is_positive_definite else 0.01
+            self.solver.factor(
+                ordering=self.ordering, reuse_analysis=reuse_analysis, pivot_tol=pivot_tol
+            )
+
+    def _solveM(self, rhs):
+        self.factor()
+        if self._transposed:
+            self.solver.mumps_instance.icntl[9] = 0
+        else:
+            self.solver.mumps_instance.icntl[9] = 1
+        sol = self.solver.solve(rhs)
+        return sol
+
+    _solve1 = _solveM

pymatsolver/direct/pardiso.py

@@ -0,0 +1,91 @@
+from pymatsolver.solvers import Base
+from pydiso.mkl_solver import MKLPardisoSolver
+from pydiso.mkl_solver import set_mkl_pardiso_threads, get_mkl_pardiso_max_threads
+
+class Pardiso(Base):
+    """
+    Pardiso Solver
+
+        https://github.com/simpeg/pydiso
+
+
+    documentation::
+
+        http://www.pardiso-project.org/
+    """
+
+    _factored = False
+
+    def __init__(self, A, **kwargs):
+        self.A = A
+        self.set_kwargs(**kwargs)
+        self.solver = MKLPardisoSolver(
+            self.A,
+            matrix_type=self._matrixType(),
+            factor=False
+        )
+
+    def _matrixType(self):
+        """
+            Set basic matrix type:
+
+            Real::
+
+                 1:  structurally symmetric
+                 2:  symmetric positive definite
+                -2:  symmetric indefinite
+                11:  nonsymmetric
+
+            Complex::
+
+                 6:  symmetric
+                 4:  hermitian positive definite
+                -4:  hermitian indefinite
+                 3:  structurally symmetric
+                13:  nonsymmetric
+
+        """
+        if self.is_real:
+            if self.is_symmetric:
+                if self.is_positive_definite:
+                    return 2
+                else:
+                    return -2
+            else:
+                return 11
+        else:
+            if self.is_symmetric:
+                return 6
+            elif self.is_hermitian:
+                if self.is_positive_definite:
+                    return 4
+                else:
+                    return -4
+            else:
+                return 13
+
+    def factor(self, A=None):
+        if A is not None:
+            self._factored = False
+            self.A = A
+        if not self._factored:
+            self.solver.refactor(self.A)
+            self._factored = True
+    def _solveM(self, rhs):
+        self.factor()
+        sol = self.solver.solve(rhs)
+        return sol
+
+    @property
+    def n_threads(self):
+        """
+        Number of threads to use for the Pardiso solver routine. This property
+        is global to all Pardiso solver objects for a single python process.
+        """
+        return get_mkl_pardiso_max_threads()
+
+    @n_threads.setter
+    def n_threads(self, n_threads):
+        set_mkl_pardiso_threads(n_threads)
+
+    _solve1 = _solveM