Wrap

SciML · Dec 22, 2023 · d569ee0 · d569ee0
1 parent 8866d03
commit d569ee0
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 5 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Sundials"
 uuid = "c3572dad-4567-51f8-b174-8c6c989267f4"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "4.22.1"
+version = "4.23.0"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"

diff --git a/src/common_interface/algorithms.jl b/src/common_interface/algorithms.jl
@@ -752,12 +752,16 @@ struct KINSOL{LinearSolver} <: SundialsNonlinearSolveAlgorithm{LinearSolver}
     jac_upper::Int
     jac_lower::Int
     userdata::Any
-end
+    prec_side::Int
+    krylov_dim::Int
+end # TODO: Add globalization options
 Base.@pure function KINSOL(;
     linear_solver = :Dense,
     jac_upper = 0,
     jac_lower = 0,
-    userdata = nothing)
+    userdata = nothing,
+    prec_side = 0,
+    krylov_dim = 0)
     if !(linear_solver in (:None,
         :Diagonal,
         :Dense,
@@ -774,7 +778,9 @@ Base.@pure function KINSOL(;
     end
     KINSOL{linear_solver}(jac_upper,
         jac_lower,
-        userdata)
+        userdata,
+        prec_side,
+        krylov_dim)
 end
 
 method_choice(alg::SundialsODEAlgorithm{Method}) where {Method} = Method

diff --git a/src/simple.jl b/src/simple.jl
@@ -71,12 +71,41 @@ function ___kinsol(f,
     if linear_solver == :Dense
         A = Sundials.SUNDenseMatrix(length(y0), length(y0))
         LS = Sundials.SUNLinSol_Dense(y0, A)
+    elseif linear_solver == :LapackDense
+        A = Sundials.SUNDenseMatrix(length(y0), length(y0))
+        LS = Sundials.SUNLinSol_LapackDense(y0, A)
     elseif linear_solver == :Band
         A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
         LS = Sundials.SUNLinSol_Band(y0, A)
+    elseif linear_solver == :LapackBand
+        A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
+        LS = Sundials.SUNLinSol_LapackBand(y0, A)
+    elseif linear_solver == :Diagonal
+        error("Diagonal linear solver not implemented")
+    elseif linear_solver == :GMRES
+        A = C_NULL
+        LS = Sundials.SUNLinSol_SPGMR(y0, alg.prec_side, alg.krylov_dim)
+    elseif linear_solver == :FGMRES
+        A = C_NULL
+        LS = Sundials.SUNLinSol_SPFGMR(y0, alg.prec_side, alg.krylov_dim)
+    elseif linear_solver == :BCG
+        A = C_NULL
+        LS = Sundials.SUNLinSol_SPBCGS(y0, alg.prec_side, alg.krylov_dim)
+    elseif linear_solver == :PCG
+        A = C_NULL
+        LS = Sundials.SUNLinSol_PCG(y0, alg.prec_side, alg.krylov_dim)
+    elseif linear_solver == :TFQMR
+        A = C_NULL
+        LS = Sundials.SUNLinSol_SPTFQMR(y0, alg.prec_side, alg.krylov_dim)
+    elseif linear_solver == :KLU
+        nnz = length(SparseArrays.nonzeros(prob.f.jac_prototype))
+        A = Sundials.SUNSparseMatrix(length(y0), length(y0), nnz, CSC_MAT)
+        LS = SUNLinSol_KLU(y0, A)
+    else
+        error("Unknown linear solver")
     end
     flag = @checkflag KINSetFuncNormTol(kmem, abstol) true
-    flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A) true
+    flag = @checkflag Sundials.KINSetLinearSolver(kmem, LS, A) true
     flag = @checkflag KINSetUserData(kmem, userfun) true
     ## Solve problem
     scale = ones(length(y0))