Mumps solver updates #39
Conversation
The Mumps test_singular test is available now. The Mumps test_1to5_T test is still failing.
Also correctly handle warnings as flags.
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.
Otherwise you get
ValueError: dimension mismatch
when the solver doesn't do the reshaping
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.
Codecov Report
@@ Coverage Diff @@
## main #39 +/- ##
==========================================
- Coverage 89.26% 87.78% -1.49%
==========================================
Files 5 5
Lines 298 311 +13
==========================================
+ Hits 266 273 +7
- Misses 32 38 +6
|
pymatsolver/mumps/_init_.py
Outdated
| 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__( |
There was a problem hiding this comment.
I think the point before was that mumps could use a factorization of A to solve systems with A.T instead of refactoring the entire matrix, as like what will happen here (This might not be the case though.)
There was a problem hiding this comment.
Nice catch - and yes, that's exactly the point here. Fixed in c848dd8.
I must have introduced this while debugging why the transpose was resulting in a Hermitian transpose for complex matrices (which is not how the pardiso or the rest of the solvers handle this case).
The following are now equivalent but AinvT1 doesn't need a refactorization:
AinvT1 = Mumps(A).T
AinvT2 = Mumps(A.T)This behavior matches all the other solvers.
I contemplated adding a .H (Hermitian transpose) operator such that:
Mumps(A).H == Mumps(A.T.conjugate())However, this presents a floating-point +0 / -0 problem that I'm not thrilled about. Given the matrix
[ 1 + 0j, 0 + 0j ]
[ 0 + 0j, 1 + 1j ]
the inverse is
[ 1 + 0j, 0 + 0j ]
[ 0 + 0j, 0.5 - 0.5j ]
Since Mumps doesn't take the Hermitian transpose itself, we'd need to emulate that behavior by conjugating the RHS, telling Mumps to solve the transpose problem, and then conjugating again.
The inverse of the matrix's conjugate transpose (i.e. Mumps(A.T.conjugate()) is
[ 1 + 0j, 0 + 0j ]
[ 0 + 0j, 0.5 + 0.5j ]
as computed by Mumps, scipy, or pardiso.
However, the inverse as computed with the emulation (i.e. Mumps(A).H) is
[ 1 - 0j, 0 - 0j ]
[ 0 - 0j, 0.5 + 0.5j ]
That's a little iffy for me. I'd rather have the solvers all behave the same way on simple cases.
There was a problem hiding this comment.
Honestly, I'm not concerned about the differences between +0 and -0 (A floating point thing), I don't think we're ever doing anything that divides by 0 which would be the only place that makes a difference.
There was a problem hiding this comment.
That has me convinced.
I went away from the .H property because the scipy.sparse API doesn't have that property and I wanted to present the same interface here. To take the Hermitian transpose of a scipy.sparse matrix A, you run A.T.conjugate(), so I added a conjugate operator here. I also added unit tests to ensure that Mumps(A).T.conjugate() gives the same results as Mumps(A.T.conjugate()).
|
What steps are you using to link this to mumps?
There are a few python interfaces already out there,
|
At the moment I'm compiling MUMPS from source using https://github.com/scivision/mumps, which lets me use the latest version of MUMPS (5.6.1) rather than the last version available on conda-forge (5.2.1). The conda-forge feedstock was last updated a year and 5 months ago. I am not using the OpenMPI-parallel version of MUMPS, relying instead on OpenMP and BLAS parallelism. Aside from a migration to With MUMPS 5.6.1 built and installed, a
I agree that a Cython interface is the way to go here. A Cython interface like pydiso's will be far easier to support than the current Fortran / I'd like to avoid completely separating the MUMPS Cython interface from pymatsolver, though. First, we can still package pymatsolver with a MUMPS Cython extension and have installation succeed even without MUMPS installed. Cython will make it very easy to package the extension with setuptools and the setuptools Second, our current MUMPS interface needs are very small in scope. It would be fairly straightforward to write a small extension that supports only the functionality we need from MUMPS -- just init, a helper to convert to MUMPS format, factor, solve, and destroy. Third, the existing Python MUMPS interfaces are dormant. The last commit to pymumps was on October 27, 2020 and its last release was on November 5, 2018. The last commit to mumpy was on November 14, 2017 and it has never had a release. Ultimately I'd like to release a version of pymatsolver with good MUMPS support, and to do so with a minimum amount of work. Lack of MUMPS support makes it impossible to have a performant SimPEG on Apple's M1 hardware. I'd like to close that gap pragmatically -- and quickly if possible. |
|
For the Apple silicon processors I was experimenting with apple’s sparse solvers in their accelerate framework. It works well for most of our systems because they are SPD, but they do not support complex matrices. https://developer.apple.com/documentation/accelerate/sparse_solvers/ |
|
It's too bad about Apple's solvers. MUMPS does really well on Apple's M1 machines though, and it's also portable to ARM servers like AWS's Graviton processors. |
|
Thanks for the discussion at last night's meeting! I'm going to move this out of draft status. As discussed we'll handle packaging separately. |
mplough-kobold
changed the title
|
Hi all, I'm taking over the development of
Our work is for now over here: https://gitlab.kwant-project.org/kwant/python-mumps, and it is very much in progress. Do you have any feedback? Is that project relevant to this package? On a separate note, I made a proposal to scipy that is similar in spirit to the main point of this package: scipy/scipy#14485. |
|
An update: python-mumps is now released on pypi and conda-forge on all platforms. |
|
@akhmerov - the python-mumps project is definitely relevant to this package. It's not the focus of my attention right now so please don't interpret my delay in responding as a lack of interest in your work. Why python-mumps is relevantThe SimPEG geophysics simulation package requires a performant sparse solver to complete simulations in a reasonable amount of time. This package provides an abstraction layer and an interface for sparse solvers so that SimPEG doesn't have to understand the specifics of every solver. The best direct solvers at this point are Intel's MKL Pardiso and MUMPS. MKL Pardiso only runs on x86-64 processors while MUMPS can also run on ARM. Performance is fairly similar. The pymatsolver package brings in MKL Pardiso support via pydiso. It currently has a very small interface to MUMPS; this PR makes that interface good enough for our purposes. However, this PR has really just served as a conversation starter; the approach is a bit of a dead end. In the future, we'd like pymatsolver to be able to use a separate interface package to talk to MUMPS in much the same manner that it uses pydiso to talk to MKL Pardiso. Your attention here is exactly what we need - thanks for making us aware of your work. Some immediate needsGetting rid of We'd also like the solver interface to be able to release the GIL when doing computations. SimPEG has some level of Dask support and we have users parallelizing using Dask right now. The inability to release the GIL means that workers time out, making simulation runs less reliable. Some attempts have been made to release the GIL in pydiso but I don't believe that the work is complete. python-mumps gapsLooking in PyPI, the python-mumps package as of 0.0.1.post1 only offers binaries for linux / x86-64. Is there a plan to offer binaries for linux / ARM? There's not really a need at this point for SimPEG to use MUMPS on x86-64 - if someone is running on that platform then the existing interface to MKL Pardiso will work fine. |
|
Thanks for the explanation.
I prioritize conda-forge over pypi because their build machinery is more standardized. As far as I understand, it arm should already be available there. |
|
Currently, |
|
Wow, that would be awesome! |
This PR contains updates that I have made to the Mumps solver.
Mumps changes:
.conjugate()method, allowing for taking the conjugate transpose viaSolver(A).T.conjugate(). The Mumps solver now supports solving the conjugate transpose problem (A^H)x = b without refactorization viaMumps(A).T.conjugate(). Unit tests are included ensure this case works.Maintenance: