Add SuperLU_MT support

[ChrisRackauckas]

It would be nice to have the option of using SuperLUMT with the sparse matrix support, since KLU is for quite weird circumstances. We'd first have to update the artifact to build with it, then add a branch to our linear solver support for it.

[ViralBShah]

@jd-lara Do you know if SuperLU is the only supported direct sparse solver for Sundials. If we could use UMFPACK, it may be trivial to enable.

[ChrisRackauckas]

KLU and SuperLU are the only two that are supported, other than things like SuperLUDist which require a different NVector type.

[jd-lara]

I had a conversation with @balos1 a few months ago about the implementation of custom linear solvers with the intention of implementing pardiso. It is not too bad to write the Sundials Wrapper for custom linear solvers.

[ChrisRackauckas]

We could use that so that way it works with the DiffEq linear solver support, and that would then allow UMFPACK.

Ref: https://docs.sciml.ai/latest/features/linear_nonlinear/

@YingboMa is currently improving this part of DiffEq a lot, so it might be good to have him in the loop here.

[ViralBShah]

We already have MKL Pardiso in BB. I am not sure if SuperLU MT receives bugfixes. How important is this feature to enable in Sundials?

[ViralBShah]

@ChrisRackauckas Are you suggesting that have Sundials call out to the sciml linsolve, and then use that to pull in any sparse direct solver?

[ChrisRackauckas]

Yes. Otherwise SuperLU is the only option left in the list that Sundials allows. So if we plan on linking in more, we can just put it to the SciML linsolve interface, and then any Julia code works there and we can call it done.

[jd-lara]

looks like there is already a wrapper for SuperLUMT
https://github.com/LLNL/sundials/tree/master/src/sunlinsol/superlumt

Building with SuperLU MT
The SuperLU MT libraries are available for download from the Lawrence Berkeley National Labo- ratory website: http://crd-legacy.lbl.gov/∼xiaoye/SuperLU/#superlu mt. sundials has been tested with SuperLU MT version 3.1. To enable SuperLU MT, set SUPERLUMT ENABLE to ON, set SUPERLUMT INCLUDE DIR to the SRC path of the SuperLU MT installation, and set the variable SUPERLUMT LIBRARY DIR to the lib path of the SuperLU MT installation. At the same time, the vari- able SUPERLUMT LIBRARIES must be set to a semi-colon separated list of other libraries SuperLU MT depends on. For example, if SuperLU MT ws build with an external blas library, then include the full path to the blas library in this list. Additionally, the variable SUPERLUMT THREAD TYPE must be set to either Pthread or OpenMP.
Do not mix thread types when building sundials solvers. If threading is enabled for sundials by having either OPENMP ENABLE or PTHREAD ENABLE set to ON then SuperLU MT should be set to use the same threading type.

[jd-lara]

Yes. Otherwise SuperLU is the only option left in the list that Sundials allows. So if we plan on linking in more, we can just put it to the SciML linsolve interface, and then any Julia code works there and we can call it done.

In order to achieve this, these methods need to be implemented
https://github.com/LLNL/sundials/blob/master/src/sunlinsol/klu/fsunlinsol_klu.h

[ChrisRackauckas]

Yes.

This is issue is that, of all of the ones that don't require a new special array type, SuperLU MT is the last one in the list to do. So we probably should do it, or just add our SciML linear operator interface to this package, in which case we'd get UMFPACK support for free from Julia Base, + MKLPardisio + etc.

[ViralBShah]

Adding SuperLU MT is much lesser effort, I suspect. Do we know if it links to 64-bit BLAS and if that combination works with Sundials? If not, no problem, because we now also have 32-bit BLAS in BB.

[jd-lara]

Adding SuperLU MT is much lesser effort, I suspect. Do we know if it links to 64-bit BLAS and if that combination works with Sundials? If not, no problem, because we now also have 32-bit BLAS in BB.

From the documentation, it seems that it depends on how SUPERLUMT is built.

At the same time, the variable SUPERLUMT LIBRARIES must be set to a semi-colon separated list of other libraries SuperLU MT depends on. For example, if SuperLU MT was build with an external blas library, then include the full path to the blas library in this list. Additionally, the variable SUPERLUMT THREAD TYPE must be set to either Pthread or OpenMP.

[ViralBShah]

@christopher-dG It seems you might be looking into this. Just checking, so that we don't end up doing the same things. It would be great to have more help maintaining this whole ecosystem of scientific software.

[christopher-dG]

Yeah I was looking into building a SuperLU_MT_jll to then be used here.

[spraharsh]

Any updates on using Julia linear solvers within sundials?
I'm also trying to figure out whether there's code anywhere that gets Pardiso to work with sundials.

[ViralBShah]

Aren't you better off just using Julia's DiffEq at that point?

[jd-lara]

Pardiso isn't compatible with the the Sundials C code at the moment. Currently we the only linear solver not supported is SuperLUMT.

[ChrisRackauckas]

Yes, Sundials has a pretty strict set of linear solvers. If you want to use Pardiso you can use it with OrdinaryDiffEq as documented on https://diffeq.sciml.ai/stable/features/linear_nonlinear/#Pre-Built-Linear-Solver-Choices . Given the performance and generality advantages of OrdinaryDiffEq, it makes a lot more sense to just do that anyways.

[jd-lara]

I will take another swing at this during the summer to improve benchmarking against some of the DiffEq solvers for stiff systems. But for now, the addition of interfaces for more linear solvers is very time-consuming and it seems that there isn't an immediate need to support more linear solvers.

[spraharsh]

I tried solvers from DiffEq, the costliest steps for my use case were LU factorizations followed by function/Jacobian calls, both of which CVODE seemed to reduce the number of, compared to the DiffEq solvers.
so it seems the best option for me is to use a performant sparse solver within CVODE (if possible a Cholesky solver since my jacobian is symmetric too).

[ChrisRackauckas]

I tried solvers from DiffEq, the costliest steps for my use case were LU factorizations followed by function/Jacobian calls, both of which CVODE seemed to reduce the number of, compared to the DiffEq solvers.

With TRBDF2 and KenCarp47? Can you share the example?

so it seems the best option for me is to use a performant sparse solver within CVODE (if possible a Cholesky solver since my jacobian is symmetric too).

LinSolveFactorize(cholesky!) should do it for the pure Julia solvers.

[ChrisRackauckas]

BTW, the latest release post shows QNDF outperforming CVODE_BDF across the board, which renders this issue somewhat moot: https://sciml.ai/news/2021/05/24/QNDF/

[spraharsh]

Thanks! will check it out. (was gonna write up an example independent of the libraries we use next weekend)