diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 3fa8607..709fbb1 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -12,17 +12,21 @@ 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 @@ -30,14 +34,24 @@ jobs: 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,7 +63,7 @@ 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 @@ -57,7 +71,7 @@ jobs: 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) diff --git a/.gitignore b/.gitignore index c2e9891..4b397e2 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,5 @@ pymatsolver.egg-info/ *.dSYM docs/_build/* docs/generated + +.idea/ diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index ead34b0..a144b9f 100644 --- a/pymatsolver/__init__.py +++ b/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' diff --git a/pymatsolver/direct.py b/pymatsolver/direct.py deleted file mode 100644 index 2262325..0000000 --- a/pymatsolver/direct.py +++ /dev/null @@ -1,95 +0,0 @@ -from pymatsolver.solvers import Base -try: - 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 -except ImportError: - pass diff --git a/pymatsolver/direct/__init__.py b/pymatsolver/direct/__init__.py new file mode 100644 index 0000000..498b9a0 --- /dev/null +++ b/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 diff --git a/pymatsolver/direct/mumps.py b/pymatsolver/direct/mumps.py new file mode 100644 index 0000000..645ef8b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pymatsolver/direct/pardiso.py b/pymatsolver/direct/pardiso.py new file mode 100644 index 0000000..43076c6 --- /dev/null +++ b/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 diff --git a/pymatsolver/iterative.py b/pymatsolver/iterative.py index d16cf9e..a815f85 100644 --- a/pymatsolver/iterative.py +++ b/pymatsolver/iterative.py @@ -1,22 +1,26 @@ import numpy as np import scipy import scipy.sparse as sp -from scipy.sparse.linalg import aslinearoperator, bicgstab +from scipy.sparse.linalg import bicgstab, cg, aslinearoperator from packaging.version import Version -from pymatsolver.solvers import Base +from .wrappers import WrapIterative +from .solvers import Base # The tol kwarg was removed from bicgstab in scipy 1.14.0. # See https://docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.sparse.linalg.bicgstab.html -TOL_ARG_NAME = "rtol" if Version(scipy.__version__) >= Version("1.14.0") else "tol" +RTOL_ARG_NAME = "rtol" if Version(scipy.__version__) >= Version("1.14.0") else "tol" +SolverCG = WrapIterative(cg, name="SolverCG") +SolverBiCG = WrapIterative(bicgstab, name="SolverBiCG") -class BicgJacobi(Base): +class BiCGJacobi(Base): """Bicg Solver with Jacobi preconditioner""" _factored = False solver = None maxiter = 1000 - tol = 1E-6 + rtol = 1E-6 + atol = 0.0 def __init__(self, A, symmetric=True): self.A = A @@ -32,16 +36,18 @@ def factor(self): self.M = aslinearoperator(Ainv) self._factored = True + @property + def _tols(self): + return {RTOL_ARG_NAME: self.rtol, 'atol': self.atol} + + def _solve1(self, rhs): self.factor() sol, info = self.solver( - self.A, - rhs, - **{ - TOL_ARG_NAME: self.tol, - "maxiter": self.maxiter, - "M": self.M, - } + self.A, rhs, + maxiter=self.maxiter, + M=self.M, + **self._tols, ) return sol @@ -49,17 +55,8 @@ def _solveM(self, rhs): self.factor() sol = [] for icol in range(rhs.shape[1]): - sol.append( - self.solver( - self.A, - rhs[:, icol].flatten(), - **{ - TOL_ARG_NAME: self.tol, - "maxiter": self.maxiter, - "M": self.M, - } - )[0] - ) + sol.append(self.solver(self.A, rhs[:, icol].flatten(), + maxiter=self.maxiter, M=self.M, **self._tols,)[0]) out = np.hstack(sol) out.shape return out @@ -68,3 +65,4 @@ def clean(self): self.M = None self.A = None self._factored = False + diff --git a/pymatsolver/mumps/Makefile b/pymatsolver/mumps/Makefile deleted file mode 100644 index e6501e4..0000000 --- a/pymatsolver/mumps/Makefile +++ /dev/null @@ -1,46 +0,0 @@ -INCMUMPS = -I/usr/include -L/usr/lib \ - -L/usr/local/opt/libevent/lib -INCMAC = -L/usr/local/opt/scotch5/lib \ - -L/usr/local/Cellar/open-mpi/2.0.1/lib \ - -L/usr/local/opt/libevent/lib \ - -L/usr/local/opt/mumps/libexec/lib \ - -I/usr/local/opt/mumps/libexec/include \ - -I/usr/local/opt/scotch5/include \ - -I/usr/local/Cellar/open-mpi/2.0.1/include - -# brew install mumps --with-scotch5 --without-mpi -# Helpful: mpicc --showme - -LIB = -lmpi -lblas -lpord -lcmumps -LIBSEQ = -lmpiseq_seq -lzmumps_seq -ldmumps_seq -lsmumps_seq -lmumps_common_seq -LIBPAR = -lmpiseq -lzmumps -ldmumps -lsmumps -lmumps_common -LIBSCOTCH = -lesmumps -lptscotch -lptscotcherrexit -lscotch -lscotcherrexit -lptesmumps -lptscotcherr -lptscotchparmetis -lscotcherr -lscotchmetis - -all: clean build build_mac - -build_mac: - gfortran -c $(INCMAC) mumps_p.f90 -fPIC - gfortran -c $(INCMAC) mumps_cmplx_p.f90 -fPIC - f2py -c mumps_interface.f90 -m MumpsInterface \ - $(INCMAC) \ - --f90flags='-fcray-pointer' \ - $(LIB) $(LIBPAR) $(LIBSCOTCH)\ - mumps_p.o mumps_cmplx_p.o - rm -f *.o *.mod - -build: - gfortran -c $(INCMUMPS) mumps_p.f90 -fPIC - gfortran -c $(INCMUMPS) mumps_cmplx_p.f90 -fPIC - f2py -c mumps_interface.f90 -m MumpsInterface \ - $(INCMUMPS) \ - --f90flags='-fcray-pointer' \ - $(LIB) $(LIBSEQ)\ - mumps_p.o mumps_cmplx_p.o - rm -f *.o *.mod - -clean: - rm -f *.o *.mod *.so - rm -rf *.dSYM - -test: - cd ..;nosetests diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py deleted file mode 100644 index 6826e25..0000000 --- a/pymatsolver/mumps/_init_.py +++ /dev/null @@ -1,165 +0,0 @@ -from __future__ import unicode_literals -from __future__ import print_function -from __future__ import division -from __future__ import absolute_import - -import gc -from pymatsolver.solvers import Base -from . import MumpsInterface as _MUMPSINT - -_mumpsErrors = { - -1: "An error occurred on processor INFO(2).", - -2: "NZ is out of range. INFO(2)=NZ.", - -3: "MUMPS was called with an invalid value for JOB. This may happen for example if the analysis (JOB=1) was not performed before the factorization (JOB=2), or the factorization was not performed before the solve (JOB=3), or the initialization phase (JOB=-1) was performed a second time on an instance not freed (JOB=-2). See description of JOB in Section 3. This error also occurs if JOB does not contain the same value on all processes on entry to MUMPS.", - -4: "Error in user-provided permutation array PERM IN at position INFO(2). This error may only occuron the host.", - -5: "Problem of REAL or COMPLEX workspace allocation of size INFO(2) during analysis.", - -6: "Matrix is singular in structure. INFO(2) holds the structural rank.", - -7: "Problem of INTEGER workspace allocation of size INFO(2) during analysis.", - -8: "Main internal integer workarray IS too small for factorization. This may happen, for example, if numerical pivoting leads to significantly more fill-in than was predicted by the analysis. The user should increase the value of ICNTL(14) before calling the factorization again (JOB=2).", - -9: "Main internal real/complex workarray S too small. If INFO(2) is positive, then the number of entries that are missing in S at the moment when the error is raised is available in INFO(2). If INFO(2) is negative, then its absolute value should be multiplied by 1 million. If an error -9 occurs, the user should increase the value of ICNTL(14) before calling the factorization (JOB=2) again, except if ICNTL(23) is provided, in which case ICNTL(23) should be increased.", - -10: "Numerically singular matrix.", - -11: "Internal real/complex workarray S too small for solution. Please contact us. If INFO(2) is positive, then the number of entries that are missing in S at the moment when the error is raised is available in INFO(2).", - -12: "Internal real/complex workarray S too small for iterative refinement. Please contact us.", - -13: "An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.", - -14: "Internal integer workarray IS too small for solution. See error INFO(1) = -8.", - -15: "Integer workarray IS too small for iterative refinement and/or error analysis. See error INFO(1) = -8.", - -16: "N is out of range. INFO(2)=N.", - -17: "The internal send buffer that was allocated dynamically by MUMPS on the processor is too small. The user should increase the value of ICNTL(14) before calling MUMPS again.", - -20: "The internal reception buffer that was allocated dynamically by MUMPS is too small. INFO(2) holds the minimum size of the reception buffer required (in bytes). The user should increase the value of ICNTL(14) before calling MUMPS again.", - -21: "Value of PAR=0 is not allowed because only one processor is available; Running node mode (the host is not a slave processor itself) requires at least two processors. The user should either set PAR to 1 or increase the number of processors.", - -22: "A pointer array is provided by the user that is either not associated, or has insufficient size, or is associated and should not be associated (for example, RHS on non-host processors).", - -23: "MPI was not initialized by the user prior to a call to MUMPS with JOB = -1.", - -24: "NELT is out of range. INFO(2)=NELT.", - -25: "A problem has occurred in the initialization of the BLACS. This may be because you are using a vendor's BLACS. Try using a BLACS version from netlib instead.", - -26: "LRHS is out of range. INFO(2)=LRHS.", - -27: "NZ RHS and IRHS PTR(NRHS+1) do not match. INFO(2) = IRHS PTR(NRHS+1).", - -28: "IRHS PTR(1) is not equal to 1. INFO(2) = IRHS PTR(1).", - -29: "LSOL loc is smaller than INFO(23). INFO(2)=LSOL loc.", - -30: "SCHUR LLD is out of range. INFO(2) = SCHUR LLD.", - -31: "A 2D block cyclic symmetric (SYM=1 or 2) Schur complement is required with the option ICNTL(19)=3, but the user has provided a process grid that does not satisfy the constraint MBLOCK=NBLOCK. INFO(2)=MBLOCK-NBLOCK.", - -32: "Incompatible values of NRHS and ICNTL(25). Either ICNTL(25) was set to -1 and NRHS is different from INFOG(28); or ICNTL(25) was set to i, 1 <= i <= INFOG(28) and NRHS is different from 1. Value of NRHS is stored in INFO(2).", - -33: "ICNTL(26) was asked for during solve phase but the Schur complement was not asked for at the analysis phase (ICNTL(19)). INFO(2)=ICNTL(26).", - -34: "LREDRHS is out of range. INFO(2)=LREDRHS.", - -35: "This error is raised when the expansion phase is called (ICNTL(26) = 2) but reduction phase (ICNTL(26)=1) was not called before. INFO(2) contains the value of ICNTL(26).", - -36: "Incompatible values of ICNTL(25) and INFOG(28). The value of ICNTL(25) is stored in INFO(2).", - -37: "Value of ICNTL(25) incompatible with some other parameter. with SYM or ICNTL(xx). If INFO(2)=0 then ICNTL(25) is incompatible with SYM: in current version, the null space basis functionality is not available for unsymmetric matrices (SYM=0). Otherwise, ICNTL(25) is incompatible with ICNTL(xx), and the index xx is stored in INFO(2).", - -38: "Parallel analysis was set (i.e., ICNTL(28)=2) but PT-SCOTCH or ParMetis were not provided.", - -39: "Incompatible values for ICNTL(28) and ICNTL(5) and/or ICNTL(19) and/or ICNTL(6). Parallel analysis is not possible in the cases where the matrix is unassembled and/or a Schur complement is requested and/or a maximum transversal is requested on the matrix.", - -40: "The matrix was indicated to be positive definite (SYM=1) by the user but a negative or null pivot was encountered during the processing of the root by ScaLAPACK. SYM=2 should be used.", - -44: "The solve phase (JOB=3) cannot be performed because the factors or part of the factors are not available. INFO(2) contains the value of ICNTL(31).", - -45: "NRHS <= 0. INFO(2) contains the value of NRHS.", - -46: "NZ RHS <= 0. This is currently not allowed in case of reduced right-hand-side (ICNTL(26)=1) and in case entries of A-1 are requested (ICNTL(30)=1). INFO(2) contains the value of NZ RHS.", - -47: "Entries of A-1 were requested during the solve phase (JOB=3, ICNTL(30)=1) but the constraint NRHS=N is not respected. The value of NRHS is provided in INFO(2).", - -48: "A-1 Incompatible values of ICNTL(30) and ICNTL(xx). xx is stored in INFO(2).", - -49: "SIZE SCHUR has an incorrect value (SIZE SCHUR < 0 or SIZE SCHUR >=N, or SIZE SCHUR was modified on the host since the analysis phase. The value of SIZE SCHUR is provided in INFO(2).", - -90: "Error in out-of-core management. See the error message returned on output unit ICNTL(1) for more information.", - +1: "Index (in IRN or JCN) out of range. Action taken by subroutine is to ignore any such entries and continue. INFO(2) is set to the number of faulty entries. Details of the first ten are printed on unit ICNTL(2).", - +2: "During error analysis the max-norm of the computed solution was found to be zero.", - +4: "User data JCN has been modified (internally) by the solver.", - +8: "Warning return from the iterative refinement routine. More than ICNTL(10) iterations are required.", -} - - -class _Pointer(object): - """Gets an int and a destroy call that gets called on garbage collection. - - There can be multiple Solvers around the place that are pointing to the same factor in memory. - - This ensures that when all references are removed, there is automatic garbage collection. - - You can always clean the Solver, and it will explicitly call the destroy method. - """ - def __init__(self, pointerINT, destroyCall): - self.INT = pointerINT - self.destroyCall = destroyCall - - def __del__(self): - self.destroyCall(self.INT) - -class Mumps(Base): - """ - - documentation:: - - http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf - - """ - - transpose = False - symmetric = False - - @property - def T(self): - newMS = self.__class__( - self.A, - symmetric=self.symmetric, - fromPointer=self.pointer - ) - newMS.transpose = not self.transpose - return newMS - - def __init__(self, A, symmetric=False, fromPointer=None): - self.A = A.tocsc() - self.symmetric = symmetric - - if fromPointer is None: - self.factor() - elif isinstance(fromPointer, _Pointer): - self.pointer = fromPointer - else: - raise Exception('Unknown pointer for construction.') - - @property - def isfactored(self): - return getattr(self, 'pointer', None) is not None - - def _funhandle(self, ftype): - """ - switches the function handle between real and complex. - - ftype in ['F','S','D'] - - means factor, solve, destroy - """ - if self.A.dtype == float: - return {'F': _MUMPSINT.factor_mumps, - 'S': _MUMPSINT.solve_mumps, - 'D': _MUMPSINT.destroy_mumps}[ftype] - elif self.A.dtype == complex: - return {'F': _MUMPSINT.factor_mumps_cmplx, - 'S': _MUMPSINT.solve_mumps_cmplx, - 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] - - def factor(self): - if self.isfactored: - return - - sym = 1 if self.symmetric else 0 - ierr, p = self._funhandle('F')( - sym, - self.A.data, - self.A.indices+1, - self.A.indptr+1 - ) - if ierr < 0: - raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) - elif ierr > 0: - print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) - - self.pointer = _Pointer(p, self._funhandle('D')) - - def _solveM(self, rhs): - self.factor() - rhs = rhs.flatten(order='F') - n = self.A.shape[0] - nrhs = rhs.size // n - T = 1 if self.transpose else 0 - sol = self._funhandle('S')(self.pointer.INT, nrhs, rhs, T) - return sol - - _solve1 = _solveM - - def clean(self): - self._funhandle('D')(self.pointer.INT) - del self.pointer - gc.collect() diff --git a/pymatsolver/mumps/mumps_cmplx_p.f90 b/pymatsolver/mumps/mumps_cmplx_p.f90 deleted file mode 100644 index 4fb056d..0000000 --- a/pymatsolver/mumps/mumps_cmplx_p.f90 +++ /dev/null @@ -1,272 +0,0 @@ - - - -module mumps_cmplx_mod - - save - ! integer,parameter:: ui_out = 10 ! unit for mumps log file - !integer,parameter:: st_unit=11 ! unit for output statistics file - - - INCLUDE 'zmumps_struc.h' - - !logical,private:: mumps_initialized = .false. - - -contains - -!--------------------------------------------------------------- - -subroutine init( mumps_par, sym ) - -implicit none - -TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. - - !mumps_par%COMM = MPI_COMM_WORLD - mumps_par%SYM = sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. - mumps_par%PAR = 1 - - mumps_par%JOB = -1 ! initialize - CALL ZMUMPS(mumps_par) - - - !if (mumps_par%MYID == 0) then - ! open(unit=ui_out, file='mumps.log', action='write') - !end if - - !mumps_initialized = .true. - -return -end subroutine init - -!--------------------------------------------------------------- - -subroutine convert_to_mumps_format( mumps_par, n, A,jA,iA, ierr ) -! A is actually a transpose. -implicit none - -TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: n ! # of rows in A -complex(kind=8),intent(in):: A(*) -integer(kind=8),intent(in):: jA(*), iA(n+1) -integer(kind=8),intent(out):: ierr - -integer nonz, i,j, j1,j2, ind, jcol, istat - - -nonz = iA(n+1) - 1 -!n = size(iA) - 1 - -mumps_par%N = n - -if (mumps_par%SYM == 0) then ! not symmetric matrix - mumps_par%NZ = nonz - allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) - if ( istat /= 0 ) then - ierr = -13 ! allocation error - return - end if - - - do i = 1, n - mumps_par%JCN( iA(i) : iA(i+1)-1 ) = i - end do ! i - - mumps_par%IRN = jA(1:nonz) - mumps_par%A = A(1:nonz) - -else ! symmetric matrix - ! Keep only the lower half of the matrix (row >= column). - - nonz = nonz/2 + n ! should be +n/2, but I'm using n just in case. - allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) - if ( istat /= 0 ) then - ierr = -13 ! allocation error - return - end if - - ind = 0 - do i = 1, n - - j1 = iA(i) - j2 = iA(i+1) - 1 - do j = j1, j2 - jcol = jA(j) - - if (i >= jcol) then - ind = ind + 1 - mumps_par%A(ind) = A(j) - mumps_par%JCN(ind) = jcol - mumps_par%IRN(ind) = i - end if - - end do ! j - end do ! i - - mumps_par%NZ = ind - !if (nonz < ind) call errorstop('nonz < ind') ! debug -end if - -ierr = 0 -return -end subroutine convert_to_mumps_format - -!--------------------------------------------------------------- - -subroutine factor_matrix( mumps_par, ierr ) - -implicit none - -TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(out):: ierr - -!mumps_par%icntl(2) = 6 ! output stream for diagnostics -mumps_par%icntl(4) = 0 ! 1 ! amount of output - - -mumps_par%icntl(2) = 0 ! ui_out ! output stream for diagnostic printing -mumps_par%icntl(3) = 0 ! ui_out ! output stream for global information - -!mumps_par%icntl(14) = 40 ! % increase in working space - -!mumps_par%icntl(11) = 0 ! return statistics -! mumps_par%icntl(7) = 5 ! ordering type -!mumps_par%cntl(1) = 0.d0 ! rel. threshold for numerical pivoting - - -mumps_factorization: do - - mumps_par%JOB = 1 ! analysis - CALL ZMUMPS(mumps_par) - - mumps_par%JOB = 2 ! factorization - CALL ZMUMPS(mumps_par) - - ierr = mumps_par%INFOG(1) - - - if ( mumps_par%INFOG(1) == -9 .or. mumps_par%INFOG(1) == -8 ) then - ! Main internal real/complex workarray S too small. - mumps_par%icntl(14) = mumps_par%icntl(14) + 10 - if ( mumps_par%MYID == 0 ) then - write(*,30) mumps_par%icntl(14) - 30 format(/'MUMPS percentage increase in the estimated working space', & - /'increased to',i4) - end if - - !else if (mumps_par%INFOG(1) == -13) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,40) - ! 40 format(/'MUMPS memory allocation error.'/) - ! end if - ! stop - - !else if (mumps_par%INFOG(1) == -10) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,45) - ! 45 format(/'MUMPS ERROR: Numerically singular matrix.'/) - ! end if - ! stop - - !else if (mumps_par%INFOG(1) == -40) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,46) - ! 46 format(/'MUMPS ERROR: matrix is not positive definite.'/) - ! end if - ! stop - - !else if ( mumps_par%INFOG(1) < 0 ) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,20) mumps_par%INFOG(1), mumps_par%INFOG(2) - ! 20 format(/'ERROR occured in MUMPS!',/,'INFOG(1), INFOG(2) ', 2i6,/) - ! end if - ! stop - - else - exit mumps_factorization ! factorization successful - end if - -end do mumps_factorization - - -if ( mumps_par%MYID == 0 ) then - !flush(ui_out) - ! Turn off mumps output. - mumps_par%icntl(2) = 0 ! output stream for diagnostic printing - mumps_par%icntl(3) = 0 ! output stream for global information -end if - - -return -end subroutine factor_matrix - -!-------------------------------------------------------------- - -subroutine solve( mumps_par, nrhs, rhs, x, transpose ) -! Solve A*x = rhs - -implicit none - -TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: nrhs ! # of right-hand-sides -complex(kind=8),intent(in):: rhs(nrhs * mumps_par%N) -complex(kind=8),intent(out),target:: x(nrhs * mumps_par%N) ! solution -logical,intent(in):: transpose ! if .true. take the transpose - - x = rhs - - mumps_par%RHS => x - - ! The following is significant only on the host cpu. - mumps_par%NRHS = nrhs ! # of right-hand-sides - mumps_par%LRHS = mumps_par%N ! size of system - - if (transpose) then - mumps_par%icntl(9) = 0 ! for solving A'x = b - else - mumps_par%icntl(9) = 1 ! for solving Ax = b - end if - - if (transpose) then - mumps_par%RHS = conjg(mumps_par%RHS) - end if - - mumps_par%JOB = 3 ! Solve the system. - CALL ZMUMPS(mumps_par) - ! At this point mumps_par%RHS (rhs) contains the solution. - - if (transpose) then - mumps_par%RHS = conjg(mumps_par%RHS) - end if - -return -end subroutine solve - -!-------------------------------------------------------------- - -subroutine destroy(mumps_par) - -implicit none -TYPE (ZMUMPS_STRUC),intent(inout):: mumps_par - -!if (.not. mumps_initialized) return - -! Destroy the instance (deallocate internal data structures) -mumps_par%JOB = -2 -CALL ZMUMPS(mumps_par) - -if (associated(mumps_par%A)) then - deallocate(mumps_par%IRN, mumps_par%JCN, mumps_par%A) -end if - -!if (mumps_par%MYID == 0) then -! close(ui_out) -!end if - -return -end subroutine destroy - - -end module mumps_cmplx_mod diff --git a/pymatsolver/mumps/mumps_interface.f90 b/pymatsolver/mumps/mumps_interface.f90 deleted file mode 100644 index 943124e..0000000 --- a/pymatsolver/mumps/mumps_interface.f90 +++ /dev/null @@ -1,192 +0,0 @@ - -!------------------------------------------------- -!------------- REAL ------------ -!------------------------------------------------- - -subroutine factor_mumps( n, nnz, sym, A, jA, iA, ierr, pm_out ) - - -use mumps_mod, only: init, convert_to_mumps_format, factor_matrix, destroy -implicit none - -INCLUDE 'dmumps_struc.h' - -integer(kind=8),intent(out):: pm_out ! mumps pointer -integer(kind=8),intent(in):: n ! # of rows in A -integer(kind=8),intent(in):: nnz ! # of non zeros in A -integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. -real(kind=8),intent(in):: A(nnz) -integer(kind=8),intent(in):: jA(nnz), iA(n+1) - -integer(kind=8),intent(out):: ierr -! integer,intent(out):: ierr ! =0 no error, < 0 error (memory, singular, etc.) - -TYPE(DMUMPS_STRUC),pointer:: pmumps_par -TYPE(DMUMPS_STRUC):: mumps_par - -pointer ( pm, mumps_par ) - -allocate(pmumps_par) -pm = loc(pmumps_par) - -call init( mumps_par, sym) -call convert_to_mumps_format(mumps_par, n, A,jA,iA, ierr ) -if (ierr < 0) return ! allocation error - -call factor_matrix(mumps_par, ierr) - -if ( ierr < 0 ) then - ! Error in factorization. - call destroy(mumps_par) - pm_out = 0 -else - pm_out = pm -end if - -return -end subroutine factor_mumps - -!------------------------------------------------- - -subroutine solve_mumps( pm_in, nrhs, lrhs, rhs, x, transpose ) -! Solve A*x = rhs - -use mumps_mod, only: solve - -implicit none - -INCLUDE 'dmumps_struc.h' - -integer(kind=8),intent(in):: pm_in ! mumps pointer -integer(kind=8),intent(in):: nrhs ! # of right-hand-sides -integer(kind=8),intent(in):: lrhs ! # length of right-hand-sides -real(kind=8),intent(in):: rhs(lrhs) ! right-hand-side -real(kind=8),intent(out):: x(lrhs) ! solution -integer(kind=8),intent(in):: transpose ! =1 for transpose - -TYPE(DMUMPS_STRUC):: mumps_par -pointer ( pm, mumps_par ) -pm = pm_in - -call solve(mumps_par, nrhs, rhs, x, (transpose==1) ) - -return -end subroutine solve_mumps - -!------------------------------------------------- - -subroutine destroy_mumps( pm_in ) -! Destroy the instance (deallocate internal data structures) - -!DIR$ ATTRIBUTES DLLEXPORT :: destroy_mumps -!DIR$ ATTRIBUTES ALIAS: 'destroy_mumps_':: destroy_mumps -use mumps_mod, only: destroy - -implicit none -INCLUDE 'dmumps_struc.h' - -integer(kind=8),intent(in):: pm_in ! mumps pointer -TYPE(DMUMPS_STRUC):: mumps_par -pointer ( pm, mumps_par ) -pm = pm_in - -call destroy(mumps_par) - -return -end subroutine destroy_mumps - -!------------------------------------------------- -!------------- COMPLEX ------------ -!------------------------------------------------- - -subroutine factor_mumps_cmplx( n, nnz, sym, A, jA, iA, ierr, pm_out ) - -use mumps_cmplx_mod, only: init, convert_to_mumps_format, factor_matrix, destroy -implicit none - -INCLUDE 'zmumps_struc.h' - -integer(kind=8),intent(out):: pm_out ! mumps pointer -integer(kind=8),intent(in):: n ! # of rows in A -integer(kind=8),intent(in):: nnz ! # of non zeros in A -integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. -complex(kind=8),intent(in):: A(nnz) -integer(kind=8),intent(in):: jA(nnz), iA(n+1) - -integer(kind=8),intent(out):: ierr -! integer,intent(out):: ierr ! =0 no error, < 0 error (memory, singular, etc.) - -TYPE(ZMUMPS_STRUC),pointer:: pmumps_par -TYPE(ZMUMPS_STRUC):: mumps_par - -pointer ( pm, mumps_par ) - -allocate(pmumps_par) -pm = loc(pmumps_par) - -call init( mumps_par, sym) -call convert_to_mumps_format(mumps_par, n, A,jA,iA, ierr ) -if (ierr < 0) return ! allocation error - -call factor_matrix(mumps_par, ierr) - -if ( ierr < 0 ) then - ! Error in factorization. - call destroy(mumps_par) - pm_out = 0 -else - pm_out = pm -end if - -return -end subroutine factor_mumps_cmplx - -!------------------------------------------------- - -subroutine solve_mumps_cmplx( pm_in, nrhs, lrhs, rhs, x, transpose ) -! Solve A*x = rhs - -use mumps_cmplx_mod, only: solve - -implicit none - -INCLUDE 'zmumps_struc.h' - -integer(kind=8),intent(in):: pm_in ! mumps pointer -integer(kind=8),intent(in):: nrhs ! # of right-hand-sides -integer(kind=8),intent(in):: lrhs ! # length of right-hand-sides -complex(kind=8),intent(in):: rhs(lrhs) ! right-hand-side -complex(kind=8),intent(out):: x(lrhs) ! solution -integer(kind=8),intent(in):: transpose ! =1 for transpose - -TYPE(ZMUMPS_STRUC):: mumps_par -pointer ( pm, mumps_par ) -pm = pm_in - -call solve(mumps_par, nrhs, rhs, x, (transpose==1) ) - -return -end subroutine solve_mumps_cmplx - -!------------------------------------------------- - -subroutine destroy_mumps_cmplx( pm_in ) -! Destroy the instance (deallocate internal data structures) - -use mumps_cmplx_mod, only: destroy - -implicit none -INCLUDE 'zmumps_struc.h' - -integer(kind=8),intent(in):: pm_in ! mumps pointer -TYPE(ZMUMPS_STRUC):: mumps_par -pointer ( pm, mumps_par ) -pm = pm_in - -call destroy(mumps_par) - -return -end subroutine destroy_mumps_cmplx - -!------------------------------------------------- - diff --git a/pymatsolver/mumps/mumps_p.f90 b/pymatsolver/mumps/mumps_p.f90 deleted file mode 100755 index 710e0c8..0000000 --- a/pymatsolver/mumps/mumps_p.f90 +++ /dev/null @@ -1,266 +0,0 @@ - - - -module mumps_mod - - save - ! integer,parameter:: ui_out = 10 ! unit for mumps log file - !integer,parameter:: st_unit=11 ! unit for output statistics file - - - INCLUDE 'dmumps_struc.h' - - !logical,private:: mumps_initialized = .false. - - -contains - -!--------------------------------------------------------------- - -subroutine init( mumps_par, sym ) - -implicit none - -TYPE (DMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. - - !mumps_par%COMM = MPI_COMM_WORLD - mumps_par%SYM = sym ! =0 unsymmetric, =1 symm. pos def, =2 general symm. - mumps_par%PAR = 1 - - mumps_par%JOB = -1 ! initialize - CALL DMUMPS(mumps_par) - - - !if (mumps_par%MYID == 0) then - ! open(unit=ui_out, file='mumps.log', action='write') - !end if - - !mumps_initialized = .true. - -return -end subroutine init - -!--------------------------------------------------------------- - -subroutine convert_to_mumps_format( mumps_par, n, A,jA,iA, ierr ) -! A is actually a transpose. -implicit none - -TYPE (DMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: n ! # of rows in A -real(kind=8),intent(in):: A(*) -integer(kind=8),intent(in):: jA(*), iA(n+1) -integer(kind=8),intent(out):: ierr - -integer nonz, i,j, j1,j2, ind, jcol, istat - - -nonz = iA(n+1) - 1 -!n = size(iA) - 1 - -mumps_par%N = n - -if (mumps_par%SYM == 0) then ! not symmetric matrix - mumps_par%NZ = nonz - allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) - if ( istat /= 0 ) then - ierr = -13 ! allocation error - return - end if - - - do i = 1, n - mumps_par%JCN( iA(i) : iA(i+1)-1 ) = i - end do ! i - - mumps_par%IRN = jA(1:nonz) - mumps_par%A = A(1:nonz) - -else ! symmetric matrix - ! Keep only the lower half of the matrix (row >= column). - - nonz = nonz/2 + n ! should be +n/2, but I'm using n just in case. - allocate( mumps_par%IRN(nonz), mumps_par%JCN(nonz), mumps_par%A(nonz), stat=istat) - if ( istat /= 0 ) then - ierr = -13 ! allocation error - return - end if - - ind = 0 - do i = 1, n - - j1 = iA(i) - j2 = iA(i+1) - 1 - do j = j1, j2 - jcol = jA(j) - - if (i >= jcol) then - ind = ind + 1 - mumps_par%A(ind) = A(j) - mumps_par%JCN(ind) = jcol - mumps_par%IRN(ind) = i - end if - - end do ! j - end do ! i - - mumps_par%NZ = ind - !if (nonz < ind) call errorstop('nonz < ind') ! debug -end if - -ierr = 0 -return -end subroutine convert_to_mumps_format - -!--------------------------------------------------------------- - -subroutine factor_matrix( mumps_par, ierr ) - -implicit none - -TYPE (DMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(out):: ierr - -!mumps_par%icntl(2) = 6 ! output stream for diagnostics -mumps_par%icntl(4) = 0 ! 1 ! amount of output - - -mumps_par%icntl(2) = 0 ! ui_out ! output stream for diagnostic printing -mumps_par%icntl(3) = 0 ! ui_out ! output stream for global information - -!mumps_par%icntl(14) = 40 ! % increase in working space - -!mumps_par%icntl(11) = 0 ! return statistics -! mumps_par%icntl(7) = 5 ! ordering type -!mumps_par%cntl(1) = 0.d0 ! rel. threshold for numerical pivoting - - -mumps_factorization: do - - mumps_par%JOB = 1 ! analysis - CALL DMUMPS(mumps_par) - - mumps_par%JOB = 2 ! factorization - CALL DMUMPS(mumps_par) - - ierr = mumps_par%INFOG(1) - - - if ( mumps_par%INFOG(1) == -9 .or. mumps_par%INFOG(1) == -8 ) then - ! Main internal real/complex workarray S too small. - mumps_par%icntl(14) = mumps_par%icntl(14) + 10 - if ( mumps_par%MYID == 0 ) then - write(*,30) mumps_par%icntl(14) - 30 format(/'MUMPS percentage increase in the estimated working space', & - /'increased to',i4) - end if - - !else if (mumps_par%INFOG(1) == -13) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,40) - ! 40 format(/'MUMPS memory allocation error.'/) - ! end if - ! stop - - !else if (mumps_par%INFOG(1) == -10) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,45) - ! 45 format(/'MUMPS ERROR: Numerically singular matrix.'/) - ! end if - ! stop - - !else if (mumps_par%INFOG(1) == -40) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,46) - ! 46 format(/'MUMPS ERROR: matrix is not positive definite.'/) - ! end if - ! stop - - !else if ( mumps_par%INFOG(1) < 0 ) then - ! if ( mumps_par%MYID == 0 ) then - ! write(*,20) mumps_par%INFOG(1), mumps_par%INFOG(2) - ! 20 format(/'ERROR occured in MUMPS!',/,'INFOG(1), INFOG(2) ', 2i6,/) - ! end if - ! stop - - else - exit mumps_factorization ! factorization successful - end if - -end do mumps_factorization - - -if ( mumps_par%MYID == 0 ) then - !flush(ui_out) - ! Turn off mumps output. - mumps_par%icntl(2) = 0 ! output stream for diagnostic printing - mumps_par%icntl(3) = 0 ! output stream for global information -end if - - -return -end subroutine factor_matrix - -!-------------------------------------------------------------- - -subroutine solve( mumps_par, nrhs, rhs, x, transpose ) -! Solve A*x = rhs - -implicit none - -TYPE (DMUMPS_STRUC),intent(inout):: mumps_par -integer(kind=8),intent(in):: nrhs ! # of right-hand-sides -real(kind=8),intent(in):: rhs(nrhs * mumps_par%N) -real(kind=8),intent(out),target:: x(nrhs * mumps_par%N) ! solution -logical,intent(in):: transpose ! if .true. take the transpose - - x = rhs - - mumps_par%RHS => x - - ! The following is significant only on the host cpu. - mumps_par%NRHS = nrhs ! # of right-hand-sides - mumps_par%LRHS = mumps_par%N ! size of system - - if (transpose) then - mumps_par%icntl(9) = 0 ! for solving A'x = b - else - mumps_par%icntl(9) = 1 ! for solving Ax = b - end if - - - mumps_par%JOB = 3 ! Solve the system. - CALL DMUMPS(mumps_par) - ! At this point mumps_par%RHS (rhs) contains the solution. - - -return -end subroutine solve - -!-------------------------------------------------------------- - -subroutine destroy(mumps_par) - -implicit none -TYPE (DMUMPS_STRUC),intent(inout):: mumps_par - -!if (.not. mumps_initialized) return - -! Destroy the instance (deallocate internal data structures) -mumps_par%JOB = -2 -CALL DMUMPS(mumps_par) - -if (associated(mumps_par%A)) then - deallocate(mumps_par%IRN, mumps_par%JCN, mumps_par%A) -end if - -!if (mumps_par%MYID == 0) then -! close(ui_out) -!end if - -return -end subroutine destroy - - -end module mumps_mod diff --git a/pymatsolver/wrappers.py b/pymatsolver/wrappers.py index e70a068..2d0ca46 100644 --- a/pymatsolver/wrappers.py +++ b/pymatsolver/wrappers.py @@ -112,7 +112,3 @@ def _solveM(self, rhs): } ) -Solver = WrapDirect(linalg.spsolve, factorize=False, name="Solver") -SolverLU = WrapDirect(linalg.splu, factorize=True, name="SolverLU") -SolverCG = WrapIterative(linalg.cg, name="SolverCG") -SolverBiCG = WrapIterative(linalg.bicgstab, name="SolverBiCG") diff --git a/setup.py b/setup.py index da46988..de98ca1 100644 --- a/setup.py +++ b/setup.py @@ -5,8 +5,7 @@ """ -from distutils.core import setup -from setuptools import find_packages +from setuptools import setup, find_packages CLASSIFIERS = [ 'Development Status :: 4 - Beta', @@ -30,19 +29,19 @@ setup( name="pymatsolver", version="0.2.0", - packages=find_packages(exclude=["*mumps", "tests"]), + packages=find_packages(exclude=["tests"]), install_requires=[ 'numpy>=1.7', 'scipy>=1.8', ], - author="Rowan Cockett", + author="SimPEG Developers", author_email="rowanc1@gmail.com", description="pymatsolver: Matrix Solvers for Python", long_description=LONG_DESCRIPTION, license="MIT", keywords="matrix solver", url="http://simpeg.xyz/", - download_url="http://github.com/rowanc1/pymatsolver", + download_url="http://github.com/simpeg/pymatsolver", classifiers=CLASSIFIERS, platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], use_2to3=False diff --git a/tests/test_BicgJacobi.py b/tests/test_BicgJacobi.py index 24e41a5..29b04ea 100644 --- a/tests/test_BicgJacobi.py +++ b/tests/test_BicgJacobi.py @@ -46,7 +46,7 @@ def test_T(self): Ainv.clean() -class TestPardisoComplex(unittest.TestCase): +class TestBicgJacobiComplex(unittest.TestCase): def setUp(self): nSize = 100 diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 1bf4afc..26d68b9 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -5,82 +5,79 @@ TOL = 1e-11 -""" -class TestMumps(unittest.TestCase): +if pymatsolver.AvailableSolvers['Mumps']: + class TestMumps(unittest.TestCase): - if pymatsolver.AvailableSolvers['Mumps']: + def setUp(self): + n = 5 + irn = np.r_[0, 1, 3, 4, 1, 0, 4, 2, 1, 2, 0, 2] + jcn = np.r_[1, 2, 2, 4, 0, 0, 1, 3, 4, 1, 2, 2] + a = np.r_[ + 3.0, -3.0, 2.0, 1.0, 3.0, 2.0, + 4.0, 2.0, 6.0, -1.0, 4.0, 1.0 + ] + rhs = np.r_[20.0, 24.0, 9.0, 6.0, 13.0] + rhs = np.c_[rhs, 10*rhs, 100*rhs] + sol = np.r_[1., 2., 3., 4., 5.] + sol = np.c_[sol, 10*sol, 100*sol] + A = sp.coo_matrix((a, (irn, jcn)), shape=(n, n)) + self.A = A + self.rhs = rhs + self.sol = sol - def setUp(self): - n = 5 - irn = np.r_[0, 1, 3, 4, 1, 0, 4, 2, 1, 2, 0, 2] - jcn = np.r_[1, 2, 2, 4, 0, 0, 1, 3, 4, 1, 2, 2] - a = np.r_[ - 3.0, -3.0, 2.0, 1.0, 3.0, 2.0, - 4.0, 2.0, 6.0, -1.0, 4.0, 1.0 - ] - rhs = np.r_[20.0, 24.0, 9.0, 6.0, 13.0] - rhs = np.c_[rhs, 10*rhs, 100*rhs] - sol = np.r_[1., 2., 3., 4., 5.] - sol = np.c_[sol, 10*sol, 100*sol] - A = sp.coo_matrix((a, (irn, jcn)), shape=(n, n)) - self.A = A - self.rhs = rhs - self.sol = sol + def test_1to5(self): + rhs = self.rhs + sol = self.sol + Ainv = pymatsolver.Mumps(self.A) + for i in range(3): + self.assertLess( + np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL + ) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - def test_1to5(self): - rhs = self.rhs - sol = self.sol - Ainv = pymatsolver.Mumps(self.A) - for i in range(3): - self.assertLess( - np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + def test_1to5_cmplx(self): + rhs = self.rhs.astype(complex) + sol = self.sol.astype(complex) + self.A = self.A.astype(complex) + Ainv = pymatsolver.Mumps(self.A) + for i in range(3): + self.assertLess( + np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL + ) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - def test_1to5_cmplx(self): - rhs = self.rhs.astype(complex) - sol = self.sol.astype(complex) - self.A = self.A.astype(complex) - Ainv = pymatsolver.Mumps(self.A) - for i in range(3): - self.assertLess( - np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + def test_1to5_T(self): + rhs = self.rhs + sol = self.sol + Ainv = pymatsolver.Mumps(self.A) + AinvT = Ainv.T + for i in range(3): + self.assertLess( + np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]), TOL + ) + self.assertLess(np.linalg.norm(AinvT.T * rhs - sol, np.inf), TOL) - def test_1to5_T(self): - rhs = self.rhs - sol = self.sol - Ainv = pymatsolver.Mumps(self.A) - AinvT = Ainv.T - for i in range(3): - self.assertLess( - np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(AinvT.T * rhs - sol, np.inf), TOL) + # def test_singular(self): + # A = sp.identity(5).tocsr() + # A[-1, -1] = 0 + # self.assertRaises(Exception, pymatsolver.Mumps, A) - # def test_singular(self): - # A = sp.identity(5).tocsr() - # A[-1, -1] = 0 - # self.assertRaises(Exception, pymatsolver.Mumps, A) + def test_multiFactorsInMem(self): + n = 100 + A = sp.rand(n, n, 0.7)+sp.identity(n) + x = np.ones((n, 10)) + rhs = A * x + solvers = [pymatsolver.Mumps(A) for _ in range(20)] - def test_multiFactorsInMem(self): - n = 100 - A = sp.rand(n, n, 0.7)+sp.identity(n) - x = np.ones((n, 10)) - rhs = A * x - solvers = [pymatsolver.Mumps(A) for _ in range(20)] + for Ainv in solvers: + self.assertLess( + np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL) + Ainv.clean() - for Ainv in solvers: - self.assertLess( - np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL) - Ainv.clean() - - for Ainv in solvers: - self.assertLess( - np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL - ) + for Ainv in solvers: + self.assertLess( + np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL + ) if __name__ == '__main__': unittest.main() -""" diff --git a/tests/test_Pardiso.py b/tests/test_Pardiso.py index da10975..adba3df 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -1,180 +1,183 @@ import unittest -from pymatsolver import Pardiso -from pydiso.mkl_solver import ( - get_mkl_pardiso_max_threads, - PardisoTypeConversionWarning -) +try: + from pymatsolver import Pardiso + from pydiso.mkl_solver import ( + get_mkl_pardiso_max_threads, + PardisoTypeConversionWarning + ) +except ImportError: + Pardiso = None import numpy as np import scipy.sparse as sp import os TOL = 1e-10 +if Pardiso: + class TestPardiso(unittest.TestCase): -class TestPardiso(unittest.TestCase): - - def setUp(self): - - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.T*A - A = A.tocsr() - np.random.seed(1) - sol = np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - - def test_refactor(self): - rhs = self.rhs - sol = self.sol - A = self.A - Ainv = Pardiso(A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - - # scale rows and collumns - D = sp.diags(np.random.rand(A.shape[0]) + 1.0) - A2 = D.T @ A @ D - - rhs2 = A2 @ sol - Ainv.factor(A2) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs2[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs2 - sol, np.inf), TOL) - - def test_T(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - - with self.assertWarns(PardisoTypeConversionWarning): - AinvT = Ainv.T - x = AinvT * rhs + def setUp(self): + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.T*A + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + self.A = A + self.rhs = rhs + self.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + for i in range(3): + self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + + def test_refactor(self): + rhs = self.rhs + sol = self.sol + A = self.A + Ainv = Pardiso(A, is_symmetric=True) for i in range(3): - self.assertLess(np.linalg.norm(x[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + + # scale rows and collumns + D = sp.diags(np.random.rand(A.shape[0]) + 1.0) + A2 = D.T @ A @ D + + rhs2 = A2 @ sol + Ainv.factor(A2) + for i in range(3): + self.assertLess(np.linalg.norm(Ainv * rhs2[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs2 - sol, np.inf), TOL) + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + + with self.assertWarns(PardisoTypeConversionWarning): + AinvT = Ainv.T + x = AinvT * rhs + + for i in range(3): + self.assertLess(np.linalg.norm(x[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + + def test_n_threads(self): + max_threads = get_mkl_pardiso_max_threads() + print(f'testing setting n_threads to 1 and {max_threads}') + Ainv = Pardiso(self.A, is_symmetric=True, n_threads=1) + self.assertEqual(Ainv.n_threads, 1) + + Ainv2 = Pardiso(self.A, is_symmetric=True, n_threads=max_threads) + self.assertEqual(Ainv2.n_threads, max_threads) + self.assertEqual(Ainv2.n_threads, Ainv.n_threads) + + Ainv.n_threads = 1 + self.assertEqual(Ainv.n_threads, 1) + self.assertEqual(Ainv2.n_threads, Ainv.n_threads) + + with self.assertRaises(TypeError): + Ainv.n_threads = "2" + + + class TestPardisoNotSymmetric(unittest.TestCase): + + def setUp(self): + + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + np.random.seed(1) + sol = np.random.rand(nSize, 5) + rhs = A.dot(sol) + + self.A = A + self.rhs = rhs + self.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True, check_accuracy=True) + self.assertRaises(Exception, lambda: Ainv * rhs) + Ainv.clean() + + Ainv = Pardiso(self.A) + for i in range(3): + self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + Ainv.clean() + + + class TestPardisoFDEM(unittest.TestCase): + + def setUp(self): + + base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') + + data = np.load(os.path.join(base_path, 'A_data.npy')) + indices = np.load(os.path.join(base_path, 'A_indices.npy')) + indptr = np.load(os.path.join(base_path, 'A_indptr.npy')) + + self.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) + self.rhs = np.load(os.path.join(base_path, 'RHS.npy')) + + def test(self): + rhs = self.rhs + Ainv = Pardiso(self.A, check_accuracy=True) + sol = Ainv * rhs + with self.assertWarns(PardisoTypeConversionWarning): + sol = Ainv * rhs.real + + + class TestPardisoComplex(unittest.TestCase): + + def setUp(self): + nSize = 100 + A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) + A.data = A.data + 1j*np.random.rand(A.nnz) + A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize) + A = A.tocsr() + + np.random.seed(1) + sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) + rhs = A.dot(sol) + + self.A = A + self.rhs = rhs + self.sol = sol - def test_n_threads(self): - max_threads = get_mkl_pardiso_max_threads() - print(f'testing setting n_threads to 1 and {max_threads}') - Ainv = Pardiso(self.A, is_symmetric=True, n_threads=1) - self.assertEqual(Ainv.n_threads, 1) - - Ainv2 = Pardiso(self.A, is_symmetric=True, n_threads=max_threads) - self.assertEqual(Ainv2.n_threads, max_threads) - self.assertEqual(Ainv2.n_threads, Ainv.n_threads) - - Ainv.n_threads = 1 - self.assertEqual(Ainv.n_threads, 1) - self.assertEqual(Ainv2.n_threads, Ainv.n_threads) - - with self.assertRaises(TypeError): - Ainv.n_threads = "2" - - -class TestPardisoNotSymmetric(unittest.TestCase): - - def setUp(self): - - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A = A + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.tocsr() - np.random.seed(1) - sol = np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True, check_accuracy=True) - self.assertRaises(Exception, lambda: Ainv * rhs) - Ainv.clean() - - Ainv = Pardiso(self.A) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - Ainv.clean() - - -class TestPardisoFDEM(unittest.TestCase): - - def setUp(self): - - base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') - - data = np.load(os.path.join(base_path, 'A_data.npy')) - indices = np.load(os.path.join(base_path, 'A_indices.npy')) - indptr = np.load(os.path.join(base_path, 'A_indptr.npy')) - - self.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) - self.rhs = np.load(os.path.join(base_path, 'RHS.npy')) - - def test(self): - rhs = self.rhs - Ainv = Pardiso(self.A, check_accuracy=True) - sol = Ainv * rhs - with self.assertWarns(PardisoTypeConversionWarning): - sol = Ainv * rhs.real - - -class TestPardisoComplex(unittest.TestCase): - - def setUp(self): - nSize = 100 - A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) - A.data = A.data + 1j*np.random.rand(A.nnz) - A = A.T.dot(A) + sp.spdiags(np.ones(nSize), 0, nSize, nSize) - A = A.tocsr() - - np.random.seed(1) - sol = np.random.rand(nSize, 5) + 1j*np.random.rand(nSize, 5) - rhs = A.dot(sol) - - self.A = A - self.rhs = rhs - self.sol = sol - - def test(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - for i in range(3): - self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) - Ainv.clean() - - def test_T(self): - rhs = self.rhs - sol = self.sol - Ainv = Pardiso(self.A, is_symmetric=True) - with self.assertWarns(PardisoTypeConversionWarning): - AinvT = Ainv.T - x = AinvT * rhs + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) for i in range(3): - self.assertLess( - np.linalg.norm(x[:, i] - sol[:, i]), TOL - ) - self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL) + self.assertLess(np.linalg.norm(Ainv * rhs - sol, np.inf), TOL) + Ainv.clean() + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Pardiso(self.A, is_symmetric=True) + with self.assertWarns(PardisoTypeConversionWarning): + AinvT = Ainv.T + x = AinvT * rhs + for i in range(3): + self.assertLess( + np.linalg.norm(x[:, i] - sol[:, i]), TOL + ) + self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) if __name__ == '__main__': unittest.main()