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.