Alternate backend for sparse solver

[lindahua]

Unlike for dense matrices where BLAS & LAPACK are de facto standard, there are multiple competing implementation (with different interfaces) for sparse solvers.

Currently, Julia is using UMFPACK. PARDISO is also widely used (partly because it has been used by MKL for solving sparse systems).

A report (perhaps outdated) here: ftp://ftp.numerical.rl.ac.uk/pub/reports/ghsRAL200505.pdf
shows that PARDISO has superior performance in general.

We may implement PARDISO (and perhaps others) as an alternative solver and provide a way to allow users to set their preferred solver.

[ViralBShah]

PARDISO is great, but it is not open source and hence that rules it out of Base. It would be great to have it as a package.

[ViralBShah]

I wonder if such wishlists should live as issues filed aganist METADATA. Just a thought.

[kmsquire]

I'd suggest that, at least currently, wishlists (and other things) are more likely to be seen by more people here (or the mailing list).

[lindahua]

Actually, what I suggested is not to incorporate PARDISO into the base, but that we should allow user to register different sparse solver globally.

Consider such a scenario, a package contributor Alice writes a data processing package that contains such a function

function process_data(A, x)
    ...
    y = A \ x;   # A is sparse, so it calls underlying sparse solver
    ...
end

Now, Bob is a user of this package, and he has a PARDISO license, and wants to use the PARDISO solver instead of UMFPACK when calling the functions in Alice's package. Then, Bob wishes to write

Base.register_sparse_solver(pardiso_solve) 
process_data(my_A, my_x);

Here, pardiso_solve is just a sparse solver function that conforms to some standardized interface. It is completely ok that this function is provided by an external package (say Pardiso.jl). But, I wish there is a seamless way that people can plugin their favorite solver, such that all functions that does sparse equation solving can change to use that solver.

[mlubin]

This will require some reorganization but it's a good idea.

There are some tricky issues like how to pass parameters. If you're interested in performance enough to switch the sparse solver, you probably also want to pass specific parameters to PARDISO; there are many different parameters to tune the solution method to the particular structure of the matrix. You may also want to separate the symbolic/numerical/solve phases. The interface will need to support all of this.

[lindahua]

What I thought was much easier than this.

A \ x can just be translated into registered_sparse_solver(A, x).

I can just register a closure like x -> pardiso_solve(A, x, what_ever_params...) if I want to bother with tuning the parameters.

[lindahua]

I think there should not be a very huge change to the base. Here is basically all what it needs to make it work:

change the implementation of \ to

function default_sparse_solve(A::SparseMatrixCSC, x::Vector)
     ... use the current implementation here ...
end

registered_sparse_solver = default_sparse_solve

\ ( a::SparseMatrixCSC, x::Vector) = registered_sparse_solver(A, x)

and add a function to register a different solver.

[mlubin]

Seems reasonable to me. @dmbates?

@JeffBezanson: will this cause problems for type inference?

[stevengj]

PaStiX is also a good package. The de-facto standard here is PETSc, however -- if our sparse package used PETSc as its backend, it would automatically gain access to a huge number of (serial and parallel) sparse-direct and iterative solvers, including UMFPACK, PaStiX, PARDISO, etcetera.

[lindahua]

I just checked. it seemed that PETSc is not working with PARDISO though.

[stevengj]

@lindahua, you're right, although PETSc supports many linear solvers in its backend, PARDISO is not one of them for some reason. (Of course, PARDISO is not free/open-source software, but PETSc already supports several proprietary back-ends. But the main priority is to support free/open-source back-ends like UMFPACK and PaSTiX, probably. Probably you could get the data out of PETSc in a format that PARDISO supports, however, with enough work; sparse-matrix formats are pretty standardized.)

[ViralBShah]

At an interface level, I think we should do a little more refactoring of the sparse solver modules and abstract types so that it is easy to use one from a package. That said, I think it is not a good idea that A \ b should transparently switch to Pardiso, when the user says using Pardiso. Every sparse solver package can implement the same interface, and you should be able to say Pardiso.solve() and so on. This makes the change of the default solver obvious to the reader, and the code more maintainable in my opinion.

Completely agree with @stevengj that we should not write any more sparse solver interfaces, and use PetSc so that we can get all of them. At least on debian systems and such, you can then get PetSc and all the other dependent solvers easily, in theory.

[lindahua]

I closed this. All relevant discussions should go to #2645