Error in MKL Pardiso for single precision matrices

[RayleighLord]

This is related to the following post post.

With the following MWE, if I try to solve a complex linear system with single precision,

using SparseArrays, MKL, Pardiso, LinearSolve

A = sprand(Float32, 100, 100, 0.1) .+ 1im * sprand(Float32, 100, 100, 0.1)
b = rand(Float32, 100)

prob = LinearProblem(A, b)
sol = LinearSolve.solve(prob, MKLPardisoFactorize())
x = sol.u

I get the error

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:37

which suggests that the single precision feature is not exposed.

[j-fu]

According to

https://github.com/JuliaSparse/Pardiso.jl/blob/a645c664fb32c09c9fd6de0700b9d8484484722b/src/Pardiso.jl#L67

it indeed seems that it is not exposed. Do you have a reference to the feature being documented in MKL?

[mipals]

As I noted in the post it does actually look like it is there for MKL (see iparm[27]. For Panua I its iparm[29], but I am pretty sure that for Panua it is only the factorization that is done in single precision, so the matrix still needs to be passed in double precision.

[KristofferC]

https://www.intel.com/content/www/us/en/developer/articles/technical/pardiso-tips.html also has it documented as:

Single and double precision: Starting version 10.2 Intel® oneAPI Math Kernel Library (oneMKL) supports computations in single precision as well as double precision. To switch between single and double precision modes you should use iparm(28). The default mode of iparm(28) is 0, which corresponds to the double precision. If iparm(28) = 1 all internal computations are made in single precision.

Although that talks about iparm(28) instead of 27...

[mipals]

I've linked to the C-library documentation. I guess the difference is just between Fortran and C indexing?

[j-fu]

Or input data type vs internal computation ?

[KristofferC]

I guess the difference is just between Fortran and C indexing?

Makes sense.

[mipals]

Super. I will have a look tomorrow when I have access to a computer with an Intel CPU.

[KristofferC]

I think this is fixed by https://github.com/JuliaSparse/Pardiso.jl/pull/108.