diff --git a/docs/src/api.md b/docs/src/api.md index 41b50fa..75e12fb 100644 --- a/docs/src/api.md +++ b/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 diff --git a/docs/src/matrix-vector.md b/docs/src/matrix-vector.md index c9ce5b0..bc148d7 100644 --- a/docs/src/matrix-vector.md +++ b/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 -`HYPREVector`, respectively. See the following sections for details: +`HYPREMatrix` and `HYPREVector` can be constructed either by assembling directly, or by +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. diff --git a/src/HYPRE.jl b/src/HYPRE.jl index 7c973f5..a88bab3 100644 --- a/src/HYPRE.jl +++ b/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") diff --git a/test/runtests.jl b/test/runtests.jl index b6baa88..1b9da60 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_assembler.jl b/test/test_assembler.jl new file mode 100644 index 0000000..3c91eb5 --- /dev/null +++ b/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 diff --git a/test/test_utils.jl b/test/test_utils.jl new file mode 100644 index 0000000..25a81ee --- /dev/null +++ b/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