Skip to content

Commit

Permalink
WIP: KSP
Browse files Browse the repository at this point in the history
  • Loading branch information
Jeremy E Kozdon committed Aug 4, 2021
1 parent 1bd216c commit 3e3e0f9
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 2 deletions.
47 changes: 45 additions & 2 deletions src/ksp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,13 +11,28 @@ Base.@kwdef mutable struct KSP{PetscLib, PetscScalar} <:
P::Union{AbstractMat, Nothing} = nothing
end

"""
KSP(A::AbstractMat, P::AbstractMat{PetscLib} = A; options...)
Create a `KSP` using the matrix `A` and preconditioner construction matrix `P`
with the `options`.
The communicator is obtained from `A` and if it has size `1` then the garbage
collector is set, otherwise the user is responsible for calling
[`destroy`](@ref).
# External Links
$(_doc_external("KSP/KSPCreate"))
$(_doc_external("KSP/KSPSetOperators"))
$(_doc_external("KSP/KSPSetFromOptions"))
"""
function KSP(
A::AbstractMat{PetscLib},
P::AbstractMat{PetscLib} = A;
kwargs...,
options...,
) where {PetscLib}
@assert initialized(PetscLib)
opts = Options(PetscLib; kwargs...)
opts = Options(PetscLib; options...)
PetscScalar = PetscLib.PetscScalar
ksp = KSP{PetscLib, PetscScalar}(opts = opts)
comm = getcomm(A)
Expand Down Expand Up @@ -70,6 +85,34 @@ function solve!(
return x
end

"""
createvecs(ksp::AbstractKSP; nright = 0, nleft = 0)
Create `nright` right and `nleft` left vectors compatible with the `ksp`.
Returned object `V` has `Tuple` members `V.right` and `V.left` containing the
vectors.
# External Links
$(_doc_external("KSP/KSPCreateVecs"))
"""
function createvecs(ksp::AbstractKSP{PetscLib}; nright = 0, nleft = 0) where PetscLib
r_right = [Ref{CVec}() for i = 1:nright]
r_left = [Ref{CVec}() for i = 1:nleft]
LibPETSc.KSPCreateVecs(PetscLib, ksp, nright, r_right, nleft, r_left)

for i = 1:nright
@show r_right[i][]
end
for i = 1:nleft
@show r_left[i][]
end

right = ntuple(i -> VecPtr(PetscLib, r_right[i][], false), nright)
left = ntuple(i -> VecPtr(PetscLib, r_left[i][], false), nleft)

(right = right, left = left)
end

function LinearAlgebra.ldiv!(x::AbstractVec, ksp::AbstractKSP, b::AbstractVec)
solve!(x, ksp, b)
end
Expand Down
5 changes: 5 additions & 0 deletions test/ksp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,11 @@ using LinearAlgebra: mul!
@test x y
end

PETSc.createvecs(ksp)
PETSc.createvecs(ksp; nright = 1)
PETSc.createvecs(ksp; nleft = 2)
PETSc.createvecs(ksp; nright = 3, nleft = 7)

# PETSc.destroy(ksp)
PETSc.destroy(A)
PETSc.destroy(x)
Expand Down

0 comments on commit 3e3e0f9

Please sign in to comment.