From 603711979e221797a282922c11ed3b453334dc6f Mon Sep 17 00:00:00 2001 From: Jeremy E Kozdon Date: Mon, 19 Jul 2021 14:09:21 -0700 Subject: [PATCH] Add back in vec.jl --- src/PETSc.jl | 2 +- src/sys.jl | 47 ++++++- src/vec.jl | 312 +++++++++++++++++++++++++++++------------------ test/runtests.jl | 2 +- test/vec.jl | 40 ++++++ 5 files changed, 280 insertions(+), 123 deletions(-) create mode 100644 test/vec.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/sys.jl b/src/sys.jl index 3118a7e0..5ec2d5f0 100644 --- a/src/sys.jl +++ b/src/sys.jl @@ -1,6 +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 @@ -10,3 +10,48 @@ Base.unsafe_convert(::Type{CPetscObject}, obj::UnionPetscTypes) = obj.ptr function Base.unsafe_convert(::Type{Ptr{CPetscObject}}, obj::UnionPetscTypes) convert(Ptr{CPetscObject}, pointer_from_objref(obj)) end + +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, + ) + + # Register PetscCommDestroy to decriment the reference count + finalizer(PetscCommDestroy, comm) + =# + + return comm +end + +#= +#XXX Not sure why this doesn't work +@for_libpetsc begin + function PetscCommDestroy( + comm::MPI.Comm + ) + @show comm.val + @chk ccall( + (:PetscCommDestroy, $libpetsc), + PetscErrorCode, + (Ptr{MPI.MPI_Comm},), + comm, + ) + return nothing + end +end +=# diff --git a/src/vec.jl b/src/vec.jl index 24d1ebbb..6ab0f1c0 100644 --- a/src/vec.jl +++ b/src/vec.jl @@ -6,157 +6,229 @@ 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) + VecSeq(petsclib, 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, +) where {PetscLib, PetscScalar} + comm = MPI.COMM_SELF + @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 assemblybegin(V::AbstractVec{$PetscScalar}) - @chk ccall((:VecAssemblyBegin, $libpetsc), PetscErrorCode, (CVec,), V) - return nothing - end - function assemblyend(V::AbstractVec{$PetscScalar}) - @chk ccall((:VecAssemblyEnd, $libpetsc), PetscErrorCode, (CVec,), V) - return nothing - end +function locallength(v::AbstractVec{PetscLib}) where {PetscLib} + petsclib = getlib(PetscLib) + PetscInt = petsclib.PetscInt + r_sz = Ref{PetscInt}() + LibPETSc.VecGetLocalSize(PetscLib, v, r_sz) + return r_sz[] +end - function ownershiprange(vec::AbstractVec{$PetscScalar}) - r_lo = Ref{$PetscInt}() - r_hi = Ref{$PetscInt}() - @chk ccall((:VecGetOwnershipRange, $libpetsc), PetscErrorCode, - (CVec, Ptr{$PetscInt}, Ptr{$PetscInt}), vec, r_lo, r_hi) - r_lo[]:(r_hi[]-$PetscInt(1)) - 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{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, getcomm(vec))) - @chk ccall((:VecView, $libpetsc), PetscErrorCode, - (CVec, CPetscViewer), - vec, viewer); - return nothing - 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, vec::AbstractVec) = _show(io, vec) +Base.show(io::IO, ::MIME"text/plain", vec::AbstractVec) = _show(io, vec) - function localsize(v::AbstractVec{$PetscScalar}) - return r_sz[] - end +""" + ownership_range(vec::AbstractVec, [base_one = true]) - function unsafe_localarray(::Type{$PetscScalar}, cv::CVec; read::Bool=true, write::Bool=true) - r_pv = Ref{Ptr{$PetscScalar}}() - if write - if read - @chk ccall((:VecGetArray, $libpetsc), PetscErrorCode, - (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv) - else - @chk ccall((:VecGetArrayWrite, $libpetsc), PetscErrorCode, - (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv) - end - else - @chk ccall((:VecGetArrayRead, $libpetsc), PetscErrorCode, - (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv) - end - r_sz = Ref{$PetscInt}() - @chk ccall((:VecGetLocalSize, $libpetsc), PetscErrorCode, - (CVec, Ptr{$PetscInt}), cv, r_sz) - v = unsafe_wrap(Array, r_pv[], r_sz[]; own = false) - - if write - if read - finalizer(v) do v - @chk ccall((:VecRestoreArray, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v))) - return nothing - end - else - finalizer(v) do v - @chk ccall((:VecRestoreArrayWrite, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v))) - return nothing - end - end - else - finalizer(v) do v - @chk ccall((:VecRestoreArrayRead, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v))) - return nothing - end - end - return v - end -end +The range of indices owned by this processor, assuming that the `vec` is laid +out with the first `n1` elements on the first processor, next `n2` elements on +the second, etc. For certain parallel layouts this range may not be well +defined. + +If the optional argument `base_one == true` then base-1 indexing is used, +otherwise base-0 index is used. +!!! note + + unlike the C function, the range returned is inclusive (`idx_first:idx_last`) + +# External Links +$(_doc_external("Vec/VecGetOwnershipRange")) """ - unsafe_localarray(PetscScalar, ptr:CVec; read=true, write=true) +function ownershiprange( + vec::AbstractVec{PetscLib}, + base_one::Bool = true, +) where {PetscLib} + petsclib = getlib(PetscLib) + PetscInt = petsclib.PetscInt + r_lo = Ref{PetscInt}() + r_hi = Ref{PetscInt}() + LibPETSc.VecGetOwnershipRange(PetscLib, vec, r_lo, r_hi) + return base_one ? ((r_lo[] + PetscInt(1)):(r_hi[])) : + ((r_lo[]):(r_hi[] - PetscInt(1))) +end -Return an `Array{PetscScalar}` containing local portion of the PETSc data. +""" + unsafe_localarray(vec::AbstractVec; read=true, write=true) -`finalize` should be called on the `Array` before the data can be used. +Return an `Array{PetscScalar}` containing local portion of the PETSc `vec` 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 +!!! note + `Base.finalize` should be called on the `Array` before the data can be used. -VecSeq(X::Vector{T}; kwargs...) where {T} = VecSeq(MPI.COMM_SELF, X; kwargs...) -AbstractVec(X::AbstractVector) = VecSeq(X) +# External Links +$(_doc_external("Vec/VecGetArray")) +$(_doc_external("Vec/VecGetArrayWrite")) +$(_doc_external("Vec/VecGetArrayRead")) +$(_doc_external("Vec/VecRestoreArray")) +$(_doc_external("Vec/VecRestoreArrayWrite")) +$(_doc_external("Vec/VecRestoreArrayRead")) +""" +function unsafe_localarray( + vec::AbstractVec{PetscLib}; + read::Bool = true, + write::Bool = true, +) where {PetscLib} + petsclib = getlib(PetscLib) + PetscScalar = petsclib.PetscScalar + r_pv = Ref{Ptr{PetscScalar}}() + if write && read + LibPETSc.VecGetArray(PetscLib, vec, r_pv) + elseif write + LibPETSc.VecGetArrayWrite(PetscLib, vec, r_pv) + elseif read + LibPETSc.VecGetArrayRead(PetscLib, vec, r_pv) + else + error("either read or write should be true") + end + sz = locallength(vec) + v = unsafe_wrap(Array, r_pv[], sz; own = false) + if write && read + finalizer(v) do v + LibPETSc.VecRestoreArray(PetscLib, vec, Ref(pointer(v))) + return nothing + end + elseif write + finalizer(v) do v + LibPETSc.VecRestoreArrayWrite(PetscLib, vec, Ref(pointer(v))) + return nothing + end + elseif read + finalizer(v) do v + LibPETSc.VecRestoreArrayRead(PetscLib, vec, Ref(pointer(v))) + return nothing + end + end + return v +end """ - ownership_range(vec::AbstractVec) + with_unsafe_localarray!( + f!, + x::AbstractVec; + read=true, + write=true, + ) -The range of indices owned by this processor, assuming that the vectors are laid out with the first n1 elements on the first processor, next n2 elements on the second, etc. For certain parallel layouts this range may not be well defined. +Convert `x` to an `Array{PetscScalar}` using [`unsafe_localarray`](@ref) and +apply the function `f!`. -Note: unlike the C function, the range returned is inclusive (`idx_first:idx_last`) +Use `read=false` if the array is write-only; `write=false` if read-only. -# External Links -$(_doc_external("Vec/VecGetOwnershipRange")) +# Examples +```julia-repl +julia> map_unsafe_localarray(x; write=true) do x + @. x .*= 2 +end + +!!! note + `Base.finalize` is automatically called on the array. """ -ownershiprange +function with_unsafe_localarray!(f!, v::AbstractVec; kwargs...) + array = unsafe_localarray(v; kwargs...) + f!(array) + Base.finalize(array) +end +#= +@for_libpetsc begin + function assemblybegin(V::AbstractVec{$PetscScalar}) + @chk ccall((:VecAssemblyBegin, $libpetsc), PetscErrorCode, (CVec,), V) + return nothing + end + function assemblyend(V::AbstractVec{$PetscScalar}) + @chk ccall((:VecAssemblyEnd, $libpetsc), PetscErrorCode, (CVec,), V) + return nothing + end +end +=# diff --git a/test/runtests.jl b/test/runtests.jl index cac13ea3..410f36ba 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ include("init.jl") include("options.jl") -# include("vec.jl") +include("vec.jl") # include("mat.jl") diff --git a/test/vec.jl b/test/vec.jl new file mode 100644 index 00000000..14201b41 --- /dev/null +++ b/test/vec.jl @@ -0,0 +1,40 @@ +using Test +using PETSc +using LinearAlgebra: norm + +@testset "VecSeq" begin + N = 10 + for petsclib in PETSc.petsclibs + PETSc.initialize(petsclib) + PetscScalar = petsclib.PetscScalar + x = rand(PetscScalar, N) + petsc_x = PETSc.VecSeq(petsclib, x) + @test length(petsc_x) == N + @test norm(petsc_x) ≈ norm(x) + + # make sure the viewer works + _stdout = stdout + (rd, wr) = redirect_stdout() + @show petsc_x + @test readline(rd) == "petsc_x = Vec Object: 1 MPI processes" + @test readline(rd) == " type: seq" + redirect_stdout(_stdout) + + _stdout = stdout + (rd, wr) = redirect_stdout() + show(stdout, "text/plain", petsc_x) + @test readline(rd) == "Vec Object: 1 MPI processes" + @test readline(rd) == " type: seq" + redirect_stdout(_stdout) + + @test PETSc.ownershiprange(petsc_x) == 1:N + @test PETSc.ownershiprange(petsc_x, false) == 0:(N - 1) + + PETSc.with_unsafe_localarray!(petsc_x) do x2 + @test x2 == x + end + + PETSc.destroy(petsc_x) + PETSc.finalize(petsclib) + end +end