How to use Pardiso with single-precision matrices

Specific Domains Numerics
1

I would like to ask if it is possible to run Pardiso from the LinearSolve library using single precision matrices. I am running

using SparseArrays, LinearSolve, MKL, Pardiso
prob = LinearProblem(A, b)
sol = LinearSolve.solve(prob, MKLPardisoFactorize())

and I get

ERROR: MethodError: no method matching pardiso(::MKLPardisoSolver, ::SparseMatrixCSC{ComplexF32, Int64}, ::Vector{ComplexF32})

Closest candidates are:
  pardiso(::Pardiso.AbstractPardisoSolver, ::SparseMatrixCSC{Tv, Ti}, ::StridedVecOrMat{Tv}) where {Ti, Tv<:Union{Float64, ComplexF64}}
   @ Pardiso ~/.julia/packages/Pardiso/z4YUw/src/Pardiso.jl:367
  pardiso(::Pardiso.AbstractPardisoSolver)
   @ Pardiso ~/.julia/packages/Pardiso/z4YUw/src/Pardiso.jl:366
  pardiso(::Pardiso.AbstractPardisoSolver, ::StridedVecOrMat{Tv}, ::SparseMatrixCSC{Tv, Ti}, ::StridedVecOrMat{Tv}) where {Ti, Tv<:Union{Float64, ComplexF64}}
   @ Pardiso ~/.julia/packages/Pardiso/z4YUw/src/Pardiso.jl:341

Stacktrace:
 [1] solve!(cache::LinearSolve.LinearCache{…}, alg::PardisoJL{…}; kwargs::@Kwargs{})
   @ LinearSolvePardisoExt ~/.julia/packages/LinearSolve/WDeMC/ext/LinearSolvePardisoExt.jl:137
 [2] solve!(cache::LinearSolve.LinearCache{…}, alg::PardisoJL{…})
   @ LinearSolvePardisoExt ~/.julia/packages/LinearSolve/WDeMC/ext/LinearSolvePardisoExt.jl:131
 [3] solve!(::LinearSolve.LinearCache{…}; kwargs::@Kwargs{})
   @ LinearSolve ~/.julia/packages/LinearSolve/WDeMC/src/common.jl:269
 [4] solve!(::LinearSolve.LinearCache{…})
   @ LinearSolve ~/.julia/packages/LinearSolve/WDeMC/src/common.jl:268
 [5] solve(::LinearProblem{…}, ::PardisoJL{…}; kwargs::@Kwargs{})
   @ LinearSolve ~/.julia/packages/LinearSolve/WDeMC/src/common.jl:265
 [6] solve(::LinearProblem{…}, ::PardisoJL{…})
   @ LinearSolve ~/.julia/packages/LinearSolve/WDeMC/src/common.jl:263
 [7] top-level scope
   @ ~/Desktop/LinearSystem_GiulioD/forced_response.jl:11
Some type information was truncated. Use `show(err)` to see complete types.
2

In the case of PanuaPardiso you still need to pass in the matrix in double precision even if the factorization is done in single precision. I thought the same was the case for MKLPardiso, but having looked at ‘iparm[27]’ in the oneMKL documentation it seems as it is not the case. As such I have a feeling that the single precision precision feature of MKLPardiso is simply not exposed.

Can you please make an issue with an MWE on Pardiso.jl?

3

Sure, I will make a MWE. By the way, what parameter do you need to pass to PanuaPardiso if you want the factorization to be done in single precision? Can you do it using LinearSolve or do I have to use the API of Pardiso.jl?

4

I do not know how to pass parameters to LinearSolve.jl. In Pardiso.jl you should change iparm by using the set_iparm function. For PanuaPardiso it should be iparm(29) that is of interest to you (see https://panua.ch/static/manual/manual.pdf) - but I am fairly certain that this only changes the factorization to single precision, meaning that the input matrix still has to be given in double precision.

5

In LinearSolve, you should be able to pass the iparm data
in the constructor:

However, this appears to be not well tested. A good MWE might be helpful here.

1 Like
6

Thank you for the answer. But I just discovered that Panua is not free and will not suit my needs. I’ll have to wait if the MKL Pardiso gets fixed for single-precision in Julia, since I tried to “hack” the library locally but I do not have the experience with it to make it work.

Thank you very much for the help, though.