From 89a118423668b2bf72834bf28b524d180f827148 Mon Sep 17 00:00:00 2001 From: Jeremy E Kozdon Date: Fri, 16 Jul 2021 14:16:09 -0700 Subject: [PATCH] Try to get Vec working --- gen/generator.toml | 1 + gen/prologue.jl | 1 + lib/petsc_library.jl | 5 +- src/LibPETSc.jl | 2 + src/PETSc.jl | 2 +- src/init.jl | 8 +-- src/lib.jl | 17 ++++++ src/sys.jl | 60 +++++++------------- src/vec.jl | 131 ++++++++++++++++++++++++++----------------- test/options.jl | 4 +- 10 files changed, 129 insertions(+), 102 deletions(-) diff --git a/gen/generator.toml b/gen/generator.toml index 5427cf43..5f58f8e8 100644 --- a/gen/generator.toml +++ b/gen/generator.toml @@ -14,6 +14,7 @@ printer_blacklist = [ "PetscOptions", "_n_PetscOptions", "PetscViewer", "_p_PetscViewer", "Vec", "_p_Vec", + "PetscObject", "_p_PetscObject", # # Remove types defined in const.jl # diff --git a/gen/prologue.jl b/gen/prologue.jl index 43678086..418af414 100644 --- a/gen/prologue.jl +++ b/gen/prologue.jl @@ -21,3 +21,4 @@ const MPI_REPLACE = MPI.MPI_REPLACE const PetscOptions = Ptr{Cvoid} const PetscViewer = Ptr{Cvoid} const Vec = Ptr{Cvoid} +const PetscObject = Ptr{Cvoid} diff --git a/lib/petsc_library.jl b/lib/petsc_library.jl index 56ba6702..daec294b 100644 --- a/lib/petsc_library.jl +++ b/lib/petsc_library.jl @@ -21,6 +21,7 @@ const MPI_REPLACE = MPI.MPI_REPLACE const PetscOptions = Ptr{Cvoid} const PetscViewer = Ptr{Cvoid} const Vec = Ptr{Cvoid} +const PetscObject = Ptr{Cvoid} const __darwin_off_t = Int64 @@ -56,10 +57,6 @@ mutable struct _p_PetscToken end const PetscToken = Ptr{_p_PetscToken} -mutable struct _p_PetscObject end - -const PetscObject = Ptr{_p_PetscObject} - const PetscObjectId = PetscInt64 const PetscObjectState = PetscInt64 diff --git a/src/LibPETSc.jl b/src/LibPETSc.jl index a1f4c8d0..4c2844f2 100644 --- a/src/LibPETSc.jl +++ b/src/LibPETSc.jl @@ -8,6 +8,8 @@ export PetscLibType export petsclibs export PetscBool export PETSC_TRUE +export UnionPetscLibType +export getlib include("const.jl") include("startup.jl") diff --git a/src/PETSc.jl b/src/PETSc.jl index c1eeb59b..e9903b51 100644 --- a/src/PETSc.jl +++ b/src/PETSc.jl @@ -22,7 +22,7 @@ end include("init.jl") include("viewer.jl") include("options.jl") -# include("vec.jl") +include("vec.jl") # include("mat.jl") # include("matshell.jl") # include("ksp.jl") diff --git a/src/init.jl b/src/init.jl index 31a08634..b1415177 100644 --- a/src/init.jl +++ b/src/init.jl @@ -6,7 +6,7 @@ Check if `petsclib` is initialized # External Links $(_doc_external("Sys/PetscInitialized")) """ -function initialized(petsclib::PetscLibType) +function initialized(petsclib) r_flag = Ref{LibPETSc.PetscBool}() LibPETSc.PetscInitialized(petsclib, r_flag) return r_flag[] == LibPETSc.PETSC_TRUE @@ -33,7 +33,7 @@ function initialize() return nothing end -function initialize(petsclib::PetscLibType) +function initialize(petsclib) if !initialized(petsclib) MPI.Initialized() || MPI.Init() LibPETSc.PetscInitializeNoArguments(petsclib) @@ -60,7 +60,7 @@ function finalize() return nothing end -function finalize(petsclib::PetscLibType) +function finalize(petsclib) if !finalized(petsclib) LibPETSc.PetscFinalize(petsclib) end @@ -75,7 +75,7 @@ Check if `petsclib` is finalized # External Links $(_doc_external("Sys/PetscFinalized")) """ -function finalized(petsclib::PetscLibType) +function finalized(petsclib) r_flag = Ref{LibPETSc.PetscBool}() LibPETSc.PetscFinalized(petsclib, r_flag) return r_flag[] == LibPETSc.PETSC_TRUE diff --git a/src/lib.jl b/src/lib.jl index 7434386f..96f94463 100644 --- a/src/lib.jl +++ b/src/lib.jl @@ -14,6 +14,19 @@ function PetscLibType{ST, IT}(petsc_library) where {ST, IT} LT = typeof(petsc_library) return PetscLibType{ST, IT, LT}(petsc_library) end +const UnionPetscLibType = Union{PetscLibType, Type{PetscLibType}} + +function Base.getproperty(petsclib::PetscLibType, name::Symbol) + if name == :PetscScalar + return scalartype(petsclib) + elseif name == :PetscReal + return realtype(petsclib) + elseif name == :PetscInt + return inttype(petsclib) + else + return getfield(petsclib, name) + end +end """ scalartype(petsclib::PetscLibType) @@ -102,3 +115,7 @@ macro for_petsc(expr) end end end + +@for_petsc begin + getlib(::Type{$PetscLib}) = $petsclib +end diff --git a/src/sys.jl b/src/sys.jl index dbc97590..654e9f79 100644 --- a/src/sys.jl +++ b/src/sys.jl @@ -1,8 +1,6 @@ const CPetscObject = Ptr{Cvoid} -const UnionPetscTypes = Union{ - Options, -} +const UnionPetscTypes = Union{Options, AbstractVec} # allows us to pass PETSc_XXX objects directly into CXXX ccall signatures Base.cconvert(::Type{CPetscObject}, obj::UnionPetscTypes) = obj @@ -13,46 +11,30 @@ function Base.unsafe_convert(::Type{Ptr{CPetscObject}}, obj::UnionPetscTypes) convert(Ptr{CPetscObject}, pointer_from_objref(obj)) end -#= -@for_libpetsc begin - function getcomm( - obj::Union{ - AbstractKSP{$PetscScalar}, - AbstractMat{$PetscScalar}, - AbstractVec{$PetscScalar}, - }, +function getcomm(obj::Union{AbstractVec{PetscLib}}) where {PetscLib} + comm = MPI.Comm() + LibPETSc.PetscObjectGetComm(PetscLib, obj, comm) + + #XXX We should really increase the petsc reference counter. + # But for for some reason the PETSc says that this communicator is + # unknown + #= + # Call the PetscCommDuplicate to increase reference count + @chk ccall( + (:PetscCommDuplicate, $libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, Ptr{MPI.MPI_Comm}, Ptr{Cvoid}), + comm, + comm, + C_NULL, ) - comm = MPI.Comm() - @chk ccall( - (:PetscObjectGetComm, $libpetsc), - PetscErrorCode, - (CPetscObject, Ptr{MPI.MPI_Comm}), - obj, - comm, - ) - - #XXX We should really increase the petsc reference counter. - # But for for some reason the PETSc says that this communicator is - # unknown - #= - # Call the PetscCommDuplicate to increase reference count - @chk ccall( - (:PetscCommDuplicate, $libpetsc), - PetscErrorCode, - (MPI.MPI_Comm, Ptr{MPI.MPI_Comm}, Ptr{Cvoid}), - comm, - comm, - C_NULL, - ) - # Register PetscCommDestroy to decriment the reference count - finalizer(PetscCommDestroy, comm) - =# + # Register PetscCommDestroy to decriment the reference count + finalizer(PetscCommDestroy, comm) + =# - return comm - end + return comm end -=# #= #XXX Not sure why this doesn't work diff --git a/src/vec.jl b/src/vec.jl index 24d1ebbb..d3025973 100644 --- a/src/vec.jl +++ b/src/vec.jl @@ -6,60 +6,99 @@ const CVec = Ptr{Cvoid} -abstract type AbstractVec{T} <: AbstractVector{T} end -scalartype(::AbstractVec{T}) where {T} = T +abstract type AbstractVec{PetscLib, PetscScalar} <: AbstractVector{PetscScalar} end + +Base.eltype( + ::Type{V}, +) where { + V <: AbstractVec{PetscLib, PetscScalar}, +} where {PetscLib, PetscScalar} = PetscScalar +Base.eltype( + v::AbstractVec{PetscLib, PetscScalar}, +) where {PetscLib, PetscScalar} = PetscScalar +Base.size(v::AbstractVec) = (length(v),) """ VecSeq(v::Vector) -A standard, sequentially-stored serial PETSc vector, wrapping the Julia vector `v`. +A standard, sequentially-stored serial PETSc vector, wrapping the Julia vector +`v`. -This reuses the array `v` as storage, and so `v` should not be `resize!`-ed or otherwise have its length modified while the PETSc object exists. +This reuses the array `v` as storage, and so `v` should not be `resize!`-ed or +otherwise have its length modified while the PETSc object exists. -This should only be need to be called for more advanced uses, for most simple usecases, users should be able to pass `Vector`s directly and have the wrapping performed automatically +This should only be need to be called for more advanced uses, for most simple +usecases, users should be able to pass `Vector`s directly and have the wrapping +performed automatically # External Links $(_doc_external("Vec/VecCreateSeqWithArray")) """ -mutable struct VecSeq{T} <: AbstractVec{T} +mutable struct VecSeq{PetscLib, PetscScalar} <: + AbstractVec{PetscLib, PetscScalar} ptr::CVec - array::Vector{T} + array::Vector{PetscScalar} +end +Base.parent(v::VecSeq) = v.array + +function VecSeq( + petsclib::PetscLib, + array::Vector{PetscScalar}; + blocksize = 1, + comm::MPI.Comm = MPI.COMM_SELF, +) where {PetscLib, PetscScalar} + @assert initialized(petsclib) + @assert PetscScalar == petsclib.PetscScalar + v = VecSeq{PetscLib, PetscScalar}(C_NULL, array) + LibPETSc.VecCreateSeqWithArray( + petsclib, + comm, + blocksize, + length(array), + array, + v, + ) + finalizer(destroy, v) + return v end -Base.eltype(::Type{V}) where {V<:AbstractVec{T}} where T = T -Base.eltype(v::AbstractVec{T}) where {T} = T -Base.size(v::AbstractVec) = (length(v),) -Base.parent(v::AbstractVec) = v.array +function destroy(v::AbstractVec{PetscLib}) where {PetscLib} + finalized(PetscLib) || LibPETSc.VecDestroy(PetscLib, v) + return nothing +end -@for_libpetsc begin - function VecSeq(comm::MPI.Comm, X::Vector{$PetscScalar}; blocksize=1) - @assert initialized($petsclib) - v = VecSeq(C_NULL, X) - @chk ccall((:VecCreateSeqWithArray, $libpetsc), PetscErrorCode, - (MPI.MPI_Comm, $PetscInt, $PetscInt, Ptr{$PetscScalar}, Ptr{CVec}), - comm, blocksize, length(X), X, v) - finalizer(destroy, v) - return v - end - function destroy(v::AbstractVec{$PetscScalar}) - finalized($petsclib) || - @chk ccall((:VecDestroy, $libpetsc), PetscErrorCode, (Ptr{CVec},), v) - return nothing - end - function Base.length(v::AbstractVec{$PetscScalar}) - r_sz = Ref{$PetscInt}() - @chk ccall((:VecGetSize, $libpetsc), PetscErrorCode, - (CVec, Ptr{$PetscInt}), v, r_sz) - return r_sz[] - end - function LinearAlgebra.norm(v::AbstractVec{$PetscScalar}, normtype::NormType=NORM_2) - r_val = Ref{$PetscReal}() - @chk ccall((:VecNorm, $libpetsc), PetscErrorCode, - (CVec, NormType, Ptr{$PetscReal}), - v, normtype,r_val) - return r_val[] - end +function Base.length(v::AbstractVec{PetscLib}) where {PetscLib} + petsclib = getlib(PetscLib) + PetscInt = petsclib.PetscInt + r_sz = Ref{PetscInt}() + LibPETSc.VecGetSize(PetscLib, v, r_sz) + return r_sz[] +end +function LinearAlgebra.norm( + v::AbstractVec{PetscLib}, + normtype::LibPETSc.NormType = LibPETSc.NORM_2, +) where {PetscLib} + petsclib = getlib(PetscLib) + PetscReal = petsclib.PetscReal + r_val = Ref{PetscReal}() + LibPETSc.VecNorm(PetscLib, v, normtype, r_val) + return r_val[] +end + + +function view( + vec::AbstractVec{PetscLib}, + viewer = LibPETSc.PETSC_VIEWER_STDOUT_(PetscLib, getcomm(vec)) + ) where PetscLib + LibPETSc.VecView(PetscLib, vec, viewer) + return nothing +end +Base.show(io::IO, opts::AbstractVec) = _show(io, opts) +Base.print_array(io::IO, opts::AbstractVec) = _show(io, opts) + +#= +@for_libpetsc begin function assemblybegin(V::AbstractVec{$PetscScalar}) @chk ccall((:VecAssemblyBegin, $libpetsc), PetscErrorCode, (CVec,), V) return nothing @@ -77,13 +116,6 @@ Base.parent(v::AbstractVec) = v.array r_lo[]:(r_hi[]-$PetscInt(1)) end - function view(vec::AbstractVec{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, getcomm(vec))) - @chk ccall((:VecView, $libpetsc), PetscErrorCode, - (CVec, CPetscViewer), - vec, viewer); - return nothing - end - function localsize(v::AbstractVec{$PetscScalar}) return r_sz[] end @@ -140,14 +172,9 @@ Use `read=false` if the array is write-only; `write=false` if read-only. """ unsafe_localarray -function Base.show(io::IO, ::MIME"text/plain", vec::AbstractVec) - _show(io, vec) -end - VecSeq(X::Vector{T}; kwargs...) where {T} = VecSeq(MPI.COMM_SELF, X; kwargs...) AbstractVec(X::AbstractVector) = VecSeq(X) - """ ownership_range(vec::AbstractVec) @@ -159,4 +186,4 @@ Note: unlike the C function, the range returned is inclusive (`idx_first:idx_las $(_doc_external("Vec/VecGetOwnershipRange")) """ ownershiprange - +=# diff --git a/test/options.jl b/test/options.jl index 886e1706..189c8bac 100644 --- a/test/options.jl +++ b/test/options.jl @@ -1,7 +1,7 @@ using Test using PETSc -@testset "options tests" begin +@testset "options" begin kw_opts = ( ksp_monitor = nothing, ksp_view = true, @@ -84,7 +84,7 @@ using PETSc end end -@testset "parse_options tests" begin +@testset "parse_options" begin @test begin julia = joinpath(Sys.BINDIR, Base.julia_exename()) run(`$(julia) --project -e "using PETSc