diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..d8737f8 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,25 @@ +# https://github.com/domluna/JuliaFormatter.jl + +# Annotates fields in a type definitions with ::Any if no type annotation is provided +annotate_untyped_fields_with_any = false + +# Transforms a short function definition to a long function definition if the short +# function definition exceeds the maximum margin. +short_to_long_function_def = true + +# Transforms a long function definition to a short function definition if the short +# function definition does not exceed the maximum margin. +long_to_short_function_def = true + +# Format code docstrings with the same options used for the code source. +format_docstrings = true + +# If true, Markdown files are also formatted. Julia code blocks will be formatted +# in addition to the Markdown being normalized. +format_markdown = true + +# If true, whitespace is added for binary operations in indices +whitespace_ops_in_indices = true + +# If true, superfluous newlines will be removed +remove_extra_newlines = true \ No newline at end of file diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 2d76bf2..3594060 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -4,7 +4,7 @@ on: push: branches: - master - tags: '*' + tags: "*" pull_request: workflow_dispatch: @@ -15,7 +15,7 @@ jobs: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@latest with: - version: '1.5' + version: "1.8" - name: Install dependencies run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' - name: Build and deploy @@ -23,4 +23,4 @@ jobs: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key DOC_DEPLOYMENT: true # Workaround to avoid finding PETSc location - run: julia --project=docs/ docs/make.jl \ No newline at end of file + run: julia --project=docs/ docs/make.jl diff --git a/Project.toml b/Project.toml index 49fb943..9ffa73c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "PetscWrap" uuid = "5be22e1c-01b5-4697-96eb-ef9ccdc854b8" authors = ["bmxam "] -version = "0.1.5" +version = "0.2.0" [deps] Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] -MPI = "0.16, 0.17, 0.18, 0.19" -julia = "1.5" +MPI = "0.20" +julia = "1.5, 1.6, 1.7, 1.8" diff --git a/README.md b/README.md index 2102bec..90b0924 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,7 @@ PetscWrap.jl is a parallel Julia wrapper for the (awesome) [PETSc](https://www.mcs.anl.gov/petsc/) library. It can be considered as a fork from the [GridapPetsc.jl](https://github.com/gridap/GridapPETSc.jl) and [Petsc.jl](https://github.com/JuliaParallel/PETSc.jl) projects : these two projects have extensively inspired this project, and some code has even been directly copied. The main differences with the two aformentionned projects are: + - parallel support : you can solve linear systems on multiple core with `mpirun -n 4 julia foo.jl`; - no dependance to other Julia packages except `MPI.jl`; - possibility to switch from one PETSc "arch" to another; @@ -13,28 +14,45 @@ The main differences with the two aformentionned projects are: Note that the primary objective of this project is to enable the wrapper of the SLEPc library through the [SlepcWrap.jl](https://github.com/bmxam/SlepcWrap.jl) project. +This project is only a wrapper to PETSc functions, the purpose is not to deliver a julia `Array` (it maybe be one day the purpose of a package `PetscArrays.jl`). + ## How to install it + You must have installed the PETSc library on your computer and set the two following environment variables : `PETSC_DIR` and `PETSC_ARCH`. At run time, PetscWrap.jl looks for the `libpetsc.so` using these environment variables and "load" the library. To install the package, use the Julia package manager: + ```Julia pkg> add PetscWrap ``` + ## Contribute + Any contribution(s) and/or remark(s) are welcome! If you need a function that is not wrapped yet but you don't think you are capable of contributing, post an issue with a minimum working example. +Conventions to be applied in future versions ("fancy" stuff is not concerned): + +- all PETSc types should have the exact same name in Julia; +- all PETSc functions should have the exact same name in julia, but without the type as a prefix, and with a lower case for the first letter. `VecSetValues` becomes `setValues`. This rule is not applied when the name conflicts with a name from `Base` (for instance `VecView` becomes `vecView` and not `view`); +- all PETSc functions must have the same number of arguments and, if possible the same names in julia, except for out-of-place arguments. +- functions arguments must all be typed. Additional functions, without type or with fewer args, can be defined if the original version is present. + ## PETSc compat. -This version of PetscWrap.jl has been tested with petsc-3.13. Complex numbers are supported. + +This version of PetscWrap.jl has been tested with petsc-3.19. Complex numbers are supported. ## How to use it -PETSc methods wrappers share the same name as their C equivalent : for instance `MatCreate` or `MatSetValue`. Furthermore, an optional "higher level" API, referred to as "fancy", is exposed : for instance `create_matrix` or `A[i,j] = v`). Note that this second way of manipulating PETSc will evolve according the package's author needs while the first one will try to follow PETSc official API. + +PETSc methods wrappers share the almost same name as their C equivalent (with the type) : for instance `MatSetValues` becomes `setValues`. Furthermore, an optional "higher level" API, referred to as "fancy", is exposed : for instance `create_matrix` or `A[i,j] = v`. Note that this second way of manipulating PETSc will evolve according the package's author needs while the first one will try to follow PETSc official API. You will find examples of use by building the documentation: `julia PetscWrap.jl/docs/make.jl`. Here is one of the examples: + ### A first demo + This example serves as a test since this project doesn't have a "testing" procedure yet. In this example, -the equation ``u'(x) = 2`` with ``u(0) = 0`` is solved on the domain ``[0,1]`` using a backward finite +the equation `u'(x) = 2` with `u(0) = 0` is solved on the domain `[0,1]` using a backward finite difference scheme. In this example, PETSc classic method names are used. For more fancy names, check the fancy version. @@ -63,36 +81,36 @@ Number of mesh points and mesh step ```julia n = 11 -Δx = 1. / (n - 1) +Δx = 1.0 / (n - 1) ``` -Create a matrix and a vector +Create a matrix and a vector (you can specify the MPI communicator if you want) ```julia -A = MatCreate() -b = VecCreate() +A = create(Mat) +b = create(Vec) ``` Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should set the number of preallocated non-zeros to increase performance. ```julia -MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) -VecSetSizes(b, PETSC_DECIDE, n) +setSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) +setSizes(b, PETSC_DECIDE, n) ``` We can then use command-line options to set our matrix/vectors. ```julia -MatSetFromOptions(A) -VecSetFromOptions(b) +setFromOptions(A) +setFromOptions(b) ``` Finish the set up ```julia -MatSetUp(A) -VecSetUp(b) +setUp(A) +setUp(b) ``` Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. @@ -101,88 +119,87 @@ As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first the fancy version of this example for a more Julia-like convention. ```julia -b_start, b_end = VecGetOwnershipRange(b) +b_start, b_end = getOwnershipRange(b) ``` Now let's build the right hand side vector. Their are various ways to do this, this is just one. ```julia -n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start... -VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc)) +n_loc = getLocalSize(b) # Note that n_loc = b_end - b_start... +setValues(b, collect(b_start:(b_end - 1)), 2 * ones(n_loc)) ``` And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. ```julia -A_start, A_end = MatGetOwnershipRange(A) -for i in A_start:A_end-1 - MatSetValues(A, [i], [i-1, i], [-1. 1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES) +A_start, A_end = getOwnershipRange(A) +for i = A_start:(A_end - 1) + setValues(A, [i], [i - 1, i], [-1.0 1.0] / Δx, INSERT_VALUES) # setValues(A, I, J, V, INSERT_VALUES) end ``` Set boundary condition (only the proc handling index `0` is acting) ```julia -(b_start == 0) && VecSetValue(b, 0, 0.) +(b_start == 0) && setValue(b, 0, 0.0) ``` Assemble matrices ```julia -MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY) -VecAssemblyBegin(b) -MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY) -VecAssemblyEnd(b) +assemblyBegin(A, MAT_FINAL_ASSEMBLY) +assemblyBegin(b) +assemblyEnd(A, MAT_FINAL_ASSEMBLY) +assemblyEnd(b) ``` At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call ```julia -MatView(A) -VecView(b) +matView(A) +vecView(b) ``` Set up the linear solver ```julia -ksp = KSPCreate() -KSPSetOperators(ksp, A, A) -KSPSetFromOptions(ksp) -KSPSetUp(ksp) +ksp = create(KSP) +setOperators(ksp, A, A) +setFromOptions(ksp) +setUp(ksp) ``` Solve the system. We first allocate the solution using the `VecDuplicate` method. ```julia -x = VecDuplicate(b) -KSPSolve(ksp, b, x) +x = duplicate(b) +solve(ksp, b, x) ``` Print the solution ```julia -VecView(x) +vecView(x) ``` Access the solution (this part is under development), getting a Julia array; and then restore it ```julia -array, ref = VecGetArray(x) # do something with array +array, ref = getArray(x) # do something with array @show array -VecRestoreArray(x, ref) +restoreArray(x, ref) ``` -Free memory +Free memory. Note that this call is faculative since, by default, +the julia GC will trigger a call to Petsc `destroy` to each object ```julia -MatDestroy(A) -VecDestroy(b) -VecDestroy(x) +destroy.((ksp, A, b, x)) ``` -Finalize Petsc +Finalize Petsc. This call is faculative : it will be triggered automatically at the end of the script. ```julia PetscFinalize() -``` \ No newline at end of file +``` diff --git a/example/linear_system.jl b/example/linear_system.jl index 8670ad3..b1f5fa6 100644 --- a/example/linear_system.jl +++ b/example/linear_system.jl @@ -22,78 +22,78 @@ PetscInitialize() # Number of mesh points and mesh step n = 11 -Δx = 1. / (n - 1) +Δx = 1.0 / (n - 1) # Create a matrix and a vector (you can specify the MPI communicator if you want) -A = MatCreate() -b = VecCreate() +A = create(Mat) +b = create(Vec) # Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should # set the number of preallocated non-zeros to increase performance. -MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) -VecSetSizes(b, PETSC_DECIDE, n) +setSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) +setSizes(b, PETSC_DECIDE, n) # We can then use command-line options to set our matrix/vectors. -MatSetFromOptions(A) -VecSetFromOptions(b) +setFromOptions(A) +setFromOptions(b) # Finish the set up -MatSetUp(A) -VecSetUp(b) +setUp(A) +setUp(b) # Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. # As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first row handled by the local processor # (starting at 0), and the last row (which is `rend-1`). This may be disturbing for non-regular PETSc users. Checkout # the fancy version of this example for a more Julia-like convention. -b_start, b_end = VecGetOwnershipRange(b) +b_start, b_end = getOwnershipRange(b) # Now let's build the right hand side vector. Their are various ways to do this, this is just one. -n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start... -VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc)) +n_loc = getLocalSize(b) # Note that n_loc = b_end - b_start... +setValues(b, collect(b_start:(b_end - 1)), 2 * ones(n_loc)) # And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. -A_start, A_end = MatGetOwnershipRange(A) -for i in A_start:A_end-1 - MatSetValues(A, [i], [i-1, i], [-1. 1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES) +A_start, A_end = getOwnershipRange(A) +for i = A_start:(A_end - 1) + setValues(A, [i], [i - 1, i], [-1.0 1.0] / Δx, INSERT_VALUES) # setValues(A, I, J, V, INSERT_VALUES) end # Set boundary condition (only the proc handling index `0` is acting) -(b_start == 0) && VecSetValue(b, 0, 0.) +(b_start == 0) && setValue(b, 0, 0.0) # Assemble matrices -MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY) -VecAssemblyBegin(b) -MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY) -VecAssemblyEnd(b) +assemblyBegin(A, MAT_FINAL_ASSEMBLY) +assemblyBegin(b) +assemblyEnd(A, MAT_FINAL_ASSEMBLY) +assemblyEnd(b) # At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call -MatView(A) -VecView(b) +matView(A) +vecView(b) # Set up the linear solver -ksp = KSPCreate() -KSPSetOperators(ksp, A, A) -KSPSetFromOptions(ksp) -KSPSetUp(ksp) +ksp = create(KSP) +setOperators(ksp, A, A) +setFromOptions(ksp) +setUp(ksp) # Solve the system. We first allocate the solution using the `VecDuplicate` method. -x = VecDuplicate(b) -KSPSolve(ksp, b, x) +x = duplicate(b) +solve(ksp, b, x) # Print the solution -VecView(x) +vecView(x) # Access the solution (this part is under development), getting a Julia array; and then restore it -array, ref = VecGetArray(x) # do something with array +array, ref = getArray(x) # do something with array @show array -VecRestoreArray(x, ref) +restoreArray(x, ref) -# Free memory -MatDestroy(A) -VecDestroy(b) -VecDestroy(x) +# Free memory. Note that this call is faculative since, by default, +# the julia GC will trigger a call to Petsc `destroy` to each object +destroy.((ksp, A, b, x)) -# Finalize Petsc +# Finalize Petsc. This call is faculative : it will +# be triggered automatically at the end of the script. PetscFinalize() -end #hide \ No newline at end of file +end #hide diff --git a/example/linear_system_fancy.jl b/example/linear_system_fancy.jl index 5cb704b..0fdc96b 100644 --- a/example/linear_system_fancy.jl +++ b/example/linear_system_fancy.jl @@ -22,19 +22,12 @@ PetscInitialize() # Number of mesh points and mesh step n = 11 -Δx = 1. / (n - 1) +Δx = 1.0 / (n - 1) -# Create a matrix of size `(n,n)` and a vector of size `(n)` -A = create_matrix(n, n) -b = create_vector(n) - -# We can then use command line options to set our matrix/vectors. -set_from_options!(A) -set_from_options!(b) - -# Finish the set up -set_up!(A) -set_up!(b) +# Create a matrix of size `(n,n)` and a vector of size `(n)`. The `autosetup` option +# triggers a call to `setFromOptions` and `setUp` +A = create_matrix(; nrows_glo = n, ncols_glo = n, autosetup = true) +b = create_vector(; nrows_glo = n, autosetup = true) # Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. # The `rstart, rend = get_range(*)` methods differ a little bit from PETSc API since the two integers it @@ -48,21 +41,19 @@ b[b_start:b_end] = 2 * ones(n_loc) # And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. A_start, A_end = get_range(A) -for i in A_start:A_end - A[i, i-1:i] = [-1. 1.] / Δx +for i = A_start:A_end + A[i, (i - 1):i] = [-1.0 1.0] / Δx end # Set boundary condition (only the proc handling index `1` is acting) -(b_start == 1) && (b[1] = 0.) +(b_start == 1) && (b[1] = 0.0) # Assemble matrice and vector assemble!(A) assemble!(b) # Set up the linear solver -ksp = create_ksp(A) -set_from_options!(ksp) -set_up!(ksp) +ksp = create_ksp(A; autosetup = true) # Solve the system x = solve(ksp, b) @@ -73,13 +64,11 @@ x = solve(ksp, b) # Convert `PetscVec` to Julia `Array` (and play with it!) array = vec2array(x) -# Free memory -destroy!(A) -destroy!(b) -destroy!(x) +# Free memory (optional, objects are garbage collected otherwise) +destroy!(A, b, x, ksp) # Note that it's also possible to build a matrix using the COO format as in `SparseArrays`: -M = create_matrix(3,3; auto_setup = true) +M = create_matrix(; nrows_glo = 3, ncols_glo = 3, autosetup = true) M_start, M_end = get_range(M) I = [1, 1, 1, 2, 3] J = [1, 3, 1, 3, 2] @@ -88,11 +77,9 @@ k = findall(x -> M_start <= x <= M_end, I) # just a trick to allow this example set_values!(M, I[k], J[k], V[k], ADD_VALUES) assemble!(M) @show M -destroy!(M) # This is very convenient in sequential since you can fill the three vectors I, J, V in your code and decide only # at the last moment if you'd like to use `SparseArrays` or `PetscMat`. -# Finalize Petsc -PetscFinalize() +destroy!(M) -end #hide \ No newline at end of file +end #hide diff --git a/src/ISLocalToGlobalMapping.jl b/src/ISLocalToGlobalMapping.jl new file mode 100644 index 0000000..70686fb --- /dev/null +++ b/src/ISLocalToGlobalMapping.jl @@ -0,0 +1,114 @@ +const CISLocalToGlobalMapping = Ptr{Cvoid} + +mutable struct ISLocalToGlobalMapping <: AbstractPetscObject + ptr::CISLocalToGlobalMapping + comm::MPI.Comm + + ISLocalToGlobalMapping(comm::MPI.Comm) = new(CISLocalToGlobalMapping(), comm) +end + +Base.unsafe_convert(::Type{CISLocalToGlobalMapping}, x::ISLocalToGlobalMapping) = x.ptr +function Base.unsafe_convert( + ::Type{Ptr{CISLocalToGlobalMapping}}, + x::ISLocalToGlobalMapping, +) + Ptr{CISLocalToGlobalMapping}(pointer_from_objref(x)) +end + +""" + create( + comm::MPI.Comm, + bs::PetscInt, + n::PetscInt, + indices::Vector{PetscInt}, + mode::PetscCopyMode; + add_finalizer = true, + ) + + create( + ::Type{ISLocalToGlobalMapping}, + comm::MPI.Comm, + indices::Vector{I}, + mode::PetscCopyMode = PETSC_COPY_VALUES, + mode::PetscCopyMode; + add_finalizer = true, + ) where {I<:Integer} + +Wrapper to `ISLocalToGlobalMappingCreate` +https://petsc.org/release/manualpages/IS/ISLocalToGlobalMappingCreate/ + +0-based indexing +""" +function create( + ::Type{ISLocalToGlobalMapping}, + comm::MPI.Comm, + bs::PetscInt, + n::PetscInt, + indices::Vector{PetscInt}, + mode::PetscCopyMode; + add_finalizer = true, +) + l2g = ISLocalToGlobalMapping(comm) + error = ccall( + (:ISLocalToGlobalMappingCreate, libpetsc), + PetscErrorCode, + ( + MPI.MPI_Comm, + PetscInt, + PetscInt, + Ptr{PetscInt}, + PetscCopyMode, + Ptr{CISLocalToGlobalMapping}, + ), + comm, + bs, + n, + indices, + mode, + l2g, + ) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, l2g) + + return l2g +end + +function create( + ::Type{ISLocalToGlobalMapping}, + comm::MPI.Comm, + indices::Vector{I}, + mode::PetscCopyMode = PETSC_COPY_VALUES; + add_finalizer = true, +) where {I<:Integer} + return create( + ISLocalToGlobalMapping, + comm, + PetscIntOne, + PetscInt(length(indices)), + PetscInt.(indices), + mode; + add_finalizer, + ) +end + +""" + destroy(mapping::ISLocalToGlobalMapping) + +Wrapper to `ISLocalToGlobalMappingDestroy` +https://petsc.org/release/manualpages/IS/ISLocalToGlobalMappingDestroy/ +""" +function destroy(mapping::ISLocalToGlobalMapping) + _is_destroyed(mapping) && return + + error = ccall( + (:ISLocalToGlobalMappingDestroy, libpetsc), + PetscErrorCode, + (Ptr{CISLocalToGlobalMapping},), + mapping, + ) + @assert iszero(error) + + _NREFS[] -= 1 +end diff --git a/src/KSP.jl b/src/KSP.jl new file mode 100644 index 0000000..7b5f29a --- /dev/null +++ b/src/KSP.jl @@ -0,0 +1,95 @@ +const CKSP = Ptr{Cvoid} + +mutable struct KSP <: AbstractPetscObject + ptr::CKSP + comm::MPI.Comm + + KSP(comm::MPI.Comm) = new(CKSP(), comm) +end + +Base.unsafe_convert(::Type{CKSP}, x::KSP) = x.ptr +Base.unsafe_convert(::Type{Ptr{CKSP}}, x::KSP) = Ptr{CKSP}(pointer_from_objref(x)) + +""" + create(::Type{KSP}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + +Wrapper for `KSPCreate` +https://petsc.org/release/manualpages/KSP/KSPCreate/ +""" +function create(::Type{KSP}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + ksp = KSP(comm) + error = + ccall((:KSPCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CKSP}), comm, ksp) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, ksp) + + return ksp +end + +""" + destroy(ksp::KSP) + +Wrapper to `KSPDestroy` +https://petsc.org/release/manualpages/KSP/KSPDestroy/ +""" +function destroy(ksp::KSP) + _is_destroyed(ksp) && return + + error = ccall((:KSPDestroy, libpetsc), PetscErrorCode, (Ptr{CKSP},), ksp) + @assert iszero(error) + + _NREFS[] -= 1 +end + +""" + setFromOptions(ksp::KSP) + +Wrapper to `KSPSetFromOptions` +https://petsc.org/release/manualpages/KSP/KSPSetFromOptions/ +""" +function setFromOptions(ksp::KSP) + error = ccall((:KSPSetFromOptions, libpetsc), PetscErrorCode, (CKSP,), ksp) + @assert iszero(error) +end + +""" + setOperators(ksp::KSP, Amat::Mat, Pmat::Mat) + +Wrapper for `KSPSetOperators`` +https://petsc.org/release/manualpages/KSP/KSPSetOperators/ +""" +function setOperators(ksp::KSP, Amat::Mat, Pmat::Mat) + error = ccall( + (:KSPSetOperators, libpetsc), + PetscErrorCode, + (CKSP, CMat, CMat), + ksp, + Amat, + Pmat, + ) + @assert iszero(error) +end + +""" + KSPSetUp(ksp::KSP) + +Wrapper to `KSPSetUp` +https://petsc.org/release/manualpages/KSP/KSPSetUp/ +""" +function setUp(ksp::KSP) + error = ccall((:KSPSetUp, libpetsc), PetscErrorCode, (CKSP,), ksp) + @assert iszero(error) +end + +""" + solve(ksp::KSP, b::Vec, x::Vec) + +Wrapper for `KSPSolve` +https://petsc.org/release/manualpages/KSP/KSPSolve/ +""" +function solve(ksp::KSP, b::Vec, x::Vec) + error = ccall((:KSPSolve, libpetsc), PetscErrorCode, (CKSP, CVec, CVec), ksp, b, x) + @assert iszero(error) +end diff --git a/src/Mat.jl b/src/Mat.jl new file mode 100644 index 0000000..3282172 --- /dev/null +++ b/src/Mat.jl @@ -0,0 +1,738 @@ +const CMat = Ptr{Cvoid} + +""" + A Petsc matrix, actually just a pointer to the actual C matrix +""" +mutable struct Mat <: AbstractPetscObject + ptr::CMat + comm::MPI.Comm + + Mat(comm::MPI.Comm) = new(CMat(), comm) +end + +Base.unsafe_convert(::Type{CMat}, x::Mat) = x.ptr +Base.unsafe_convert(::Type{Ptr{CMat}}, x::Mat) = Ptr{CMat}(pointer_from_objref(x)) + +""" + assemblyBegin(mat::Mat, type::MatAssemblyType) + +Wrapper to `MatAssemblyBegin` +https://petsc.org/release/manualpages/Mat/MatAssemblyBegin/ +""" +function assemblyBegin(mat::Mat, type::MatAssemblyType) + error = ccall( + (:MatAssemblyBegin, libpetsc), + PetscErrorCode, + (CMat, MatAssemblyType), + mat, + type, + ) + @assert iszero(error) +end + +""" + assemblyEnd(mat::Mat, type::MatAssemblyType) + +Wrapper to `MatAssemblyEnd` +https://petsc.org/release/manualpages/Mat/MatAssemblyEnd/ +""" +function assemblyEnd(mat::Mat, type::MatAssemblyType) + error = ccall( + (:MatAssemblyEnd, libpetsc), + PetscErrorCode, + (CMat, MatAssemblyType), + mat, + type, + ) + @assert iszero(error) +end + +""" + compositeAddMat(mat::Mat, smat::Mat) + +Wrapper to `MatCompositeAddMat` +https://petsc.org/release/manualpages/Mat/MatCompositeAddMat/ +""" +function compositeAddMat(mat::Mat, smat::Mat) + error = ccall((:MatCompositeAddMat, libpetsc), PetscErrorCode, (CMat, CMat), mat, smat) + @assert iszero(error) +end + +""" + create(::Type{Mat}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + +Wrapper to `MatCreate` +https://petsc.org/release/manualpages/Mat/MatCreate/ +""" +function create(::Type{Mat}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + mat = Mat(comm) + error = + ccall((:MatCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CMat}), comm, mat) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, mat) + + return mat +end + +""" + createDense( + comm::MPI.Comm, + m::PetscInt, + n::PetscInt, + M::PetscInt, + N::PetscInt; + add_finalizer = true, + ) + createDense( + comm::MPI.Comm, + m::Integer = PETSC_DECIDE, + n::Integer = PETSC_DECIDE, + M::Integer = PETSC_DECIDE, + N::Integer = PETSC_DECIDE; + add_finalizer = true, + ) + +Wrapper to `MatCreateDense` +https://petsc.org/release/manualpages/Mat/MatCreateDense/ + +Last argument `data` is not supported yet (NULL is passed). +""" +function createDense( + comm::MPI.Comm, + m::PetscInt, + n::PetscInt, + M::PetscInt, + N::PetscInt; + add_finalizer = true, +) + mat = Mat(comm) + error = ccall( + (:MatCreateDense, libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, PetscInt, PetscInt, PetscInt, PetscInt, Ptr{PetscScalar}, Ptr{CMat}), + comm, + m, + n, + M, + N, + C_NULL, + mat, + ) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, mat) + + return mat +end + +function createDense( + comm::MPI.Comm, + m::Integer = PETSC_DECIDE, + n::Integer = PETSC_DECIDE, + M::Integer = PETSC_DECIDE, + N::Integer = PETSC_DECIDE; + add_finalizer = true, +) + return createDense( + comm, + PetscInt(m), + PetscInt(n), + PetscInt(M), + PetscInt(N); + add_finalizer, + ) +end + +""" + createVecs(mat::Mat, vecr::Vec, veci::Vec) + createVecs(mat::Mat; add_finalizer = true) + +Wrapper to `MatCreateVecs` +https://petsc.org/release/manualpages/Mat/MatCreateVecs/ +""" +function createVecs(mat::Mat, right::Vec, left::Vec) + error = ccall( + (:MatCreateVecs, libpetsc), + PetscErrorCode, + (CMat, Ptr{CVec}, Ptr{CVec}), + mat, + right, + left, + ) + @assert iszero(error) +end + +function createVecs(mat::Mat; add_finalizer = true) + right = Vec(mat.comm) + left = Vec(mat.comm) + createVecs(mat, right, left) + + _NREFS[] += 1 + _NREFS[] += 1 + if add_finalizer + finalizer(destroy, right) + finalizer(destroy, left) + end + + return right, left +end + +""" + destroy(A::Mat) + +Wrapper to `MatDestroy` +https://petsc.org/release/manualpages/Mat/MatDestroy/ +""" +function destroy(A::Mat) + _is_destroyed(A) && return + + error = ccall((:MatDestroy, libpetsc), PetscErrorCode, (Ptr{CMat},), A) + @assert iszero(error) + + _NREFS[] -= 1 +end + +""" + duplicate(mat::Mat, op::MatDuplicateOption; add_finalizer = true) + +Wrapper for `MatDuplicate` +https://petsc.org/release/manualpages/Mat/MatDuplicate/ +""" +function duplicate(mat::Mat, op::MatDuplicateOption; add_finalizer = true) + M = Mat(mat.comm) + error = ccall( + (:MatDuplicate, libpetsc), + PetscErrorCode, + (CMat, MatDuplicateOption, Ptr{CMat}), + mat, + op, + M, + ) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, M) + + return M +end + +""" + getLocalSize(mat::Mat) + +Wrapper to `MatGetLocalSize` +https://petsc.org/release/manualpages/Mat/MatGetLocalSize/ + +Return the number of local rows and cols of the matrix (i.e on the processor). +""" +function getLocalSize(mat::Mat) + m = Ref{PetscInt}(0) + n = Ref{PetscInt}(0) + + error = ccall( + (:MatGetLocalSize, libpetsc), + PetscErrorCode, + (CMat, Ref{PetscInt}, Ref{PetscInt}), + mat, + m, + n, + ) + @assert iszero(error) + + return m[], n[] +end + +""" + getOwnershipRange(mat::Mat) + +Wrapper to `MatGetOwnershipRange` +https://petsc.org/release/manualpages/Mat/MatGetOwnershipRange/ + +The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor. + +# Warning + +`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row +handled by the local processor. +""" +function getOwnershipRange(mat::Mat) + rstart = Ref{PetscInt}(0) + rend = Ref{PetscInt}(0) + + error = ccall( + (:MatGetOwnershipRange, libpetsc), + PetscErrorCode, + (CMat, Ref{PetscInt}, Ref{PetscInt}), + mat, + rstart, + rend, + ) + @assert iszero(error) + + return rstart[], rend[] +end + +""" + getOwnershipRangeColumn(mat::Mat) + +Wrapper to `MatGetOwnershipRangeColumn` +https://petsc.org/release/manualpages/Mat/MatGetOwnershipRangeColumn/ + +The result `(cstart, cend)` is a Tuple indicating the columns handled by the local processor. + +# Warning + +`PETSc` indexing starts at zero (so `cstart` may be zero) and `cend-1` is the last column +handled by the local processor. +""" +function getOwnershipRangeColumn(mat::Mat) + cstart = Ref{PetscInt}(0) + cend = Ref{PetscInt}(0) + + error = ccall( + (:MatGetOwnershipRangeColumn, libpetsc), + PetscErrorCode, + (CMat, Ref{PetscInt}, Ref{PetscInt}), + mat, + cstart, + cend, + ) + @assert iszero(error) + + return cstart[], cend[] +end + +""" + getSize(mat::Mat) + +Wrapper to `MatGetSize` +https://petsc.org/release/manualpages/Mat/MatGetSize/ + +Return the number of rows and cols of the matrix (global number). +""" +function getSize(mat::Mat) + nrows_glo = Ref{PetscInt}(0) + ncols_glo = Ref{PetscInt}(0) + + error = ccall( + (:MatGetSize, libpetsc), + PetscErrorCode, + (CMat, Ref{PetscInt}, Ref{PetscInt}), + mat, + nrows_glo, + ncols_glo, + ) + @assert iszero(error) + + return nrows_glo[], ncols_glo[] +end + +""" + getType(mat::Mat) + +Wrapper to `MatGetType` +https://petsc.org/release/manualpages/Mat/MatType/ + +Return the matrix type as a string. See matrix types here: +https://petsc.org/release/manualpages/Mat/MatType.html +""" +function getType(mat::Mat) + type = Ref{Cstring}() + + error = ccall((:MatGetType, libpetsc), PetscErrorCode, (CMat, Ptr{Cstring}), mat, type) + @assert iszero(error) + + return unsafe_string(type[]) +end + +""" + MPIAIJSetPreallocation( + B::Mat, + d_nz::PetscInt, + d_nnz::Vector{PetscInt}, + o_nz::PetscInt, + o_nnz::Vector{PetscInt}, + ) + MPIAIJSetPreallocation(mat::Mat, dnz::PetscInt, onz::PetscInt) + +Wrapper to `MatMPIAIJSetPreallocation` +https://petsc.org/release//manualpages/Mat/MatMPIAIJSetPreallocation/ + +# Warning + +`dnnz` and `onnz` not tested yet. + +# Arguments + + - `dnz::PetscInt`: number of nonzeros per row in DIAGONAL portion of local submatrix (same value is used for all local rows) + - `dnnz::Vector{PetscInt}`: array containing the number of nonzeros in the various rows of the DIAGONAL portion of the local + submatrix (possibly different for each row) or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero + structure. The size of this array is equal to the number of local rows, i.e 'm'. For matrices that will be factored, you + must leave room for (and set) the diagonal entry even if it is zero. + - `onz::PetscInt`: number of nonzeros per row in the OFF-DIAGONAL portion of local submatrix (same value is used for all local rows). + - `onnz::Vector{PetscInt}`: array containing the number of nonzeros in the various rows of the OFF-DIAGONAL portion of the local + submatrix (possibly different for each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero structure. + The size of this array is equal to the number of local rows, i.e 'm'. +""" +function MPIAIJSetPreallocation( + B::Mat, + d_nz::PetscInt, + d_nnz::Vector{PetscInt}, + o_nz::PetscInt, + o_nnz::Vector{PetscInt}, +) + error = ccall( + (:MatMPIAIJSetPreallocation, libpetsc), + PetscErrorCode, + (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}), + B, + d_nz, + d_nnz, + o_nz, + o_nnz, + ) + @assert iszero(error) +end + +function MPIAIJSetPreallocation(mat::Mat, dnz::PetscInt, onz::PetscInt) + error = ccall( + (:MatMPIAIJSetPreallocation, libpetsc), + PetscErrorCode, + (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}), + mat, + dnz, + C_NULL, + onz, + C_NULL, + ) + @assert iszero(error) +end + +""" + mult(mat::Mat, x::Vec, y::Vec) + +Wrapper to `MatMult` +https://petsc.org/release/manualpages/Mat/MatMult/ + +Compute `y = Ax` +""" +function mult(mat::Mat, x::Vec, y::Vec) + error = ccall((:MatMult, libpetsc), PetscErrorCode, (CMat, CVec, CVec), mat, x, y) + @assert iszero(error) +end + +""" + multAdd(A::Mat, v1::Vec, v2::Vec, v3::Vec) + +Wrapper to `MatMultAdd` +https://petsc.org/release/manualpages/Mat/MatMultAdd/ + +Compute `v3 = v2 + A * v1`. +""" +function multAdd(A::Mat, v1::Vec, v2::Vec, v3::Vec) + error = ccall( + (:MatMultAdd, libpetsc), + PetscErrorCode, + (CMat, CVec, CVec, CVec), + A, + v1, + v2, + v3, + ) + @assert iszero(error) +end + +""" + SeqAIJSetPreallocation(mat::Mat, nz::PetscInt, nnz::Vector{PetscInt}) + SeqAIJSetPreallocation(mat::Mat, nz::PetscInt) + +Wrapper to `MatSeqAIJSetPreallocation` +https://petsc.org/release/manualpages/Mat/MatSeqAIJSetPreallocation/ +""" +function SeqAIJSetPreallocation(mat::Mat, nz::PetscInt, nnz::Vector{PetscInt}) + error = ccall( + (:MatSeqAIJSetPreallocation, libpetsc), + PetscErrorCode, + (CMat, PetscInt, Ptr{PetscInt}), + mat, + nz, + nnz, + ) + @assert iszero(error) +end + +function SeqAIJSetPreallocation(mat::Mat, nz::PetscInt) + error = ccall( + (:MatSeqAIJSetPreallocation, libpetsc), + PetscErrorCode, + (CMat, PetscInt, Ptr{PetscInt}), + mat, + nz, + C_NULL, + ) + @assert iszero(error) +end + +""" + setFromOptions(mat::Mat) + +Wrapper to `MatSetFromOptions` +https://petsc.org/release/manualpages/Mat/MatSetFromOptions/ +""" +function setFromOptions(mat::Mat) + error = ccall((:MatSetFromOptions, libpetsc), PetscErrorCode, (CMat,), mat) + @assert iszero(error) +end + +""" + setLocalToGlobalMapping( + x::Mat, + rmapping::ISLocalToGlobalMapping, + cmapping::ISLocalToGlobalMapping, + ) + +Wrapper to `MatSetLocalToGlobalMapping` +https://petsc.org/release/manualpages/Mat/MatSetLocalToGlobalMapping/ +""" +function setLocalToGlobalMapping( + x::Mat, + rmapping::ISLocalToGlobalMapping, + cmapping::ISLocalToGlobalMapping, +) + error = ccall( + (:MatSetLocalToGlobalMapping, libpetsc), + PetscErrorCode, + (CMat, CISLocalToGlobalMapping, CISLocalToGlobalMapping), + x, + rmapping, + cmapping, + ) + @assert iszero(error) +end + +""" + setOption(mat::Mat, op::MatOption, flg::PetscBool) + +Wrapper for `MatSetOption` +https://petsc.org/release/manualpages/Mat/MatSetOption/ +""" +function setOption(mat::Mat, op::MatOption, flg::PetscBool) + error = ccall( + (:MatSetOption, libpetsc), + PetscErrorCode, + (CMat, MatOption, Cuchar), + mat, + op, + flg, + ) + @assert iszero(error) +end + +setOption(mat::Mat, op::MatOption, flg::Bool) = setOption(mat, op, bool2petsc(flg)) + +""" + setSizes(mat::Mat, m::PetscInt, n::PetscInt, M::PetscInt, N::PetscInt) + setSizes( + mat::Mat, + nrows_loc::Integer, + ncols_loc::Integer, + nrows_glo::Integer, + ncols_glo::Integer, + ) + +Wrapper to `MatSetSizes` +https://petsc.org/release/manualpages/Mat/MatSetSizes/ +""" +function setSizes(mat::Mat, m::PetscInt, n::PetscInt, M::PetscInt, N::PetscInt) + error = ccall( + (:MatSetSizes, libpetsc), + PetscErrorCode, + (CMat, PetscInt, PetscInt, PetscInt, PetscInt), + mat, + m, + n, + M, + N, + ) + @assert iszero(error) +end + +function setSizes( + mat::Mat, + nrows_loc::Integer, + ncols_loc::Integer, + nrows_glo::Integer, + ncols_glo::Integer, +) + setSizes( + mat, + PetscInt(nrows_loc), + PetscInt(ncols_loc), + PetscInt(nrows_glo), + PetscInt(ncols_glo), + ) +end + +""" + setType(mat::Mat, type::String) + +Wrapper for `MatSetType` +https://petsc.org/release/manualpages/Mat/MatSetType/ + +Values for `type` alors available here: +https://petsc.org/release/manualpages/Mat/MatType.html#MatType +""" +function setType(mat::Mat, type::String) + error = ccall((:MatSetType, libpetsc), PetscErrorCode, (CMat, Cstring), mat, type) + @assert iszero(error) +end + +""" + setUp(mat::Mat) + +Wrapper to `MatSetUp` +https://petsc.org/release/manualpages/Mat/MatSetUp/ +""" +function setUp(mat::Mat) + error = ccall((:MatSetUp, libpetsc), PetscErrorCode, (CMat,), mat) + @assert iszero(error) +end + +# Avoid allocating an array of size 1 for each call to MatSetValue +const _ivec = zeros(PetscInt, 1) +const _jvec = zeros(PetscInt, 1) +const _vvec = zeros(PetscScalar, 1) + +""" + setValue( + m::Mat, + row::PetscInt, + col::PetscInt, + value::PetscScalar, + mode::InsertMode, + ) + setValue(m::Mat, row::Integer, col::Integer, value::Number, mode::InsertMode) + +Wrapper to `MatSetValue` +https://petsc.org/release/manualpages/Mat/MatSetValue/ + +Indexing starts at 0 (as in PETSc). + +# Implementation + +For an unknow reason, calling PETSc.MatSetValue leads to an "undefined symbol: MatSetValue" error. +So this wrapper directly call MatSetValues (anyway, this is what is done in PETSc...) +""" +function setValue( + m::Mat, + row::PetscInt, + col::PetscInt, + value::PetscScalar, + mode::InsertMode, +) + # Convert to arrays + _ivec[1] = row + _jvec[1] = col + _vvec[1] = value + + setValues(m, PetscIntOne, _ivec, PetscIntOne, _jvec, _vvec, mode) + #MatSetValues(mat, PetscIntOne, [i], PetscIntOne, [j], [v], mode) +end + +function setValue(m::Mat, row::Integer, col::Integer, value::Number, mode::InsertMode) + setValue(m, PetscInt(row), PetscInt(col), PetscScalar(value), mode) +end + +""" + setValues( + mat::Mat, + m::PetscInt, + idxm::Vector{PetscInt}, + n::PetscInt, + idxn::Vector{PetscInt}, + V::Array{PetscScalar}, + addv::InsertMode, + ) + setValues( + mat::Mat, + I::Vector{PetscInt}, + J::Vector{PetscInt}, + V::Array{PetscScalar}, + mode::InsertMode, + ) + setValues(mat::Mat, I, J, V, mode::InsertMode) + +Wrapper to `MatSetValues` +https://petsc.org/release/manualpages/Mat/MatSetValues/ + +Indexing starts at 0 (as in PETSc) +""" +function setValues( + mat::Mat, + m::PetscInt, + idxm::Vector{PetscInt}, + n::PetscInt, + idxn::Vector{PetscInt}, + V::Array{PetscScalar}, + addv::InsertMode, +) + error = ccall( + (:MatSetValues, libpetsc), + PetscErrorCode, + ( + CMat, + PetscInt, + Ptr{PetscInt}, + PetscInt, + Ptr{PetscInt}, + Ptr{PetscScalar}, + InsertMode, + ), + mat, + m, + idxm, + n, + idxn, + V, + addv, + ) + @assert iszero(error) +end + +function setValues( + mat::Mat, + I::Vector{PetscInt}, + J::Vector{PetscInt}, + V::Array{PetscScalar}, + mode::InsertMode, +) + setValues(mat, PetscInt(length(I)), I, PetscInt(length(J)), J, V, mode) +end + +function setValues(mat::Mat, I, J, V, mode::InsertMode) + setValues(mat, PetscInt.(I), PetscInt.(collect(J)), PetscScalar.(V), mode) +end + +""" + matView(mat::Mat, viewer::PetscViewer = StdWorld()) + +Wrapper to `MatView` +https://petsc.org/release/manualpages/Mat/MatView/ +""" +function matView(mat::Mat, viewer::PetscViewer = StdWorld(mat.comm)) + error = ccall((:MatView, libpetsc), PetscErrorCode, (CMat, CViewer), mat, viewer) + @assert iszero(error) +end + +""" + zeroEntries(mat::Mat) + +Wrapper to `MatZeroEntries` +https://petsc.org/release/manualpages/Mat/MatZeroEntries/ +""" +function zeroEntries(mat::Mat) + error = ccall((:MatZeroEntries, libpetsc), PetscErrorCode, (CMat,), mat) + @assert iszero(error) +end diff --git a/src/PetscViewer.jl b/src/PetscViewer.jl new file mode 100644 index 0000000..92dc50e --- /dev/null +++ b/src/PetscViewer.jl @@ -0,0 +1,204 @@ +# DEV NOTE : using `finalizer` for PetscViewer failed because for some reason I don't understand, +# calling `PetscViewerDestroy` often fails (even outside the context of `finalizer`) + +const CViewer = Ptr{Cvoid} +mutable struct PetscViewer <: AbstractPetscObject + ptr::CViewer + comm::MPI.Comm + + PetscViewer(comm::MPI.Comm) = new(CViewer(), comm) +end + +Base.unsafe_convert(::Type{CViewer}, x::PetscViewer) = x.ptr +function Base.unsafe_convert(::Type{Ptr{CViewer}}, x::PetscViewer) + Ptr{CViewer}(pointer_from_objref(x)) +end + +""" + ASCIIGetStdout(comm::MPI.Comm) + +Wrapper for `PetscViewerASCIIGetStdout` +https://petsc.org/release/manualpages/Viewer/PetscViewerASCIIGetStdout/ +""" +function ASCIIGetStdout(comm::MPI.Comm = MPI.COMM_WORLD) + viewer = PetscViewer(comm) + error = ccall( + (:PetscViewerASCIIGetStdout, libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, Ptr{CViewer}), + comm, + viewer, + ) + @assert iszero(error) + return viewer +end + +const StdWorld = ASCIIGetStdout + +""" + ASCIIOpen(comm::MPI.Comm, name) + +Wrapper for `PetscViewerASCIIOpen` +https://petsc.org/release/manualpages/Viewer/PetscViewerASCIIOpen/ +""" +function ASCIIOpen(comm::MPI.Comm, name::String) + viewer = PetscViewer(comm) + error = ccall( + (:PetscViewerASCIIOpen, libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, Cstring, Ptr{CViewer}), + comm, + name, + viewer, + ) + @assert iszero(error) + return viewer +end + +""" + create(::Type{PetscViewer}, comm::MPI.Comm = MPI.COMM_WORLD) + +Wrapper for `PetscViewerCreate` +https://petsc.org/release/manualpages/Viewer/PetscViewerCreate/ +""" +function create(::Type{PetscViewer}, comm::MPI.Comm = MPI.COMM_WORLD) + viewer = PetscViewer(comm) + error = ccall( + (:PetscViewerCreate, libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, Ptr{CViewer}), + comm, + viewer, + ) + @assert iszero(error) + return viewer +end + +""" + destroy(viewer::PetscViewer) + +Wrapper for `PetscViewerDestroy` +https://petsc.org/release/manualpages/Viewer/PetscViewerDestroy/ + +Warning : from what I understand, all viewers must not be destroyed explicitely using `PetscViewerDestroy`. +""" +function destroy(viewer::PetscViewer) + _is_destroyed(viewer) && return + + error = ccall((:PetscViewerDestroy, libpetsc), PetscErrorCode, (Ptr{CViewer},), viewer) + @assert iszero(error) +end + +""" + fileSetMode(viewer::PetscViewer, mode::PetscFileMode) + +Wrapper for `PetscViewerFileSetMode` +https://petsc.org/release/manualpages/Viewer/PetscViewerFileSetMode/ +""" +function fileSetMode(viewer::PetscViewer, mode::PetscFileMode) + error = ccall( + (:PetscViewerFileSetMode, libpetsc), + PetscErrorCode, + (CViewer, PetscFileMode), + viewer, + mode, + ) + @assert iszero(error) +end + +""" + fileSetName(viewer::PetscViewer, name::String) + +Wrapper for `PetscViewerFileSetName` +https://petsc.org/release/manualpages/Viewer/PetscViewerFileSetName/ +""" +function fileSetName(viewer::PetscViewer, name::String) + error = ccall( + (:PetscViewerFileSetName, libpetsc), + PetscErrorCode, + (CViewer, Cstring), + viewer, + name, + ) + @assert iszero(error) +end + +""" + HDF5Open(comm::MPI.Comm, name::String, type::PetscFileMode) + +Wrapper for `PetscViewerHDF5Open` +https://petsc.org/release/manualpages/Viewer/PetscViewerHDF5Open/ +""" +function HDF5Open(comm::MPI.Comm, name::String, type::PetscFileMode) + viewer = PetscViewer(comm) + error = ccall( + (:PetscViewerHDF5Open, libpetsc), + PetscErrorCode, + (MPI.MPI_Comm, Cstring, PetscFileMode, Ptr{CViewer}), + comm, + name, + type, + viewer, + ) + @assert iszero(error) + return viewer +end + +""" + popFormat(viewer::PetscViewer) + +Wrapper for `PetscViewerPopFormat` +https://petsc.org/release/manualpages/Viewer/PetscViewerPopFormat/ +""" +function popFormat(viewer::PetscViewer) + error = ccall((:PetscViewerPopFormat, libpetsc), PetscErrorCode, (CViewer,), viewer) + @assert iszero(error) +end + +""" + pushFormat(viewer::PetscViewer, format::PetscViewerFormat) + +Wrapper for `PetscViewerPushFormat` +https://petsc.org/release/manualpages/Viewer/PetscViewerPushFormat/ +""" +function pushFormat(viewer::PetscViewer, format::PetscViewerFormat) + error = ccall( + (:PetscViewerPushFormat, libpetsc), + PetscErrorCode, + (CViewer, PetscViewerFormat), + viewer, + format, + ) + @assert iszero(error) +end + +""" + setType(viewer::PetscViewer, type::String) + +Wrapper for `PetscViewerSetType` +https://petsc.org/release/manualpages/Viewer/PetscViewerSetType/ + +Values for `type` alors available here: +https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Viewer/PetscViewerType.html#PetscViewerType +""" +function setType(viewer::PetscViewer, type::String) + error = ccall( + (:PetscViewerSetType, libpetsc), + PetscErrorCode, + (CViewer, Cstring), + viewer, + type, + ) + @assert iszero(error) +end + +""" + viewerView(v::PetscViewer, viewer::PetscViewer) + +Wrapper to `PetscViewerView` +""" +function viewerView(v::PetscViewer, viewer::PetscViewer) + error = + ccall((:PetscViewerView, libpetsc), PetscErrorCode, (CViewer, CViewer), v, viewer) + @assert iszero(error) +end diff --git a/src/PetscWrap.jl b/src/PetscWrap.jl index 2dc9b20..13b017d 100644 --- a/src/PetscWrap.jl +++ b/src/PetscWrap.jl @@ -7,16 +7,22 @@ module PetscWrap using Libdl using MPI +using LinearAlgebra + +# GridapPETSc trick to track allocs/deallocs +const _NREFS = Ref(0) include("const_arch_ind.jl") -export PetscErrorCode, PETSC_DECIDE +export PetscErrorCode, PETSC_DECIDE # export all items of some enums for item in Iterators.flatten(( instances(InsertMode), instances(MatAssemblyType), + instances(MatOption), + instances(MatDuplicateOption), instances(PetscViewerFormat), instances(PetscFileMode), - )) +)) @eval export $(Symbol(item)) end @@ -26,108 +32,96 @@ export show_petsc_path include("const_arch_dep.jl") export PetscReal, PetscScalar, PetscInt, PetscIntOne +include("common.jl") + include("init.jl") -export PetscInitialize, PetscFinalize - -include("viewer.jl") -export PetscViewer, CViewer, - PetscViewerASCIIOpen, - PetscViewerCreate, - PetscViewerDestroy, - PetscViewerFileSetMode, - PetscViewerFileSetName, - PetscViewerHDF5Open, - PetscViewerPopFormat, - PetscViewerPushFormat, - PetscViewerSetType, - PetscViewerStdWorld, - PetscViewerView - -include("vec.jl") -export PetscVec, CVec, - VecAssemble, - VecAssemblyBegin, - VecAssemblyEnd, - VecCreate, - VecDestroy, - VecDuplicate, - VecGetArray, - VecGetLocalSize, - VecGetOwnershipRange, - VecGetSize, - VecRestoreArray, - VecSetFromOptions, - VecSetSizes, - VecSetUp, - VecSetValue, - VecSetValues, - VecView - -include("mat.jl") -export PetscMat, CMat, - MatAssemble, - MatAssemblyBegin, - MatAssemblyEnd, - MatCreate, - MatCreateComposite, - MatCreateDense, - MatCreateVecs, - MatDestroy, - MatGetLocalSize, - MatGetOwnershipRange, - MatGetOwnershipRangeColumn, - MatGetSize, - MatGetType, - MatMPIAIJSetPreallocation, - MatMult, - MatSeqAIJSetPreallocation, - MatSetFromOptions, - MatSetSizes, - MatSetUp, - MatSetValue, - MatSetValues, - MatView - -include("ksp.jl") -export PetscKSP, CKSP, - KSPCreate, - KSPDestroy, - KSPSetFromOptions, - KSPSetOperators, - KSPSetUp, - KSPSolve +export PetscInitialize, PetscInitialized, PetscFinalize + +include("PetscViewer.jl") +export PetscViewer, + CViewer, + ASCIIOpen, + create, + destroy, + fileSetMode, + fileSetName, + HDF5Open, + popFormat, + pushFormat, + setType, + stdWorld, + viewerView + +include("ISLocalToGlobalMapping.jl") +include("Vec.jl") +export Vec, + CVec, + assemblyBegin, + assemblyEnd, + copy, + duplicate, + getArray, + getLocalSize, + getOwnershipRange, + getSize, + getValues, + restoreArray, + scale, + setFromOptions, + setSizes, + setUp, + setValue, + setValueLocal, + setValues, + setValuesLocal, + vecView + +include("Mat.jl") +export Mat, + CMat, + createComposite, + createDense, + createVecs, + getOwnershipRangeColumn, + getType, + MPIAIJSetPreallocation, + mult, + multAdd, + SeqAIJSetPreallocation, + matView, + zeroEntries + +include("KSP.jl") +export KSP, CKSP, setOperators, solve # fancy +include("fancy/common.jl") include("fancy/viewer.jl") -export destroy!, - push_format!, - set_mode!, - set_name!, - set_type! +export destroy!, push_format!, set_mode!, set_name!, set_type! include("fancy/vec.jl") -export assemble!, - create_vector, - destroy!, - duplicate, - get_range, - get_urange, - set_local_size!, set_global_size!, - set_from_options!, - set_up!, - vec2array, - vec2file +export assemble!, + create_vector, + destroy!, + duplicate, + get_range, + get_urange, + set_local_size!, + set_local_to_global!, + set_global_size!, + set_from_options!, + set_up!, + set_value!, + set_value_local!, + set_values!, + set_values_local!, + vec2array, + vec2file include("fancy/mat.jl") -export create_composite_add, - create_matrix, - mat2file, - preallocate!, - set_values! +export create_composite_add, create_matrix, mat2file, preallocate! include("fancy/ksp.jl") -export create_ksp, - solve, solve!, - set_operators! +export create_ksp, solve, solve!, set_operators! -end \ No newline at end of file +end diff --git a/src/Vec.jl b/src/Vec.jl new file mode 100644 index 0000000..f3e6caf --- /dev/null +++ b/src/Vec.jl @@ -0,0 +1,504 @@ +const CVec = Ptr{Cvoid} + +mutable struct Vec <: AbstractPetscObject + ptr::CVec + comm::MPI.Comm + + Vec(comm::MPI.Comm) = new(CVec(), comm) +end + +Base.unsafe_convert(::Type{CVec}, x::Vec) = x.ptr +Base.unsafe_convert(::Type{Ptr{CVec}}, x::Vec) = Ptr{CVec}(pointer_from_objref(x)) + +""" + assemblyBegin(vec::Vec) + +Wrapper to `VecAssemblyBegin` +https://petsc.org/release/manualpages/Vec/VecAssemblyBegin/ +""" +function assemblyBegin(vec::Vec) + error = ccall((:VecAssemblyBegin, libpetsc), PetscErrorCode, (CVec,), vec) + @assert iszero(error) +end + +""" + assemblyEnd(vec::Vec) + +Wrapper to `VecAssemblyEnd` +https://petsc.org/release/manualpages/Vec/VecAssemblyEnd/ +""" +function assemblyEnd(vec::Vec) + error = ccall((:VecAssemblyEnd, libpetsc), PetscErrorCode, (CVec,), vec) + @assert iszero(error) +end + +""" + copy(x::Vec, y::Vec) + copy(x::Vec; add_finalizer = true) + +Wrapper to `VecCopy` +https://petsc.org/release/manualpages/Vec/VecCopy/ + +Not really sure if I should extend Base here... +""" +function Base.copy(x::Vec, y::Vec) + error = ccall((:VecCopy, libpetsc), PetscErrorCode, (CVec, CVec), x, y) + @assert iszero(error) +end + +function Base.copy(x::Vec; add_finalizer = true) + y = duplicate(x) + copy(x, y) + _NREFS[] += 1 + add_finalizer && finalizer(destroy, y) + return y +end + +""" + create(::Type{Vec}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + +Wrapper to `VecCreate` +https://petsc.org/release/manualpages/Vec/VecCreate/ +""" +function create(::Type{Vec}, comm::MPI.Comm = MPI.COMM_WORLD; add_finalizer = true) + vec = Vec(comm) + error = + ccall((:VecCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CVec}), comm, vec) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, vec) + + return vec +end + +""" + destroy(v::Vec) + +Wrapper to `VecDestroy` +https://petsc.org/release/manualpages/Vec/VecDestroy/ +""" +function destroy(v::Vec) + _is_destroyed(v) && return + + error = ccall((:VecDestroy, libpetsc), PetscErrorCode, (Ptr{CVec},), v) + @assert iszero(error) + + _NREFS[] -= 1 +end + +function LinearAlgebra.dot(x::Vec, y::Vec) + val = Ref{PetscScalar}() + error = ccall( + (:VecDot, libpetsc), + PetscErrorCode, + (CVec, CVec, Ref{PetscScalar}), + x, + y, + val, + ) + @assert iszero(error) + + return val[] +end + +""" + duplicate(v::Vec; add_finalizer = true) + +Wrapper for `VecDuplicate` +https://petsc.org/release/manualpages/Vec/VecDuplicate/ +""" +function duplicate(v::Vec; add_finalizer = true) + newv = Vec(v.comm) + error = ccall((:VecDuplicate, libpetsc), PetscErrorCode, (CVec, Ptr{CVec}), v, newv) + @assert iszero(error) + + _NREFS[] += 1 + add_finalizer && finalizer(destroy, newv) + + return newv +end + +""" + getArray(x::Vec, own::Bool = false) + +Wrapper for `VecGetArray` +https://petsc.org/release/manualpages/Vec/VecGetArray/ + +# Warning + +I am not confortable at all with memory management, both on the C side and on the Julia side. Use +this at you own risk. + +According to Julia documentation, `own` optionally specifies whether Julia should take ownership +of the memory, calling free on the pointer when the array is no longer referenced." +""" +function getArray(x::Vec, own::Bool = false) + # Get array pointer + array_ref = Ref{Ptr{PetscScalar}}() + error = ccall( + (:VecGetArray, libpetsc), + PetscErrorCode, + (CVec, Ref{Ptr{PetscScalar}}), + x, + array_ref, + ) + @assert iszero(error) + + # Get array size + rstart, rend = getOwnershipRange(x) + n = rend - rstart # this is not `rend - rstart + 1` because of VecGetOwnershipRange convention + + array = unsafe_wrap(Array, array_ref[], n; own = own) + + return array, array_ref +end + +""" + getLocalSize(vec::Vec) + +Wrapper for `VecGetLocalSize` +https://petsc.org/release/manualpages/Vec/VecGetLocalSize/ +""" +function getLocalSize(x::Vec) + n = Ref{PetscInt}() + + error = ccall((:VecGetLocalSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), x, n) + @assert iszero(error) + + return n[] +end + +""" + getOwnershipRange(x::Vec) + +Wrapper to `VecGetOwnershipRange` +https://petsc.org/release/manualpages/Vec/VecGetOwnershipRange/ + +The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor. + +# Warning + +`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row +handled by the local processor. +""" +function getOwnershipRange(x::Vec) + rstart = Ref{PetscInt}(0) + rend = Ref{PetscInt}(0) + + error = ccall( + (:VecGetOwnershipRange, libpetsc), + PetscErrorCode, + (CVec, Ref{PetscInt}, Ref{PetscInt}), + x, + rstart, + rend, + ) + @assert iszero(error) + + return rstart[], rend[] +end + +""" + getSize(x::Vec) + +Wrapper for `VecGetSize` +https://petsc.org/release/manualpages/Vec/VecGetSize/ +""" +function getSize(x::Vec) + n = Ref{PetscInt}() + + error = ccall((:VecGetSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), x, n) + @assert iszero(error) + + return n[] +end + +""" + getValues(x::Vec, ni::PetscInt, ix::Vector{PetscInt}) + getValues(x::Vec, ix::Vector{PetscInt}) + getValues(x::Vec, ix::Vector{I}) where {I<:Integer} + +Wrapper for `VecGetValues` +https://petsc.org/release/manualpages/Vec/VecGetValues/ +""" +function getValues(x::Vec, ni::PetscInt, ix::Vector{PetscInt}) + y = zeros(PetscScalar, ni) + + error = ccall( + (:VecGetValues, libpetsc), + PetscErrorCode, + (CVec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}), + x, + ni, + ix, + y, + ) + @assert iszero(error) + + return y +end + +getValues(x::Vec, ix::Vector{PetscInt}) = getValues(x, PetscInt(length(ix)), ix) + +getValues(x::Vec, ix::Vector{I}) where {I<:Integer} = getValues(x, PetscInt.(ix)) + +""" + restoreArray(x::Vec, array_ref) + +Wrapper for `VecRestoreArray`. `array_ref` is obtained from `VecGetArray` +https://petsc.org/release/manualpages/Vec/VecRestoreArray/ +""" +function restoreArray(x::Vec, array_ref) + error = ccall( + (:VecRestoreArray, libpetsc), + PetscErrorCode, + (CVec, Ref{Ptr{PetscScalar}}), + x, + array_ref, + ) + @assert iszero(error) +end + +""" + scale(x::Vec, alpha::PetscScalar) + scale(x::Vec, alpha::Number) + +Wrapper for `VecScale` +https://petsc.org/release/manualpages/Vec/VecScale/ +""" +function scale(x::Vec, alpha::PetscScalar) + error = ccall((:VecScale, libpetsc), PetscErrorCode, (CVec, PetscScalar), x, alpha) + @assert iszero(error) +end + +scale(x::Vec, alpha::Number) = scale(x, PetscScalar(alpha)) + +""" + setFromOptions(vec::Vec) + +Wrapper to `VecSetFromOptions` +https://petsc.org/release/manualpages/Vec/VecSetFromOptions/ +""" +function setFromOptions(vec::Vec) + error = ccall((:VecSetFromOptions, libpetsc), PetscErrorCode, (CVec,), vec) + @assert iszero(error) +end + +""" + setLocalToGlobalMapping(x::Vec, mapping::ISLocalToGlobalMapping) + +Wrapper to `VecSetLocalToGlobalMapping` +https://petsc.org/release/manualpages/Vec/VecSetLocalToGlobalMapping/ + +0-based indexing +""" +function setLocalToGlobalMapping(x::Vec, mapping::ISLocalToGlobalMapping) + error = ccall( + (:VecSetLocalToGlobalMapping, libpetsc), + PetscErrorCode, + (CVec, CISLocalToGlobalMapping), + x, + mapping, + ) + @assert iszero(error) +end + +""" + setSizes(v::Vec, n::PetscInt, N::PetscInt) + setSizes(v::Vec, n::Integer, N::Integer) + +Wrapper to `VecSetSizes` +https://petsc.org/release/manualpages/Vec/VecSetSizes/ +""" +function setSizes(v::Vec, n::PetscInt, N::PetscInt) + nr_loc = PetscInt(n) + nr_glo = PetscInt(N) + error = ccall( + (:VecSetSizes, libpetsc), + PetscErrorCode, + (CVec, PetscInt, PetscInt), + v, + nr_loc, + nr_glo, + ) + @assert iszero(error) +end + +setSizes(v::Vec, n::Integer, N::Integer) = setSizes(v, PetscInt(n), PetscInt(N)) + +""" + setUp(vec::Vec) + +Wrapper to `VecSetUp` +https://petsc.org/release/manualpages/Vec/VecSetUp/ +""" +function setUp(v::Vec) + error = ccall((:VecSetUp, libpetsc), PetscErrorCode, (CVec,), v) + @assert iszero(error) +end + +""" + setValue(v::Vec, row::PetscInt, value::PetscScalar, mode::InsertMode) + setValue(v::Vec, row, value, mode::InsertMode = INSERT_VALUES) + +Wrapper to `setValue`. Indexing starts at 0 (as in PETSc). +https://petsc.org/release/manualpages/Vec/VecSetValue/ + +# Implementation + +For an unknow reason, calling PETSc.VecSetValue leads to an "undefined symbol: VecSetValue" error. +So this wrapper directly call VecSetValues (anyway, this is what is done in PETSc...) +""" +function setValue(v::Vec, row::PetscInt, value::PetscScalar, mode::InsertMode) + setValues(v, PetscIntOne, [row], [value], mode) +end + +function setValue(v::Vec, row, value, mode::InsertMode = INSERT_VALUES) + setValue(v, PetscInt(row), PetscScalar(value), mode) +end + +""" + setValues( + x::Vec, + ni::PetscInt, + ix::Vector{PetscInt}, + y::Vector{PetscScalar}, + iora::InsertMode, + ) + + setValues( + x::Vec, + I::Vector{PetscInt}, + V::Vector{PetscScalar}, + mode::InsertMode = INSERT_VALUES, + ) + setValues(x::Vec, I, V, mode::InsertMode = INSERT_VALUES) + +Wrapper to `VecSetValues`. Indexing starts at 0 (as in PETSc) +https://petsc.org/release/manualpages/Vec/VecSetValues/ +""" +function setValues( + x::Vec, + ni::PetscInt, + ix::Vector{PetscInt}, + y::Vector{PetscScalar}, + iora::InsertMode, +) + error = ccall( + (:VecSetValues, libpetsc), + PetscErrorCode, + (CVec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), + x, + ni, + ix, + y, + iora, + ) + @assert iszero(error) +end + +function setValues( + x::Vec, + I::Vector{PetscInt}, + V::Vector{PetscScalar}, + mode::InsertMode = INSERT_VALUES, +) + setValues(x, PetscInt(length(I)), I, V, mode) +end + +function setValues(x::Vec, I, V, mode::InsertMode = INSERT_VALUES) + setValues(x, PetscInt.(I), PetscScalar.(V), mode) +end + +""" + setValueLocal(v::Vec, row::PetscInt, value::PetscScalar, mode::InsertMode) + setValueLocal(v::Vec, row, value, mode::InsertMode = INSERT_VALUES) + +Wrapper to `setValueLocal`. Indexing starts at 0 (as in PETSc). +https://petsc.org/release/manualpages/Vec/VecSetValueLocal/ + +# Implementation + +For an unknow reason, calling PETSc.VecSetValue leads to an "undefined symbol: VecSetValue" error. +So this wrapper directly call VecSetValues (anyway, this is what is done in PETSc...) +""" +function setValueLocal(v::Vec, row::PetscInt, value::PetscScalar, mode::InsertMode) + setValuesLocal(v, PetscIntOne, [row], [value], mode) +end + +function setValueLocal(v::Vec, row, value, mode::InsertMode = INSERT_VALUES) + setValueLocal(v, PetscInt(row), PetscScalar(value), mode) +end + +""" + setValuesLocal( + x::Vec, + ni::PetscInt, + ix::Vector{PetscInt}, + y::Vector{PetscScalar}, + iora::InsertMode, + ) + + setValuesLocal( + x::Vec, + I::Vector{PetscInt}, + V::Vector{PetscScalar}, + mode::InsertMode = INSERT_VALUES, + ) + setValues(x::Vec, I, V, mode::InsertMode = INSERT_VALUES) + +Wrapper to `VecSetValues`. Indexing starts at 0 (as in PETSc) +https://petsc.org/release/manualpages/Vec/VecSetValuesLocal/ +""" +function setValuesLocal( + x::Vec, + ni::PetscInt, + ix::Vector{PetscInt}, + y::Vector{PetscScalar}, + iora::InsertMode, +) + error = ccall( + (:VecSetValuesLocal, libpetsc), + PetscErrorCode, + (CVec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), + x, + ni, + ix, + y, + iora, + ) + @assert iszero(error) +end + +function setValuesLocal( + x::Vec, + I::Vector{PetscInt}, + V::Vector{PetscScalar}, + mode::InsertMode = INSERT_VALUES, +) + setValuesLocal(x, PetscInt(length(I)), I, V, mode) +end + +function setValuesLocal(x::Vec, I, V, mode::InsertMode = INSERT_VALUES) + setValuesLocal(x, PetscInt.(I), PetscScalar.(V), mode) +end + +function Base.sum(x::Vec) + s = Ref{PetscScalar}() + error = ccall((:VecSum, libpetsc), PetscErrorCode, (CVec, Ref{PetscScalar}), x, s) + @assert iszero(error) + + return s[] +end + +""" + vecView(vec::Vec, viewer::PetscViewer = StdWorld()) + +Wrapper to `VecView` +https://petsc.org/release/manualpages/Vec/VecView/ +""" +function vecView(vec::Vec, viewer::PetscViewer = StdWorld(vec.comm)) + error = ccall((:VecView, libpetsc), PetscErrorCode, (CVec, CViewer), vec, viewer) + @assert iszero(error) +end diff --git a/src/common.jl b/src/common.jl new file mode 100644 index 0000000..92df969 --- /dev/null +++ b/src/common.jl @@ -0,0 +1,12 @@ +abstract type AbstractPetscObject end +_get_ptr(obj::AbstractPetscObject) = obj.ptr + +const CPetscObject = Ptr{Cvoid} + +function objectRegisterDestroy(obj::AbstractPetscObject) + error = + ccall((:PetscObjectRegisterDestroy, libpetsc), PetscErrorCode, (CPetscObject,), obj) + @assert iszero(error) +end + +_is_destroyed(obj::AbstractPetscObject) = _get_ptr(obj) == C_NULL diff --git a/src/const_arch_dep.jl b/src/const_arch_dep.jl index 7d174f2..dd5d9e2 100644 --- a/src/const_arch_dep.jl +++ b/src/const_arch_dep.jl @@ -7,9 +7,14 @@ function DataTypeFromString(name::AbstractString) dtype_ref = Ref{PetscDataType}() found_ref = Ref{PetscBool}() - ccall((:PetscDataTypeFromString, libpetsc), PetscErrorCode, - (Cstring, Ptr{PetscDataType}, Ptr{PetscBool}), - name, dtype_ref, found_ref) + ccall( + (:PetscDataTypeFromString, libpetsc), + PetscErrorCode, + (Cstring, Ptr{PetscDataType}, Ptr{PetscBool}), + name, + dtype_ref, + found_ref, + ) @assert found_ref[] == PETSC_TRUE return dtype_ref[] end @@ -19,9 +24,13 @@ end """ function PetscDataTypeGetSize(dtype::PetscDataType) datasize_ref = Ref{Csize_t}() - ccall((:PetscDataTypeGetSize, libpetsc), PetscErrorCode, - (PetscDataType, Ptr{Csize_t}), - dtype, datasize_ref) + ccall( + (:PetscDataTypeGetSize, libpetsc), + PetscErrorCode, + (PetscDataType, Ptr{Csize_t}), + dtype, + datasize_ref, + ) return datasize_ref[] end @@ -76,4 +85,4 @@ const PetscReal = PetscReal2Type() const PetscScalar = PetscScalar2Type() const PetscInt = PetscInt2Type() -const PetscIntOne = PetscInt(1) # Integer `1` with the type of PetscInt, usefull to go back and forth with julia/petsc indexing \ No newline at end of file +const PetscIntOne = PetscInt(1) # Integer `1` with the type of PetscInt, usefull to go back and forth with julia/petsc indexing diff --git a/src/const_arch_ind.jl b/src/const_arch_ind.jl index 2c5ea75..eab3cad 100644 --- a/src/const_arch_ind.jl +++ b/src/const_arch_ind.jl @@ -4,7 +4,8 @@ const PetscErrorCode = Cint const PETSC_DEFAULT = -2 -const PETSC_DECIDE = -1 +const PETSC_DECIDE = -1 +const PETSC_DETERMINE = PETSC_DECIDE """ Macro to automatically export all items of an enum @@ -14,34 +15,41 @@ macro exported_enum(name, args...) @enum($name, $(args...)) export $name $([:(export $arg) for arg in args]...) - end) + end) end - @enum PetscBool PETSC_FALSE PETSC_TRUE +Base.Bool(x::PetscBool) = x == PETSC_TRUE +bool2petsc(x::Bool) = x ? PETSC_TRUE : PETSC_FALSE + +@enum PetscCopyMode begin + PETSC_COPY_VALUES + PETSC_OWN_POINTER + PETSC_USE_POINTER +end @enum PetscDataType begin - PETSC_DATATYPE_UNKNOWN = 0 - PETSC_DOUBLE = 1 - PETSC_COMPLEX = 2 - PETSC_LONG = 3 - PETSC_SHORT = 4 - PETSC_FLOAT = 5 - PETSC_CHAR = 6 - PETSC_BIT_LOGICAL = 7 - PETSC_ENUM = 8 - PETSC_BOOL = 9 - PETSC_FLOAT128 = 10 - PETSC_OBJECT = 11 - PETSC_FUNCTION = 12 - PETSC_STRING = 13 - PETSC___FP16 = 14 - PETSC_STRUCT = 15 - PETSC_INT = 16 - PETSC_INT64 = 17 + PETSC_DATATYPE_UNKNOWN = 0 + PETSC_DOUBLE = 1 + PETSC_COMPLEX = 2 + PETSC_LONG = 3 + PETSC_SHORT = 4 + PETSC_FLOAT = 5 + PETSC_CHAR = 6 + PETSC_BIT_LOGICAL = 7 + PETSC_ENUM = 8 + PETSC_BOOL = 9 + PETSC_FLOAT128 = 10 + PETSC_OBJECT = 11 + PETSC_FUNCTION = 12 + PETSC_STRING = 13 + PETSC___FP16 = 14 + PETSC_STRUCT = 15 + PETSC_INT = 16 + PETSC_INT64 = 17 end -const Petsc64bitInt = Int64 +const Petsc64bitInt = Int64 const PetscLogDouble = Cdouble @enum InsertMode begin @@ -70,12 +78,12 @@ end end @enum MatFactorType begin - MAT_FACTOR_NONE = 0 - MAT_FACTOR_LU = 1 + MAT_FACTOR_NONE = 0 + MAT_FACTOR_LU = 1 MAT_FACTOR_CHOLESKY = 2 - MAT_FACTOR_ILU = 3 - MAT_FACTOR_ICC = 4 - MAT_FACTOR_ILUDT = 5 + MAT_FACTOR_ILU = 3 + MAT_FACTOR_ICC = 4 + MAT_FACTOR_ILUDT = 5 end @enum MatOption begin @@ -109,6 +117,12 @@ end MAT_OPTION_MAX = 25 end +@enum MatDuplicateOption begin + MAT_DO_NOT_COPY_VALUES + MAT_COPY_VALUES + MAT_SHARE_NONZERO_PATTERN +end + @enum PetscViewerFormat begin PETSC_VIEWER_DEFAULT PETSC_VIEWER_ASCII_MATLAB @@ -159,15 +173,15 @@ end end @enum MatOperation begin - MATOP_SET_VALUES=0 - MATOP_GET_ROW=1 - MATOP_RESTORE_ROW=2 - MATOP_MULT=3 - MATOP_MULT_ADD=4 - MATOP_MULT_TRANSPOSE=5 - MATOP_MULT_TRANSPOSE_ADD=6 - MATOP_SOLVE=7 - MATOP_SOLVE_ADD=8 - MATOP_SOLVE_TRANSPOSE=9 - MATOP_SOLVE_TRANSPOSE_ADD=10 -end \ No newline at end of file + MATOP_SET_VALUES = 0 + MATOP_GET_ROW = 1 + MATOP_RESTORE_ROW = 2 + MATOP_MULT = 3 + MATOP_MULT_ADD = 4 + MATOP_MULT_TRANSPOSE = 5 + MATOP_MULT_TRANSPOSE_ADD = 6 + MATOP_SOLVE = 7 + MATOP_SOLVE_ADD = 8 + MATOP_SOLVE_TRANSPOSE = 9 + MATOP_SOLVE_TRANSPOSE_ADD = 10 +end diff --git a/src/fancy/common.jl b/src/fancy/common.jl new file mode 100644 index 0000000..30d1d4b --- /dev/null +++ b/src/fancy/common.jl @@ -0,0 +1,3 @@ +destroy!(obj::Vararg{<:AbstractPetscObject,N}) where {N} = destroy.(obj) +set_up!(obj::Vararg{<:AbstractPetscObject,N}) where {N} = setUp.(obj) +set_from_options!(obj::Vararg{<:AbstractPetscObject,N}) where {N} = setFromOptions.(obj) diff --git a/src/fancy/ksp.jl b/src/fancy/ksp.jl index a7b3d9f..0577f65 100644 --- a/src/fancy/ksp.jl +++ b/src/fancy/ksp.jl @@ -1,34 +1,25 @@ -create_ksp() = KSPCreate() - -function create_ksp(A::PetscMat) - ksp = KSPCreate() - KSPSetOperators(ksp, A, A) - return ksp -end - -function create_ksp(Amat::PetscMat, Pmat::PetscMat) - ksp = KSPCreate() - KSPSetOperators(ksp, Amat, Pmat) +function create_ksp(Amat::Mat, Pmat::Mat; autosetup = false) + ksp = create(KSP) + setOperators(ksp, Amat, Pmat) + + if autosetup + set_from_options!(ksp) + set_up!(ksp) + end return ksp end -set_operators!(ksp::PetscKSP, Amat::PetscMat) = KSPSetOperators(ksp, Amat, Amat) -set_operators!(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat) = KSPSetOperators(ksp, Amat, Pmat) +create_ksp(A::Mat; autosetup = false) = create_ksp(A, A; autosetup) +set_operators!(ksp::KSP, Amat::Mat) = setOperators(ksp, Amat, Amat) +set_operators!(ksp::KSP, Amat::Mat, Pmat::Mat) = setOperators(ksp, Amat, Pmat) -solve!(ksp::PetscKSP, b::PetscVec, x::PetscVec) = KSPSolve(ksp, b, x) +solve!(ksp::KSP, b::Vec, x::Vec) = solve(ksp, b, x) -function solve(ksp::PetscKSP, b::PetscVec) - x = VecDuplicate(b) - KSPSolve(ksp, b, x) +function solve(ksp::KSP, b::Vec) + x = duplicate(b) + solve(ksp, b, x) return x end -set_up!(ksp::PetscKSP) = KSPSetUp(ksp) - -set_from_options!(ksp::PetscKSP) = KSPSetFromOptions(ksp) - - -Base.show(::IO, ksp::PetscKSP) = KSPView(ksp) - -destroy!(ksp::PetscKSP) = KSPDestroy(ksp) \ No newline at end of file +Base.show(::IO, ksp::KSP) = KSPView(ksp) diff --git a/src/fancy/mat.jl b/src/fancy/mat.jl index 227ce04..3e2af46 100644 --- a/src/fancy/mat.jl +++ b/src/fancy/mat.jl @@ -2,30 +2,47 @@ `row` and `col` must be in [1,size(mat)], i.e indexing starts at 1 (Julia). # Implementation + For some unkwnown reason, calling `MatSetValue` fails. """ -function Base.setindex!(mat::PetscMat, value::Number, row::Integer, col::Integer) - MatSetValues(mat, PetscInt[row - 1], PetscInt[col - 1], PetscScalar[value], INSERT_VALUES) +function Base.setindex!(mat::Mat, value::Number, row::Integer, col::Integer) + setValues(mat, PetscInt[row - 1], PetscInt[col - 1], PetscScalar[value], INSERT_VALUES) end # This is stupid but I don't know how to do better yet -Base.setindex!(mat::PetscMat, values, row::Integer, cols) = MatSetValues(mat, [row-1], collect(cols) .- 1, values, INSERT_VALUES) -Base.setindex!(mat::PetscMat, values, rows, col::Integer) = MatSetValues(mat, collect(rows) .- 1, [col-1], values, INSERT_VALUES) +function Base.setindex!(mat::Mat, values, row::Integer, cols) + setValues(mat, [row - 1], collect(cols) .- 1, values, INSERT_VALUES) +end +function Base.setindex!(mat::Mat, values, rows, col::Integer) + setValues(mat, collect(rows) .- 1, [col - 1], values, INSERT_VALUES) +end -Base.ndims(::Type{PetscMat}) = 2 +Base.ndims(::Type{Mat}) = 2 """ - create_matrix(nrows, ncols, nrows_loc = PETSC_DECIDE, ncols_loc = PETSC_DECIDE; auto_setup = false) + create_matrix( + comm::MPI.Comm = MPI.COMM_WORLD; + nrows_loc = PETSC_DECIDE, + ncols_loc = PETSC_DECIDE, + nrows_glo = PETSC_DECIDE, + ncols_glo = PETSC_DECIDE, + autosetup = false, + ) -Create a `PetscMat` matrix of global size `(nrows, ncols)`. - -Use `auto_setup = true` to immediatly call `set_from_options!` and `set_up!`. +Use `autosetup = true` to immediatly call `set_from_options!` and `set_up!`. """ -function create_matrix(nrows, ncols, nrows_loc = PETSC_DECIDE, ncols_loc = PETSC_DECIDE; auto_setup = false, comm::MPI.Comm = MPI.COMM_WORLD) - mat = MatCreate() - MatSetSizes(mat::PetscMat, nrows_loc, ncols_loc, nrows, ncols) +function create_matrix( + comm::MPI.Comm = MPI.COMM_WORLD; + nrows_loc = PETSC_DECIDE, + ncols_loc = PETSC_DECIDE, + nrows_glo = PETSC_DETERMINE, + ncols_glo = PETSC_DETERMINE, + autosetup = false, +) + mat = create(Mat, comm) + setSizes(mat, nrows_loc, ncols_loc, nrows_glo, ncols_glo) - if (auto_setup) + if (autosetup) set_from_options!(mat) set_up!(mat) end @@ -36,27 +53,45 @@ end Wrapper to `MatCreateComposite` using the "alternative construction" from the PETSc documentation. """ function create_composite_add(matrices) - N, M = MatGetSize(matrices[1]) - n, m = MatGetLocalSize(matrices[1]) - mat = create_matrix(N, M, n, m; auto_setup = false, comm = matrices[1].comm) - MatSetType(mat, "composite") - for m in matrices - MatCompositeAddMat(mat, m) + M, N = getSize(matrices[1]) + m, n = getLocalSize(matrices[1]) + mat = create_matrix( + matrices[1].comm; + nrows_loc = m, + ncols_loc = n, + nrows_glo = M, + ncols_glo = N, + autosetup = false, + ) + setType(mat, "composite") + for _mat in matrices + compositeAddMat(mat, _mat) end assemble!(mat) return mat end +function set_global_size!(mat::Mat, nrows, ncols) + setSizes(mat, PETSC_DECIDE, PETSC_DECIDE, nrows, ncols) +end -set_global_size!(mat::PetscMat, nrows, ncols) = MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, nrows, ncols) -set_local_size!(mat::PetscMat, nrows, ncols) = MatSetSizes(mat, nrows, ncols, PETSC_DECIDE, PETSC_DECIDE) - -set_from_options!(mat::PetscMat) = MatSetFromOptions(mat) +function set_local_to_global!( + mat::Mat, + rlid2gid::Vector{I}, + clid2gid::Vector{I}, +) where {I<:Integer} + rmapping = create(ISLocalToGlobalMapping, rlid2gid) + cmapping = create(ISLocalToGlobalMapping, clid2gid) + setLocalToGlobalMapping(mat, rmapping, cmapping) + destroy.(rmapping, cmapping) +end -set_up!(mat::PetscMat) = MatSetUp(mat) +function set_local_size!(mat::Mat, nrows, ncols) + setSizes(mat, nrows, ncols, PETSC_DECIDE, PETSC_DECIDE) +end """ - get_range(mat::PetscMat) + get_range(mat::Mat) Wrapper to `MatGetOwnershipRange` @@ -64,71 +99,115 @@ However, the result `(rstart, rend)` is such that `mat[rstart:rend]` are the row This is different from the default `PETSc.MatGetOwnershipRange` result where the indexing starts at zero and where `rend-1` is last row handled by the local processor. """ -function get_range(mat::PetscMat) - rstart, rend = MatGetOwnershipRange(mat) +function get_range(mat::Mat) + rstart, rend = getOwnershipRange(mat) return (rstart + 1, rend) end """ - get_urange(mat::PetscMat) + get_urange(mat::Mat) Provide a `UnitRange` from the method `get_range`. """ -function get_urange(mat::PetscMat) - rstart, rend = MatGetOwnershipRange(mat) - return rstart+1:rend +function get_urange(mat::Mat) + rstart, rend = getOwnershipRange(mat) + return (rstart + 1):rend end """ Wrapper to `MatAssemblyBegin` and `MatAssemblyEnd` successively. """ -function assemble!(mat::PetscMat, type::MatAssemblyType = MAT_FINAL_ASSEMBLY) - MatAssemblyBegin(mat, type) - MatAssemblyEnd(mat, type) +function assemble!(mat::Mat, type::MatAssemblyType = MAT_FINAL_ASSEMBLY) + assemblyBegin(mat, type) + assemblyEnd(mat, type) +end + +""" + set_value!(mat::Mat, I, J, V, mode = ADD_VALUES) + +Set value of `mat` +`mat[i, j] = v`. + +1-based indexing +""" +function set_value!(mat::Mat, i::PetscInt, j::PetscInt, v::PetscScalar, mode = ADD_VALUES) + setValue(mat, i - 1, j - 1, v, mode) +end +function set_value!(mat, i, j, v, mode = ADD_VALUES) + set_value!(mat, PetscInt(i), PetscInt(j), PetscScalar(v), mode) end """ - set_values!(mat::PetscMat, I, J, V, mode = ADD_VALUES) + set_values!(mat::Mat, I, J, V, mode = ADD_VALUES) Set values of `mat` in `SparseArrays` fashion : using COO format: `mat[I[k], J[k]] = V[k]`. """ -function set_values!(mat::PetscMat, I::Vector{PetscInt}, J::Vector{PetscInt}, V::Vector{PetscScalar}, mode = ADD_VALUES) - for (i,j,v) in zip(I, J, V) - MatSetValue(mat, i - PetscIntOne, j - PetscIntOne, v, mode) +function set_values!( + mat::Mat, + I::Vector{PetscInt}, + J::Vector{PetscInt}, + V::Vector{PetscScalar}, + mode = ADD_VALUES, +) + for (i, j, v) in zip(I, J, V) + setValue(mat, i - PetscIntOne, j - PetscIntOne, v, mode) end end -set_values!(mat, I, J, V, mode = ADD_VALUES) = set_values!(mat, PetscInt.(I), PetscInt.(J), PetscScalar.(V), mode) +function set_values!(mat, I, J, V, mode = ADD_VALUES) + set_values!(mat, PetscInt.(I), PetscInt.(J), PetscScalar.(V), mode) +end # Warning : cannot use Vector{Integer} because `[1, 2] isa Vector{Integer}` is `false` -_preallocate!(mat::PetscMat, dnz::Integer, onz::Integer, ::Val{:mpiaij}) = MatMPIAIJSetPreallocation(mat, PetscInt(dnz), PetscInt(onz)) -_preallocate!(mat::PetscMat, d_nnz::Vector{I}, o_nnz::Vector{I}, ::Val{:mpiaij}) where I = MatMPIAIJSetPreallocation(mat, PetscInt(0), PetscInt.(d_nnz), PetscInt(0), PetscInt.(o_nnz)) -_preallocate!(mat::PetscMat, nz::Integer, ::Integer, ::Val{:seqaij}) = MatSeqAIJSetPreallocation(mat, PetscInt(nz)) -_preallocate!(mat::PetscMat, nnz::Vector{I}, ::Vector{I}, ::Val{:seqaij}) where I = MatSeqAIJSetPreallocation(mat, PetscInt(0), PetscInt.(nnz)) +function _preallocate!(mat::Mat, dnz::Integer, onz::Integer, ::Val{:mpiaij}) + MPIAIJSetPreallocation(mat, PetscInt(dnz), PetscInt(onz)) +end +function _preallocate!( + mat::Mat, + d_nnz::Vector{I}, + o_nnz::Vector{I}, + ::Val{:mpiaij}, +) where {I} + MPIAIJSetPreallocation( + mat, + PetscInt(0), + PetscInt.(d_nnz), + PetscInt(0), + PetscInt.(o_nnz), + ) +end +function _preallocate!(mat::Mat, nz::Integer, ::Integer, ::Val{:seqaij}) + SeqAIJSetPreallocation(mat, PetscInt(nz)) +end +function _preallocate!(mat::Mat, nnz::Vector{I}, ::Vector{I}, ::Val{:seqaij}) where {I} + SeqAIJSetPreallocation(mat, PetscInt(0), PetscInt.(nnz)) +end """ - preallocate!(mat::PetscMat, dnz, onz, warn::Bool = true) + preallocate!(mat::Mat, dnz, onz, warn::Bool = true) Dispatch preallocation according matrix type (seq or mpiaij for instance). TODO: should use kwargs. """ -function preallocate!(mat::PetscMat, dnz, onz, warn::Bool = true) - _preallocate!(mat, dnz, onz, Val(Symbol(MatGetType(mat)))) - MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR, warn) +function preallocate!(mat::Mat, dnz, onz, warn::Bool = true) + _preallocate!(mat, dnz, onz, Val(Symbol(getType(mat)))) + setOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR, warn) end - """ - mat2file(mat::PetscMat, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") + mat2file(mat::Mat, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") -Write a PetscMat to a file. +Write a Mat to a file. """ -function mat2file(mat::PetscMat, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") +function mat2file( + mat::Mat, + filename::String, + format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, + type::String = "ascii", +) viewer = PetscViewer(mat.comm, filename, format, type) - MatView(mat, viewer) + matView(mat, viewer) destroy!(viewer) end -Base.show(::IO, mat::PetscMat) = MatView(mat) - -destroy!(mat::PetscMat) = MatDestroy(mat) \ No newline at end of file +Base.show(::IO, mat::Mat) = matView(mat) diff --git a/src/fancy/vec.jl b/src/fancy/vec.jl index f059e6b..54f2989 100644 --- a/src/fancy/vec.jl +++ b/src/fancy/vec.jl @@ -1,14 +1,36 @@ -""" - create_vector(nrows, nrows_loc = PETSC_DECIDE) +function Base.:*(x::Vec, alpha::Number) + y = VecCopy(x) + scale!(y, alpha) + return y +end +Base.:*(alpha::Number, vec::Vec) = vec * alpha -Create a `PetscVec` vector of global size `(nrows)`. -""" -function create_vector(nrows, nrows_loc = PETSC_DECIDE; auto_setup = false, comm::MPI.Comm = MPI.COMM_WORLD) - vec = VecCreate(comm) - VecSetSizes(vec::PetscVec, nrows_loc, nrows) +function assemble!(vec::Vec) + assemblyBegin(vec) + assemblyEnd(vec) +end - if (auto_setup) +""" + create_vector( + comm::MPI.Comm = MPI.COMM_WORLD; + nrows_loc = PETSC_DECIDE, + nrows_glo = PETSC_DECIDE, + autosetup = false, + ) + +Create a `Vec` vector of global size `(nrows_glo)`. +""" +function create_vector( + comm::MPI.Comm = MPI.COMM_WORLD; + nrows_loc = PETSC_DECIDE, + nrows_glo = PETSC_DECIDE, + autosetup = false, +) + vec = create(Vec, comm) + setSizes(vec, nrows_loc, nrows_glo) + + if autosetup set_from_options!(vec) set_up!(vec) end @@ -16,102 +38,148 @@ function create_vector(nrows, nrows_loc = PETSC_DECIDE; auto_setup = false, comm return vec end -Base.ndims(::Type{PetscVec}) = 1 +duplicate(vec::Vec, n::Int) = ntuple(i -> duplicate(vec), n) """ - Base.setindex!(vec::PetscVec, value::Number, row::Integer) + get_range(vec::Vec) -`row` must be in [1,size(vec)], i.e indexing starts at 1 (Julia). +Wrapper to `VecGetOwnershipRange` -# Implementation -For some unkwnown reason, calling `VecSetValue` fails. +However, the result `(rstart, rend)` is such that `vec[rstart:rend]` are the rows handled by the local processor. +This is different from the default `PETSc.VecGetOwnershipRange` result where the indexing starts at zero and where +`rend-1` is last row handled by the local processor. """ -function Base.setindex!(vec::PetscVec, value::Number, row::Integer) - VecSetValues(vec, PetscInt[row .- 1], PetscScalar[value], INSERT_VALUES) +function get_range(vec::Vec) + rstart, rend = getOwnershipRange(vec) + return (rstart + 1, rend) end +""" + get_urange(vec::Vec) -# This is stupid but I don't know how to do better yet -Base.setindex!(vec::PetscVec, values, rows) = VecSetValues(vec, collect(rows .- 1), values, INSERT_VALUES) +Provide a `UnitRange` from the method `get_range`. +""" +function get_urange(vec::Vec) + rstart, rend = getOwnershipRange(vec) + return (rstart + 1):rend +end +Base.ndims(::Type{Vec}) = 1 -set_global_size!(vec::PetscVec, nrows) = VecSetSizes(vec, PETSC_DECIDE, nrows) -set_local_size!(vec::PetscVec, nrows) = VecSetSizes(vec, nrows, PETSC_DECIDE) +const scale! = scale -set_from_options!(vec::PetscVec) = VecSetFromOptions(vec) +""" + Base.setindex!(vec::Vec, value::Number, row::Integer) + +`row` must be in [1,size(vec)], i.e indexing starts at 1 (Julia). -set_up!(vec::PetscVec) = VecSetUp(vec) +# Implementation +For some unkwnown reason, calling `VecSetValue` fails. """ - get_range(vec::PetscVec) +function Base.setindex!(vec::Vec, value::Number, row::Integer) + setValues(vec, PetscInt[row .- 1], PetscScalar[value], INSERT_VALUES) +end -Wrapper to `VecGetOwnershipRange` +# This is stupid but I don't know how to do better yet +function Base.setindex!(vec::Vec, values, rows) + setValues(vec, collect(rows .- 1), values, INSERT_VALUES) +end + +set_global_size!(vec::Vec, nrows) = setSizes(vec, PETSC_DECIDE, nrows) -However, the result `(rstart, rend)` is such that `vec[rstart:rend]` are the rows handled by the local processor. -This is different from the default `PETSc.VecGetOwnershipRange` result where the indexing starts at zero and where -`rend-1` is last row handled by the local processor. """ -function get_range(vec::PetscVec) - rstart, rend = VecGetOwnershipRange(vec) - return (rstart + 1, rend) +Wrapper to `ISSetLocalToGlobalMapping` + +1-based indexing +""" +function set_local_to_global!(vec::Vec, lid2gid::Vector{I}) where {I<:Integer} + mapping = create(ISLocalToGlobalMapping, vec.comm, lid2gid .- 1) + setLocalToGlobalMapping(vec, mapping) + destroy(mapping) end +set_local_size!(vec::Vec, nrows) = setSizes(vec, nrows, PETSC_DECIDE) + """ - get_urange(vec::PetscVec) + set_value!(vec::Vec, row, value, mode::InsertMode = INSERT_VALUES) -Provide a `UnitRange` from the method `get_range`. +Wrapper to `VecSetValue`, using julia 1-based indexing. """ -function get_urange(vec::PetscVec) - rstart, rend = VecGetOwnershipRange(vec) - return rstart+1:rend +function set_value!(vec::Vec, row, value, mode::InsertMode = INSERT_VALUES) + setValue(vec, row - 1, value, mode) end +""" + set_values!(vec::Vec, rows, values, mode::InsertMode = INSERT_VALUES) + set_values!(vec::Vec, values) -# Discutable choice -Base.length(vec::PetscVec) = VecGetLocalSize(vec) -Base.size(vec::PetscVec) = (length(vec),) +Wrapper to `VecSetValues`, using julia 1-based indexing. +""" +function set_values!(vec::Vec, rows, values, mode::InsertMode = INSERT_VALUES) + setValues(vec, rows .- 1, values, mode) +end -function assemble!(vec::PetscVec) - VecAssemblyBegin(vec) - VecAssemblyEnd(vec) +function set_values!(vec::Vec, values) + setValues(vec, collect(get_urange(vec) .- 1), values, INSERT_VALUES) end -duplicate(vec::PetscVec) = VecDuplicate(vec) +""" + set_value_local!(vec::Vec, row, value, mode::InsertMode = INSERT_VALUES) -set_values!(vec::PetscVec, values) = VecSetValues(vec, collect(get_urange(vec) .- 1), values, INSERT_VALUES) +Wrapper to `VecSetValueLocal`, using julia 1-based indexing. +""" +function set_value_local!(vec::Vec, row, value, mode::InsertMode = INSERT_VALUES) + setValueLocal(vec, row - 1, value, mode) +end """ -Wrapper to `VecSetValues`, using julia 1-based indexing. + set_values_local!(vec::Vec, rows, values, mode::InsertMode = INSERT_VALUES) + set_values_local!(vec::Vec, values) + +Wrapper to `VecSetValuesLocal`, using julia 1-based indexing. """ -function set_values!(vec::PetscVec, rows::Vector{PetscInt}, values::Vector{PetscScalar}, mode::InsertMode = INSERT_VALUES) - VecSetValues(vec, rows .- PetscIntOne, values, mode) +function set_values_local!(vec::Vec, rows, values, mode::InsertMode = INSERT_VALUES) + setValuesLocal(vec, rows .- 1, values, mode) +end + +function set_values_local!(vec::Vec, values) + # Since we are using the "local" numbering, we set element "1" to "nloc" + setValuesLocal(vec, collect(1:getLocalSize(vec)) .- 1, values, INSERT_VALUES) end +Base.show(::IO, vec::Vec) = vecView(vec) + """ - vec2array(vec::PetscVec) + vec2array(vec::Vec) -Convert a `PetscVec` into a Julia `Array`. Allocation is involved in the process since the `PetscVec` +Convert a `Vec` into a Julia `Array`. Allocation is involved in the process since the `Vec` allocated by PETSC is copied into a freshly allocated array. If you prefer not to allocate memory, use `VectGetArray` and `VecRestoreArray` """ -function vec2array(vec::PetscVec) - arrayFromC, array_ref = VecGetArray(vec) +function vec2array(vec::Vec) + arrayFromC, array_ref = getArray(vec) array = copy(arrayFromC) - VecRestoreArray(vec, array_ref) + restoreArray(vec, array_ref) return array end """ - vec2file(vec::PetscVec, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") + vec2file(vec::Vec, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") -Write a PetscVec to a file. +Write a Vec to a file. """ -function vec2file(vec::PetscVec, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii") +function vec2file( + vec::Vec, + filename::String, + format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, + type::String = "ascii", +) viewer = PetscViewer(vec.comm, filename, format, type) - VecView(vec, viewer) + vecView(vec, viewer) destroy!(viewer) end -Base.show(::IO, vec::PetscVec) = VecView(vec) - -destroy!(vec::PetscVec) = VecDestroy(vec) \ No newline at end of file +# Discutable choice(s) +Base.length(vec::Vec) = getLocalSize(vec) +Base.size(vec::Vec) = (length(vec),) diff --git a/src/fancy/viewer.jl b/src/fancy/viewer.jl index 3659e8c..676a520 100644 --- a/src/fancy/viewer.jl +++ b/src/fancy/viewer.jl @@ -1,20 +1,30 @@ -const push_format! = PetscViewerPushFormat -const set_mode! = PetscViewerFileSetMode -const set_name! = PetscViewerFileSetName -const set_type! = PetscViewerSetType +const push_format! = pushFormat +const set_mode! = fileSetMode +const set_name! = fileSetName +const set_type! = setType """ - PetscViewer(comm::MPI.Comm, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii", mode::PetscFileMode = FILE_MODE_WRITE) + PetscViewer( + comm::MPI.Comm, + filename::String, + format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, + type::String = "ascii", + mode::PetscFileMode = FILE_MODE_WRITE, + ) Constructor for a PetscViewer intended to read/write a matrix or a vector with the supplied type and format. """ -function PetscViewer(comm::MPI.Comm, filename::String, format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, type::String = "ascii", mode::PetscFileMode = FILE_MODE_WRITE) - viewer = PetscViewerCreate(comm) +function PetscViewer( + comm::MPI.Comm, + filename::String, + format::PetscViewerFormat = PETSC_VIEWER_ASCII_CSV, + type::String = "ascii", + mode::PetscFileMode = FILE_MODE_WRITE, +) + viewer = create(PetscViewer, comm) set_type!(viewer, type) - set_mode!(viewer, FILE_MODE_WRITE) + set_mode!(viewer, mode) push_format!(viewer, format) set_name!(viewer, filename) return viewer end - -destroy!(viewer::PetscViewer) = PetscViewerDestroy(viewer) \ No newline at end of file diff --git a/src/init.jl b/src/init.jl index 7764b09..4595f5c 100644 --- a/src/init.jl +++ b/src/init.jl @@ -2,34 +2,49 @@ Wrapper to `PetscInitializeNoPointers`. Initialize PETCs with arguments. # Implementation + I don't know if I am supposed to use PetscInt or not... """ -function PetscInitialize(args::Vector{String}, filename::String, help::String) +function PetscInitialize( + args::Vector{String}, + filename::String, + help::String; + finalize_atexit = true, +) + MPI.Initialized() || MPI.Init() + args2 = ["julia"; args] nargs = Cint(length(args2)) - error = ccall( (:PetscInitializeNoPointers, libpetsc), - PetscErrorCode, - (Cint, - Ptr{Ptr{UInt8}}, - Cstring, - Cstring), - nargs, args2, filename, help + error = ccall( + (:PetscInitializeNoPointers, libpetsc), + PetscErrorCode, + (Cint, Ptr{Ptr{UInt8}}, Cstring, Cstring), + nargs, + args2, + filename, + help, ) @assert iszero(error) + + finalize_atexit && atexit(PetscFinalize) end """ Initialize PETSc with arguments stored in a vector of string """ -PetscInitialize(args::Vector{String}) = PetscInitialize(args, "", "") +function PetscInitialize(args::Vector{String}; finalize_atexit = true) + PetscInitialize(args, "", ""; finalize_atexit) +end """ Initialize PETSc with arguments concatenated in a unique string. """ -PetscInitialize(args::String) = PetscInitialize(convert(Vector{String}, split(args)), "", "") +function PetscInitialize(args::String; finalize_atexit = true) + PetscInitialize(convert(Vector{String}, split(args)), "", ""; finalize_atexit) +end """ - PetscInitialize(cmd_line_args::Bool = true) + PetscInitialize(cmd_line_args::Bool = true; finalize_atexit = true) Initialize PETSc. @@ -39,12 +54,16 @@ arguments for PETSc (leading to a call to `PetscInitializeNoPointers`). Otherwise, if `cmd_line_args == false`, initialize PETSc without arguments (leading to a call to `PetscInitializeNoArguments`). """ -function PetscInitialize(cmd_line_args::Bool = true) +function PetscInitialize(cmd_line_args::Bool = true; finalize_atexit = true) if (cmd_line_args) - PetscInitialize(ARGS) + PetscInitialize(ARGS; finalize_atexit) else + MPI.Initialized() || MPI.Init() + error = ccall((:PetscInitializeNoArguments, libpetsc), PetscErrorCode, ()) @assert iszero(error) + + finalize_atexit && atexit(PetscFinalize) end end @@ -52,6 +71,54 @@ end Wrapper to PetscFinalize """ function PetscFinalize() - error = ccall( (:PetscFinalize, libpetsc), PetscErrorCode, ()) + PetscFinalized() && return + + GC.gc() + if _NREFS[] != 0 + @warn "$(_NREFS[]) objects still not finalized before calling PetscWrap.Finalize()" + end + + error = ccall((:PetscFinalize, libpetsc), PetscErrorCode, ()) + @assert iszero(error) + + if _NREFS[] != 0 + @warn "$(_NREFS[]) objects still not finalized after calling PetscWrap.Finalize()" + end + + _NREFS[] = 0 +end + +""" + Wrapper to PetscInitialized +""" +function PetscInitialized() + isInitialized = Ref{PetscBool}() + error = ccall( + (:PetscInitialized, libpetsc), + PetscErrorCode, + (Ref{PetscBool},), + isInitialized, + ) @assert iszero(error) -end \ No newline at end of file + return Bool(isInitialized[]) +end + +""" + Wrapper to PetscFinalized + +Cmd line options: +-options_view - Calls PetscOptionsView() +-options_left - Prints unused options that remain in the database +-objects_dump [all] - Prints list of objects allocated by the user that have not been freed, the option all cause all outstanding objects to be listed +-mpidump - Calls PetscMPIDump() +-malloc_dump - Calls PetscMallocDump(), displays all memory allocated that has not been freed +-malloc_info - Prints total memory usage +-malloc_view - Prints list of all memory allocated and where +""" +function PetscFinalized() + isFinalized = Ref{PetscBool}() + error = + ccall((:PetscFinalized, libpetsc), PetscErrorCode, (Ref{PetscBool},), isFinalized) + @assert iszero(error) + return Bool(isFinalized[]) +end diff --git a/src/ksp.jl b/src/ksp.jl deleted file mode 100644 index 05c5013..0000000 --- a/src/ksp.jl +++ /dev/null @@ -1,73 +0,0 @@ -const CKSP = Ptr{Cvoid} - -struct PetscKSP - ptr::Ref{CKSP} - comm::MPI.Comm - - PetscKSP(comm::MPI.Comm) = new(Ref{CKSP}(), comm) -end - -# allows us to pass PetscKSP objects directly into CKSP ccall signatures -Base.cconvert(::Type{CKSP}, ksp::PetscKSP) = ksp.ptr[] - -""" - KSPCreate(comm::MPI.Comm = MPI.COMM_WORLD) - -Wrapper for `KSPCreate` -""" -function KSPCreate(comm::MPI.Comm = MPI.COMM_WORLD) - ksp = PetscKSP(comm) - error = ccall((:KSPCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CKSP}), comm, ksp.ptr) - @assert iszero(error) - return ksp -end - -""" - KSPDestroy(ksp::PetscKSP) - -Wrapper to `KSPDestroy` -""" -function KSPDestroy(ksp::PetscKSP) - error = ccall((:KSPDestroy, libpetsc), PetscErrorCode, (Ptr{CKSP},), ksp.ptr) - @assert iszero(error) -end - -""" - KSPSetFromOptions(ksp::PetscKSP) - -Wrapper to KSPSetFromOptions -""" -function KSPSetFromOptions(ksp::PetscKSP) - error = ccall((:KSPSetFromOptions, libpetsc), PetscErrorCode, (CKSP,), ksp) - @assert iszero(error) -end - -""" - KSPSetOperators(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat) - -Wrapper for KSPSetOperators -""" -function KSPSetOperators(ksp::PetscKSP, Amat::PetscMat, Pmat::PetscMat) - error = ccall((:KSPSetOperators, libpetsc), PetscErrorCode, (CKSP, CMat, CMat), ksp, Amat, Pmat) - @assert iszero(error) -end - -""" - KSPSetUp(ksp::PetscKSP) - -Wrapper to KSPSetUp -""" -function KSPSetUp(ksp::PetscKSP) - error = ccall((:KSPSetUp, libpetsc), PetscErrorCode, (CKSP,), ksp) - @assert iszero(error) -end - -""" - KSPSolve(ksp::PetscKSP, b::PetscVec, x::PetscVec) - -Wrapper for KSPSolve -""" -function KSPSolve(ksp::PetscKSP, b::PetscVec, x::PetscVec) - error = ccall((:KSPSolve, libpetsc), PetscErrorCode, (CKSP, CVec, CVec), ksp, b, x) - @assert iszero(error) -end \ No newline at end of file diff --git a/src/load.jl b/src/load.jl index be99133..9b45f17 100644 --- a/src/load.jl +++ b/src/load.jl @@ -2,32 +2,36 @@ Function from GridapPETSC to find PETSc lib location. """ function get_petsc_location() - PETSC_DIR = haskey(ENV,"PETSC_DIR") ? ENV["PETSC_DIR"] : "/usr/lib/petsc" - PETSC_ARCH = haskey(ENV,"PETSC_ARCH") ? ENV["PETSC_ARCH"] : "" + PETSC_DIR = haskey(ENV, "PETSC_DIR") ? ENV["PETSC_DIR"] : "/usr/lib/petsc" + PETSC_ARCH = haskey(ENV, "PETSC_ARCH") ? ENV["PETSC_ARCH"] : "" PETSC_LIB = "" # Check PETSC_DIR exists if isdir(PETSC_DIR) # Define default paths - PETSC_LIB_DIR = joinpath(PETSC_DIR,PETSC_ARCH,"lib") + PETSC_LIB_DIR = joinpath(PETSC_DIR, PETSC_ARCH, "lib") # Check PETSC_LIB (.../libpetsc.so or .../libpetsc_real.so file) exists - if isfile(joinpath(PETSC_LIB_DIR,"libpetsc.so")) - PETSC_LIB = joinpath(PETSC_LIB_DIR,"libpetsc.so") - elseif isfile(joinpath(PETSC_LIB_DIR,"libpetsc_real.so")) - PETSC_LIB = joinpath(PETSC_LIB_DIR,"libpetsc_real.so") + if isfile(joinpath(PETSC_LIB_DIR, "libpetsc.so")) + PETSC_LIB = joinpath(PETSC_LIB_DIR, "libpetsc.so") + elseif isfile(joinpath(PETSC_LIB_DIR, "libpetsc_real.so")) + PETSC_LIB = joinpath(PETSC_LIB_DIR, "libpetsc_real.so") end end # PETSc lib not found - if(length(PETSC_LIB) == 0) + if (length(PETSC_LIB) == 0) # Workaround for automerging on RegistryCI or doc deployment - if(haskey(ENV,"JULIA_REGISTRYCI_AUTOMERGE") || haskey(ENV, "DOC_DEPLOYMENT")) + if (haskey(ENV, "JULIA_REGISTRYCI_AUTOMERGE") || haskey(ENV, "DOC_DEPLOYMENT")) @warn "Setting fictive PETSc path because of RegistryCI or Doc deployment" PETSC_LIB = "JULIA_REGISTRYCI_AUTOMERGE" else - throw(ErrorException("PETSc shared library (libpetsc.so) not found. Please check that PETSC_DIR and PETSC_ARCH env. variables are set.")) + throw( + ErrorException( + "PETSc shared library (libpetsc.so) not found. Please check that PETSC_DIR and PETSC_ARCH env. variables are set.", + ), + ) end end @@ -37,4 +41,4 @@ end # Absolute path to libpetsc.so const libpetsc = get_petsc_location() -show_petsc_path() = return libpetsc \ No newline at end of file +show_petsc_path() = return libpetsc diff --git a/src/mat.jl b/src/mat.jl deleted file mode 100644 index ddf5526..0000000 --- a/src/mat.jl +++ /dev/null @@ -1,362 +0,0 @@ -const CMat = Ptr{Cvoid} - -""" - A Petsc matrix, actually just a pointer to the actual C matrix -""" -struct PetscMat - ptr::Ref{CMat} - comm::MPI.Comm - - PetscMat(comm::MPI.Comm) = new(Ref{CMat}(), comm) -end - - -# allows us to pass PetscMat objects directly into CMat ccall signatures -Base.cconvert(::Type{CMat}, mat::PetscMat) = mat.ptr[] - -""" - Wrapper to `MatAssemblyBegin` -""" -function MatAssemblyBegin(mat::PetscMat, type::MatAssemblyType) - error = ccall((:MatAssemblyBegin, libpetsc), PetscErrorCode, (CMat, MatAssemblyType), mat, type) - @assert iszero(error) -end - -""" - Wrapper to `MatAssemblyEnd` -""" -function MatAssemblyEnd(mat::PetscMat, type::MatAssemblyType) - error = ccall((:MatAssemblyEnd, libpetsc), PetscErrorCode, (CMat, MatAssemblyType), mat, type) - @assert iszero(error) -end - - -""" - MatCompositeAddMat(mat::PetscMat, smat::PetscMat) - -Wrapper to `MatCompositeAddMat`. -""" -function MatCompositeAddMat(mat::PetscMat, smat::PetscMat) - error = ccall((:MatCompositeAddMat, libpetsc), PetscErrorCode, (CMat, CMat), mat, smat) - @assert iszero(error) -end - -""" - MatCreate(comm::MPI.Comm, mat::PetscMat) - -Wrapper to `MatCreate` -""" -function MatCreate(comm::MPI.Comm = MPI.COMM_WORLD) - mat = PetscMat(comm) - error = ccall((:MatCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CMat}), comm, mat.ptr) - @assert iszero(error) - return mat -end - -""" - MatCreateDense(comm::MPI.Comm, m::PetscInt, n::PetscInt, M::PetscInt, N::PetscInt) - -Wrapper to `MatCreateDense`. Last argument `data` is not supported yet (NULL is passed). -""" -function MatCreateDense(comm::MPI.Comm, m::PetscInt, n::PetscInt, M::PetscInt, N::PetscInt) - mat = PetscMat(comm) - error = ccall((:MatCreateDense, libpetsc), PetscErrorCode, - (MPI.MPI_Comm, PetscInt, PetscInt, PetscInt, PetscInt, Ptr{PetscScalar}, Ptr{CMat}), - comm, m, n, M, N, C_NULL, mat.ptr) - @assert iszero(error) - return mat -end - -function MatCreateDense(comm::MPI.Comm, - m::Integer = PETSC_DECIDE, - n::Integer = PETSC_DECIDE, - M::Integer = PETSC_DECIDE, - N::Integer = PETSC_DECIDE) - - return MatCreateDense(comm, PetscInt(m), PetscInt(n), PetscInt(M), PetscInt(N)) -end - -""" - Wrapper to MatCreateVecs -""" -function MatCreateVecs(mat::PetscMat, vecr::PetscVec, veci::PetscVec) - error = ccall((:MatCreateVecs, libpetsc), PetscErrorCode, (CMat, Ptr{CVec}, Ptr{CVec}), mat, vecr.ptr, veci.ptr) - @assert iszero(error) -end - -function MatCreateVecs(mat::PetscMat) - vecr = PetscVec(mat.comm); veci = PetscVec(mat.comm) - MatCreateVecs(mat, vecr, veci) - return vecr, veci -end - -""" - MatDestroy(mat::PetscMat) - -Wrapper to MatDestroy -""" -function MatDestroy(mat::PetscMat) - error = ccall((:MatDestroy, libpetsc), PetscErrorCode, (Ptr{CMat},), mat.ptr) - @assert iszero(error) -end - -""" - MatGetLocalSize(mat::PetscMat) - -Wrapper to `MatGetLocalSize`. Return the number of local rows and cols of the matrix (i.e on the processor). -""" -function MatGetLocalSize(mat::PetscMat) - nrows_loc = Ref{PetscInt}(0) - ncols_loc = Ref{PetscInt}(0) - - error = ccall((:MatGetLocalSize, libpetsc), PetscErrorCode, (CMat, Ref{PetscInt}, Ref{PetscInt}), mat, nrows_loc, ncols_loc) - @assert iszero(error) - - return nrows_loc[], ncols_loc[] -end - -""" - MatGetOwnershipRange(mat::PetscMat) - -Wrapper to `MatGetOwnershipRange` - -The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor. - -# Warning -`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row -handled by the local processor. -""" -function MatGetOwnershipRange(mat::PetscMat) - rstart = Ref{PetscInt}(0) - rend = Ref{PetscInt}(0) - - error = ccall((:MatGetOwnershipRange, libpetsc), PetscErrorCode, (CMat, Ref{PetscInt}, Ref{PetscInt}), mat, rstart, rend) - @assert iszero(error) - - return rstart[], rend[] -end - -""" - MatGetOwnershipRangeColumn(mat::PetscMat) - -Wrapper to `MatGetOwnershipRangeColumn` - -The result `(cstart, cend)` is a Tuple indicating the columns handled by the local processor. - -# Warning -`PETSc` indexing starts at zero (so `cstart` may be zero) and `cend-1` is the last column -handled by the local processor. -""" -function MatGetOwnershipRangeColumn(mat::PetscMat) - cstart = Ref{PetscInt}(0) - cend = Ref{PetscInt}(0) - - error = ccall((:MatGetOwnershipRangeColumn, libpetsc), PetscErrorCode, (CMat, Ref{PetscInt}, Ref{PetscInt}), mat, cstart, cend) - @assert iszero(error) - - return cstart[], cend[] -end - -""" - MatGetSize(mat::PetscMat) - -Wrapper to `MatGetSize`. Return the number of rows and cols of the matrix (global number). -""" -function MatGetSize(mat::PetscMat) - nrows_glo = Ref{PetscInt}(0) - ncols_glo = Ref{PetscInt}(0) - - error = ccall((:MatGetSize, libpetsc), PetscErrorCode, (CMat, Ref{PetscInt}, Ref{PetscInt}), mat, nrows_glo, ncols_glo) - @assert iszero(error) - - return nrows_glo[], ncols_glo[] -end - -""" - MatGetType(mat::PetscMat) - -Wrapper to `MatGetType`. Return the matrix type as a string. See matrix types here: -https://petsc.org/release/docs/manualpages/Mat/MatType.html#MatType -""" -function MatGetType(mat::PetscMat) - type = Ref{Cstring}() - - error = ccall((:MatGetType, libpetsc), PetscErrorCode, (CMat, Ptr{Cstring}), mat, type) - @assert iszero(error) - - return unsafe_string(type[]) -end - - -""" - MatMult(mat::PetscMat, x::PetscVec, y::PetscVec) - -Wrapper to `MatMult` -""" -function MatMult(mat::PetscMat, x::PetscVec, y::PetscVec) - error = ccall((:MatMult, libpetsc), PetscErrorCode, (CMat, CVec, CVec), mat, x, y) - @assert iszero(error) -end - -""" - MatSetFromOptions(mat::PetscMat) - -Wrapper to MatSetFromOptions -""" -function MatSetFromOptions(mat::PetscMat) - error = ccall((:MatSetFromOptions, libpetsc), PetscErrorCode, (CMat,), mat) - @assert iszero(error) -end - -""" - MatSetSizes(mat::PetscMat, nrows_loc::PetscInt, ncols_loc::PetscInt, nrows_glo::PetscInt, ncols_glo::PetscInt) - -Wrapper to MatSetSizes -""" -function MatSetSizes(mat::PetscMat, nrows_loc::PetscInt, ncols_loc::PetscInt, nrows_glo::PetscInt, ncols_glo::PetscInt) - error = ccall((:MatSetSizes, libpetsc), - PetscErrorCode, - (CMat, PetscInt, PetscInt, PetscInt, PetscInt), - mat, nrows_loc, ncols_loc, nrows_glo, ncols_glo - ) - @assert iszero(error) -end - -function MatSetSizes(mat::PetscMat, nrows_loc::Integer, ncols_loc::Integer, nrows_glo::Integer, ncols_glo::Integer) - MatSetSizes(mat, PetscInt(nrows_loc), PetscInt(ncols_loc), PetscInt(nrows_glo), PetscInt(ncols_glo)) -end - -""" - MatSetType(mat::PetscMat, type::String) - -Wrapper for `MatSetType`. Values for `type` alors available here: -https://petsc.org/release/docs/manualpages/Mat/MatType.html#MatType -""" -function MatSetType(mat::PetscMat, type::String) - error = ccall((:MatSetType, libpetsc), PetscErrorCode, (CMat, Cstring), mat, type) - @assert iszero(error) -end - - -""" - MatSetUp(mat::PetscMat) - -Wrapper to MatSetUp -""" -function MatSetUp(mat::PetscMat) - error = ccall((:MatSetUp, libpetsc), PetscErrorCode, (CMat,), mat) - @assert iszero(error) -end - -""" - MatMPIAIJSetPreallocation(mat::PetscMat, dnz::PetscInt, dnnz::Vector{PetscInt}, onz::PetscInt, onnz::Vector{PetscInt}) - -Wrapper to `MatMPIAIJSetPreallocation`. `dnnz` and `onnz` not tested yet. - -# Arguments -- `dnz::PetscInt`: number of nonzeros per row in DIAGONAL portion of local submatrix (same value is used for all local rows) -- `dnnz::Vector{PetscInt}`: array containing the number of nonzeros in the various rows of the DIAGONAL portion of the local - submatrix (possibly different for each row) or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero - structure. The size of this array is equal to the number of local rows, i.e 'm'. For matrices that will be factored, you - must leave room for (and set) the diagonal entry even if it is zero. -- `onz::PetscInt`: number of nonzeros per row in the OFF-DIAGONAL portion of local submatrix (same value is used for all local rows). -- `onnz::Vector{PetscInt}`: array containing the number of nonzeros in the various rows of the OFF-DIAGONAL portion of the local - submatrix (possibly different for each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero structure. - The size of this array is equal to the number of local rows, i.e 'm'. -""" -function MatMPIAIJSetPreallocation(mat::PetscMat, dnz::PetscInt, dnnz::Vector{PetscInt}, onz::PetscInt, onnz::Vector{PetscInt}) - error = ccall((:MatMPIAIJSetPreallocation, libpetsc), PetscErrorCode, - (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}), - mat, dnz, dnnz, onz, onnz) - @assert iszero(error) -end - -function MatMPIAIJSetPreallocation(mat::PetscMat, dnz::PetscInt, onz::PetscInt) - error = ccall((:MatMPIAIJSetPreallocation, libpetsc), PetscErrorCode, - (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}), - mat, dnz, C_NULL, onz, C_NULL) - @assert iszero(error) -end - -""" - MatSeqAIJSetPreallocation(mat::PetscMat, nz::PetscInt, nnz::Vector{PetscInt}) - -Wrapper to `MatSeqAIJSetPreallocation`. -""" -function MatSeqAIJSetPreallocation(mat::PetscMat, nz::PetscInt, nnz::Vector{PetscInt}) - error = ccall((:MatSeqAIJSetPreallocation, libpetsc), PetscErrorCode, - (CMat, PetscInt, Ptr{PetscInt}), - mat, nz, nnz) - @assert iszero(error) -end - -function MatSeqAIJSetPreallocation(mat::PetscMat, nz::PetscInt) - error = ccall((:MatSeqAIJSetPreallocation, libpetsc), PetscErrorCode, - (CMat, PetscInt, Ptr{PetscInt}), - mat, nz, C_NULL) - @assert iszero(error) -end - -""" -Wrapper for `MatSetOption` -""" -function MatSetOption(mat::PetscMat, option::MatOption, value::Bool) - error = ccall((:MatSetOption, libpetsc), PetscErrorCode, (CMat, MatOption, Cuchar), mat, option, value) - @assert iszero(error) -end - -# Avoid allocating an array of size 1 for each call to MatSetValue -const _ivec = zeros(PetscInt,1) -const _jvec = zeros(PetscInt,1) -const _vvec = zeros(PetscScalar,1) - -""" -MatSetValue(mat::PetscMat, i::PetscInt, j::PetscInt, v::PetscScalar, mode::InsertMode) - -Wrapper to `MatSetValue`. Indexing starts at 0 (as in PETSc). - -# Implementation -For an unknow reason, calling PETSc.MatSetValue leads to an "undefined symbol: MatSetValue" error. -So this wrapper directly call MatSetValues (anyway, this is what is done in PETSc...) -""" -function MatSetValue(mat::PetscMat, i::PetscInt, j::PetscInt, v::PetscScalar, mode::InsertMode) - # Convert to arrays - _ivec[1] = i; _jvec[1] = j; _vvec[1] = v - - MatSetValues(mat, PetscIntOne, _ivec, PetscIntOne, _jvec, _vvec, mode) - #MatSetValues(mat, PetscIntOne, [i], PetscIntOne, [j], [v], mode) -end - -function MatSetValue(mat::PetscMat, i::Integer, j::Integer, v::Number, mode::InsertMode) - MatSetValue(mat, PetscInt(i), PetscInt(j), PetscScalar(v), mode) -end - -""" -MatSetValues(mat::PetscMat, nI::PetscInt, I::Vector{PetscInt}, nJ::PetscInt, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode) - -Wrapper to `MatSetValues`. Indexing starts at 0 (as in PETSc) -""" -function MatSetValues(mat::PetscMat, nI::PetscInt, I::Vector{PetscInt}, nJ::PetscInt, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode) - error = ccall((:MatSetValues, libpetsc), PetscErrorCode, - (CMat, PetscInt, Ptr{PetscInt}, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), - mat, nI, I, nJ, J, V, mode) - @assert iszero(error) -end - -function MatSetValues(mat::PetscMat, I::Vector{PetscInt}, J::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode) - MatSetValues(mat, PetscInt(length(I)), I, PetscInt(length(J)), J, V, mode) -end - -function MatSetValues(mat::PetscMat, I, J, V, mode::InsertMode) - MatSetValues(mat, PetscInt.(I), PetscInt.(collect(J)), PetscScalar.(V), mode) -end - -""" - MatView(mat::PetscMat, viewer::PetscViewer = PetscViewerStdWorld()) - -Wrapper to `MatView` -""" -function MatView(mat::PetscMat, viewer::PetscViewer = PetscViewerStdWorld()) - error = ccall((:MatView, libpetsc), PetscErrorCode, (CMat, CViewer), mat, viewer) - @assert iszero(error) -end \ No newline at end of file diff --git a/src/vec.jl b/src/vec.jl deleted file mode 100644 index fbac20f..0000000 --- a/src/vec.jl +++ /dev/null @@ -1,240 +0,0 @@ -const CVec = Ptr{Cvoid} - -struct PetscVec - ptr::Ref{CVec} - comm::MPI.Comm - - PetscVec(comm::MPI.Comm) = new(Ref{Ptr{Cvoid}}(), comm) -end - -# allows us to pass PetscVec objects directly into CVec ccall signatures -Base.cconvert(::Type{CVec}, vec::PetscVec) = vec.ptr[] - -""" - VecCreate(comm::MPI.Comm, vec::PetscVec) - -Wrapper to VecCreate -""" -function VecCreate(comm::MPI.Comm = MPI.COMM_WORLD) - vec = PetscVec(comm) - error = ccall((:VecCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CVec}), comm, vec.ptr) - @assert iszero(error) - return vec -end - -""" - VecSetValue(vec::PetscVec, i::PetscInt, v::PetscScalar, mode::InsertMode = INSERT_VALUES) - -Wrapper to `VecSetValue`. Indexing starts at 0 (as in PETSc). - -# Implementation -For an unknow reason, calling PETSc.VecSetValue leads to an "undefined symbol: VecSetValue" error. -So this wrapper directly call VecSetValues (anyway, this is what is done in PETSc...) -""" -function VecSetValue(vec::PetscVec, i::PetscInt, v::PetscScalar, mode::InsertMode = INSERT_VALUES) - VecSetValues(vec, PetscIntOne, [i], [v], mode) -end - -function VecSetValue(vec::PetscVec, i, v, mode::InsertMode = INSERT_VALUES) - VecSetValue(vec, PetscInt(i), PetscScalar(v), mode) -end - -""" - VecSetValues(vec::PetscVec, nI::PetscInt, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES) - -Wrapper to `VecSetValues`. Indexing starts at 0 (as in PETSc) -""" -function VecSetValues(vec::PetscVec, nI::PetscInt, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES) - error = ccall((:VecSetValues, libpetsc), PetscErrorCode, - (CVec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode), - vec, nI, I, V, mode) - @assert iszero(error) -end - -function VecSetValues(vec::PetscVec, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES) - VecSetValues(vec, PetscInt(length(I)), I, V, mode) -end - -function VecSetValues(vec::PetscVec, I, V, mode::InsertMode = INSERT_VALUES) - VecSetValues(vec, PetscInt.(I), PetscScalar.(V), mode) -end - - -""" - VecSetSizes(vec::PetscVec, nrows_loc, nrows_glo) - -Wrapper to VecSetSizes -""" -function VecSetSizes(vec::PetscVec, nrows_loc, nrows_glo) - nr_loc = PetscInt(nrows_loc) - nr_glo = PetscInt(nrows_glo) - error = ccall((:VecSetSizes, libpetsc), - PetscErrorCode, - (CVec, PetscInt, PetscInt), - vec, nr_loc, nr_glo - ) - @assert iszero(error) -end - -""" - VecSetFromOptions(vec::PetscVec) - -Wrapper to VecSetFromOptions -""" -function VecSetFromOptions(vec::PetscVec) - error = ccall((:VecSetFromOptions, libpetsc), PetscErrorCode, (CVec,), vec) - @assert iszero(error) -end - -""" - VecSetUp(vec::PetscVec) - -Wrapper to VecSetUp -""" -function VecSetUp(vec::PetscVec) - error = ccall((:VecSetUp, libpetsc), PetscErrorCode, (CVec,), vec) - @assert iszero(error) -end - -""" - VecGetOwnershipRange(vec::PetscVec) - -Wrapper to `VecGetOwnershipRange` - -The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor. - -# Warning -`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row -handled by the local processor. -""" -function VecGetOwnershipRange(vec::PetscVec) - rstart = Ref{PetscInt}(0) - rend = Ref{PetscInt}(0) - - error = ccall((:VecGetOwnershipRange, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}, Ref{PetscInt}), vec, rstart, rend) - @assert iszero(error) - - return rstart[], rend[] -end - -""" - VecGetSize(vec::PetscVec) - -Wrapper for VecGetSize -""" -function VecGetSize(vec::PetscVec) - n = Ref{PetscInt}() - - error = ccall((:VecGetSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), vec, n) - @assert iszero(error) - - return n[] -end - -""" - VecGetLocalSize(vec::PetscVec) - -Wrapper for VecGetLocalSize -""" -function VecGetLocalSize(vec::PetscVec) - n = Ref{PetscInt}() - - error = ccall((:VecGetLocalSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), vec, n) - @assert iszero(error) - - return n[] -end - -""" - VecAssemblyBegin(vec::PetscVec) - -Wrapper to VecAssemblyBegin -""" -function VecAssemblyBegin(vec::PetscVec) - error = ccall((:VecAssemblyBegin, libpetsc), PetscErrorCode, (CVec,), vec) - @assert iszero(error) -end - -""" - VecAssemblyEnd(vec::PetscVec) - -Wrapper to VecAssemblyEnd -""" -function VecAssemblyEnd(vec::PetscVec) - error = ccall((:VecAssemblyEnd, libpetsc), PetscErrorCode, (CVec,), vec) - @assert iszero(error) -end - -""" - VecDuplicate(vec::PetscVec) - -Wrapper for VecDuplicate, except that it returns the new vector instead of taking it as an input. -""" -function VecDuplicate(vec::PetscVec) - x = PetscVec(vec.comm) - error = ccall((:VecDuplicate, libpetsc), PetscErrorCode, (CVec, Ptr{CVec}), vec, x.ptr) - @assert iszero(error) - - return x -end - -""" - VecGetArray(vec::PetscVec, own = false) - -Wrapper for VecGetArray. - -# Warning -I am not confortable at all with memory management, both on the C side and on the Julia side. Use -this at you own risk. - -According to Julia documentation, `own` optionally specifies whether Julia should take ownership -of the memory, calling free on the pointer when the array is no longer referenced." - -""" -function VecGetArray(vec::PetscVec, own = false) - # Get array pointer - array_ref = Ref{Ptr{PetscScalar}}() - error = ccall((:VecGetArray, libpetsc), PetscErrorCode, (CVec, Ref{Ptr{PetscScalar}}), vec, array_ref) - @assert iszero(error) - - # Get array size - rstart, rend = VecGetOwnershipRange(vec) - n = rend - rstart # this is not `rend - rstart + 1` because of VecGetOwnershipRange convention - - array = unsafe_wrap(Array, array_ref[], n; own = own) - - return array, array_ref -end - -""" - VecRestoreArray(vec::PetscVec, array_ref) - -Wrapper for VecRestoreArray. `array_ref` is obtained from `VecGetArray` -""" -function VecRestoreArray(vec::PetscVec, array_ref) - error = ccall((:VecRestoreArray, libpetsc), PetscErrorCode, (CVec, Ref{Ptr{PetscScalar}}), vec, array_ref) - @assert iszero(error) -end - -""" - VecView(vec::PetscVec, viewer::PetscViewer = PetscViewerStdWorld()) - -Wrapper to `VecView`. -""" -function VecView(vec::PetscVec, viewer::PetscViewer = PetscViewerStdWorld()) - error = ccall( (:VecView, libpetsc), PetscErrorCode, (CVec, CViewer), vec, viewer); - @assert iszero(error) -end - -""" - VecDestroy(vec::PetscVec) - -Wrapper to VecDestroy -""" -function VecDestroy(vec::PetscVec) - error = ccall( (:VecDestroy, libpetsc), - PetscErrorCode, - (Ptr{CVec},), - vec.ptr) - @assert iszero(error) -end \ No newline at end of file diff --git a/src/viewer.jl b/src/viewer.jl deleted file mode 100644 index 02662fb..0000000 --- a/src/viewer.jl +++ /dev/null @@ -1,137 +0,0 @@ -const CViewer = Ptr{Cvoid} -struct PetscViewer - ptr::Ref{CViewer} - comm::MPI.Comm - - PetscViewer(comm::MPI.Comm) = new(Ref{CViewer}(), comm) -end - -# allows us to pass PetscViewer objects directly into CViewer ccall signatures -Base.cconvert(::Type{CViewer}, viewer::PetscViewer) = viewer.ptr[] - -""" - PetscViewerASCIIGetStdout(comm::MPI.Comm = MPI.COMM_WORLD) - -Wrapper for `PetscViewerASCIIGetStdout` -""" -function PetscViewerASCIIGetStdout(comm::MPI.Comm = MPI.COMM_WORLD) - viewer = PetscViewer(comm) - error = ccall((:PetscViewerASCIIGetStdout, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr) - @assert iszero(error) - return viewer -end - -const PetscViewerStdWorld = PetscViewerASCIIGetStdout - -""" - PetscViewerASCIIOpen(comm::MPI.Comm, filename) - -Wrapper for `PetscViewerASCIIOpen` -""" -function PetscViewerASCIIOpen(comm::MPI.Comm, filename) - viewer = PetscViewer(comm) - error = ccall((:PetscViewerASCIIOpen, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Cstring, Ptr{CViewer}), comm, filename, viewer.ptr) - @assert iszero(error) - return viewer -end - -""" - PetscViewerCreate(comm::MPI.Comm = MPI.COMM_WORLD) - -Wrapper for `PetscViewerCreate` -""" -function PetscViewerCreate(comm::MPI.Comm = MPI.COMM_WORLD) - viewer = PetscViewer(comm) - error = ccall((:PetscViewerCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CViewer}), comm, viewer.ptr) - @assert iszero(error) - return viewer -end - -""" - PetscViewerDestroy(viewer::PetscViewer) - -Wrapper for `PetscViewerDestroy` - -Warning : from what I understand, all viewers must not be destroyed explicitely using `PetscViewerDestroy`. - -""" -function PetscViewerDestroy(viewer::PetscViewer) - error = ccall((:PetscViewerDestroy, libpetsc), PetscErrorCode, (Ptr{CViewer}, ), viewer.ptr) - @assert iszero(error) -end - -""" - PetscViewerPopFormat(viewer::PetscViewer) - -Wrapper for `PetscViewerPopFormat` -""" -function PetscViewerPopFormat(viewer::PetscViewer) - error = ccall((:PetscViewerPopFormat, libpetsc), PetscErrorCode, (CViewer,), viewer) - @assert iszero(error) -end - -""" - PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat) - -Wrapper for `PetscViewerPushFormat` -""" -function PetscViewerPushFormat(viewer::PetscViewer, format::PetscViewerFormat) - error = ccall((:PetscViewerPushFormat, libpetsc), PetscErrorCode, (CViewer, PetscViewerFormat), viewer, format) - @assert iszero(error) -end - -""" - PetscViewerFileSetMode(viewer::PetscViewer, mode::PetscFileMode = FILE_MODE_WRITE) - -Wrapper for `PetscViewerFileSetMode` -""" -function PetscViewerFileSetMode(viewer::PetscViewer, mode::PetscFileMode = FILE_MODE_WRITE) - error = ccall((:PetscViewerFileSetMode, libpetsc), PetscErrorCode, (CViewer, PetscFileMode), viewer, mode) - @assert iszero(error) -end - - -""" - PetscViewerFileSetName(viewer::PetscViewer, filename) - -Wrapper for `PetscViewerFileSetName` -""" -function PetscViewerFileSetName(viewer::PetscViewer, filename::String) - error = ccall((:PetscViewerFileSetName, libpetsc), PetscErrorCode, (CViewer, Cstring), viewer, filename) - @assert iszero(error) -end - -""" - PetscViewerHDF5Open(comm::MPI.Comm, filename::String, type::PetscFileMode) - -Wrapper for `PetscViewerHDF5Open` -""" -function PetscViewerHDF5Open(comm::MPI.Comm, filename::String, type::PetscFileMode) - viewer = PetscViewer(comm) - error = ccall((:PetscViewerHDF5Open, libpetsc), PetscErrorCode, - (MPI.MPI_Comm, Cstring, PetscFileMode, Ptr{CViewer}), - comm, filename, type, viewer.ptr) - @assert iszero(error) - return viewer -end - -""" - PetscViewerSetType(viewer::PetscViewer, type::String) - -Wrapper for `PetscViewerSetType`. Values for `type` alors available here: -https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Viewer/PetscViewerType.html#PetscViewerType -""" -function PetscViewerSetType(viewer::PetscViewer, type::String) - error = ccall((:PetscViewerSetType, libpetsc), PetscErrorCode, (CViewer, Cstring), viewer, type) - @assert iszero(error) -end - -""" - PetscViewerView(v::PetscViewer, viewer::PetscViewer = PetscViewerStdWorld()) - -Wrapper to `PetscViewerView` -""" -function PetscViewerView(v::PetscViewer, viewer::PetscViewer = PetscViewerStdWorld()) - error = ccall((:PetscViewerView, libpetsc), PetscErrorCode, (CViewer, CViewer), v, viewer) - @assert iszero(error) -end \ No newline at end of file diff --git a/test/linear_system.jl b/test/linear_system.jl index f6c76e7..e6d64d3 100644 --- a/test/linear_system.jl +++ b/test/linear_system.jl @@ -1,92 +1,86 @@ @testset "linear system" begin -# Only on one processor... - -# Import package -using PetscWrap - -# Initialize PETSc. Command line arguments passed to Julia are parsed by PETSc. Alternatively, you can -# also provide "command line arguments by defining them in a string, for instance -# `PetscInitialize("-ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always")` or by providing each argument in -# separate strings : `PetscInitialize(["-ksp_monitor_short", "-ksp_gmres_cgs_refinement_type", "refine_always")` -PetscInitialize() - -# Number of mesh points and mesh step -n = 11 -Δx = 1. / (n - 1) - -# Create a matrix and a vector -A = MatCreate() -b = VecCreate() - -# Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should -# set the number of preallocated non-zeros to increase performance. -MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) -VecSetSizes(b, PETSC_DECIDE, n) - -# We can then use command-line options to set our matrix/vectors. -MatSetFromOptions(A) -VecSetFromOptions(b) - -# Finish the set up -MatSetUp(A) -VecSetUp(b) - -# Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. -# As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first row handled by the local processor -# (starting at 0), and the last row (which is `rend-1`). This may be disturbing for non-regular PETSc users. Checkout -# the fancy version of this example for a more Julia-like convention. -b_start, b_end = VecGetOwnershipRange(b) - -# Now let's build the right hand side vector. Their are various ways to do this, this is just one. -n_loc = VecGetLocalSize(b) # Note that n_loc = b_end - b_start... -VecSetValues(b, collect(b_start:b_end-1), 2 * ones(n_loc)) - -# And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. -A_start, A_end = MatGetOwnershipRange(A) -for i in A_start:A_end-1 - MatSetValues(A, [i], [i-1, i], [-1. 1.] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES) -end - -# Set boundary condition (only the proc handling index `0` is acting) -(b_start == 0) && VecSetValue(b, 0, 0.) - -# Assemble matrices -MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY) -VecAssemblyBegin(b) -MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY) -VecAssemblyEnd(b) - -# At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call -MatView(A) -VecView(b) - -# Set up the linear solver -ksp = KSPCreate() -KSPSetOperators(ksp, A, A) -KSPSetFromOptions(ksp) -KSPSetUp(ksp) - -# Solve the system. We first allocate the solution using the `VecDuplicate` method. -x = VecDuplicate(b) -KSPSolve(ksp, b, x) - -# Print the solution -VecView(x) - -# Access the solution (this part is under development), getting a Julia array; and then restore it -array, ref = VecGetArray(x) # do something with array -@test isapprox(array, range(0., 2.; length = n)) -VecRestoreArray(x, ref) - -# Free memory -MatDestroy(A) -VecDestroy(b) -VecDestroy(x) - -# Finalize Petsc -PetscFinalize() - -# Reach this point? -@test true - + # Number of mesh points and mesh step + n = 11 + Δx = 1.0 / (n - 1) + + println("before create Mat") + + # Create a matrix and a vector + A = create(Mat) + println("before create Vec") + b = create(Vec) + + println("before setSizes") + + # Set the size of the different objects, leaving PETSC to decide how to distribute. Note that we should + # set the number of preallocated non-zeros to increase performance. + setSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n) + setSizes(b, PETSC_DECIDE, n) + + # We can then use command-line options to set our matrix/vectors. + setFromOptions(A) + setFromOptions(b) + + # Finish the set up + setUp(A) + setUp(b) + + println("before getOwnershipRange") + + # Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. + # As in PETSc, the `rstart, rend = *GetOwnershipRange` methods indicate the first row handled by the local processor + # (starting at 0), and the last row (which is `rend-1`). This may be disturbing for non-regular PETSc users. Checkout + # the fancy version of this example for a more Julia-like convention. + b_start, b_end = getOwnershipRange(b) + + # Now let's build the right hand side vector. Their are various ways to do this, this is just one. + n_loc = getLocalSize(b) # Note that n_loc = b_end - b_start... + setValues(b, collect(b_start:(b_end - 1)), 2 * ones(n_loc)) + + # And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. + A_start, A_end = getOwnershipRange(A) + for i = A_start:(A_end - 1) + setValues(A, [i], [i - 1, i], [-1.0 1.0] / Δx, INSERT_VALUES) # MatSetValues(A, I, J, V, INSERT_VALUES) + end + + # Set boundary condition (only the proc handling index `0` is acting) + (b_start == 0) && setValue(b, 0, 0.0) + + # Assemble matrices + assemblyBegin(A, MAT_FINAL_ASSEMBLY) + assemblyBegin(b) + assemblyEnd(A, MAT_FINAL_ASSEMBLY) + assemblyEnd(b) + + println("before matview") + # At this point, you can inspect `A` or `b` using a viewer (stdout by default), simply call + matView(A) + vecView(b) + + # Set up the linear solver + ksp = create(KSP) + setOperators(ksp, A, A) + setFromOptions(ksp) + setUp(ksp) + + # Solve the system. We first allocate the solution using the `VecDuplicate` method. + x = duplicate(b) + solve(ksp, b, x) + + println("before vecview") + # Print the solution + vecView(x) + + # Access the solution (this part is under development), getting a Julia array; and then restore it + array, ref = getArray(x) # do something with array + @test isapprox(array, range(0.0, 2.0; length = n)) + restoreArray(x, ref) + + # Free memory + destroy(A) + destroy(b) + destroy(x) + + # Reach this point? + @test true end diff --git a/test/linear_system_fancy.jl b/test/linear_system_fancy.jl index 206564a..4499917 100644 --- a/test/linear_system_fancy.jl +++ b/test/linear_system_fancy.jl @@ -1,89 +1,67 @@ @testset "linear system fancy" begin -# Only on one processor... - -# Initialize PETSc. Command line arguments passed to Julia are parsed by PETSc. Alternatively, you can -# also provide "command line arguments by defining them in a string, for instance -# `PetscInitialize("-ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always")` or by providing each argument in -# separate strings : `PetscInitialize(["-ksp_monitor_short", "-ksp_gmres_cgs_refinement_type", "refine_always")` -PetscInitialize() - -# Number of mesh points and mesh step -n = 11 -Δx = 1. / (n - 1) - -# Create a matrix of size `(n,n)` and a vector of size `(n)` -A = create_matrix(n, n) -b = create_vector(n) - -# We can then use command line options to set our matrix/vectors. -set_from_options!(A) -set_from_options!(b) - -# Finish the set up -set_up!(A) -set_up!(b) - -# Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. -# The `rstart, rend = get_range(*)` methods differ a little bit from PETSc API since the two integers it -# returns are the effective Julia range of rows handled by the local processor. If `n` is the total -# number of rows, then `rstart ∈ [1,n]` and `rend` is the last row handled by the local processor. -b_start, b_end = get_range(b) - -# Now let's build the right hand side vector. Their are various ways to do this, this is just one. -n_loc = length(b) ## Note that n_loc = b_end - b_start + 1... -b[b_start:b_end] = 2 * ones(n_loc) - -# And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. -A_start, A_end = get_range(A) -for i in A_start:A_end - A[i, i-1:i] = [-1. 1.] / Δx -end - -# Set boundary condition (only the proc handling index `1` is acting) -(b_start == 1) && (b[1] = 0.) - -# Assemble matrice and vector -assemble!(A) -assemble!(b) - -# Set up the linear solver -ksp = create_ksp(A) -set_from_options!(ksp) -set_up!(ksp) - -# Solve the system -x = solve(ksp, b) - -# Print the solution -@show x - -# Convert to Julia array -array = vec2array(x) -@test isapprox(array, range(0., 2.; length = n)) - -# Free memory -destroy!(A) -destroy!(b) -destroy!(x) - -# Note that it's also possible to build a matrix using the COO format as in `SparseArrays`: -M = create_matrix(3,3; auto_setup = true) -M_start, M_end = get_range(M) -I = [1, 1, 1, 2, 3] -J = [1, 3, 1, 3, 2] -V = [1, 2, 3, 4, 5] -k = findall(x -> M_start <= x <= M_end, I) # just a stupid trick to allow this example to run in parallel -set_values!(M, I[k], J[k], V[k], ADD_VALUES) -assemble!(M) -@show M -destroy!(M) -# This is very convenient in sequential since you can fill the three vectors I, J, V in your code and decide only -# at the last moment if you'd like to use `SparseArrays` or `PetscMat`. - -# Finalize Petsc -PetscFinalize() - -# Reach this point? -@test true - + # Only on one processor... + + # Number of mesh points and mesh step + n = 11 + Δx = 1.0 / (n - 1) + + # Create a matrix of size `(n,n)` and a vector of size `(n)` + A = create_matrix(; nrows_glo = n, ncols_glo = n, autosetup = true) + b = create_vector(; nrows_glo = n, autosetup = true) + + # Let's build the right hand side vector. We first get the range of rows of `b` handled by the local processor. + # The `rstart, rend = get_range(*)` methods differ a little bit from PETSc API since the two integers it + # returns are the effective Julia range of rows handled by the local processor. If `n` is the total + # number of rows, then `rstart ∈ [1,n]` and `rend` is the last row handled by the local processor. + b_start, b_end = get_range(b) + + # Now let's build the right hand side vector. Their are various ways to do this, this is just one. + n_loc = length(b) ## Note that n_loc = b_end - b_start + 1... + b[b_start:b_end] = 2 * ones(n_loc) + + # And here is the differentiation matrix. Rembember that PETSc.MatSetValues simply ignores negatives rows indices. + A_start, A_end = get_range(A) + for i = A_start:A_end + A[i, (i - 1):i] = [-1.0 1.0] / Δx + end + + # Set boundary condition (only the proc handling index `1` is acting) + (b_start == 1) && (b[1] = 0.0) + + # Assemble matrice and vector + assemble!(A) + assemble!(b) + + # Set up the linear solver + ksp = create_ksp(A; autosetup = true) + + # Solve the system + x = solve(ksp, b) + + # Print the solution + @show x + + # Convert to Julia array + array = vec2array(x) + @test isapprox(array, range(0.0, 2.0; length = n)) + + # Free memory + destroy!(A, b, x) + + # Note that it's also possible to build a matrix using the COO format as in `SparseArrays`: + M = create_matrix(; nrows_glo = 3, ncols_glo = 3, autosetup = true) + M_start, M_end = get_range(M) + I = [1, 1, 1, 2, 3] + J = [1, 3, 1, 3, 2] + V = [1, 2, 3, 4, 5] + k = findall(x -> M_start <= x <= M_end, I) # just a stupid trick to allow this example to run in parallel + set_values!(M, I[k], J[k], V[k], ADD_VALUES) + assemble!(M) + @show M + destroy!(M) + # This is very convenient in sequential since you can fill the three vectors I, J, V in your code and decide only + # at the last moment if you'd like to use `SparseArrays` or `PetscMat`. + + # Reach this point? + @test true end diff --git a/test/misc.jl b/test/misc.jl index d7826ad..3f71bdf 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -1,52 +1,74 @@ @testset "composite" begin + # Create two matrices + A = create_matrix(; nrows_glo = 2, ncols_glo = 2, autosetup = true) + B = create_matrix(; nrows_glo = 2, ncols_glo = 2, autosetup = true) -PetscInitialize() - -# Create two matrices -A = create_matrix(2, 2; auto_setup = true) -B = create_matrix(2, 2; auto_setup = true) - -# Fill - -# A -# 1 - -# 3 4 -A[1,1] = 1 -A[2,1] = 3 -A[2,2] = 4 - -# B -# 4 3 -# - 1 -B[1,1] = 4 -B[1,2] = 3 -B[2,2] = 1 - -# Assemble matrices -assemble!.((A,B)) - -# Create composite mat from A and B -# C -# 5 3 -# 3 5 -C = create_composite_add([A, B]) - -# Create vectors to check C (see below) -x1 = create_vector(2; auto_setup = true) -x2 = create_vector(2; auto_setup = true) -y = create_vector(2; auto_setup = true) -x1[1] = 1.; x1[2] = 0. -x2[1] = 0.; x2[2] = 1. -assemble!.((x1,x2,y)) - -# Check result by multiplying with test vectors -MatMult(C, x1, y) -@test vec2array(y) == [5., 3.] - -MatMult(C, x2, y) -@test vec2array(y) == [3., 5.] - -# Free memory -destroy!.((A, B, C, x1, x2, y)) -PetscFinalize() + # Fill + + # A + # 1 - + # 3 4 + A[1, 1] = 1 + A[2, 1] = 3 + A[2, 2] = 4 + + # B + # 4 3 + # - 1 + B[1, 1] = 4 + B[1, 2] = 3 + B[2, 2] = 1 + + # Assemble matrices + assemble!.((A, B)) + + # Create composite mat from A and B + # C + # 5 3 + # 3 5 + C = create_composite_add([A, B]) + + # Create vectors to check C (see below) + x1 = create_vector(; nrows_glo = 2, autosetup = true) + x2 = create_vector(; nrows_glo = 2, autosetup = true) + y = create_vector(; nrows_glo = 2, autosetup = true) + x1[1] = 1.0 + x1[2] = 0.0 + x2[1] = 0.0 + x2[2] = 1.0 + assemble!.((x1, x2, y)) + + # Check result by multiplying with test vectors + mult(C, x1, y) + @test vec2array(y) == [5.0, 3.0] + + mult(C, x2, y) + @test vec2array(y) == [3.0, 5.0] + + # Free memory + destroy!(A, B, C, x1, x2, y) +end + +@testset "mapping" begin + n = 10 + x = create_vector(; nrows_glo = n, autosetup = true) + y = create_vector(; nrows_glo = n, autosetup = true) + l2g = reverse(collect(1:n)) + set_local_to_global!(x, l2g) + + set_value_local!(x, 1, 1.0) + set_value!(y, 1, 33) + set_value!(y, 10, 22) + + assemble!.((x, y)) + + # @show dot(x,y) + @test dot(x, x) == 1.0 + @test dot(x, y) == 22.0 + @test sum(x) == 1.0 + + _x = vec2array(x) + @test _x[end] == 1.0 + + destroy!(x, y) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index b6d8d19..9d6c8f8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,20 +1,37 @@ using PetscWrap +using MPI using Test - +using LinearAlgebra # from : # https://discourse.julialang.org/t/what-general-purpose-commands-do-you-usually-end-up-adding-to-your-projects/4889 @generated function compare_struct(x, y) if !isempty(fieldnames(x)) && x == y - mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x)) + mapreduce(n -> :(x.$n == y.$n), (a, b) -> :($a && $b), fieldnames(x)) else :(x == y) end end +""" +Custom way to "include" a file to print infos. +""" +function custom_include(path) + filename = split(path, "/")[end] + print("Running test file " * filename * "...") + include(path) + println("done.") +end + +MPI.Init() +PetscInitialize() +error("Bug : tests are not running anymore when using `test`") @testset "PetscWrap.jl" begin - include("./linear_system.jl") - include("./linear_system_fancy.jl") - include("./misc.jl") + custom_include("./linear_system.jl") + custom_include("./linear_system_fancy.jl") + custom_include("./misc.jl") end + +PetscFinalize() +MPI.Finalize()