From 3f19a2fdcb0ff53f779b430c216409a337aa88b0 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Sun, 29 Sep 2024 03:01:30 +0200 Subject: [PATCH] Move SparseArrays support to an extension --- .github/workflows/Check.yml | 5 ++- CHANGELOG.md | 9 +++- Project.toml | 9 ++-- ext/HYPRESparseArrays.jl | 86 +++++++++++++++++++++++++++++++++++++ src/HYPRE.jl | 66 +++------------------------- src/solvers.jl | 20 --------- 6 files changed, 107 insertions(+), 88 deletions(-) create mode 100644 ext/HYPRESparseArrays.jl diff --git a/.github/workflows/Check.yml b/.github/workflows/Check.yml index df4123c..732e133 100644 --- a/.github/workflows/Check.yml +++ b/.github/workflows/Check.yml @@ -34,19 +34,20 @@ jobs: Pkg.add([ PackageSpec(name = "ExplicitImports", version = "1.9"), PackageSpec(name = "PartitionedArrays"), + PackageSpec(name = "SparseArrays"), PackageSpec(name = "SparseMatricesCSR"), ]) - name: ExplicitImports.jl code checks shell: julia --project=@explicit-imports {0} run: | - using HYPRE, ExplicitImports, PartitionedArrays, SparseMatricesCSR + using HYPRE, ExplicitImports, PartitionedArrays, SparseArrays, SparseMatricesCSR # Check HYPRE check_no_implicit_imports(HYPRE) check_no_stale_explicit_imports(HYPRE) check_all_qualified_accesses_via_owners(HYPRE) check_no_self_qualified_accesses(HYPRE) # Check extension modules - for ext in (:HYPREPartitionedArrays, :HYPRESparseMatricesCSR) + for ext in (:HYPREPartitionedArrays, :HYPRESparseArrays, :HYPRESparseMatricesCSR) extmod = Base.get_extension(HYPRE, ext) if extmod !== nothing check_no_implicit_imports(extmod) diff --git a/CHANGELOG.md b/CHANGELOG.md index ed9a4b2..634d1ed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,8 +11,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 0.5.x. ([#17], [#18]) - CEnum.jl dependency upgraded to release series 0.5.x (release series 0.4.x still allowed). ([#17], [#18]) - - PartitionedArrays.jl support is now moved to a package extension. ([#23]) - - SparseMatricesCSR.jl support is now moved to a package extension. ([#24]) + - PartitionedArrays.jl support (`PSparseMatrix`, `PVector`) is now moved to a package + extension. ([#23]) + - SparseMatricesCSR.jl support (`SparseMatrixCSR`) is now moved to a package extension. + ([#24]) + - SparseArrays.jl support (`SparseMatrixCSC`) is now moved to a package extension. ([#25]) ## [v1.5.0] - 2023-05-26 ### Changed @@ -93,3 +96,5 @@ Initial release of HYPRE.jl. [#17]: https://github.com/fredrikekre/HYPRE.jl/issues/17 [#18]: https://github.com/fredrikekre/HYPRE.jl/issues/18 [#23]: https://github.com/fredrikekre/HYPRE.jl/issues/23 +[#24]: https://github.com/fredrikekre/HYPRE.jl/issues/24 +[#25]: https://github.com/fredrikekre/HYPRE.jl/issues/25 diff --git a/Project.toml b/Project.toml index 629dc3b..f368e49 100644 --- a/Project.toml +++ b/Project.toml @@ -13,11 +13,13 @@ SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" [weakdeps] PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" [extensions] -HYPREPartitionedArrays = ["PartitionedArrays", "SparseMatricesCSR"] -HYPRESparseMatricesCSR = "SparseMatricesCSR" +HYPREPartitionedArrays = ["PartitionedArrays", "SparseArrays", "SparseMatricesCSR"] +HYPRESparseArrays = "SparseArrays" +HYPRESparseMatricesCSR = ["SparseArrays", "SparseMatricesCSR"] [compat] CEnum = "0.4, 0.5" @@ -29,8 +31,9 @@ julia = "1.6" [extras] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["LinearAlgebra", "PartitionedArrays", "SparseMatricesCSR", "Test"] +test = ["LinearAlgebra", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR", "Test"] diff --git a/ext/HYPRESparseArrays.jl b/ext/HYPRESparseArrays.jl new file mode 100644 index 0000000..e09ffbe --- /dev/null +++ b/ext/HYPRESparseArrays.jl @@ -0,0 +1,86 @@ +module HYPRESparseArrays + +using HYPRE.LibHYPRE: @check, HYPRE_BigInt, HYPRE_Complex, HYPRE_Int +using HYPRE: + HYPRE, HYPREMatrix, HYPRESolver, HYPREVector, HYPRE_IJMatrixSetValues, Internals +using MPI: MPI +using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals + +################################## +# SparseMatrixCSC -> HYPREMatrix # +################################## + +function Internals.to_hypre_data(A::SparseMatrixCSC, ilower, iupper) + Internals.check_n_rows(A, ilower, iupper) + nnz = SparseArrays.nnz(A) + A_rows = rowvals(A) + A_vals = nonzeros(A) + + # Initialize the data buffers HYPRE wants + nrows = HYPRE_Int(iupper - ilower + 1) # Total number of rows + ncols = zeros(HYPRE_Int, nrows) # Number of colums for each row + rows = collect(HYPRE_BigInt, ilower:iupper) # The row indices + cols = Vector{HYPRE_BigInt}(undef, nnz) # The column indices + values = Vector{HYPRE_Complex}(undef, nnz) # The values + + # First pass to count nnz per row + @inbounds for j in 1:size(A, 2) + for i in nzrange(A, j) + row = A_rows[i] + ncols[row] += 1 + end + end + + # Keep track of the last index used for every row + lastinds = zeros(Int, nrows) + cumsum!((@view lastinds[2:end]), (@view ncols[1:end-1])) + + # Second pass to populate the output + @inbounds for j in 1:size(A, 2) + for i in nzrange(A, j) + row = A_rows[i] + k = lastinds[row] += 1 + val = A_vals[i] + cols[k] = j + values[k] = val + end + end + @assert nrows == length(ncols) == length(rows) + return nrows, ncols, rows, cols, values +end + +# Note: keep in sync with the SparseMatrixCSR method +function HYPRE.HYPREMatrix(comm::MPI.Comm, B::SparseMatrixCSC, ilower, iupper) + A = HYPREMatrix(comm, ilower, iupper) + nrows, ncols, rows, cols, values = Internals.to_hypre_data(B, ilower, iupper) + @check HYPRE_IJMatrixSetValues(A, nrows, ncols, rows, cols, values) + Internals.assemble_matrix(A) + return A +end + +# Note: keep in sync with the SparseMatrixCSC method +function HYPRE.HYPREMatrix(B::SparseMatrixCSC, ilower=1, iupper=size(B, 1)) + return HYPREMatrix(MPI.COMM_SELF, B, ilower, iupper) +end + + +#################################### +# SparseMatrixCSC solver interface # +#################################### + +# Note: keep in sync with the SparseMatrixCSR method +function HYPRE.solve(solver::HYPRESolver, A::SparseMatrixCSC, b::Vector) + hypre_x = HYPRE.solve(solver, HYPREMatrix(A), HYPREVector(b)) + x = copy!(similar(b, HYPRE_Complex), hypre_x) + return x +end + +# Note: keep in sync with the SparseMatrixCSR method +function HYPRE.solve!(solver::HYPRESolver, x::Vector, A::SparseMatrixCSC, b::Vector) + hypre_x = HYPREVector(x) + HYPRE.solve!(solver, hypre_x, HYPREMatrix(A), HYPREVector(b)) + copy!(x, hypre_x) + return x +end + +end # module HYPRESparseMatricesCSR diff --git a/src/HYPRE.jl b/src/HYPRE.jl index 7e8e105..07f376b 100644 --- a/src/HYPRE.jl +++ b/src/HYPRE.jl @@ -3,7 +3,6 @@ module HYPRE using MPI: MPI -using SparseArrays: SparseArrays, SparseMatrixCSC, nonzeros, nzrange, rowvals export HYPREMatrix, HYPREVector @@ -183,9 +182,10 @@ function Base.zero(b::HYPREVector) return x end -################################## -# SparseMatrixCSC -> HYPREMatrix # -################################## + +######################### +# Vector -> HYPREVector # +######################### function Internals.check_n_rows(A, ilower, iupper) if size(A, 1) != (iupper - ilower + 1) @@ -193,63 +193,6 @@ function Internals.check_n_rows(A, ilower, iupper) end end -function Internals.to_hypre_data(A::SparseMatrixCSC, ilower, iupper) - Internals.check_n_rows(A, ilower, iupper) - nnz = SparseArrays.nnz(A) - A_rows = rowvals(A) - A_vals = nonzeros(A) - - # Initialize the data buffers HYPRE wants - nrows = HYPRE_Int(iupper - ilower + 1) # Total number of rows - ncols = zeros(HYPRE_Int, nrows) # Number of colums for each row - rows = collect(HYPRE_BigInt, ilower:iupper) # The row indices - cols = Vector{HYPRE_BigInt}(undef, nnz) # The column indices - values = Vector{HYPRE_Complex}(undef, nnz) # The values - - # First pass to count nnz per row - @inbounds for j in 1:size(A, 2) - for i in nzrange(A, j) - row = A_rows[i] - ncols[row] += 1 - end - end - - # Keep track of the last index used for every row - lastinds = zeros(Int, nrows) - cumsum!((@view lastinds[2:end]), (@view ncols[1:end-1])) - - # Second pass to populate the output - @inbounds for j in 1:size(A, 2) - for i in nzrange(A, j) - row = A_rows[i] - k = lastinds[row] += 1 - val = A_vals[i] - cols[k] = j - values[k] = val - end - end - @assert nrows == length(ncols) == length(rows) - return nrows, ncols, rows, cols, values -end - -# Note: keep in sync with the SparseMatrixCSR method -function HYPREMatrix(comm::MPI.Comm, B::SparseMatrixCSC, ilower, iupper) - A = HYPREMatrix(comm, ilower, iupper) - nrows, ncols, rows, cols, values = Internals.to_hypre_data(B, ilower, iupper) - @check HYPRE_IJMatrixSetValues(A, nrows, ncols, rows, cols, values) - Internals.assemble_matrix(A) - return A -end - -# Note: keep in sync with the SparseMatrixCSC method -function HYPREMatrix(B::SparseMatrixCSC, ilower=1, iupper=size(B, 1)) - return HYPREMatrix(MPI.COMM_SELF, B, ilower, iupper) -end - -######################### -# Vector -> HYPREVector # -######################### - function Internals.to_hypre_data(x::Vector, ilower, iupper) Internals.check_n_rows(x, ilower, iupper) indices = collect(HYPRE_BigInt, ilower:iupper) @@ -475,6 +418,7 @@ include("solver_options.jl") # Compatibility for Julias that doesn't support package extensions if !(isdefined(Base, :get_extension)) include("../ext/HYPREPartitionedArrays.jl") + include("../ext/HYPRESparseArrays.jl") include("../ext/HYPRESparseMatricesCSR.jl") end diff --git a/src/solvers.jl b/src/solvers.jl index c242707..be472fa 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -51,26 +51,6 @@ See also [`solve`](@ref). solve!(pcg::HYPRESolver, x::HYPREVector, A::HYPREMatrix, ::HYPREVector) -#################################### -# SparseMatrixCSC solver interface # -#################################### - -# Note: keep in sync with the SparseMatrixCSR method -function solve(solver::HYPRESolver, A::SparseMatrixCSC, b::Vector) - hypre_x = solve(solver, HYPREMatrix(A), HYPREVector(b)) - x = copy!(similar(b, HYPRE_Complex), hypre_x) - return x -end - -# Note: keep in sync with the SparseMatrixCSR method -function solve!(solver::HYPRESolver, x::Vector, A::SparseMatrixCSC, b::Vector) - hypre_x = HYPREVector(x) - solve!(solver, hypre_x, HYPREMatrix(A), HYPREVector(b)) - copy!(x, hypre_x) - return x -end - - ##################################### ## Concrete solver implementations ## #####################################