introduce optional iparm etc parameters for pardiso constructors

Author: j-fu

src/factorizations.jl

@@ -117,7 +117,7 @@ LinearAlgebra.ldiv!(fact::AbstractFactorization, v)=ldiv!(fact.fact,v)
 
 macro makefrommatrix(fact)
     return quote
-        $(esc(fact))(A::ExtendableSparseMatrix{Tv,Ti}; kwargs...) where {Tv,Ti} = factorize!($(esc(fact))(;valuetype=Tv, indextype=Ti),A; kwargs...)
+        $(esc(fact))(A::ExtendableSparseMatrix{Tv,Ti}; kwargs...) where {Tv,Ti} = factorize!($(esc(fact))(;valuetype=Tv, indextype=Ti,  kwargs...),A)
         $(esc(fact))(A::SparseMatrixCSC{Tv,Ti}; kwargs...) where {Tv,Ti} = $(esc(fact))(ExtendableSparseMatrix(A); kwargs...)
     end
 end

src/pardiso_lu.jl

@@ -6,31 +6,37 @@ mutable struct PardisoLU{Tv, Ti} <: AbstractPardisoLU{Tv,Ti}
     phash::UInt64
 end
 
-function PardisoLU{Tv,Ti}() where {Tv,Ti}
+function PardisoLU{Tv,Ti}(;iparm=nothing,dparm=nothing,matrixtype=nothing) where {Tv,Ti}
     fact=PardisoLU{Tv,Ti}(nothing,Pardiso.PardisoSolver(),0)
-    set_default_matrixtype!(fact)
+    default_initialize!(fact,iparm,dparm,matrixtype)
 end
 
 """
 ```
-PardisoLU(;valuetype=Float64, indextype=Int64)
-PardisoLU(matrix)
+PardisoLU(;valuetype=Float64, 
+           indextype=Int64,
+           iparm::Vector, 
+           dparm::Vector, 
+           matrixtype=Int)
+
+PardisoLU(matrix; iparm,dparm,matrixtype)
 ```
 
 LU factorization based on pardiso. For using this, you need to issue `using Pardiso`
 and have the pardiso library from  [pardiso-project.org](https://pardiso-project.org) 
 [installed](https://github.com/JuliaSparse/Pardiso.jl#pardiso-60).
 
-For (setting Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme], 
-you can access the `PardisoSolver` e.g. like
+The optional keyword arguments `matrixtype`, `iparm`  and `dparm` are 
+(Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme].
+
+Forsetting them, one can also access the `PardisoSolver` e.g. like
 ```
 using Pardiso
 plu=PardisoLU()
 Pardiso.set_iparm!(plu.ps,5,13.0)
 ```
 """
-PardisoLU(;valuetype::Type=Float64, indextype::Type=Int64)=PardisoLU{valuetype,indextype}()
-
+PardisoLU(;valuetype::Type=Float64, indextype::Type=Int64, kwargs...)=PardisoLU{valuetype,indextype}(;kwargs...)
 
 
 #############################################################################################
@@ -40,39 +46,60 @@ mutable struct MKLPardisoLU{Tv, Ti} <: AbstractPardisoLU{Tv,Ti}
     phash::UInt64
 end
 
-function MKLPardisoLU{Tv,Ti}() where {Tv,Ti}
+function MKLPardisoLU{Tv,Ti}(;iparm=nothing,matrixtype=nothing) where {Tv,Ti}
     fact=MKLPardisoLU{Tv,Ti}(nothing,Pardiso.MKLPardisoSolver(),0)
-    set_default_matrixtype!(fact)
+    default_initialize!(fact, iparm,nothing,matrixtype)
 end
 
 
 """
 ```
-MKLPardisoLU(;valuetype=Float64, indextype=Int64)
-MKLPardisoLU(matrix)
+MKLPardisoLU(;valuetype=Float64, 
+           indextype=Int64,
+           iparm::Vector, 
+           matrixtype=Int)
+
+MKLPardisoLU(matrix; iparm, matrixtype)
 ```
 
 LU factorization based on pardiso. For using this, you need to issue `using Pardiso`.
 This version  uses the early 2000's fork in Intel's MKL library.
 
+The optional keyword arguments `matrixtype` and `iparm`  are  
+(Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme].
 
-For (setting Pardiso internal parameters)[https://github.com/JuliaSparse/Pardiso.jl#readme], 
-you can access the `PardisoSolver` e.g. like
+For setting them you can also access the `PardisoSolver` e.g. like
 ```
 using Pardiso
 plu=MKLPardisoLU()
 Pardiso.set_iparm!(plu.ps,5,13.0)
 ```
 """
-MKLPardisoLU(;valuetype::Type=Float64, indextype::Type=Int64)=MKLPardisoLU{valuetype,indextype}()
+MKLPardisoLU(;valuetype::Type=Float64, indextype::Type=Int64,kwargs...)=MKLPardisoLU{valuetype,indextype}(;kwargs...)
 
 
 ##########################################################################################
-function set_default_matrixtype!(fact::AbstractPardisoLU{Tv,Ti}) where {Tv, Ti}
-    if Tv<:Complex
-        Pardiso.set_matrixtype!(fact.ps,Pardiso.COMPLEX_NONSYM)
+function default_initialize!(fact::AbstractPardisoLU{Tv,Ti}, iparm,dparm,matrixtype) where {Tv, Ti}
+    if !isnothing(matrixtype)
+        my_matrixtype=matrixtype
+    elseif Tv<:Complex
+        my_matrixtype=Pardiso.COMPLEX_NONSYM
     else
-        Pardiso.set_matrixtype!(fact.ps,Pardiso.REAL_NONSYM)
+        my_matrixtype=Pardiso.REAL_NONSYM
+    end
+    
+    Pardiso.set_matrixtype!(fact.ps,Pardiso.REAL_NONSYM)
+
+    if !isnothing(iparm)
+        for i=1:min(length(iparm),length(fact.ps.iparm))
+            Pardiso.set_iparm!(fact.ps,i,iparm[i])
+        end
+    end
+
+    if !isnothing(dparm)
+        for i=1:min(length(dparm),length(fact.ps.dparm))
+            Pardiso.set_dparm!(fact.ps,i,dparm[i])
+        end
     end
     fact
 end
@@ -82,7 +109,6 @@ function update!(lufact::AbstractPardisoLU{Tv,Ti}) where {Tv, Ti}
     flush!(lufact.A)
     Acsc=lufact.A.cscmatrix
     if lufact.phash!=lufact.A.phash 
-        Pardiso.pardisoinit(ps)
         Pardiso.set_phase!(ps, Pardiso.RELEASE_ALL)
         Pardiso.pardiso(ps, Tv[], Acsc, Tv[])
         Pardiso.set_phase!(ps, Pardiso.ANALYSIS_NUM_FACT)

test/runtests.jl

@@ -362,7 +362,7 @@ end
     @test test_lu2(1000,1,1,lufac=CholeskyFactorization())
 end
 
-@testset "pardiso" begin
+@testset "mkl-pardiso" begin
     @test test_lu1(10,10,10,lufac=MKLPardisoLU())
     @test test_lu1(25,40,1,lufac=MKLPardisoLU())
     @test test_lu1(1000,1,1,lufac=MKLPardisoLU())
@@ -373,3 +373,15 @@ end
 end
 
 
+if Pardiso.PARDISO_LOADED[]
+    @testset "pardiso" begin
+        @test test_lu1(10,10,10,lufac=PardisoLU())
+        @test test_lu1(25,40,1,lufac=PardisoLU())
+        @test test_lu1(1000,1,1,lufac=PardisoLU())
+        
+        @test test_lu2(10,10,10,lufac=PardisoLU())
+        @test test_lu2(25,40,1,lufac=PardisoLU())
+        @test test_lu2(1000,1,1,lufac=PardisoLU())
+    end
+end
+