Merge pull request #43 from simpeg/mumps

Add mumps solver

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

@@ -12,32 +12,46 @@ on:
 
 jobs:
   build-and-test:
-    name: Testing (${{ matrix.python-version }}, ${{ matrix.os }})
+    name: Testing (${{ matrix.python-version }} on ${{ matrix.os }}) with ${{ matrix.solver }}.
     runs-on: ${{ matrix.os }}
     defaults:
       run:
         shell: bash -l {0}
     strategy:
+      fail-fast: false
       matrix:
         # NOTE: macOS-13 is the last Intel runner.
-        # When we move past that version we'll need to deal with Apple Silicon, likely using MUMPS.
-        os: [ubuntu-latest, windows-latest, macOS-13]
-        python-version: ["3.8", "3.9", "3.10", "3.11"]
+        os: [ubuntu-latest, windows-latest, macOS-13, macOS-latest]
+        solver: [mumps, pardiso]
+        python-version : ['3.10', '3.11', '3.12']
+        exclude:
+          - os: macOS-latest
+            solver: pardiso
 
     steps:
     - uses: actions/checkout@v4
     - name: Setup Conda
       uses: conda-incubator/setup-miniconda@v3
       with:
         auto-update-conda: true
-        activate-environment: dev
+        channels: conda-forge
         python-version: ${{ matrix.python-version }}
-    - name: Install Env
+    - name: Install Base Env
       run: |
         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 pip numpy scipy pytest pytest-cov
+
+    - name: Install MKL solver interface
+      if: ${{ matrix.solver == 'pardiso' }}
+      run:
+        conda install --quiet --yes "pydiso>=0.1"
+
+    - name: Install MUMPS solver interface
+      if: ${{ matrix.solver == 'mumps' }}
+      run:
+        conda install --quiet --yes python-mumps
 
     - name: Install Our Package
       run: |
@@ -49,15 +63,15 @@ jobs:
         pytest --cov-config=.coveragerc --cov-report=xml --cov=pymatsolver -s -v
 
     - name: Test Documentation
-      if: ${{ matrix.os == 'ubuntu-latest' }} and {{ matrix.python-version == '3.8' }}
+      if: ${{ (matrix.os == 'ubuntu-latest') && (matrix.python-version == '3.11') }}
       run: |
         pip install -r requirements_docs.txt
         cd docs
         make html
         cd ..
 
     - name: Upload coverage
-      if: ${{ matrix.os == 'ubuntu-latest' }} and {{ matrix.python-version == '3.8' }}
+      if: ${{ (matrix.os == 'ubuntu-latest') && (matrix.python-version  == '3.11') }}
       uses: codecov/codecov-action@v4
       with:
         verbose: true # optional (default = false)

.gitignore

@@ -8,3 +8,5 @@ pymatsolver.egg-info/
 *.dSYM
 docs/_build/*
 docs/generated
+
+.idea/

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,71 @@
+from pymatsolver.solvers import Base
+from mumps import Context
+
+class Mumps(Base):
+    """
+    Mumps solver
+    """
+    _transposed = False
+    ordering = ''
+
+    def __init__(self, A, **kwargs):
+        self.set_kwargs(**kwargs)
+        self.solver = Context()
+        self._set_A(A)
+        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_obj = Mumps.__new__(Mumps)
+        trans_obj.A = self.A
+        trans_obj.solver = self.solver
+        trans_obj.is_symmetric = self.is_symmetric
+        trans_obj.is_positive_definite = self.is_positive_definite
+        trans_obj.ordering = self.ordering
+        trans_obj._transposed = not self._transposed
+        return trans_obj
+
+    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