1
$\begingroup$

Does anyone know a library that implements sparse direct solver working in already open omp parallel region? The only library that I know that works with this requirement is Pardiso7.2 worth ~8K USD for a single-user license.

Intel MKL Pardiso and MUMPS solvers assume the execution environment is single-threaded, so when calling the interface inside parallel programs it requires protection with omp single region. Activating nested parallelism does not help in this context and iterative solvers are not useful because of the high cost of synchronization (norms and dot products computation) for a kernel that is going to be called in the orders of ~1M times.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

At the risk of self-promotion, I can point you to MyraMath, which is a package I authored in this space out of my own/similar frustrations. [Since you mention cost as a driver, MyraMath is open sourced under GPL/copyleft but closed source licenses require a fee].

I am assuming your intent is to invoke multiple sparse factorizations simultaneously, but each one is perhaps not big enough to "throw the whole machine at it"? (ie something like N problem instances with P cores applied to each, more or less?). MyraMath does support this kind of workflow, the basic idea is that all parallel operations in MyraMath are accomplished using task graphs (think TBB's flowgraph) and MyraMath provides an an API to "compose" multiple taskgraphs into one, then solve that "composed" problem in parallel. (See this tutorial to get a bit more elaboration).

All that said, this is not a workflow that gets a lot exercise, so if you do decide to pursue it, I would encourage thorough testing/profiling to make sure you are actually coming out ahead by organizing the calculation like this.

PS: Mine is certainly not the only game in town, you might also search for Clique/Elemental which might (?) expose similar functionality (IIRC it's a task graph based code, too). In fact, I must mention that I regularly use MKL's Pardiso from multiple threads (with no locking mechanisms in my application) with no ill effects. However I would have to double check if my application runs those instances with multiple (MKL) threads .. of that I am not sure. Either way, I would encourage you to take a deeper look there.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Thanks for your answer! Yes I intend to use it to solve one problem at the time with nThreads already running, so the interface is expected to receive shared memory regions and support execution call inside a omp parallel region. I will test MyraMath shortly I expect it fit well inside our problems. Could you clarify if MyraMath support to call the symbolic, numeric and solve phase within this constrain? PS. MKL Pardiso is based on Pardiso6.x which has a student free license. Unfortunately it does not support calls inside already open parallel regions. $\endgroup$
    – Cesar Conopoima
    Jan 28, 2022 at 13:37
  • $\begingroup$ The reordering routines in MyraMath are pretty anemic, I would steer you towards using METIS (or something similar) to compute a fill-reducing permutation a-priori, then provide it along with your matrix. Even with a permutation, there is still a symbolic phase in MyraMath, it's sequential but reasonably quick and should be thread safe. The solve phase is similar to factorization, in that you can bake it into a task graph and perform it alongside others. That said, backsolution involves a lot of data movement and (in my experience) doesn't show great gains when run in parallel. But YMMV. $\endgroup$
    – rchilton1980
    Jan 28, 2022 at 15:35
  • $\begingroup$ I would try suitesparse, which has LLT, LU, and QR direct sparse solvers. The degree of parallelism is a mixed bag: LLT has CUDA functionality, and the QR uses TBB. Even in single thread applications, I have found these to be faster than some parallel solvers. It also comes packaged with COLAMD and METIS. For many applications, the quality of the re-ordering is just as important as the implementation of the solver. people.engr.tamu.edu/davis/suitesparse.html $\endgroup$
    – Charlie S
    Jan 28, 2022 at 16:01
  • $\begingroup$ @Charlie S: Probably better off as a distinct answer. $\endgroup$
    – rchilton1980
    Jan 28, 2022 at 16:16
  • $\begingroup$ @rchilton1980 I understand. Pardiso also uses metis to generate fill-reducing permutation. Is ok to perform the symbolic setup single threaded because it is the less called kernel on the simulation process. Numeric setup and solve phases are critical because the matrix entries change continuously in a time dependent form. Any performance improvement would be welcomed even if the scalability is limited. $\endgroup$
    – Cesar Conopoima
    Jan 28, 2022 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.