simpeg/pymatsolver

Solve matrix equations in python.

pymatsolver's Issues

issues with python setup.py install

Hello,

first, thank you for making this package available. I managed to get the tests running in the pymatsolver directory without problems, using the instructions in the README,rst (not using anaconda). There are, however, some issues, when trying to use the python setup.py install command.

1. This does not build/install the MUMPS wrappers. (I could make and copy them manually, but...).
2. The tests directory is installed into .../python/dist-packages directly (on ubuntu, this is site-packages) and not into .../python/dist-packages/pymatsolver.

`is_summetric` with `SolverLU`

Copying over from Slack:

Dieter:

Alfredo Sanchez brought a problem to my attention, and I think it is a bug: It is about the following example: https://docs.simpeg.xyz/content/tutorials/05-dcr/plot_inv_2_dcr2d.html Running it with PARDISO works fine. However, running it with SolverLU results in:

Lindsey:

hmm, yes, I think there is an inconsistency between the kwargs that splu / pardiso expect (the is_symmetric flag in Pardiso changes how the matrix is factored, whereas it is not an option in splu https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html)
This is something we should fix in the pymatsolver layer

Sparse RHS

Mumps should take a sparse vector for the RHS.

https://github.com/JuliaSparse/MUMPS.jl/blob/master/src/MUMPS_4.10.0/SCILAB/examples/sparseRHS_example.sce

Improve warning when not enough memory (output is all zero)

Issue carried over from discussion on Slack with @lheagy, @grosenkj (also pinging @rowanc1).

Using Pardiso (it might or might not affect other solvers as well, I don't know), my output was all zeros (array of correct size, all with 0+0j). No error message printed, everything run smooth as usual. (You only note it when matplotlib complains about plotting it on log-scale "Data has no positive values, and therefore cannot be log-scaled" 😄)

Turns out that there was not enough memory. Making the model smaller resolved it. However, while running the notebook the memory usage never went up, so it didn't occur to me that it could be a memory-related issue. I assume Pardiso checks how much memory is available with the model size, and doesn't even start if there is not sufficient.

The obvious solution might be to check output of Pardiso and pass it along. However, that might not be that straight forward or the output might be cryptic. In that case one could make a check at the end if the result is all 0's, and print a warning (related with some hints, as mkl-installation or mesh-size/memory-requirements).

pymatsolver docs

The pymatsolver docs do not appear to be building? (https://readthedocs.org/pymatsolver). @rowanc1: is this a leftover from moving it to the SimPEG org?

Transfer to the SimPEG organization.

This should allow us to collaborate on things.

@jcapriot @sgkang @bsmithyman

PyLops and spgl1

Would spgl1 be something for pymatsolver? Might be less memory-intensive than PARDISO; we would have to check runtime.

BUG? Pydiso not a requirement

Not sure if this is a bug or by design, but pydiso is not a requirement in the setup and therefore needs a separate install of pydiso to use the Paradiso solver. It would be nice to install pymatsolver and get the Paradiso solver with it.

PARDISO compatibility issue with MKL 2021

It seems that the current wrapper for PARDISO doesn't work for MKL 2021.

instantiate solvers without setting a system matrix?

Motivation

In developing a Simulation class in SimPEG (simpeg/simpeg#672), I want to be able to readily serialize and deserialize an instance of a Simulation, so we are injecting properties throughout, which gives us serialization, deserialization and validation (thanks @fwkoch!!). Right now, the solvers are attached as the simulation as a class, and the solver_opts are attached to the simulation as a dict. We then instantiate a solver as needed in SimPEG.

State the problem

There are a couple potential snags with this approach

• if your solver_opts are not valid, we won't find out until we try and use the solver
• It is not readily clear how to serialize the solver class (a Class property is one option seequent/properties#163, we could also think about serializing / deserializing the name of the class as a string)

Another approach?

We could instantiate the solver with its options and then call it to create Ainv. This would solve the serialization problem and also allow immediate validation of the solver options on its creation.

import pymatsolver
solver = pymatsolver.Pardiso(**opts)
Ainv = solver(A)

so the base class might look like:

class Base(properties.HasProperties):

    check_accuracy = properties.Bool(
        "check the accuracy of the solve?",
        default = False
    )

    accuracy_tol = properties.Float(
        "tolerance on the accuracy of the solver",
        default=1e-6
    )

    def __init__(self, **kwargs):
        super(Base, self).__init__(**kwargs)

    def __call__(self, A):
        self.A = A.tocsr()

There are a couple impacts of this within SimPEG. The one that comes to mind is:

• right now in the inversion, if a solver is not set for the data misfit, we use the same one as on the simulation, here we would likely want to make a copy of the instance, rather than re-use the same instance (so that factors are not cleared)

tests failing / not running

It looks like we need to go in and update our github workflows for treating as seen in #44

Use newer, non-mkl, version of PARDISO

To follow up on the Slack-discussion and to not forget about it.

For more info: https://www.pardiso-project.org

From the website:

Important: Please note that the Intel MKL version of PARDISO is based on our version from 2006 and that a lot of new features and improvements of PARDISO are not available in the Intel MKL library.

In their examples this means a speed-up factor of anything between 1.7 and 26, see the comparison charts in the feature-list.

MKL bug for reordering.

https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/657934

Error running setup.py

When running setup.py, finish with:
error: file '[...]/pymatsolver/pymatsolver/Mumps/MumpsInterface.so' does not exist

move to SimPEG org?

Since we are moving towards having SimPEG rely on pymatsolver for all of the solvers (see simpeg/simpeg#672), would it make sense to bring this repo over to the SimPEG org? We can make sure the docs are linked, pypi deploys automated, etc.