haasad / PyPardiso

Python interface to the Intel MKL Pardiso library to solve large sparse linear systems of equations
BSD 3-Clause "New" or "Revised" License
135 stars 20 forks source link

Pardiso 6.0 support? #10

Closed ma-sadeghi closed 3 years ago

ma-sadeghi commented 5 years ago

Hi @haasad, I wonder if you have any plans for supporting the latest PARDISO solver (v6). On their website, they show it gives a huge performance boost over the one that comes with MKL. I wonder if that's even possible (due to license issues perhaps?).

Cheers, Amin

haasad commented 5 years ago

Hi @ma-sadeghi, I don't have any plans currently of supporting version 6. The advantage of pypardiso is the ease of use; you can conda install pypardiso on all major OSs and it works out of the box. This is unfortunately not possible with the Pardiso Solver 6.0.

It should however be possible to use the wrapper code of pypardiso to call v6. It used to work with v5, so I assume this is still possible. I'm not in academia anymore, so I don't have a copy of the v6, but I would be very interested if you want to give it a try.

In the best case scenario this only requires changing of one line of code: https://github.com/haasad/PyPardisoProject/blob/f666ea4718b32fa1359e5ca94bedac710b09a428/pypardiso/pardiso_wrapper.py#L63-L70

Replace the library in the ctypes.CDLL call corresponding to your operating system with the full path to your copy of the Pardiso v6 library.

If this works and you can confirm this performance boost, I'm absolutely willing to add an option to pypardiso to use Pardiso v6.

Cheers, Adrian

ma-sadeghi commented 4 years ago

Hi @haasad,

Sorry for the very late reply (exactly a full year)! I was very busy finishing up grad school.

Anyways, I tried what you said, but I get this error

>>> from pypardiso import spsolve
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/amin/Code/PyPardisoProject/pypardiso/__init__.py", line 4, in <module>
    from .scipy_aliases import spsolve, factorized
  File "/home/amin/Code/PyPardisoProject/pypardiso/scipy_aliases.py", line 9, in <module>
    pypardiso_solver = PyPardisoSolver()
  File "/home/amin/Code/PyPardisoProject/pypardiso/pardiso_wrapper.py", line 70, in __init__
    self.libmkl = ctypes.CDLL(lib)
  File "/home/amin/anaconda3/envs/pmeal/lib/python3.7/ctypes/__init__.py", line 364, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: /home/amin/Code/PyPardisoProject/libpardiso600-GNU720-X86-64.so: undefined symbol: sgetrf_

I have no clue where this comes from. I'd really appreciate it if you could shed some light.

Many thanks, Amin

PS. If I can get this to work, I'll try to compile a comprehensive performance comparison between Pardiso v6 and MKL Pardiso, so you can decide whether you want to add the option to use v6.

Update: I realized that I was using OpenBLAS, rather than MKL for numpy/scipy. After reinstalling w/ MKL, I got the following error:

>>> import pypardiso
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/amin/Code/PyPardisoProject/pypardiso/__init__.py", line 4, in <module>
    from .scipy_aliases import spsolve, factorized
  File "/home/amin/Code/PyPardisoProject/pypardiso/scipy_aliases.py", line 9, in <module>
    pypardiso_solver = PyPardisoSolver()
  File "/home/amin/Code/PyPardisoProject/pypardiso/pardiso_wrapper.py", line 70, in __init__
    self.libmkl = ctypes.CDLL(lib)
  File "/home/amin/anaconda3/envs/pmeal/lib/python3.7/ctypes/__init__.py", line 364, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: /home/amin/Code/PyPardisoProject/libpardiso600-GNU720-X86-64.so: undefined symbol: _gfortran_st_close
victorminden commented 4 years ago

Hi @ma-sadeghi

I have also been trying to wrap PARDISO 6+ and faced some similar issues as your last errors. I have not quite gotten full results yet, but I wanted to mention for posterity that I think there are a a few things here that must be done still to get PARDISO 6 working here.

First, note that PARDISO 6 has an extra control parameter vector dparm and I don't think you can leave it out, so the argtypes / callsite should be modified accordingly.

Second, I had crazy trouble getting PARDISO 6 to play nicely with the other shared libraries like gfortran (which it seems like you are seeing in your error too). I had to copy some code to explicitly search for and load those as part of initialization:

        # Load the numerical libraries required by PARDISO.
        _shared_libs = ["lapack", "blas", "omp", "gfortran"]
        for lib in _shared_libs:
            # Fetch the proper name of the dependency
            libname = ctypes.util.find_library(lib)
            # Load the dependency and make it globally accessible
            ctypes.CDLL(libname, mode=ctypes.RTLD_GLOBAL)

I don't know if this helps or if you are even still facing these issues, but I got excited at the prospect of calling PARDISO 6 easily from python and saw this issue was open, so I thought I'd chime in.

ma-sadeghi commented 4 years ago

@victorminden Thanks for posting this. After all, I think PARDISO 6 was not worth it. I compared its speed against MKL PARDISO, and it's almost 2X slower. I have no idea how the authors achieved the amazing speedups reported in their website. I also found a YouTube video that reported the same issue (PARDISO 6 being slower than MKL PARDISO). Have you tried it yourself?

victorminden commented 4 years ago

Thanks for the note @ma-sadeghi. have not done much with either PARDISO 6 or MKL PARDISO but was similarly attracted by the crazy nice benchmark results on the website. I will look more critically at those comparisons before investing too much time in reviving this issue / wrapping it anew -- glad to hear your experience, thanks.

haasad commented 3 years ago

Closing this issue, since the apparently version 6.0 doesn't provide any demonstrable speedups over MKL pardiso, at least for the type of problems pypardiso is currently used for (ie. mainly https://github.com/brightway-lca and https://github.com/PMEAL/OpenPNM that I know of). This was already the case for version 5.0 when I initially created this library.