Matrix free support

.github/workflows/Documentation.yml

@@ -10,17 +10,23 @@ on:
 
 jobs:
   build:
+    permissions:
+      actions: write
+      contents: write
+      pull-requests: read
+      statuses: write
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@latest
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
-          version: "1.8"
+          version: '1'
+      - uses: julia-actions/cache@v2
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
         env:
-          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
-          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
           DOC_DEPLOYMENT: true # Workaround to avoid finding PETSc location
         run: julia --project=docs/ docs/make.jl

Project.toml

@@ -12,4 +12,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 MPI = "0.20"
-julia = "1.5, 1.6, 1.7, 1.8"
+julia = "1.5"

docs/make.jl

@@ -1,46 +1,46 @@
-push!(LOAD_PATH,"../src/")
+push!(LOAD_PATH, "../src/")
 
 using PetscWrap
 using Documenter
 using Literate
 
 # Generate examples
-example_src = joinpath(@__DIR__,"..","example")
-example_dir = joinpath(@__DIR__,"src","example")
-Sys.rm(example_dir; recursive=true, force=true)
-Literate.markdown(joinpath(example_src, "linear_system.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
-Literate.markdown(joinpath(example_src, "linear_system_fancy.jl"), example_dir; documenter = false, execute = false) # documenter = false to avoid Documenter to execute cells
+example_src = joinpath(@__DIR__, "..", "example")
+example_dir = joinpath(@__DIR__, "src", "example")
+Sys.rm(example_dir; recursive = true, force = true)
+Literate.markdown(
+    joinpath(example_src, "linear_system.jl"),
+    example_dir;
+    documenter = false,
+    execute = false,
+) # documenter = false to avoid Documenter to execute cells
+Literate.markdown(
+    joinpath(example_src, "linear_system_fancy.jl"),
+    example_dir;
+    documenter = false,
+    execute = false,
+) # documenter = false to avoid Documenter to execute cells
 
 makedocs(;
-    modules=[PetscWrap],
-    authors="bmxam",
-    sitename="PetscWrap.jl",
-    clean=true,doctest=false,
-    format=Documenter.HTML(;
-        prettyurls=get(ENV, "CI", "false") == "true",
-        canonical="https://github.com/bmxam/PetscWrap.jl",
-        assets=String[],
+    modules = [PetscWrap],
+    authors = "bmxam",
+    sitename = "PetscWrap.jl",
+    clean = true,
+    doctest = false,
+    checkdocs = :none,
+    format = Documenter.HTML(;
+        prettyurls = get(ENV, "CI", "false") == "true",
+        canonical = "https://github.com/bmxam/PetscWrap.jl",
+        assets = String[],
     ),
-    pages=[
+    pages = [
         "Home" => "index.md",
-        "Examples" => Any[
-            "example/linear_system.md",
-        ],
-        "Fancy examples" => Any[
-            "example/linear_system_fancy.md",
-        ],
-        "API Reference" => Any[
-            "api/init.md",
-            "api/viewer.md",
-            "api/vec.md",
-            "api/mat.md",
-            "api/ksp.md",
-        ],
+        "Examples" => Any["example/linear_system.md",],
+        "Fancy examples" => Any["example/linear_system_fancy.md",],
+        "API Reference" =>
+            Any["api/init.md", "api/viewer.md", "api/vec.md", "api/mat.md", "api/ksp.md"],
         "API fancy" => "api/fancy/fancy.md",
     ],
 )
 
-deploydocs(;
-    repo="github.com/bmxam/PetscWrap.jl.git",
-    push_preview = true
-)
+deploydocs(; repo = "github.com/bmxam/PetscWrap.jl.git", push_preview = true)

example/matrix_free.jl

@@ -0,0 +1,133 @@
+module Example #hide
+using MPI
+using MPIUtils
+using PetscWrap
+using LinearMaps
+using HauntedArrays
+using HauntedArrays2PetscWrap
+
+const disable_gc = false
+
+disable_gc && GC.enable(false)
+
+MPI.Initialized() || MPI.Init()
+
+comm = MPI.COMM_WORLD
+rank = MPI.Comm_rank(comm)
+np = MPI.Comm_size(comm)
+
+const lx = 1.0
+const nx = 10 ÷ np # on each process
+
+mypart = rank + 1
+Δx = lx / (np * nx - 1)
+
+# The "mesh" partitioning
+lid2gid, lid2part = HauntedArrays.generate_1d_partitioning(nx, mypart, np)
+
+x_h = HauntedVector(comm, lid2gid, lid2part; cacheType = PetscCache)
+y_h = similar(x_h)
+
+f!(y::HauntedVector, x::HauntedVector) =
+    for li in eachindex(x)
+        gi = local_to_global(x, li)
+        if gi == 1
+            y[li] = x[li]
+        else
+            for lj in eachindex(x)
+                gj = local_to_global(x, lj)
+                if gj == gi - 1
+                    y[li] = (x[li] - x[lj]) / Δx
+                end
+            end
+        end
+    end
+
+function f!(y, x)
+    println("calling custom mult operator...")
+
+    # Update HauntedVectors values `x_h` with Petsc values `x`
+    println("updating HauntedVector...")
+    HauntedArrays2PetscWrap.update!(x_h, x)
+    HauntedArrays.update_ghosts!(x_h)
+
+    # Apply the operator to obtain a result as a HauntedVector (`y_h`)
+    println("Applying operator...")
+    f!(y_h, x_h)
+
+    # Update Petsc values `y` with the HauntedVector values (`y_h`)
+    println("Updating Petsc vector...")
+    HauntedArrays2PetscWrap.update!(y, y_h)
+
+    # MPI.Barrier(PetscWrap._get_comm(y))
+    println("done ! ")
+end
+
+# Create a matrix and a vector
+b = create(Vec, comm)
+
+A = PetscWrap.createShell(comm, Int32(nx), Int32(nx))
+cfunc = PetscWrap.set_shell_mul!(A, f!)
+
+# function _f!(::CMat, x::CVec, y::CVec)::Cint
+#     comm = PetscWrap._get_comm(A)
+#     f!(Vec(comm, y), Vec(comm, x))
+#     return PetscErrorCode(0) # return "success" (mandatory)
+# end
+# shell_mul_c = @cfunction($_f!, Cint, (CMat, CVec, CVec))
+# PetscWrap.shellSetOperation(A, PetscWrap.MATOP_MULT, shell_mul_c)
+
+setSizes(b, PETSC_DECIDE, nx * np)
+
+# We can then use command-line options to set our matrix/vectors.
+setFromOptions(A)
+setFromOptions(b)
+
+# Finish the set up
+setUp(A)
+setUp(b)
+
+b_start, b_end = getOwnershipRange(b)
+n_loc = getLocalSize(b) # Note that n_loc = b_end - b_start...
+setValues(b, collect(b_start:(b_end - 1)), 2 * ones(n_loc))
+
+# 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)
+
+# 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
+println("setting up")
+ksp = create(KSP, comm)
+setOperators(ksp, A, A)
+setFromOptions(ksp)
+setUp(ksp)
+
+# Solve the system. We first allocate the solution using the `VecDuplicate` method.
+x = duplicate(b)
+println("solving")
+solve(ksp, b, x)
+
+# Print the solution
+println("viewing")
+vecView(x)
+
+@only_root @show [2 * (i - 1) * Δx for i = 1:(np * nx)]
+
+# 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))
+
+# x_h .= 0.0
+# y_h .= 0.0
+disable_gc && GC.enable(true)
+
+end #hide

src/KSP.jl

@@ -4,7 +4,7 @@ mutable struct KSP <: AbstractPetscObject
     ptr::CKSP
     comm::MPI.Comm
 
-    KSP(comm::MPI.Comm) = new(CKSP(), comm)
+    KSP(comm::MPI.Comm, ptr = CKSP()) = new(ptr, comm)
 end
 
 Base.unsafe_convert(::Type{CKSP}, x::KSP) = x.ptr

src/Mat.jl

@@ -7,7 +7,7 @@ mutable struct Mat <: AbstractPetscObject
     ptr::CMat
     comm::MPI.Comm
 
-    Mat(comm::MPI.Comm) = new(CMat(), comm)
+    Mat(comm::MPI.Comm, ptr = CMat()) = new(ptr, comm)
 end
 
 Base.unsafe_convert(::Type{CMat}, x::Mat) = x.ptr
@@ -79,18 +79,10 @@ 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;
+        m::PetscInt = PETSC_DECIDE,
+        n::PetscInt = PETSC_DECIDE,
+        M::PetscInt = PETSC_DECIDE,
+        N::PetscInt = PETSC_DECIDE;
         add_finalizer = true,
     )
 
@@ -101,10 +93,10 @@ Last argument `data` is not supported yet (NULL is passed).
 """
 function createDense(
     comm::MPI.Comm,
-    m::PetscInt,
-    n::PetscInt,
-    M::PetscInt,
-    N::PetscInt;
+    m::PetscInt = PETSC_DECIDE,
+    n::PetscInt = PETSC_DECIDE,
+    M::PetscInt = PETSC_DECIDE,
+    N::PetscInt = PETSC_DECIDE;
     add_finalizer = true,
 )
     mat = Mat(comm)
@@ -128,22 +120,48 @@ function createDense(
     return mat
 end
 
-function createDense(
+"""
+    createShell(
+        comm::MPI.Comm,
+        m::PetscInt,
+        n::PetscInt,
+        M::PetscInt = PETSC_DETERMINE,
+        N::PetscInt = PETSC_DETERMINE;
+        add_finalizer = true,
+    )
+
+Wrapper to `MatCreateShell`
+https://petsc.org/release/manualpages/Mat/MatCreateShell/
+
+Last argument `ctx` is not supported yet (NULL is passed).
+"""
+function createShell(
     comm::MPI.Comm,
-    m::Integer = PETSC_DECIDE,
-    n::Integer = PETSC_DECIDE,
-    M::Integer = PETSC_DECIDE,
-    N::Integer = PETSC_DECIDE;
+    m::PetscInt,
+    n::PetscInt,
+    M::PetscInt = PETSC_DETERMINE,
+    N::PetscInt = PETSC_DETERMINE;
     add_finalizer = true,
 )
-    return createDense(
+    mat = Mat(comm)
+    error = ccall(
+        (:MatCreateShell, libpetsc),
+        PetscErrorCode,
+        (MPI.MPI_Comm, PetscInt, PetscInt, PetscInt, PetscInt, Ptr{Cvoid}, Ptr{CMat}),
         comm,
-        PetscInt(m),
-        PetscInt(n),
-        PetscInt(M),
-        PetscInt(N);
-        add_finalizer,
+        m,
+        n,
+        M,
+        N,
+        C_NULL,
+        mat,
     )
+    @assert iszero(error)
+
+    _NREFS[] += 1
+    add_finalizer && finalizer(destroy, mat)
+
+    return mat
 end
 
 """
@@ -770,6 +788,27 @@ function setValuesCOO(A::Mat, coo_v::Vector{PetscScalar}, imode::InsertMode)
     @assert iszero(error)
 end
 
+"""
+    shellSetOperation(mat::Mat, op::MatOperation, g)
+
+`g` must have been built with `@cfunction`. For a more convenient API,
+see `src/fancy/mat.jl->set_shell_mul!`.
+
+Wrapper to `MatShellSetOperation`
+https://petsc.org/release/manualpages/Mat/MatShellSetOperation/
+"""
+function shellSetOperation(mat::Mat, op::MatOperation, g)
+    error = ccall(
+        (:MatShellSetOperation, libpetsc),
+        PetscErrorCode,
+        (CMat, MatOperation, Ptr{Cvoid}),
+        mat,
+        op,
+        g,
+    )
+    @assert iszero(error)
+end
+
 """
     matView(mat::Mat, viewer::PetscViewer = StdWorld())