diff --git a/Project.toml b/Project.toml
index 834a6b246..5056f3ff1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -96,7 +96,7 @@ MPI = "0.20"
 Markdown = "1.10"
 Metal = "1"
 MultiFloats = "1"
-Pardiso = "0.5"
+Pardiso = "0.5.7"
 Pkg = "1"
 PrecompileTools = "1.2"
 Preferences = "1.4"
diff --git a/ext/LinearSolvePardisoExt.jl b/ext/LinearSolvePardisoExt.jl
index 98db67a09..0b4cfbb11 100644
--- a/ext/LinearSolvePardisoExt.jl
+++ b/ext/LinearSolvePardisoExt.jl
@@ -22,26 +22,44 @@ function LinearSolve.init_cacheval(alg::PardisoJL,
         reltol,
         verbose::Bool,
         assumptions::LinearSolve.OperatorAssumptions)
-    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm = alg
+    @unpack nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor = alg
     A = convert(AbstractMatrix, A)
 
+    if isnothing(vendor)
+        if Pardiso.panua_is_available()
+            vendor = :Panua
+        else
+            vendor = :MKL
+        end
+    end
+
     transposed_iparm = 1
-    solver = if Pardiso.PARDISO_LOADED[]
-        solver = Pardiso.PardisoSolver()
-        Pardiso.pardisoinit(solver)
-        solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
+    solver = if vendor == :MKL
+        solver = if Pardiso.mkl_is_available()
+            solver = Pardiso.MKLPardisoSolver()
+            Pardiso.pardisoinit(solver)
+            nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
 
-        solver
-    else
-        solver = Pardiso.MKLPardisoSolver()
-        Pardiso.pardisoinit(solver)
-        nprocs !== nothing && Pardiso.set_nprocs!(solver, nprocs)
+            # for mkl 1 means conjugated and 2 means transposed.
+            # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
+            transposed_iparm = 2
 
-        # for mkl 1 means conjugated and 2 means transposed.
-        # https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-0/pardiso-iparm-parameter.html#IPARM37
-        transposed_iparm = 2
+            solver
+        else
+            error("MKL Pardiso is not available. On MacOSX, possibly, try Panua Pardiso.")
+        end
+    elseif vendor == :Panua
+        solver = if Pardiso.panua_is_available()
+            solver = Pardiso.PardisoSolver()
+            Pardiso.pardisoinit(solver)
+            solver_type !== nothing && Pardiso.set_solver!(solver, solver_type)
 
-        solver
+            solver
+        else
+            error("Panua Pardiso is not available.")
+        end
+    else
+        error("Pardiso vendor must be either `:MKL` or `:Panua`")
     end
 
     if matrix_type !== nothing
@@ -113,7 +131,6 @@ end
 function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs...)
     @unpack A, b, u = cache
     A = convert(AbstractMatrix, A)
-
     if cache.isfresh
         phase = alg.cache_analysis ? Pardiso.NUM_FACT : Pardiso.ANALYSIS_NUM_FACT
         Pardiso.set_phase!(cache.cacheval, phase)
diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index 5df86c925..fcec925d8 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -253,6 +253,7 @@ export SimpleGMRES
 export HYPREAlgorithm
 export CudaOffloadFactorization
 export MKLPardisoFactorize, MKLPardisoIterate
+export PanuaPardisoFactorize, PanuaPardisoIterate
 export PardisoJL
 export MKLLUFactorization
 export AppleAccelerateLUFactorization
diff --git a/src/extension_algs.jl b/src/extension_algs.jl
index d60f6895b..c5e2db955 100644
--- a/src/extension_algs.jl
+++ b/src/extension_algs.jl
@@ -108,7 +108,7 @@ All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
-MKLPardisoFactorize(; kwargs...) = PardisoJL(; solver_type = 0, kwargs...)
+MKLPardisoFactorize(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 0, kwargs...)
 
 """
 ```julia
@@ -136,19 +136,18 @@ All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
-MKLPardisoIterate(; kwargs...) = PardisoJL(; solver_type = 1, kwargs...)
+MKLPardisoIterate(; kwargs...) = PardisoJL(; vendor = :MKL, solver_type = 1, kwargs...)
 
 """
 ```julia
-PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
-    solver_type = nothing,
+PanuaPardisoFactorize(; nprocs::Union{Int, Nothing} = nothing,
     matrix_type = nothing,
     cache_analysis = false,
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
 ```
 
-A generic method using MKL Pardiso. Specifying `solver_type` is required.
+A sparse factorization method using Panua Pardiso.
 
 !!! note
 
@@ -165,6 +164,59 @@ All values default to `nothing` and the solver internally determines the values
 given the input types, and these keyword arguments are only for overriding the
 default handling process. This should not be required by most users.
 """
+PanuaPardisoFactorize(; kwargs...) = PardisoJL(;
+    vendor = :Panua, solver_type = 0, kwargs...)
+
+"""
+```julia
+PanuaPardisoIterate(; nprocs::Union{Int, Nothing} = nothing,
+    matrix_type = nothing,
+    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing)
+```
+
+A mixed factorization+iterative method using Panua Pardiso.
+
+!!! note
+
+    Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso`
+
+## Keyword Arguments
+
+For the definition of the keyword arguments, see the Pardiso.jl documentation.
+All values default to `nothing` and the solver internally determines the values
+given the input types, and these keyword arguments are only for overriding the
+default handling process. This should not be required by most users.
+"""
+PanuaPardisoIterate(; kwargs...) = PardisoJL(; vendor = :Panua, solver_type = 1, kwargs...)
+
+"""
+```julia
+PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
+    solver_type = nothing,
+    matrix_type = nothing,
+    iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+    vendor::Union{Symbol, Nothing} = nothing
+)
+```
+
+A generic method using  Pardiso. Specifying `solver_type` is required.
+
+!!! note
+
+    Using this solver requires adding the package Pardiso.jl, i.e. `using Pardiso`
+
+## Keyword Arguments
+
+The `vendor` keyword allows to choose between Panua pardiso  (former pardiso-project.org; `vendor=:Panua`)
+and  MKL Pardiso (`vendor=:MKL`). If `vendor==nothing`, Panua pardiso is preferred over MKL Pardiso.
+
+For the definition of the other keyword arguments, see the Pardiso.jl documentation.
+All values default to `nothing` and the solver internally determines the values
+given the input types, and these keyword arguments are only for overriding the
+default handling process. This should not be required by most users.
+"""
 struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
     nprocs::Union{Int, Nothing}
     solver_type::T1
@@ -172,13 +224,15 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
     cache_analysis::Bool
     iparm::Union{Vector{Tuple{Int, Int}}, Nothing}
     dparm::Union{Vector{Tuple{Int, Int}}, Nothing}
+    vendor::Union{Symbol, Nothing}
 
     function PardisoJL(; nprocs::Union{Int, Nothing} = nothing,
             solver_type = nothing,
             matrix_type = nothing,
             cache_analysis = false,
             iparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
-            dparm::Union{Vector{Tuple{Int, Float64}}, Nothing} = nothing)
+            dparm::Union{Vector{Tuple{Int, Int}}, Nothing} = nothing,
+            vendor::Union{Symbol, Nothing} = nothing)
         ext = Base.get_extension(@__MODULE__, :LinearSolvePardisoExt)
         if ext === nothing
             error("PardisoJL requires that Pardiso is loaded, i.e. `using Pardiso`")
@@ -188,7 +242,7 @@ struct PardisoJL{T1, T2} <: LinearSolve.SciMLLinearSolveAlgorithm
             @assert T1 <: Union{Int, Nothing, ext.Pardiso.Solver}
             @assert T2 <: Union{Int, Nothing, ext.Pardiso.MatrixType}
             return new{T1, T2}(
-                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm)
+                nprocs, solver_type, matrix_type, cache_analysis, iparm, dparm, vendor)
         end
     end
 end
diff --git a/test/pardiso/pardiso.jl b/test/pardiso/pardiso.jl
index f05cca53b..d244c4e33 100644
--- a/test/pardiso/pardiso.jl
+++ b/test/pardiso/pardiso.jl
@@ -14,20 +14,38 @@ e = ones(n)
 e2 = ones(n - 1)
 A2 = spdiagm(-1 => im * e2, 0 => lambda * e, 1 => -im * e2)
 b2 = rand(n) + im * zeros(n)
-prob2 = LinearProblem(A2, b2)
-
 cache_kwargs = (; abstol = 1e-8, reltol = 1e-8, maxiter = 30)
 
-for alg in (PardisoJL(), MKLPardisoFactorize(), MKLPardisoIterate())
-     u = solve(prob1, alg; cache_kwargs...).u
-     @test A1 * u ≈ b1
+prob2 = LinearProblem(A2, b2)
+
+algs=LinearSolve.SciMLLinearSolveAlgorithm[PardisoJL()]
+solvers=Pardiso.AbstractPardisoSolver[]
+extended_algs=LinearSolve.SciMLLinearSolveAlgorithm[PardisoJL()]
 
-     u = solve(prob2, alg; cache_kwargs...).u
-     @test eltype(u) <: Complex
-     @test A2 * u ≈ b2
+if Pardiso.mkl_is_available()
+    push!(algs,MKLPardisoFactorize())
+    push!(solvers,Pardiso.MKLPardisoSolver())
+    extended_algs=vcat(extended_algs,[MKLPardisoFactorize(), MKLPardisoIterate()])
+    @info "Testing MKL Pardiso"
+end
+    
+if Pardiso.panua_is_available()
+    push!(algs,PanuaPardisoFactorize())
+    push!(solvers,Pardiso.PardisoSolver())
+    extended_algs=vcat(extended_algs,[PanuaPardisoFactorize(), PanuaPardisoIterate()])
+    @info "Testing Panua Pardiso"
+end    
+
+
+for alg in extended_algs
+    u = solve(prob1, alg; cache_kwargs...).u
+    @test A1 * u ≈ b1
+
+    u = solve(prob2, alg; cache_kwargs...).u
+    @test eltype(u) <: Complex
+    @test A2 * u ≈ b2
 end
 
-return
 
 
 Random.seed!(10)
@@ -37,7 +55,7 @@ b1 = rand(n);
 b2 = rand(n);
 prob = LinearProblem(copy(A), copy(b1))
 
-prob = LinearProblem(copy(A), copy(b1))
+
 linsolve = init(prob, UMFPACKFactorization())
 sol11 = solve(linsolve)
 linsolve = LinearSolve.set_b(sol11.cache, copy(b2))
@@ -45,16 +63,19 @@ sol12 = solve(linsolve)
 linsolve = LinearSolve.set_A(sol12.cache, copy(A2))
 sol13 = solve(linsolve)
 
-linsolve = init(prob, MKLPardisoFactorize())
-sol31 = solve(linsolve)
-linsolve = LinearSolve.set_b(sol31.cache, copy(b2))
-sol32 = solve(linsolve)
-linsolve = LinearSolve.set_A(sol32.cache, copy(A2))
-sol33 = solve(linsolve)
 
-@test sol11.u ≈ sol31.u
-@test sol12.u ≈ sol32.u
-@test sol13.u ≈ sol33.u
+
+for alg in algs
+    linsolve = init(prob, alg)
+    sol31 = solve(linsolve)
+    linsolve = LinearSolve.set_b(sol31.cache, copy(b2))
+    sol32 = solve(linsolve)
+    linsolve = LinearSolve.set_A(sol32.cache, copy(A2))
+    sol33 = solve(linsolve)
+    @test sol11.u ≈ sol31.u
+    @test sol12.u ≈ sol32.u
+    @test sol13.u ≈ sol33.u
+end
 
 
 # Test for problem from #497
@@ -67,87 +88,92 @@ function makeA()
     return(A)
 end
 
-A=makeA()
-u0=fill(0.1,size(A,2))
-linprob = LinearProblem(A, A*u0)
-u = LinearSolve.solve(linprob, PardisoJL())
-@test norm(u-u0) < 1.0e-14
+
+for alg in algs
+    A=makeA()
+    u0=fill(0.1,size(A,2))
+    linprob = LinearProblem(A, A*u0)
+    u = LinearSolve.solve(linprob, alg)
+    @test norm(u-u0) < 1.0e-14
+end
 
 
 
-# Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver
-solver = Pardiso.MKLPardisoSolver()
-iparm = [
-    (1, 1),
-    (2, 2),
-    (3, 0),
-    (4, 0),
-    (5, 0),
-    (6, 0),
-    (7, 0),
-    (8, 20),
-    (9, 0),
-    (10, 13),
-    (11, 1),
-    (12, 1),
-    (13, 1),
-    (14, 0),
-    (15, 0),
-    (16, 0),
-    (17, 0),
-    (18, -1),
-    (19, -1),
-    (20, 0),
-    (21, 0),
-    (22, 0),
-    (23, 0),
-    (24, 10),
-    (25, 0),
-    (26, 0),
-    (27, 1),
-    (28, 0),
-    (29, 0),
-    (30, 0),
-    (31, 0),
-    (32, 0),
-    (33, 0),
-    (34, 0),
-    (35, 0),
-    (36, 0),
-    (37, 0),
-    (38, 0),
-    (39, 0),
-    (40, 0),
-    (41, 0),
-    (42, 0),
-    (43, 0),
-    (44, 0),
-    (45, 0),
-    (46, 0),
-    (47, 0),
-    (48, 0),
-    (49, 0),
-    (50, 0),
-    (51, 0),
-    (52, 0),
-    (53, 0),
-    (54, 0),
-    (55, 0),
-    (56, 0),
-    (57, 0),
-    (58, 0),
-    (59, 0),
-    (60, 0),
-    (61, 0),
-    (62, 0),
-    (63, 0),
-    (64, 0)
-]
-
-for i in iparm
-    Pardiso.set_iparm!(solver, i...)
-end
 
-for i in Base.OneTo(length(iparm))
-    @test Pardiso.get_iparm(solver, i) == iparm[i][2]
+# Testing and demonstrating Pardiso.set_iparm! for MKLPardisoSolver
+for solver in solvers
+    iparm = [
+        (1, 1),
+        (2, 2),
+        (3, 0),
+        (4, 0),
+        (5, 0),
+        (6, 0),
+        (7, 0),
+        (8, 20),
+        (9, 0),
+        (10, 13),
+        (11, 1),
+        (12, 1),
+        (13, 1),
+        (14, 0),
+        (15, 0),
+        (16, 0),
+        (17, 0),
+        (18, -1),
+        (19, -1),
+        (20, 0),
+        (21, 0),
+        (22, 0),
+        (23, 0),
+        (24, 10),
+        (25, 0),
+        (26, 0),
+        (27, 1),
+        (28, 0),
+        (29, 0),
+        (30, 0),
+        (31, 0),
+        (32, 0),
+        (33, 0),
+        (34, 0),
+        (35, 0),
+        (36, 0),
+        (37, 0),
+        (38, 0),
+        (39, 0),
+        (40, 0),
+        (41, 0),
+        (42, 0),
+        (43, 0),
+        (44, 0),
+        (45, 0),
+        (46, 0),
+        (47, 0),
+        (48, 0),
+        (49, 0),
+        (50, 0),
+        (51, 0),
+        (52, 0),
+        (53, 0),
+        (54, 0),
+        (55, 0),
+        (56, 0),
+        (57, 0),
+        (58, 0),
+        (59, 0),
+        (60, 0),
+        (61, 0),
+        (62, 0),
+        (63, 0),
+        (64, 0)
+    ]
+
+    for i in iparm
+        Pardiso.set_iparm!(solver, i...)
+    end
+    
+    for i in Base.OneTo(length(iparm))
+        @test Pardiso.get_iparm(solver, i) == iparm[i][2]
+    end
 end