From 6126a45f2c4b9a056deacd3254366322ef7dfc4a Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 12:48:41 -0400 Subject: [PATCH 01/23] Uncomment MUMPS import and MUMPS tests --- pymatsolver/__init__.py | 40 ++++++++++++++++++++-------------------- tests/test_Mumps.py | 2 -- 2 files changed, 20 insertions(+), 22 deletions(-) diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index ead34b0..d4564be 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -71,31 +71,31 @@ AvailableSolvers['Pardiso'] = True PardisoSolver = Pardiso # backwards compatibility except ImportError: - SolverHelp['Pardiso'] = """Pardiso is not working + SolverHelp['Pardiso'] = """Pardiso is not working. Ensure that you have pydiso installed, which may also require Python 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 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 +""" __version__ = '0.2.0' __author__ = 'Rowan Cockett' diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 1bf4afc..79b5533 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -5,7 +5,6 @@ TOL = 1e-11 -""" class TestMumps(unittest.TestCase): if pymatsolver.AvailableSolvers['Mumps']: @@ -83,4 +82,3 @@ def test_multiFactorsInMem(self): if __name__ == '__main__': unittest.main() -""" From 12aeda2f0d1ce8ef16fe5ee35159463348b3262e Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 12:49:29 -0400 Subject: [PATCH 02/23] Distinguish between MUMPS general symmetric and symmetric positive definite --- pymatsolver/mumps/_init_.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 6826e25..a1d4ed0 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -60,7 +60,7 @@ } -class _Pointer(object): +class _Pointer: """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. @@ -81,26 +81,27 @@ class Mumps(Base): documentation:: - http://mumps.enseeiht.fr/doc/userguide_4.10.0.pdf + https://mumps-solver.org/doc/userguide_5.6.1.pdf """ transpose = False - symmetric = False @property def T(self): newMS = self.__class__( self.A, - symmetric=self.symmetric, + symmetric=self.is_symmetric, + positive_definite=self.is_positive_definite, fromPointer=self.pointer ) newMS.transpose = not self.transpose return newMS - def __init__(self, A, symmetric=False, fromPointer=None): + def __init__(self, A, is_symmetric=False, is_positive_definite=False, fromPointer=None): self.A = A.tocsc() - self.symmetric = symmetric + self.is_symmetric = symmetric + self.is_positive_definite = is_positive_definite if fromPointer is None: self.factor() @@ -113,6 +114,14 @@ def __init__(self, A, symmetric=False, fromPointer=None): def isfactored(self): return getattr(self, 'pointer', None) is not None + @property + def _matrix_type(self): + if self.is_symmetric: + if self.is_positive_definite: + return 1 # symmetric, positive definite + return 2 # general symmetric + return 0 # general symmetric + def _funhandle(self, ftype): """ switches the function handle between real and complex. @@ -134,9 +143,8 @@ def factor(self): if self.isfactored: return - sym = 1 if self.symmetric else 0 ierr, p = self._funhandle('F')( - sym, + self._matrix_type, self.A.data, self.A.indices+1, self.A.indptr+1 From afcef6289c1f9f713c698cc33b13c572d15d791d Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 15:14:49 -0400 Subject: [PATCH 03/23] Switch from nose to pytest --- .gitignore | 1 + Makefile | 7 ++++--- README.rst | 2 +- tests/test_Basic.py | 4 ---- tests/test_BicgJacobi.py | 4 ---- tests/test_Mumps.py | 3 --- tests/test_Pardiso.py | 3 --- tests/test_Scipy.py | 4 ---- tests/test_Triangle.py | 3 --- 9 files changed, 6 insertions(+), 25 deletions(-) diff --git a/.gitignore b/.gitignore index c2e9891..8538e31 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.pyc +.coverage *.so build/ dist/ diff --git a/Makefile b/Makefile index 08a96fe..fdaa6c0 100644 --- a/Makefile +++ b/Makefile @@ -4,8 +4,9 @@ build: python setup.py build_ext --inplace coverage: - nosetests --logging-level=INFO --with-coverage --cover-package=pymatsolver --cover-html - open cover/index.html + coverage run --source pymatsolver -m pytest + coverage report -m + coverage html lint: pylint --output-format=html pymatsolver > pylint.html @@ -14,7 +15,7 @@ graphs: pyreverse -my -A -o pdf -p pymatsolver pymatsolver/**.py pymatsolver/**/**.py tests: - nosetests --logging-level=INFO + pytest docs: cd docs;make html diff --git a/README.rst b/README.rst index b321ee4..53c4e5c 100644 --- a/README.rst +++ b/README.rst @@ -64,7 +64,7 @@ From a clean install on Ubuntu: ./miniconda.sh -b export PATH=/root/anaconda/bin:/root/miniconda/bin:$PATH conda update --yes conda - conda install --yes numpy scipy matplotlib cython ipython nose + conda install --yes numpy scipy matplotlib cython ipython pytest coverage git clone https://github.com/rowanc1/pymatsolver.git cd pymatsolver diff --git a/tests/test_Basic.py b/tests/test_Basic.py index adcd4c5..ac3cbc6 100644 --- a/tests/test_Basic.py +++ b/tests/test_Basic.py @@ -24,7 +24,3 @@ def test_DiagonalSolver(self): self.assertLess(np.linalg.norm(sol-X, np.inf), TOL) self.assertLess(np.linalg.norm(sol[:, 0]-x, np.inf), TOL) - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_BicgJacobi.py b/tests/test_BicgJacobi.py index 24e41a5..2735f78 100644 --- a/tests/test_BicgJacobi.py +++ b/tests/test_BicgJacobi.py @@ -84,7 +84,3 @@ def test_T(self): self.A.T*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) self.assertLess(err, TOL) Ainv.clean() - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 79b5533..cf87ce4 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -79,6 +79,3 @@ def test_multiFactorsInMem(self): 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..83387a0 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -175,6 +175,3 @@ def test_T(self): np.linalg.norm(x[:, i] - sol[:, i]), TOL ) self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_Scipy.py b/tests/test_Scipy.py index eb77135..dd4a9de 100644 --- a/tests/test_Scipy.py +++ b/tests/test_Scipy.py @@ -81,7 +81,3 @@ def test_iterative_cg_1(self): def test_iterative_cg_M(self): self.assertLess(dotest(SolverCG, True), TOLI) - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/test_Triangle.py b/tests/test_Triangle.py index 241d1e0..2b1177f 100644 --- a/tests/test_Triangle.py +++ b/tests/test_Triangle.py @@ -29,6 +29,3 @@ def test_directLower_1(self): x = AUinv * self.rhsU[:, 0] self.assertLess(np.linalg.norm(self.sol-X, np.inf), TOL) self.assertLess(np.linalg.norm(self.sol[:, 0]-x, np.inf), TOL) - -if __name__ == '__main__': - unittest.main() From 0c4617f31c1ca0bcbc690515c0d9b18ebbf55fc5 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 15:22:44 -0400 Subject: [PATCH 04/23] Rework skipping of Mumps and Pardiso tests --- tests/test_Mumps.py | 30 +++- tests/test_Pardiso.py | 347 ++++++++++++++++++++++-------------------- 2 files changed, 203 insertions(+), 174 deletions(-) diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index cf87ce4..069c286 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -1,13 +1,21 @@ import unittest +import warnings + import numpy as np import scipy.sparse as sp -import pymatsolver + +try: + from pymatsolver import Mumps + should_run = True +except ImportError: + should_run = False TOL = 1e-11 -class TestMumps(unittest.TestCase): - if pymatsolver.AvailableSolvers['Mumps']: +if should_run: + + class TestMumps(unittest.TestCase): def setUp(self): n = 5 @@ -29,7 +37,7 @@ def setUp(self): def test_1to5(self): rhs = self.rhs sol = self.sol - Ainv = pymatsolver.Mumps(self.A) + Ainv = Mumps(self.A) for i in range(3): self.assertLess( np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL @@ -40,7 +48,7 @@ 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) + Ainv = Mumps(self.A) for i in range(3): self.assertLess( np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]), TOL @@ -50,7 +58,7 @@ def test_1to5_cmplx(self): def test_1to5_T(self): rhs = self.rhs sol = self.sol - Ainv = pymatsolver.Mumps(self.A) + Ainv = Mumps(self.A) AinvT = Ainv.T for i in range(3): self.assertLess( @@ -61,14 +69,14 @@ def test_1to5_T(self): # def test_singular(self): # A = sp.identity(5).tocsr() # A[-1, -1] = 0 - # self.assertRaises(Exception, pymatsolver.Mumps, A) + # self.assertRaises(Exception, 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)] + solvers = [Mumps(A) for _ in range(20)] for Ainv in solvers: self.assertLess( @@ -79,3 +87,9 @@ def test_multiFactorsInMem(self): self.assertLess( np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL ) + +else: + + def test_mumps_not_available(): + """Run only when Mumps is not available.""" + warnings.warn("NOTE: Mumps is not available so we're skipping its tests.") diff --git a/tests/test_Pardiso.py b/tests/test_Pardiso.py index 83387a0..f598b08 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -1,177 +1,192 @@ +import os import unittest -from pymatsolver import Pardiso -from pydiso.mkl_solver import ( - get_mkl_pardiso_max_threads, - PardisoTypeConversionWarning -) +import warnings + import numpy as np import scipy.sparse as sp -import os + +try: + from pymatsolver import Pardiso + from pydiso.mkl_solver import ( + get_mkl_pardiso_max_threads, + PardisoTypeConversionWarning + ) + should_run = True +except ImportError: + should_run = False + TOL = 1e-10 +if should_run: + + 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) -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 + # 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(x[:, i] - sol[:, i]), TOL) - self.assertLess(np.linalg.norm(x - sol, np.inf), TOL) + 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) + +else: + + def test_pardiso_not_available(): + """Run only when Pardiso is not available.""" + warnings.warn("NOTE: Pardiso is not available so we're skipping its tests.") From a7a767061d719a4f80cfe69110f5703304c319d7 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 15:26:44 -0400 Subject: [PATCH 05/23] more nose --- pymatsolver/mumps/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymatsolver/mumps/Makefile b/pymatsolver/mumps/Makefile index e6501e4..c29b41a 100644 --- a/pymatsolver/mumps/Makefile +++ b/pymatsolver/mumps/Makefile @@ -43,4 +43,4 @@ clean: rm -rf *.dSYM test: - cd ..;nosetests + cd ../..;pytest From 645bfde8ad547dd7c0f047f6e3d6aee879c41638 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 15:27:47 -0400 Subject: [PATCH 06/23] Use ImportError instead of exception to detect Mumps availability --- pymatsolver/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymatsolver/__init__.py b/pymatsolver/__init__.py index d4564be..a155254 100644 --- a/pymatsolver/__init__.py +++ b/pymatsolver/__init__.py @@ -81,7 +81,7 @@ from pymatsolver.mumps import Mumps AvailableSolvers['Mumps'] = True MumpsSolver = Mumps # backwards compatibility -except Exception: +except ImportError: SolverHelp['Mumps'] = """Mumps is not working. Ensure that you have Mumps installed, and know where the path to it is. From e0e410e151e133e4f6f63b500cf1e020a42bfc4b Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 16:38:50 -0400 Subject: [PATCH 07/23] Refactor unit tests - now pytest format The Mumps test_singular test is available now. The Mumps test_1to5_T test is still failing. --- tests/test_Basic.py | 13 +++--- tests/test_BicgJacobi.py | 31 ++++++------- tests/test_Mumps.py | 47 +++++++++----------- tests/test_Pardiso.py | 95 +++++++++++++++++++++------------------- tests/test_Scipy.py | 25 ++++------- tests/test_Triangle.py | 22 +++++----- 6 files changed, 111 insertions(+), 122 deletions(-) diff --git a/tests/test_Basic.py b/tests/test_Basic.py index ac3cbc6..bd01a6d 100644 --- a/tests/test_Basic.py +++ b/tests/test_Basic.py @@ -1,4 +1,4 @@ -import unittest +import pytest import numpy as np import scipy.sparse as sp from pymatsolver import Diagonal @@ -6,10 +6,9 @@ TOL = 1e-12 -class TestBasic(unittest.TestCase): +class TestBasic: def test_DiagonalSolver(self): - A = sp.identity(5)*2.0 rhs = np.c_[np.arange(1, 6), np.arange(2, 11, 2)] X = Diagonal(A) * rhs @@ -17,10 +16,10 @@ def test_DiagonalSolver(self): sol = rhs/2.0 - with self.assertRaises(TypeError): + with pytest.raises(TypeError): Diagonal(A, check_accuracy=np.array([1, 2, 3])) - with self.assertRaises(TypeError): + with pytest.raises(TypeError): Diagonal(A, accuracy_tol=0) - self.assertLess(np.linalg.norm(sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(sol[:, 0]-x, np.inf), TOL) + assert np.linalg.norm(sol-X, np.inf) < TOL + assert np.linalg.norm(sol[:, 0]-x, np.inf) < TOL diff --git a/tests/test_BicgJacobi.py b/tests/test_BicgJacobi.py index 2735f78..cd01b3f 100644 --- a/tests/test_BicgJacobi.py +++ b/tests/test_BicgJacobi.py @@ -1,4 +1,3 @@ -import unittest from pymatsolver import BicgJacobi import numpy as np import scipy.sparse as sp @@ -6,9 +5,10 @@ TOL = 1e-6 -class TestBicgJacobi(unittest.TestCase): +class TestBicgJacobi: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) @@ -19,9 +19,9 @@ def setUp(self): sol = np.random.rand(nSize, 4) rhs = A.dot(sol) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.sol = sol def test(self): rhs = self.rhs @@ -30,7 +30,7 @@ def test(self): for i in range(3): err = np.linalg.norm( self.A*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() def test_T(self): @@ -42,13 +42,14 @@ def test_T(self): for i in range(3): err = np.linalg.norm( self.A.T*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() -class TestPardisoComplex(unittest.TestCase): +class TestBicgJacobiComplex: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) A.data = A.data + 1j*np.random.rand(A.nnz) @@ -59,9 +60,9 @@ def setUp(self): 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 + cls.A = A + cls.rhs = rhs + cls.sol = sol def test(self): rhs = self.rhs @@ -70,7 +71,7 @@ def test(self): for i in range(3): err = np.linalg.norm( self.A*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() def test_T(self): @@ -82,5 +83,5 @@ def test_T(self): for i in range(3): err = np.linalg.norm( self.A.T*solb[:, i] - rhs[:, i]) / np.linalg.norm(rhs[:, i]) - self.assertLess(err, TOL) + assert err < TOL Ainv.clean() diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 069c286..ab3d37f 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -1,7 +1,7 @@ -import unittest import warnings import numpy as np +import pytest import scipy.sparse as sp try: @@ -15,9 +15,10 @@ if should_run: - class TestMumps(unittest.TestCase): + class TestMumps: - def setUp(self): + @classmethod + def setup_class(cls): 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] @@ -30,19 +31,17 @@ def setUp(self): 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 + cls.A = A + cls.rhs = rhs + cls.sol = sol def test_1to5(self): rhs = self.rhs sol = self.sol Ainv = 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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL def test_1to5_cmplx(self): rhs = self.rhs.astype(complex) @@ -50,10 +49,8 @@ def test_1to5_cmplx(self): self.A = self.A.astype(complex) Ainv = 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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL def test_1to5_T(self): rhs = self.rhs @@ -61,15 +58,14 @@ def test_1to5_T(self): Ainv = 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) + assert np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]) < TOL + assert 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, Mumps, A) + def test_singular(self): + A = sp.identity(5).tocsr() + A[-1, -1] = 0 + with pytest.raises(Exception): + Mumps(A) def test_multiFactorsInMem(self): n = 100 @@ -79,14 +75,11 @@ def test_multiFactorsInMem(self): solvers = [Mumps(A) for _ in range(20)] for Ainv in solvers: - self.assertLess( - np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs), TOL) + assert 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 - ) + assert np.linalg.norm(Ainv * rhs - x)/np.linalg.norm(rhs) < TOL else: diff --git a/tests/test_Pardiso.py b/tests/test_Pardiso.py index f598b08..a321895 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -1,8 +1,8 @@ import os -import unittest import warnings import numpy as np +import pytest import scipy.sparse as sp try: @@ -20,9 +20,10 @@ if should_run: - class TestPardiso(unittest.TestCase): + class TestPardiso: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) @@ -33,17 +34,17 @@ def setUp(self): sol = np.random.rand(nSize, 5) rhs = A.dot(sol) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL def test_refactor(self): rhs = self.rhs @@ -51,8 +52,8 @@ def test_refactor(self): 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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL # scale rows and collumns D = sp.diags(np.random.rand(A.shape[0]) + 1.0) @@ -61,43 +62,44 @@ def test_refactor(self): 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) + assert np.linalg.norm(Ainv * rhs2[:, i] - sol[:, i]) < TOL + assert 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): + with pytest.warns(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) + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert 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) + assert 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) + assert Ainv2.n_threads == max_threads + assert Ainv2.n_threads == Ainv.n_threads Ainv.n_threads = 1 - self.assertEqual(Ainv.n_threads, 1) - self.assertEqual(Ainv2.n_threads, Ainv.n_threads) + assert Ainv.n_threads == 1 + assert Ainv2.n_threads == Ainv.n_threads - with self.assertRaises(TypeError): + with pytest.raises(TypeError): Ainv.n_threads = "2" - class TestPardisoNotSymmetric(unittest.TestCase): + class TestPardisoNotSymmetric: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) @@ -107,27 +109,29 @@ def setUp(self): sol = np.random.rand(nSize, 5) rhs = A.dot(sol) - self.A = A - self.rhs = rhs - self.sol = sol + cls.A = A + cls.rhs = rhs + cls.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) + with pytest.raises(Exception): + 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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL Ainv.clean() - class TestPardisoFDEM(unittest.TestCase): + class TestPardisoFDEM: - def setUp(self): + @classmethod + def setup_class(cls): base_path = os.path.join(os.path.split(os.path.abspath(__file__))[0], 'fdem') @@ -135,20 +139,21 @@ def setUp(self): 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')) + cls.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) + cls.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): + with pytest.warns(PardisoTypeConversionWarning): sol = Ainv * rhs.real - class TestPardisoComplex(unittest.TestCase): + class TestPardisoComplex: - def setUp(self): + @classmethod + def setup_class(cls): nSize = 100 A = sp.rand(nSize, nSize, 0.05, format='csr', random_state=100) A.data = A.data + 1j*np.random.rand(A.nnz) @@ -159,31 +164,29 @@ def setUp(self): 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 + cls.A = A + cls.rhs = rhs + cls.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) + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert 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): + with pytest.warns(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) + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL else: diff --git a/tests/test_Scipy.py b/tests/test_Scipy.py index dd4a9de..8969ba7 100644 --- a/tests/test_Scipy.py +++ b/tests/test_Scipy.py @@ -1,4 +1,3 @@ -import unittest from pymatsolver import Solver, Diagonal, SolverCG, SolverLU import scipy.sparse as sp import numpy as np @@ -50,34 +49,28 @@ def dotest(MYSOLVER, multi=False, A=None, **solverOpts): return np.linalg.norm(e-x, np.inf) -class TestSolver(unittest.TestCase): +class TestSolver: def test_direct_spsolve_1(self): - self.assertLess(dotest(Solver, False), TOLD) + assert dotest(Solver, False) < TOLD def test_direct_spsolve_M(self): - self.assertLess(dotest(Solver, True), TOLD) + assert dotest(Solver, True) < TOLD def test_direct_splu_1(self): - self.assertLess(dotest(SolverLU, False), TOLD) + assert dotest(SolverLU, False) < TOLD def test_direct_splu_M(self): - self.assertLess(dotest(SolverLU, True), TOLD) + assert dotest(SolverLU, True) < TOLD def test_iterative_diag_1(self): - self.assertLess(dotest( - Diagonal, False, - A=sp.diags(np.random.rand(10)+1.0) - ), TOLI) + assert dotest(Diagonal, False, A=sp.diags(np.random.rand(10)+1.0)) < TOLI def test_iterative_diag_M(self): - self.assertLess(dotest( - Diagonal, True, - A=sp.diags(np.random.rand(10)+1.0) - ), TOLI) + assert dotest(Diagonal, True, A=sp.diags(np.random.rand(10)+1.0)) < TOLI def test_iterative_cg_1(self): - self.assertLess(dotest(SolverCG, False), TOLI) + assert dotest(SolverCG, False) < TOLI def test_iterative_cg_M(self): - self.assertLess(dotest(SolverCG, True), TOLI) + assert dotest(SolverCG, True) < TOLI diff --git a/tests/test_Triangle.py b/tests/test_Triangle.py index 2b1177f..98ee279 100644 --- a/tests/test_Triangle.py +++ b/tests/test_Triangle.py @@ -1,4 +1,3 @@ -import unittest import numpy as np import scipy.sparse as sp import pymatsolver @@ -6,26 +5,27 @@ TOL = 1e-12 -class TestTriangle(unittest.TestCase): +class TestTriangle: - def setUp(self): + @classmethod + def setup_class(cls): n = 50 nrhs = 20 - self.A = sp.rand(n, n, 0.4) + sp.identity(n) - self.sol = np.ones((n, nrhs)) - self.rhsU = sp.triu(self.A) * self.sol - self.rhsL = sp.tril(self.A) * self.sol + cls.A = sp.rand(n, n, 0.4) + sp.identity(n) + cls.sol = np.ones((n, nrhs)) + cls.rhsU = sp.triu(cls.A) * cls.sol + cls.rhsL = sp.tril(cls.A) * cls.sol def test_directLower(self): ALinv = pymatsolver.Forward(sp.tril(self.A)) X = ALinv * self.rhsL x = ALinv * self.rhsL[:, 0] - self.assertLess(np.linalg.norm(self.sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(self.sol[:, 0]-x, np.inf), TOL) + assert np.linalg.norm(self.sol-X, np.inf) < TOL + assert np.linalg.norm(self.sol[:, 0]-x, np.inf) < TOL def test_directLower_1(self): AUinv = pymatsolver.Backward(sp.triu(self.A)) X = AUinv * self.rhsU x = AUinv * self.rhsU[:, 0] - self.assertLess(np.linalg.norm(self.sol-X, np.inf), TOL) - self.assertLess(np.linalg.norm(self.sol[:, 0]-x, np.inf), TOL) + assert np.linalg.norm(self.sol-X, np.inf) < TOL + assert np.linalg.norm(self.sol[:, 0]-x, np.inf) < TOL From 8f94db02b42c80f0baa258c2d53701febf1f5f6c Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 16:43:29 -0400 Subject: [PATCH 08/23] Gross Mumps workaround for repeated transpose operation --- tests/test_Mumps.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index ab3d37f..73c33b1 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -57,8 +57,10 @@ def test_1to5_T(self): sol = self.sol Ainv = Mumps(self.A) AinvT = Ainv.T + AinvTT = AinvT.T for i in range(3): - assert np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]) < TOL + # WARNING: This test fails if AinvT.T is supplied instead of AinvTT + assert np.linalg.norm(AinvTT * rhs[:, i] - sol[:, i]) < TOL assert np.linalg.norm(AinvT.T * rhs - sol, np.inf) < TOL def test_singular(self): From b3d44b2b0781b8733edb3cfe5c2f9e8d20cf4b0a Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 16:51:46 -0400 Subject: [PATCH 09/23] Clean up Mumps changes --- pymatsolver/mumps/_init_.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index a1d4ed0..224e18a 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -4,6 +4,7 @@ from __future__ import absolute_import import gc +import warnings from pymatsolver.solvers import Base from . import MumpsInterface as _MUMPSINT @@ -91,27 +92,27 @@ class Mumps(Base): def T(self): newMS = self.__class__( self.A, - symmetric=self.is_symmetric, - positive_definite=self.is_positive_definite, - fromPointer=self.pointer + is_symmetric=self.is_symmetric, + is_positive_definite=self.is_positive_definite, + from_pointer=self.pointer ) newMS.transpose = not self.transpose return newMS - def __init__(self, A, is_symmetric=False, is_positive_definite=False, fromPointer=None): + def __init__(self, A, is_symmetric=False, is_positive_definite=False, from_pointer=None): self.A = A.tocsc() - self.is_symmetric = symmetric + self.is_symmetric = is_symmetric self.is_positive_definite = is_positive_definite - if fromPointer is None: + if from_pointer is None: self.factor() - elif isinstance(fromPointer, _Pointer): - self.pointer = fromPointer + elif isinstance(from_pointer, _Pointer): + self.pointer = from_pointer else: raise Exception('Unknown pointer for construction.') @property - def isfactored(self): + def is_factored(self): return getattr(self, 'pointer', None) is not None @property @@ -120,7 +121,7 @@ def _matrix_type(self): if self.is_positive_definite: return 1 # symmetric, positive definite return 2 # general symmetric - return 0 # general symmetric + return 0 # unsymmetric def _funhandle(self, ftype): """ @@ -140,7 +141,7 @@ def _funhandle(self, ftype): 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] def factor(self): - if self.isfactored: + if self.is_factored: return ierr, p = self._funhandle('F')( @@ -152,7 +153,7 @@ def factor(self): if ierr < 0: raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) elif ierr > 0: - print("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + warnings.warn("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) self.pointer = _Pointer(p, self._funhandle('D')) From 6fcd2a7b62228daa452c4eef8520667d06301590 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 20:43:02 -0400 Subject: [PATCH 10/23] Handle unknown MUMPS errors and warnings Also correctly handle warnings as flags. --- pymatsolver/mumps/_init_.py | 131 +++++++++++++++++++++--------------- 1 file changed, 78 insertions(+), 53 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 224e18a..12a31e1 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -8,57 +8,82 @@ 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.", -} + +def _mumps_message_from_exit_code(ierr): + """Convert a MUMPS routine's exit code to a message. + + See section 8 of the MUMPS user guide for a full catalog of these messages. + """ + error_messages = { + -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.", + } + + warning_messages = { + 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.", + } + + if ierr == 0: + return "No error occurred." + elif ierr < 0: + return error_messages.get(ierr, "An unknown error occurred.") + else: + # The warning values are actually bitwise flags. + messages = [] + for bit in warning_messages.keys(): + if ierr & bit != 0: + messages.append(warning_messages[bit]) + + if not messages: + messages.append("An unknown warning occurred.") + + return "\n".join(messages) class _Pointer: @@ -151,9 +176,9 @@ def factor(self): self.A.indptr+1 ) if ierr < 0: - raise Exception("Mumps Exception [{}] - {}".format(ierr, _mumpsErrors[ierr])) + raise Exception(f"Mumps Exception [{ierr}] - {_mumps_message_from_exit_code(ierr)}") elif ierr > 0: - warnings.warn("Mumps Warning [{}] - {}".format(ierr, _mumpsErrors[ierr])) + warnings.warn(f"Mumps Warning [{ierr}] - {_mumps_message_from_exit_code(ierr)}") self.pointer = _Pointer(p, self._funhandle('D')) From 24c9935f9dabd461006299d32e4f2956f7acb59c Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Tue, 15 Aug 2023 20:52:35 -0400 Subject: [PATCH 11/23] New errors from MUMPS 5.6.1 --- pymatsolver/mumps/_init_.py | 55 +++++++++++++++++++++++++++---------- 1 file changed, 40 insertions(+), 15 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 12a31e1..c63d4b2 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -16,24 +16,26 @@ def _mumps_message_from_exit_code(ierr): """ error_messages = { -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.", + -2: "NNZ/NZ, NNZ loc/NZ loc or P NNZ loc/P NZ loc are out of range. INFO(2)=NNZ/NZ, NNZ loc/NZ loc or P NNZ loc/P NZ loc.", + -3: "MUMPS was called with an invalid value for JOB. This may happen if the analysis (JOB=1) was not performed (or failed) before the factorization (JOB=2), or the factorization was not performed (or failed) 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 4. This error also occurs if JOB does not contain the same value on all processes on entry to MUMPS. INFO(2) is then set to the local value of JOB.", + -4: "Error in user-provided permutation array PERM IN at position INFO(2). This error may only occur on the host.", + -5: 'Problem of real workspace allocation of size INFO(2) during analysis. The unit for INFO(2) is the number of real values (single precision for SMUMPS/CMUMPS, double precision for DMUMPS/ZMUMPS), in the Fortran "ALLOCATE" statement that did not succeed. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.', -6: "Matrix is singular in structure. INFO(2) holds the structural rank.", - -7: "Problem of INTEGER workspace allocation of size INFO(2) during analysis.", + -7: "Problem of integer workspace allocation of size INFO(2) during analysis. The unit for INFO(2) is the number of integer values that MUMPS tried to allocate in the Fortran ALLOCATE statement that did not succeed. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.", -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).", + -9: "The main internal real/complex workarray S is 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 LWK USER is provided LWK USER should be increased.", + -10: "Numerically singular matrix. INFO(2) holds the number of eliminated pivots.", + -11: "Internal real/complex workarray S or LWK USER too small for solution. If INFO(2) is positive, then the number of entries that are missing in S/LWK USER at the moment when the error is raised is available in INFO(2). If the numerical phases are out-of-core and LWK USER is provided for the solution phase and is smaller than the value provided for the factorization, it should be increased by at least INFO(2). In other cases, please contact us.", -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.", + -13: "Problem of workspace allocation of size INFO(2) during the factorization or solve steps. The size that the package tried to allocate with a Fortran ALLOCATE statement 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. In general, the unit for INFO(2) is the number of scalar entries of the type of the input matrix (real, complex, single or double precision).", -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.", + -18: "The blocking size for multiple RHS (ICNTL(27)) is too large and may lead to an integer overflow. This error may only occurs for very large matrices with large values of ICNTL(27) (e.g., several thousands). INFO(2) provides an estimate of the maximum value of ICNTL(27) that should be used.", + -19: "The maximum allowed size of working memory ICNTL(23) is too small to run the factorization phase and should be increased. If INFO(2) is positive, then the number of entries that are missing 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.", + -20: "The internal reception buffer that was allocated dynamically by MUMPS is too small. Normally, this error is raised on the sender side when detecting that the message to be sent is too large for the reception buffer on the receiver. 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 MUMPS in host-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.", @@ -45,21 +47,44 @@ def _mumps_message_from_exit_code(ierr): -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).", + -33: "ICNTL(26) was asked for during solve phase (or during the factorization - see ICNTL(32)) 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).", + -35: "This error is raised when the expansion phase is called (ICNTL(26) = 2) but reduction phase (ICNTL(26)=1) was not called before. This error also occurs in case the reduction phase (ICNTL(26) = 1) is asked for at the solution phase (JOB=3) but the forward elimination was already performed during the factorization phase (JOB=2 and ICNTL(32)=1). 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).", + -37: "Value of ICNTL(25) incompatible with some other parameter. If ICNTL(25) is incompatible with ICNTL(xx), 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.", + -41: "Incompatible value of LWK USER between factorization and solution phases. This error may only occur when the factorization is in-core (ICNTL(22)=1), in which case both the contents of WK USER and LWK USER should be passed unchanged between the factorization (JOB=2) and solution (JOB=3) phases.", + -42: "ICNTL(32) was set to 1 (forward during factorization), but the value of NRHS on the host processor is incorrect: either the value of NRHS provided at analysis is negative or zero, or the value provided at factorization or solve is different from the value provided at analysis. INFO(2) holds the value of id%NRHS that was provided at analysis.", + -43: "Incompatible values of ICNTL(32) and ICNTL(xx). The index xx is stored in INFO(2).", -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).", + -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).", + -50: "An error occurred while computing the fill-reducing ordering during the analysis phase. This commonly happens when an (external) ordering tool returns an error code or a wrong result.", + -51: "An external ordering (Metis/ParMetis, SCOTCH/PT-SCOTCH, PORD), with 32-bit default integers, is invoked to processing a graph of size larger than 231 - 1. INFO(2) holds the size required to store the graph as a number of integer values; it is negative and its absolute value should be multiplied by 1 million.", + -52: "When default Fortran integers are 64 bit (e.g. Fortran compiler flag -i8 -fdefault-integer-8 or something equivalent depending on your compiler) then external ordering libraries (Metis/ParMetis, SCOTCH/PT-SCOTCH, PORD) should also have 64-bit default integers. INFO(2) = 1, 2, 3 means that respectively Metis/ParMetis, SCOTCH/PT-SCOTCH or PORD were invoked and were not generated with 64-bit default integers.", + -53: "Internal error that could be due to inconsistent input data between two consecutive calls.", + -54: "The analysis phase (JOB=1) was called with ICNTL(35)=0 but the factorization phase was called with ICNTL(35)=1, 2 or 3. In order to perform the factorization with BLR compression, please perform the analysis phase again using ICNTL(35)=1, 2 or 3 (see the documentation of ICNTL(35)).", + -55: "During a call to MUMPS including the solve phase with distributed right-hand side, LRHS loc was detected to be smaller than Nloc RHS. INFO(2)=LRHS loc.", + -56: "During a call to MUMPS including the solve phase with distributed right-hand side and distributed solution, RHS loc and SOL loc point to the same workarray but LRHS loc < LSOL loc. INFO(2)=LRHS loc.", + -57: "During a call to MUMPS analysis phase with a block format (ICNTL(15) != 0), an error in the interface provided by the user was detected. INFO(2) holds additional information about the issue.", + -70: "During a call to MUMPS with JOB=7, the file specified to save the current instance, as derived from SAVE DIR and/or SAVE PREFIX, already exists. Before saving an instance into this file, it should be first suppressed (see JOB=-3). Otherwise, a different file should be specified by changing the values of SAVE DIR and/or SAVE PREFIX.", + -71: "An error has occured during the creation of one of the files needed to save MUMPS data (JOB=7).", + -72: "Error while saving data (JOB=7); a write operation did not succeed (e.g., disk full, I/O error, . . . ). INFO(2) is the size that should have been written during that operation. If INFO(2) is negative, then its absolute value should be multiplied by 1 million.", + -73: "During a call to MUMPS with JOB=8, one parameter of the current instance is not compatible with the corresponding one in the saved instance.", + -74: "The file resulting from the setting of SAVE DIR and SAVE PREFIX could not be opened for restoring data (JOB=8). INFO(2) is the rank of the process (in the communicator COMM) on which the error was detected.", + -75: "Error while restoring data (JOB=8); a read operation did not succeed (e.g., end of file reached, I/O error, . . . ). INFO(2) is the size still to be read. 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.", + -76: "Error while deleting the files (JOB=-3); some files to be erased were not found or could not be suppressed. INFO(2) is the rank of the process (in the communicator COMM) on which the error was detected.", + -77: "Neither SAVE DIR nor the environment variable MUMPS SAVE DIR are defined.", + -78: "Problem of workspace allocation during the restore step. The size still to be allocated 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.", + -79: "MUMPS could not find a Fortran file unit to perform I/O's.", + -89: "Internal error during SCOTCH kway-partitioning in SCOTCHFGRAPHPART. METIS package should be made available to MUMPS.", -90: "Error in out-of-core management. See the error message returned on output unit ICNTL(1) for more information.", + -800: "Temporary error associated to the current MUMPS release, subject to change or disappearance in the future. If INFO(2)=5, then this error is due to the fact that the elemental matrix format (ICNTL(5)=1) is currently incompatible with a BLR factorization (ICNTL(35) != 0).", } warning_messages = { From 17e51261e2073bd5ed5db602e6e608d015e9175b Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Wed, 16 Aug 2023 14:48:21 -0400 Subject: [PATCH 12/23] Fix misplaced call to Mumps destroy function There was an issue where if you copied a matrix and the copy passed out of scope, the original matrix could no longer be used. This issue was the reason that the test_1to5_T test was failing. --- pymatsolver/mumps/_init_.py | 35 +++++++++++++++++++---------------- tests/test_Mumps.py | 4 +--- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index c63d4b2..56b4906 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -112,20 +112,26 @@ def _mumps_message_from_exit_code(ierr): class _Pointer: - """Gets an int and a destroy call that gets called on garbage collection. + """Store a pointer with a destroy function that gets called on garbage collection. - There can be multiple Solvers around the place that are pointing to the same factor in memory. + Multiple Mumps solver instances can point to the same factor in memory. + This class ensures that the factor is garbage collected only when all references are removed. - 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. + You can always clean the Solver, and it will explicitly call the destroy method. """ - def __init__(self, pointerINT, destroyCall): + def __init__(self, pointerINT, dtype): self.INT = pointerINT - self.destroyCall = destroyCall + # The destroy function should only be called from this class upon deallocation. + if dtype == float: + self._destroy_func = _MUMPSINT.destroy_mumps + elif dtype == complex: + self._destroy_func = _MUMPSINT.destroy_mumps_cmplx + else: + raise ValueError(f"Attempted to use an invalid data type ({dtype})") def __del__(self): - self.destroyCall(self.INT) + self._destroy_func(self.INT) + class Mumps(Base): """ @@ -177,18 +183,16 @@ def _funhandle(self, ftype): """ switches the function handle between real and complex. - ftype in ['F','S','D'] + ftype in ['F','S'] - means factor, solve, destroy + means factor, solve """ if self.A.dtype == float: return {'F': _MUMPSINT.factor_mumps, - 'S': _MUMPSINT.solve_mumps, - 'D': _MUMPSINT.destroy_mumps}[ftype] + 'S': _MUMPSINT.solve_mumps}[ftype] elif self.A.dtype == complex: return {'F': _MUMPSINT.factor_mumps_cmplx, - 'S': _MUMPSINT.solve_mumps_cmplx, - 'D': _MUMPSINT.destroy_mumps_cmplx}[ftype] + 'S': _MUMPSINT.solve_mumps_cmplx}[ftype] def factor(self): if self.is_factored: @@ -205,7 +209,7 @@ def factor(self): elif ierr > 0: warnings.warn(f"Mumps Warning [{ierr}] - {_mumps_message_from_exit_code(ierr)}") - self.pointer = _Pointer(p, self._funhandle('D')) + self.pointer = _Pointer(p, self.A.dtype) def _solveM(self, rhs): self.factor() @@ -219,6 +223,5 @@ def _solveM(self, rhs): _solve1 = _solveM def clean(self): - self._funhandle('D')(self.pointer.INT) del self.pointer gc.collect() diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 73c33b1..ab3d37f 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -57,10 +57,8 @@ def test_1to5_T(self): sol = self.sol Ainv = Mumps(self.A) AinvT = Ainv.T - AinvTT = AinvT.T for i in range(3): - # WARNING: This test fails if AinvT.T is supplied instead of AinvTT - assert np.linalg.norm(AinvTT * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(AinvT.T * rhs[:, i] - sol[:, i]) < TOL assert np.linalg.norm(AinvT.T * rhs - sol, np.inf) < TOL def test_singular(self): From 539c330ebe29a039316108ed3d41546eef6c7b94 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Wed, 16 Aug 2023 15:13:14 -0400 Subject: [PATCH 13/23] Make function handle logic easier to follow --- pymatsolver/mumps/_init_.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 56b4906..cfb8a32 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -179,26 +179,26 @@ def _matrix_type(self): return 2 # general symmetric return 0 # unsymmetric - def _funhandle(self, ftype): - """ - switches the function handle between real and complex. - - ftype in ['F','S'] - - means factor, solve - """ + def _funhandle(self, function): + """Switch the function handle between real and complex.""" if self.A.dtype == float: - return {'F': _MUMPSINT.factor_mumps, - 'S': _MUMPSINT.solve_mumps}[ftype] + return { + 'factor': _MUMPSINT.factor_mumps, + 'solve': _MUMPSINT.solve_mumps + }[function] elif self.A.dtype == complex: - return {'F': _MUMPSINT.factor_mumps_cmplx, - 'S': _MUMPSINT.solve_mumps_cmplx}[ftype] + return { + 'factor': _MUMPSINT.factor_mumps_cmplx, + 'solve': _MUMPSINT.solve_mumps_cmplx + }[function] + + raise ValueError(f"Attempted to use an invalid data type ({self.A.dtype})") def factor(self): if self.is_factored: return - ierr, p = self._funhandle('F')( + ierr, p = self._funhandle('factor')( self._matrix_type, self.A.data, self.A.indices+1, @@ -217,7 +217,7 @@ def _solveM(self, rhs): 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) + sol = self._funhandle('solve')(self.pointer.INT, nrhs, rhs, T) return sol _solve1 = _solveM From 7cec0a1054d65aa6dcc45e825905a65979d2128d Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Wed, 16 Aug 2023 15:58:03 -0400 Subject: [PATCH 14/23] Update GitHub repo URL --- README.rst | 6 +++--- setup.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/README.rst b/README.rst index 53c4e5c..74557b3 100644 --- a/README.rst +++ b/README.rst @@ -9,7 +9,7 @@ pymatsolver :target: https://github.com/simpeg/pymatsolver/blob/master/LICENSE :alt: MIT license. -.. image:: https://img.shields.io/travis/rowanc1/pymatsolver.svg +.. image:: https://img.shields.io/travis/simpeg/pymatsolver.svg :target: https://travis-ci.org/simpeg/pymatsolver :alt: Travis CI build status @@ -66,7 +66,7 @@ From a clean install on Ubuntu: conda update --yes conda conda install --yes numpy scipy matplotlib cython ipython pytest coverage - git clone https://github.com/rowanc1/pymatsolver.git + git clone https://github.com/simpeg/pymatsolver.git cd pymatsolver make mumps @@ -79,7 +79,7 @@ This assumes that you have Brew and some python installed (numpy, scipy): brew install mumps --with-scotch5 --without-mpi - git clone https://github.com/rowanc1/pymatsolver.git + git clone https://github.com/simpeg/pymatsolver.git cd pymatsolver make mumps_mac diff --git a/setup.py b/setup.py index f7681ef..8fc2c98 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,7 @@ 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 From b22147adc079150c4ac97b6e296c3f93a15c2ea7 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Wed, 16 Aug 2023 16:22:30 -0400 Subject: [PATCH 15/23] Mumps is_factored is private --- pymatsolver/mumps/_init_.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index cfb8a32..822e069 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -168,7 +168,7 @@ def __init__(self, A, is_symmetric=False, is_positive_definite=False, from_point raise Exception('Unknown pointer for construction.') @property - def is_factored(self): + def _is_factored(self): return getattr(self, 'pointer', None) is not None @property @@ -195,7 +195,7 @@ def _funhandle(self, function): raise ValueError(f"Attempted to use an invalid data type ({self.A.dtype})") def factor(self): - if self.is_factored: + if self._is_factored: return ierr, p = self._funhandle('factor')( From 26f3cb8dec14e63b211f815df08fdc6ad24c2402 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Thu, 17 Aug 2023 13:32:47 -0400 Subject: [PATCH 16/23] Fix "collumns" typo --- tests/test_Pardiso.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_Pardiso.py b/tests/test_Pardiso.py index a321895..53d5e3f 100644 --- a/tests/test_Pardiso.py +++ b/tests/test_Pardiso.py @@ -55,7 +55,7 @@ def test_refactor(self): assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - # scale rows and collumns + # scale rows and columns D = sp.diags(np.random.rand(A.shape[0]) + 1.0) A2 = D.T @ A @ D From 67351f814e65992635357e174dfbbde3fe897365 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Thu, 17 Aug 2023 13:33:13 -0400 Subject: [PATCH 17/23] Update shape, then check accuracy Otherwise you get ValueError: dimension mismatch when the solver doesn't do the reshaping --- pymatsolver/solvers.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/pymatsolver/solvers.py b/pymatsolver/solvers.py index 890e508..0bf7111 100644 --- a/pymatsolver/solvers.py +++ b/pymatsolver/solvers.py @@ -87,13 +87,16 @@ def _solve(self, rhs): else: x = self._solveM(rhs) + if nrhs == 1: + x = x.flatten() + elif nrhs > 1: + x = x.reshape((n, nrhs), order='F') + if self.check_accuracy: self._compute_accuracy(rhs, x) - if nrhs == 1: - return x.flatten() - elif nrhs > 1: - return x.reshape((n, nrhs), order='F') + return x + def clean(self): pass From 47ef682a317d8535d201fd5794912a462caa79c8 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Thu, 17 Aug 2023 13:34:22 -0400 Subject: [PATCH 18/23] Full attribute parity and result parity between Pardiso and Mumps It's now possible to use the same solver_opts values in SimPEG for both the Pardiso and MUMPS solvers. We now run all the Pardiso test cases for Mumps as well. Running these test cases uncovered a Mumps issue where matrix transposes didn't commute when operating on complex matrices. If you asked the Mumps solver: Ainv = Mumps(A) AinvT_1 = Mumps(A).T AinvT_2 = Mumps(A.T) then the two inverse transposes weren't equal when A is a complex matrix. For small matrices, you can look at the inverse transposes by solving this trivial system: AinvT_1 * np.identity(A.size) != AinvT_2 * np.identity(A.size) This issue occurred because pymatsolver treated the complex case as the conjugate transpose -- but we weren't asking for the conjugate transpose; we were explicitly asking for just the transpose. There was no issue with Mumps itself, just a few extra lines of Fortran code in the Mumps interface. Now that those lines are removed, the Mumps solver and the Pardiso solver operate identically with these transposes. --- pymatsolver/mumps/_init_.py | 20 ++-- pymatsolver/mumps/mumps_cmplx_p.f90 | 6 -- tests/test_Mumps.py | 137 +++++++++++++++++++++++++++- 3 files changed, 145 insertions(+), 18 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 822e069..98a40f8 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -146,19 +146,23 @@ class Mumps(Base): @property def T(self): + properties_with_setters = [] + for a in dir(self.__class__): + attr = getattr(self.__class__, a) + if hasattr(attr, "fset") and attr.fset is not None: + properties_with_setters.append(a) + kwargs = {attr: getattr(self, attr) for attr in properties_with_setters} + newMS = self.__class__( - self.A, - is_symmetric=self.is_symmetric, - is_positive_definite=self.is_positive_definite, - from_pointer=self.pointer + self.A.T, + **kwargs, ) - newMS.transpose = not self.transpose return newMS - def __init__(self, A, is_symmetric=False, is_positive_definite=False, from_pointer=None): + def __init__(self, A, from_pointer=None, **kwargs): self.A = A.tocsc() - self.is_symmetric = is_symmetric - self.is_positive_definite = is_positive_definite + # As of 5.6.1, MUMPS has no optimizations for Hermitian matrices + self.set_kwargs(**kwargs, ignore="is_hermitian") if from_pointer is None: self.factor() diff --git a/pymatsolver/mumps/mumps_cmplx_p.f90 b/pymatsolver/mumps/mumps_cmplx_p.f90 index 4fb056d..bdd0502 100644 --- a/pymatsolver/mumps/mumps_cmplx_p.f90 +++ b/pymatsolver/mumps/mumps_cmplx_p.f90 @@ -229,17 +229,11 @@ subroutine solve( mumps_par, nrhs, rhs, x, transpose ) 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 diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index ab3d37f..8fa19bc 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -1,3 +1,4 @@ +import os import warnings import numpy as np @@ -17,6 +18,134 @@ class TestMumps: + @classmethod + def setup_class(cls): + + 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) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + AinvT = Ainv.T + x = AinvT * rhs + + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + + class TestMumpsNotSymmetric: + + @classmethod + def setup_class(cls): + + 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) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True, check_accuracy=True) + with pytest.raises(Exception): + Ainv * rhs + Ainv.clean() + + Ainv = Mumps(self.A, check_accuracy=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + + + class TestMumpsFDEM: + + @classmethod + def setup_class(cls): + + 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')) + + cls.A = sp.csr_matrix((data, indices, indptr), shape=(13872, 13872)) + cls.rhs = np.load(os.path.join(base_path, 'RHS.npy')) + + def test(self): + rhs = self.rhs + Ainv = Mumps(self.A, check_accuracy=True) + sol = Ainv * rhs + sol = Ainv * rhs.real + + + class TestMumpsComplex: + + @classmethod + def setup_class(cls): + 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) + + cls.A = A + cls.rhs = rhs + cls.sol = sol + + def test(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + for i in range(3): + assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL + Ainv.clean() + + def test_T(self): + rhs = self.rhs + sol = self.sol + Ainv = Mumps(self.A, is_symmetric=True) + AinvT = Ainv.T + x = AinvT * rhs + for i in range(3): + assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL + assert np.linalg.norm(x - sol, np.inf) < TOL + + + class TestMumps1to5: + @classmethod def setup_class(cls): n = 5 @@ -35,7 +164,7 @@ def setup_class(cls): cls.rhs = rhs cls.sol = sol - def test_1to5(self): + def test(self): rhs = self.rhs sol = self.sol Ainv = Mumps(self.A) @@ -43,7 +172,7 @@ def test_1to5(self): assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - def test_1to5_cmplx(self): + def test_cmplx(self): rhs = self.rhs.astype(complex) sol = self.sol.astype(complex) self.A = self.A.astype(complex) @@ -52,7 +181,7 @@ def test_1to5_cmplx(self): assert np.linalg.norm(Ainv * rhs[:, i] - sol[:, i]) < TOL assert np.linalg.norm(Ainv * rhs - sol, np.inf) < TOL - def test_1to5_T(self): + def test_T(self): rhs = self.rhs sol = self.sol Ainv = Mumps(self.A) @@ -67,7 +196,7 @@ def test_singular(self): with pytest.raises(Exception): Mumps(A) - def test_multiFactorsInMem(self): + def test_multi_factors_in_mem(self): n = 100 A = sp.rand(n, n, 0.7)+sp.identity(n) x = np.ones((n, 10)) From 8f5cb6e5ef2e8032393bf4c9a80f6f11f5d81d4d Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Thu, 17 Aug 2023 15:57:50 -0400 Subject: [PATCH 19/23] Migrate away from distutils, which will be removed in Python 3.12 --- setup.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 8fc2c98..977e738 100644 --- a/setup.py +++ b/setup.py @@ -5,8 +5,7 @@ """ -from distutils.core import setup -from setuptools import find_packages +from setuptools import find_packages, setup CLASSIFIERS = [ 'Development Status :: 4 - Beta', From 3d6be06f601f66f0469899aee5ca4ffe49dd8ec6 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Fri, 18 Aug 2023 09:57:15 -0400 Subject: [PATCH 20/23] chmod -x mumps_p.f90 to match other fortran files --- pymatsolver/mumps/mumps_p.f90 | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 pymatsolver/mumps/mumps_p.f90 diff --git a/pymatsolver/mumps/mumps_p.f90 b/pymatsolver/mumps/mumps_p.f90 old mode 100755 new mode 100644 From 483660403fbae9cbc4725e4c4f02ed7610730eb6 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Fri, 18 Aug 2023 10:09:58 -0400 Subject: [PATCH 21/23] Make coverage like in CI --- .gitignore | 1 + Makefile | 4 +--- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 8538e31..4f76999 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ *.pyc .coverage +coverage.xml *.so build/ dist/ diff --git a/Makefile b/Makefile index fdaa6c0..88a3dd5 100644 --- a/Makefile +++ b/Makefile @@ -4,9 +4,7 @@ build: python setup.py build_ext --inplace coverage: - coverage run --source pymatsolver -m pytest - coverage report -m - coverage html + pytest --cov-config=.coveragerc --cov-report=xml --cov=pymatsolver -s -v lint: pylint --output-format=html pymatsolver > pylint.html From c848dd8427f743f8d2fd9290ddd1c0ec4f008ac8 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Fri, 18 Aug 2023 14:29:51 -0400 Subject: [PATCH 22/23] Fix transpose issue identified by jcapriot --- pymatsolver/mumps/_init_.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 98a40f8..7e29537 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -146,6 +146,10 @@ class Mumps(Base): @property def T(self): + """Transpose operator. + + Allows solving A^T * x = b without needing to factor again. + """ properties_with_setters = [] for a in dir(self.__class__): attr = getattr(self.__class__, a) @@ -154,9 +158,11 @@ def T(self): kwargs = {attr: getattr(self, attr) for attr in properties_with_setters} newMS = self.__class__( - self.A.T, + self.A, + from_pointer=self.pointer, **kwargs, ) + newMS.transpose = not self.transpose return newMS def __init__(self, A, from_pointer=None, **kwargs): From cde8692e69f9b9679a42a2288e4abe375b259da4 Mon Sep 17 00:00:00 2001 From: Matthew Plough Date: Mon, 21 Aug 2023 10:18:20 -0400 Subject: [PATCH 23/23] Use same conjugate API as scipy.sparse --- pymatsolver/mumps/_init_.py | 34 +++++++++++++++++++++++++------- pymatsolver/solvers.py | 12 ++++++++++++ tests/test_Mumps.py | 39 +++++++++++++++++++++++++++++++++++++ 3 files changed, 78 insertions(+), 7 deletions(-) diff --git a/pymatsolver/mumps/_init_.py b/pymatsolver/mumps/_init_.py index 7e29537..524f888 100644 --- a/pymatsolver/mumps/_init_.py +++ b/pymatsolver/mumps/_init_.py @@ -5,6 +5,9 @@ import gc import warnings + +import numpy as np + from pymatsolver.solvers import Base from . import MumpsInterface as _MUMPSINT @@ -142,7 +145,8 @@ class Mumps(Base): """ - transpose = False + _transpose = False + _conjugate = False @property def T(self): @@ -150,6 +154,16 @@ def T(self): Allows solving A^T * x = b without needing to factor again. """ + return self._make_copy(transpose=True) + + def conjugate(self): + """Conjugate operator. + + Allows solving \\bar(A) * x = b without needing to factor again. + """ + return self._make_copy(conjugate=True) + + def _make_copy(self, *, transpose=False, conjugate=False): properties_with_setters = [] for a in dir(self.__class__): attr = getattr(self.__class__, a) @@ -157,13 +171,14 @@ def T(self): properties_with_setters.append(a) kwargs = {attr: getattr(self, attr) for attr in properties_with_setters} - newMS = self.__class__( + copy = self.__class__( self.A, from_pointer=self.pointer, **kwargs, ) - newMS.transpose = not self.transpose - return newMS + copy._transpose = (not self._transpose) if transpose else self._transpose + copy._conjugate = (not self._conjugate) if conjugate else self._conjugate + return copy def __init__(self, A, from_pointer=None, **kwargs): self.A = A.tocsc() @@ -226,9 +241,14 @@ def _solveM(self, rhs): rhs = rhs.flatten(order='F') n = self.A.shape[0] nrhs = rhs.size // n - T = 1 if self.transpose else 0 - sol = self._funhandle('solve')(self.pointer.INT, nrhs, rhs, T) - return sol + T = 1 if self._transpose else 0 + sol = self._funhandle('solve')( + self.pointer.INT, + nrhs, + np.conjugate(rhs) if self._conjugate else rhs, + T + ) + return np.conjugate(sol) if self._conjugate else sol _solve1 = _solveM diff --git a/pymatsolver/solvers.py b/pymatsolver/solvers.py index 0bf7111..1111300 100644 --- a/pymatsolver/solvers.py +++ b/pymatsolver/solvers.py @@ -65,6 +65,18 @@ def T(self): newS = self._transposeClass(self.A.T) return newS + @property + def conjugate(self): + "The element-wise complex conjugate for this class" + if self._transposeClass is None: + raise Exception( + 'The complex conjugate for the {} class is not possible.'.format( + self.__name__ + ) + ) + newS = self._transposeClass(self.A.conjugate()) + return newS + def _compute_accuracy(self, rhs, x): nrm = np.linalg.norm(np.ravel(self.A*x - rhs), np.inf) nrm_rhs = np.linalg.norm(np.ravel(rhs), np.inf) diff --git a/tests/test_Mumps.py b/tests/test_Mumps.py index 8fa19bc..b3c505a 100644 --- a/tests/test_Mumps.py +++ b/tests/test_Mumps.py @@ -143,6 +143,45 @@ def test_T(self): assert np.linalg.norm(x[:, i] - sol[:, i]) < TOL assert np.linalg.norm(x - sol, np.inf) < TOL + def test_T_congruent(self): + rhs = self.rhs + sol = self.sol + + AinvT1 = Mumps(self.A, is_symmetric=True).T + AinvT2 = Mumps(self.A.T, is_symmetric=True) + x1 = AinvT2 * rhs + x2 = AinvT1 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + + def test_conjugate_congruent(self): + rhs = self.rhs + sol = self.sol + + Ainv_conj1 = Mumps(self.A).conjugate() + Ainv_conj2 = Mumps(self.A.conjugate()) + x1 = Ainv_conj1 * rhs + x2 = Ainv_conj2 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + + def test_complex_conjugate_congruent(self): + rhs = self.rhs + sol = self.sol + + Ainv_conj1 = Mumps(self.A).T.conjugate() + Ainv_conj2 = Mumps(self.A.T.conjugate()) + x1 = Ainv_conj1 * rhs + x2 = Ainv_conj2 * rhs + + for i in range(3): + assert np.linalg.norm(x1[:, i] - x2[:, i]) < TOL + assert np.linalg.norm(x1 - x2, np.inf) < TOL + class TestMumps1to5: