Add assembler interface to assemble HYPRE(Vector|Matrix) directly (#5)

This patch adds an assembler interface such that it is possible to
directly assemble HYPRE(Vector|Matrix) without going through a sparse
matrix datastructure in Julia. This is done as follows:

 1. Create a new (empty) matrix/vector using the constructor.
 2. Create an assembler and initialize the assembly using
    HYPRE.start_assemble!.
 3. Assemble all contributions using HYPRE.assemble!.
 4. Finalize the assembly using HYPRE.finish_assemble!.
    HYPRE.start_assemble!

The assembler caches some buffers that are (re)used by every call to
HYPRE.assemble! so this should be efficient. All MPI communication
should happen in the finalization step.

docs/src/api.md

@@ -6,6 +6,14 @@
 HYPRE.Init
 ```
 
+## Matrix/vector creation
+
+```@docs
+HYPRE.start_assemble!
+HYPRE.assemble!
+HYPRE.finish_assemble!
+```
+
 ## Solvers and preconditioners
 
 ```@docs

docs/src/matrix-vector.md

@@ -5,15 +5,59 @@ datastructures. Specifically it uses the [IJ System
 Interface](https://hypre.readthedocs.io/en/latest/api-int-ij.html) which can be used for
 general sparse matrices.
 
-HYPRE.jl defines conversion methods from standard Julia datastructures to `HYPREMatrix` and
+`HYPREMatrix` and `HYPREVector` can be constructed either by assembling directly, or by
-`HYPREVector`, respectively. See the following sections for details:
+first assembling into a Julia datastructure and the converting it. These various methods are
+outlined in the following sections:
 
 ```@contents
-Pages = ["hypre-matrix-vector.md"]
+Pages = ["matrix-vector.md"]
 Depth = 2:2
 ```
 
-## PartitionedArrays.jl (multi-process)
+
+## Direct assembly (multi-/single-process)
+
+Creating `HYPREMatrix` and/or `HYPREVector` directly is possible by first creating an
+assembler which is used to add all individual contributions to the matrix/vector. The
+required steps are:
+
+  1. Create a new matrix and/or vector using the constructor.
+  2. Create an assembler and initialize the assembling procedure using
+     [`HYPRE.start_assemble!`](@ref).
+  3. Assemble all non-zero contributions (e.g. element matrix/vector in a finite element
+     simulation) using [`HYPRE.assemble!`](@ref).
+  4. Finalize the assembly using [`HYPRE.finish_assemble!`](@ref).
+
+After these steps the matrix and vector are ready to pass to the solver. In case of multiple
+consecutive solves with the same sparsity pattern (e.g. multiple Newton steps, multiple time
+steps, ...) it is possible to reuse the same matrix by simply skipping the first step above.
+
+**Example pseudocode**
+
+```julia
+# MPI communicator
+comm = MPI.COMM_WORLD # MPI.COMM_SELF for single-process setups
+
+# Create empty matrix and vector -- this process owns rows ilower to iupper
+A = HYPREMatrix(comm, ilower, iupper)
+b = HYPREVector(comm, ilower, iupper)
+
+# Create assembler
+assembler = HYPRE.start_assemble!(A, b)
+
+# Assemble contributions from all elements owned by this process
+for element in owned_elements
+    Ae, be = compute_element_contribution(...)
+    global_indices = get_global_indices(...)
+    HYPRE.assemble!(assembler, global_indices, Ae, be)
+end
+
+# Finalize the assembly
+A, b = HYPRE.finish_assemble!(assembler)
+```
+
+
+## Create from PartitionedArrays.jl (multi-process)
 
 HYPRE.jl integrates seemlessly with `PSparseMatrix` and `PVector` from the
 [PartitionedArrays.jl](https://github.com/fverdugo/PartitionedArrays.jl) package. These can
@@ -71,7 +115,7 @@ copy!(x, x_h)
 ```
 
 
-## `SparseMatrixCSC` / `SparseMatrixCSR` (single-process)
+## Create from `SparseMatrixCSC` / `SparseMatrixCSR` (single-process)
 
 HYPRE.jl also support working directly with `SparseMatrixCSC` (from the
 [SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl) standard library) and
@@ -100,10 +144,3 @@ x = solve(solver, A, b)
 x = zeros(length(b))
 solve!(solver, x, A, b)
 ```
-
-
-## `SparseMatrixCSC` / `SparseMatrixCSR` (multi-process)
-
-!!! warning
-    This interface isn't finalized yet and is therefore not documented since it
-    is subject to change.

src/HYPRE.jl

@@ -524,7 +524,174 @@ function Base.copy!(dst::HYPREVector, src::PVector{HYPRE_Complex})
     return dst
 end
 
-# Solver interface
+
+####################
+## HYPREAssembler ##
+####################
+
+struct HYPREMatrixAssembler
+    A::HYPREMatrix
+    ncols::Vector{HYPRE_Int}
+    rows::Vector{HYPRE_BigInt}
+    cols::Vector{HYPRE_BigInt}
+    values::Vector{HYPRE_Complex}
+end
+
+struct HYPREVectorAssembler
+    b::HYPREVector
+    indices::Vector{HYPRE_BigInt}
+    values::Vector{HYPRE_Complex}
+end
+
+struct HYPREAssembler
+    A::HYPREMatrixAssembler
+    b::HYPREVectorAssembler
+end
+
+"""
+    HYPRE.start_assemble!(A::HYPREMatrix)                 -> HYPREMatrixAssembler
+    HYPRE.start_assemble!(b::HYPREVector)                 -> HYPREVectorAssembler
+    HYPRE.start_assemble!(A::HYPREMatrix, b::HYPREVector) -> HYPREAssembler
+
+Initialize a new assembly for matrix `A`, vector `b`, or for both. This zeroes out any
+previous data in the arrays. Return a `HYPREAssembler` with allocated data buffers needed to
+perform the assembly efficiently.
+
+See also: [`HYPRE.assemble!`](@ref), [`HYPRE.finish_assemble!`](@ref).
+"""
+start_assemble!
+
+function start_assemble!(A::HYPREMatrix)
+    if A.parmatrix != C_NULL
+        # This matrix have been assembled before, reset to 0
+        @check HYPRE_IJMatrixSetConstantValues(A.ijmatrix, 0)
+    end
+    @check HYPRE_IJMatrixInitialize(A.ijmatrix)
+    return HYPREMatrixAssembler(A, HYPRE_Int[], HYPRE_BigInt[], HYPRE_BigInt[], HYPRE_Complex[])
+end
+
+function start_assemble!(b::HYPREVector)
+    if b.parvector != C_NULL
+        # This vector have been assembled before, reset to 0
+        # See https://github.com/hypre-space/hypre/pull/689
+        # @check HYPRE_IJVectorSetConstantValues(b.ijvector, 0)
+    end
+    @check HYPRE_IJVectorInitialize(b.ijvector)
+    if b.parvector != C_NULL
+        nvalues = HYPRE_Int(b.iupper - b.ilower + 1)
+        indices = collect(HYPRE_BigInt, b.ilower:b.iupper)
+        values = zeros(HYPRE_Complex, nvalues)
+        @check HYPRE_IJVectorSetValues(b.ijvector, nvalues, indices, values)
+        # TODO: Do I need to assemble here?
+    end
+    return HYPREVectorAssembler(b, HYPRE_BigInt[], HYPRE_Complex[])
+end
+
+function start_assemble!(A::HYPREMatrix, b::HYPREVector)
+    return HYPREAssembler(start_assemble!(A), start_assemble!(b))
+end
+
+"""
+    HYPRE.assemble!(A::HYPREMatrixAssembler, ij, a::Matrix)
+    HYPRE.assemble!(A::HYPREVectorAssembler, ij, b::Vector)
+    HYPRE.assemble!(A::HYPREAssembler,       ij, a::Matrix, b::Vector)
+
+Assemble (by adding) matrix contribution `a`, vector contribution `b`, into the underlying
+array(s) of the assembler at global row and column indices `ij`.
+
+This is roughly equivalent to:
+```julia
+# A.A::HYPREMatrix
+A.A[ij, ij] += a
+
+# A.b::HYPREVector
+A.b[ij] += b
+```
+
+See also: [`HYPRE.start_assemble!`](@ref), [`HYPRE.finish_assemble!`](@ref).
+"""
+assemble!
+
+function assemble!(A::HYPREMatrixAssembler, ij::Vector, a::Matrix)
+    nrows, ncols, rows, cols, values = Internals.to_hypre_data(A, a, ij, ij)
+    @check HYPRE_IJMatrixAddToValues(A.A.ijmatrix, nrows, ncols, rows, cols, values)
+    return A
+end
+
+function assemble!(A::HYPREVectorAssembler, ij::Vector, a::Vector)
+    nvalues, indices, values = Internals.to_hypre_data(A, a, ij)
+    @check HYPRE_IJVectorAddToValues(A.b.ijvector, nvalues, indices, values)
+    return A
+end
+
+function assemble!(A::HYPREAssembler, ij::Vector, a::Matrix, b::Vector)
+    assemble!(A.A, ij, a)
+    assemble!(A.b, ij, b)
+    return A
+end
+
+function Internals.to_hypre_data(A::HYPREMatrixAssembler, a::Matrix, I::Vector, J::Vector)
+    size(a, 1) == length(I) || error("mismatching number of rows")
+    size(a, 2) == length(J) || error("mismatch number of cols")
+    nrows = HYPRE_Int(length(I))
+    # Resize cache vectors
+    ncols = resize!(A.ncols, nrows)
+    rows = resize!(A.rows, length(a))
+    cols = resize!(A.cols, length(a))
+    values = resize!(A.values, length(a))
+    # Fill vectors
+    ncols = fill!(ncols, HYPRE_Int(length(J)))
+    copyto!(rows, I)
+    idx = 0
+    for i in 1:length(I), j in 1:length(J)
+        idx += 1
+        cols[idx] = J[j]
+        values[idx] = a[i, j]
+    end
+    @assert idx == length(a)
+    return nrows, ncols, rows, cols, values
+end
+
+function Internals.to_hypre_data(A::HYPREVectorAssembler, b::Vector, I::Vector)
+    length(b) == length(I) || error("mismatching number of entries")
+    nvalues = HYPRE_Int(length(I))
+    # Resize cache vectors
+    indices = resize!(A.indices, nvalues)
+    values = resize!(A.values, nvalues)
+    # Fill vectors
+    copyto!(indices, I)
+    copyto!(values, b)
+    return nvalues, indices, values
+end
+
+"""
+    HYPRE.finish_assemble!(A::HYPREMatrixAssembler)
+    HYPRE.finish_assemble!(A::HYPREVectorAssembler)
+    HYPRE.finish_assemble!(A::HYPREAssembler)
+
+Finish the assembly. This synchronizes the data between processors.
+"""
+finish_assemble!
+
+function finish_assemble!(A::HYPREMatrixAssembler)
+    Internals.assemble_matrix(A.A)
+    return A.A
+end
+
+function finish_assemble!(A::HYPREVectorAssembler)
+    Internals.assemble_vector(A.b)
+    return A.b
+end
+
+function finish_assemble!(A::HYPREAssembler)
+    return finish_assemble!(A.A), finish_assemble!(A.b)
+end
+
+
+######################
+## Solver interface ##
+######################
+
 include("solvers.jl")
 include("solver_options.jl")
 

test/runtests.jl

@@ -10,6 +10,8 @@ using SparseArrays
 using SparseMatricesCSR
 using Test
 
+include("test_utils.jl")
+
 # Init HYPRE and MPI
 HYPRE.Init()
 
@@ -296,6 +298,62 @@ end
     @test tomain(pbc) == tomain(pb)
 end
 
+@testset "HYPRE(Matrix|Vector)?Assembler" begin
+    comm = MPI.COMM_WORLD
+    # Assembly HYPREMatrix from ::Matrix
+    A = HYPREMatrix(comm, 1, 3)
+    AM = zeros(3, 3)
+    for i in 1:2
+        assembler = HYPRE.start_assemble!(A)
+        fill!(AM, 0)
+        for idx in ([1, 2], [3, 1])
+            a = rand(2, 2)
+            HYPRE.assemble!(assembler, idx, a)
+            AM[idx, idx] += a
+        end
+        f = HYPRE.finish_assemble!(assembler)
+        @test f === A
+        @test getindex_debug(A, 1:3, 1:3) == AM
+    end
+    # Assembly HYPREVector from ::Vector
+    b = HYPREVector(comm, 1, 3)
+    bv = zeros(3)
+    for i in 1:2
+        assembler = HYPRE.start_assemble!(b)
+        fill!(bv, 0)
+        for idx in ([1, 2], [3, 1])
+            c = rand(2)
+            HYPRE.assemble!(assembler, idx, c)
+            bv[idx] += c
+        end
+        f = HYPRE.finish_assemble!(assembler)
+        @test f === b
+        @test getindex_debug(b, 1:3) == bv
+    end
+    # Assembly HYPREMatrix/HYPREVector from ::Array
+    A = HYPREMatrix(comm, 1, 3)
+    AM = zeros(3, 3)
+    b = HYPREVector(comm, 1, 3)
+    bv = zeros(3)
+    for i in 1:2
+        assembler = HYPRE.start_assemble!(A, b)
+        fill!(AM, 0)
+        fill!(bv, 0)
+        for idx in ([1, 2], [3, 1])
+            a = rand(2, 2)
+            c = rand(2)
+            HYPRE.assemble!(assembler, idx, a, c)
+            AM[idx, idx] += a
+            bv[idx] += c
+        end
+        F, f = HYPRE.finish_assemble!(assembler)
+        @test F === A
+        @test f === b
+        @test getindex_debug(A, 1:3, 1:3) == AM
+        @test getindex_debug(b, 1:3) == bv
+    end
+end
+
 @testset "BiCGSTAB" begin
     # Solver constructor and options
     @test_throws(
@@ -654,3 +712,15 @@ end
     x = HYPRE.solve(pcg, A, b)
     @test x ≈ A \ b atol=tol
 end
+
+@testset "MPI execution" begin
+    testfiles = joinpath.(@__DIR__, [
+        "test_assembler.jl",
+    ])
+    for file in testfiles
+        mpiexec() do mpi
+            r = run(ignorestatus(`$(mpi) -n 2 $(Base.julia_cmd()) $(file)`))
+            @test r.exitcode == 0
+        end
+    end
+end

test/test_assembler.jl

@@ -0,0 +1,115 @@
+# SPDX-License-Identifier: MIT
+
+using HYPRE
+using MPI
+using Test
+
+MPI.Init()
+HYPRE.Init()
+
+include("test_utils.jl")
+
+comm = MPI.COMM_WORLD
+comm_rank = MPI.Comm_rank(comm)
+comm_size = MPI.Comm_size(comm)
+
+if comm_size != 2
+    error("Must run with 2 ranks.")
+end
+
+if comm_rank == 0
+    ilower = 1
+    iupper = 10
+    N = 2:10
+else
+    ilower = 11
+    iupper = 20
+    N = 11:19
+end
+
+function values_and_indices(n)
+    idx = [n - 1, n, n + 1]
+    a = Float64[
+          n -2n  -n
+        -2n   n -2n
+         -n -2n   n
+    ]
+    b = Float64[n, n/2, n/3]
+    return idx, a, b
+end
+
+##########################
+## HYPREMatrixAssembler ##
+##########################
+
+# Dense local matrix
+
+A = HYPREMatrix(comm, ilower, iupper)
+AM = zeros(20, 20)
+for i in 1:2
+    assembler = HYPRE.start_assemble!(A)
+    fill!(AM, 0)
+    for n in N
+        idx, a, _ = values_and_indices(n)
+        HYPRE.assemble!(assembler, idx, a)
+        AM[idx, idx] += a
+    end
+    f = HYPRE.finish_assemble!(assembler)
+    @test f === A
+    MPI.Allreduce!(AM, +, comm)
+    @test getindex_debug(A, ilower:iupper, 1:20) == AM[ilower:iupper, 1:20]
+    MPI.Barrier(comm)
+end
+
+##########################
+## HYPREVectorAssembler ##
+##########################
+
+# Dense local vector
+
+b = HYPREVector(comm, ilower, iupper)
+bv = zeros(20)
+for i in 1:2
+    assembler = HYPRE.start_assemble!(b)
+    fill!(bv, 0)
+    for n in N
+        idx, _, a = values_and_indices(n)
+        HYPRE.assemble!(assembler, idx, a)
+        bv[idx] += a
+    end
+    f = HYPRE.finish_assemble!(assembler)
+    @test f === b
+    MPI.Allreduce!(bv, +, comm)
+    @test getindex_debug(b, ilower:iupper) == bv[ilower:iupper]
+    MPI.Barrier(comm)
+end
+
+####################
+## HYPREAssembler ##
+####################
+
+# Dense local arrays
+
+A = HYPREMatrix(comm, ilower, iupper)
+AM = zeros(20, 20)
+b = HYPREVector(comm, ilower, iupper)
+bv = zeros(20)
+for i in 1:2
+    assembler = HYPRE.start_assemble!(A, b)
+    fill!(AM, 0)
+    fill!(bv, 0)
+    for n in N
+        idx, a, c = values_and_indices(n)
+        HYPRE.assemble!(assembler, idx, a, c)
+        AM[idx, idx] += a
+        bv[idx] += c
+    end
+    F, f = HYPRE.finish_assemble!(assembler)
+    @test F === A
+    @test f === b
+    MPI.Allreduce!(AM, +, comm)
+    MPI.Allreduce!(bv, +, comm)
+    @test getindex_debug(A, ilower:iupper, 1:20) == AM[ilower:iupper, 1:20]
+    @test getindex_debug(b, ilower:iupper) == bv[ilower:iupper]
+    MPI.Barrier(comm)
+end

test/test_utils.jl

@@ -0,0 +1,23 @@
+# SPDX-License-Identifier: MIT
+
+using HYPRE
+using HYPRE.LibHYPRE
+using HYPRE.LibHYPRE: @check
+
+function getindex_debug(A::HYPREMatrix, i::AbstractVector, j::AbstractVector)
+    nrows = HYPRE_Int(length(i))
+    ncols = fill(HYPRE_Int(length(j)), length(i))
+    rows = convert(Vector{HYPRE_BigInt}, i)
+    cols = convert(Vector{HYPRE_BigInt}, repeat(j, length(i)))
+    values = Vector{HYPRE_Complex}(undef, length(i) * length(j))
+    @check HYPRE_IJMatrixGetValues(A.ijmatrix, nrows, ncols, rows, cols, values)
+    return permutedims(reshape(values, (length(j), length(i))))
+end
+
+function getindex_debug(b::HYPREVector, i::AbstractVector)
+    nvalues = HYPRE_Int(length(i))
+    indices = convert(Vector{HYPRE_BigInt}, i)
+    values = Vector{HYPRE_Complex}(undef, length(i))
+    @check HYPRE_IJVectorGetValues(b.ijvector, nvalues, indices, values)
+    return values
+end