From b757195b18c82cf4ba02bc4186831b351eed3522 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 11 Mar 2026 22:38:27 +1100 Subject: [PATCH 1/3] Minor bugfixes in BlockPartitionedArrays --- NEWS.md | 7 + src/Algebra.jl | 2 +- src/BlockPartitionedArrays.jl | 2 +- test/BlockPartitionedArraysTests.jl | 169 +++++++----------- test/TestApp/src/TestApp.jl | 1 + test/mpi/runtests_np4_body.jl | 1 + .../BlockSparseMatrixAssemblersTests.jl | 4 + test/sequential/runtests.jl | 4 +- 8 files changed, 81 insertions(+), 109 deletions(-) diff --git a/NEWS.md b/NEWS.md index d791a38a..cacf25a5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,13 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Fixed + +- Fixed `BlockPMatrix{V}(::UndefInitializer, rows, cols)` constructor dropping the `cols` argument, causing a `MethodError` at runtime. +- Fixed `local_views(::BlockPMatrix, rows, cols)` indexing 1D block-range vectors with a 2D `CartesianIndex`, causing `BoundsError` for any multi-field problem with ≥2 fields. + ## [0.4.11] - 2026-02-20 ### Added diff --git a/src/Algebra.jl b/src/Algebra.jl index 772cea1a..8964fc88 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -387,7 +387,7 @@ end function local_views(a::BlockPMatrix,new_rows::BlockPRange,new_cols::BlockPRange) vals = map(CartesianIndices(blocksize(a))) do I - local_views(blocks(a)[I],blocks(new_rows)[I],blocks(new_cols)[I]) + local_views(blocks(a)[I],blocks(new_rows)[I[1]],blocks(new_cols)[I[2]]) end |> to_parray_of_arrays return map(mortar,vals) end diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl index 41e3f6dc..e0adeb7d 100644 --- a/src/BlockPartitionedArrays.jl +++ b/src/BlockPartitionedArrays.jl @@ -94,7 +94,7 @@ function BlockPMatrix{V}(::UndefInitializer,rows::BlockPRange,cols::BlockPRange) c = block_cols[I[2]] PSparseMatrix{V}(undef,partition(r),partition(c)) end - return BlockPMatrix(vals,rows) + return BlockPMatrix(vals,rows,cols) end # AbstractArray API diff --git a/test/BlockPartitionedArraysTests.jl b/test/BlockPartitionedArraysTests.jl index 3e810085..e55aa6cf 100644 --- a/test/BlockPartitionedArraysTests.jl +++ b/test/BlockPartitionedArraysTests.jl @@ -1,112 +1,71 @@ +module BlockPartitionedArraysTests + +using Test, LinearAlgebra, BlockArrays, SparseArrays -using Test using Gridap -using PartitionedArrays using GridapDistributed -using BlockArrays -using SparseArrays -using LinearAlgebra - +using PartitionedArrays using GridapDistributed: BlockPArray, BlockPVector, BlockPMatrix, BlockPRange -ranks = with_debug() do distribute - distribute(LinearIndices((2,))) +function main(distribute, parts) + ranks = distribute(LinearIndices((prod(parts),))) + + n_global = 4*prod(parts) + indices = uniform_partition(ranks, n_global, true) + + pr = PRange(indices) + block_range = BlockPRange([pr, pr]) + + # BlockPRange API + @test blocklength(block_range) == 2 + @test blocksize(block_range) == (2,) + + # BlockPVector: constructor, fill!, similar, arithmetic + v = BlockPVector{Vector{Float64}}(undef, block_range) + fill!(v, 1.0) + v2 = similar(v, block_range) + fill!(v2, 2.0) + + v2 = v .+ 1.0 + v2 = v .- 1.0 + v2 = v .* 2.0 + v2 = v ./ 2.0 + copy!(v2, v) + + @test norm(v) ≈ sqrt(2*n_global) + + # PartitionedArrays vector API + consistent!(v) |> wait + assemble!(v) |> wait + partition(v) + local_values(v) + own_values(v) + ghost_values(v) + + # local_views on vector with and without new range + local_views(v) + local_views(v, block_range) + + # BlockPMatrix: constructor (exercises BUG3 fix: cols must not be dropped) + m = BlockPMatrix{SparseMatrixCSC{Float64,Int64}}(undef, block_range, block_range) + @test blocksize(m) == (2,2) + + m2 = similar(m, (block_range, block_range)) + @test blocksize(m2) == (2,2) + + # PartitionedArrays matrix API + assemble!(m) |> wait + LinearAlgebra.fillstored!(m, 0.0) + partition(m) + own_ghost_values(m) + ghost_own_values(m) + + # local_views on matrix (exercises BUG5 fix: correct I[1]/I[2] indexing) + local_views(m) + local_views(m, block_range, block_range) + + maximum(abs, v) + minimum(abs, v) end -indices = map(ranks) do r - if r == 1 - own_gids = [1,2,3,4,5] - ghost_gids = [6,7] - ghost_owners = [2,2] - else - own_gids = [6,7,8,9,10] - ghost_gids = [5] - ghost_owners = [1] - end - own_idx = OwnIndices(10,r,own_gids) - ghost_idx = GhostIndices(10,ghost_gids,ghost_owners) - OwnAndGhostIndices(own_idx,ghost_idx) -end - -block_range = BlockPRange([PRange(indices),PRange(indices)]) - -_v = PVector{OwnAndGhostVectors{Vector{Float64}}}(undef,indices) -v = BlockPArray([_v,_v],(block_range,)) -fill!(v,1.0) - -_m = map(CartesianIndices((2,2))) do I - i,j = I[1],I[2] - local_mats = map(ranks,indices) do r, ind - n = local_length(ind) - if i==j && r == 1 - SparseMatrixCSC(n,n,Int[1,3,5,7,9,10,11,13],Int[1,2,2,3,3,4,4,5,5,6,6,7],fill(1.0,12)) - elseif i==j && r == 2 - SparseMatrixCSC(n,n,Int[1,2,4,6,8,10,11],Int[1,1,2,2,3,3,4,4,5,6],fill(1.0,10)) - else - SparseMatrixCSC(n,n,fill(Int(1),n+1),Int[],Float64.([])) - end - end - PSparseMatrix(local_mats,indices,indices) -end; -m = BlockPArray(_m,(block_range,block_range)) - -x = similar(_v) -mul!(x,_m[1,1],_v) - -# BlockPRange - -@test blocklength(block_range) == 2 -@test blocksize(block_range) == (2,) - -# AbstractArray API - -__v = similar(v,block_range) -__m = similar(m,(block_range,block_range)) -fill!(__v,1.0) - -__v = __v .+ 1.0 -__v = __v .- 1.0 -__v = __v .* 1.0 -__v = __v ./ 1.0 - -# PartitionedArrays API - -consistent!(__v) |> wait -t = assemble!(__v) -assemble!(__m) |> wait -fetch(t); - -PartitionedArrays.to_trivial_partition(m) - -partition(v) -partition(m) -local_values(v) -own_values(v) -ghost_values(v) -own_ghost_values(m) -ghost_own_values(m) - -# LinearAlgebra API -fill!(v,1.0) -x = similar(v) -mul!(x,m,v) -consistent!(x) |> wait - -@test dot(v,x) ≈ 36 -norm(v) -copy!(x,v) - -LinearAlgebra.fillstored!(__m,1.0) - -__v = BlockPVector{Vector{Float64}}(undef,block_range) -#__m = BlockPMatrix{SparseMatrixCSC{Float64,Int64}}(undef,block_range,block_range) - -maximum(abs,v) -minimum(abs,v) - -# GridapDistributed API -v_parts = local_views(v) -m_parts = local_views(m) - -v_parts = local_views(v,block_range) -m_parts = local_views(m,block_range,block_range) +end # module diff --git a/test/TestApp/src/TestApp.jl b/test/TestApp/src/TestApp.jl index 56a30701..a129128d 100644 --- a/test/TestApp/src/TestApp.jl +++ b/test/TestApp/src/TestApp.jl @@ -18,6 +18,7 @@ module TestApp include("../../AdaptivityMultiFieldTests.jl") include("../../PolytopalCoarseningTests.jl") include("../../BlockSparseMatrixAssemblersTests.jl") + include("../../BlockPartitionedArraysTests.jl") include("../../VisualizationTests.jl") include("../../AutodiffTests.jl") include("../../ConstantFESpacesTests.jl") diff --git a/test/mpi/runtests_np4_body.jl b/test/mpi/runtests_np4_body.jl index eca6f4c9..c46bea12 100644 --- a/test/mpi/runtests_np4_body.jl +++ b/test/mpi/runtests_np4_body.jl @@ -53,6 +53,7 @@ function all_tests(distribute,parts) end TestApp.BlockSparseMatrixAssemblersTests.main(distribute,parts) + TestApp.BlockPartitionedArraysTests.main(distribute,parts) PArrays.toc!(t,"BlockSparseMatrixAssemblers") if prod(parts) == 4 diff --git a/test/sequential/BlockSparseMatrixAssemblersTests.jl b/test/sequential/BlockSparseMatrixAssemblersTests.jl index 9c60126d..b74fd0b8 100644 --- a/test/sequential/BlockSparseMatrixAssemblersTests.jl +++ b/test/sequential/BlockSparseMatrixAssemblersTests.jl @@ -1,17 +1,21 @@ module BlockSparseMatrixAssemblersTestsSeq using PartitionedArrays include("../BlockSparseMatrixAssemblersTests.jl") +include("../BlockPartitionedArraysTests.jl") with_debug() do distribute BlockSparseMatrixAssemblersTests.main(distribute,(2,2)) + BlockPartitionedArraysTests.main(distribute,(2,2)) end with_debug() do distribute BlockSparseMatrixAssemblersTests.main(distribute,(2,1)) + BlockPartitionedArraysTests.main(distribute,(2,1)) end with_debug() do distribute BlockSparseMatrixAssemblersTests.main(distribute,(1,2)) + BlockPartitionedArraysTests.main(distribute,(1,2)) end end # module \ No newline at end of file diff --git a/test/sequential/runtests.jl b/test/sequential/runtests.jl index 7dabeb8d..981cef0f 100644 --- a/test/sequential/runtests.jl +++ b/test/sequential/runtests.jl @@ -38,8 +38,8 @@ end @time @testset "StokesHdivDGTests.jl" begin include("StokesHdivDGTests.jl") end -@time @testset "BlockSparseMatrixAssemblers" begin - include("BlockSparseMatrixAssemblersTests.jl") +@time @testset "BlockSparseMatrixAssemblers" begin + include("BlockSparseMatrixAssemblersTests.jl") end @time @testset "AdaptivityTests" begin From 5366d9afb214e4b6111502ca749995a4b2899f02 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 11 Mar 2026 22:40:10 +1100 Subject: [PATCH 2/3] UPdate news --- NEWS.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index cacf25a5..7e401ec8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,8 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed -- Fixed `BlockPMatrix{V}(::UndefInitializer, rows, cols)` constructor dropping the `cols` argument, causing a `MethodError` at runtime. -- Fixed `local_views(::BlockPMatrix, rows, cols)` indexing 1D block-range vectors with a 2D `CartesianIndex`, causing `BoundsError` for any multi-field problem with ≥2 fields. +- Fixed `BlockPMatrix{V}(::UndefInitializer, rows, cols)` constructor dropping the `cols` argument, causing a `MethodError` at runtime. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194). +- Fixed `local_views(::BlockPMatrix, rows, cols)` indexing 1D block-range vectors with a 2D `CartesianIndex`, causing `BoundsError` for any multi-field problem with ≥2 fields. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194). ## [0.4.11] - 2026-02-20 From f76fc84df871e808c8a9fb0122e3d9413ae6a497 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 11 Mar 2026 23:31:06 +1100 Subject: [PATCH 3/3] Fixed issue in mul! for block arrays --- NEWS.md | 1 + src/BlockPartitionedArrays.jl | 12 ++++++++---- test/BlockPartitionedArraysTests.jl | 9 +++++++++ test/BlockSparseMatrixAssemblersTests.jl | 2 +- 4 files changed, 19 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7e401ec8..6c979727 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fixed `BlockPMatrix{V}(::UndefInitializer, rows, cols)` constructor dropping the `cols` argument, causing a `MethodError` at runtime. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194). - Fixed `local_views(::BlockPMatrix, rows, cols)` indexing 1D block-range vectors with a 2D `CartesianIndex`, causing `BoundsError` for any multi-field problem with ≥2 fields. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194). +- Fixed `mul!(y::BlockPVector, A::BlockPMatrix, x::BlockPVector, α, β)` computing `α*β*(A*x)` instead of `α*(A*x) + β*y`; the 3-arg `mul!` was also updated to correctly zero `y` before accumulating. Since PR[#194](https://github.com/gridap/GridapDistributed.jl/pull/194). ## [0.4.11] - 2026-02-20 diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl index e0adeb7d..eb918d63 100644 --- a/src/BlockPartitionedArrays.jl +++ b/src/BlockPartitionedArrays.jl @@ -329,19 +329,23 @@ end function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector) o = one(eltype(A)) - mul!(y,A,x,o,o) + z = zero(eltype(A)) + mul!(y,A,x,o,z) end function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector,α::Number,β::Number) yb, Ab, xb = blocks(y), blocks(A), blocks(x) - z = zero(eltype(y)) o = one(eltype(A)) + z = zero(eltype(y)) for i in 1:blocksize(A,1) - fill!(yb[i],z) + if iszero(β) + fill!(yb[i],z) + else + rmul!(yb[i],β) + end for j in 1:blocksize(A,2) mul!(yb[i],Ab[i,j],xb[j],α,o) end - rmul!(yb[i],β) end return y end diff --git a/test/BlockPartitionedArraysTests.jl b/test/BlockPartitionedArraysTests.jl index e55aa6cf..ec7bd60e 100644 --- a/test/BlockPartitionedArraysTests.jl +++ b/test/BlockPartitionedArraysTests.jl @@ -66,6 +66,15 @@ function main(distribute, parts) maximum(abs, v) minimum(abs, v) + + # BUG4 fix: mul!(y,A,x,α,β) must compute α*(A*x) + β*y, not α*β*(A*x) + # With A=0: y ← α*(0*x) + β*y_old = β*y_old + fill!(v, 1.0) + v4 = similar(v, block_range) + fill!(v4, 3.0) + LinearAlgebra.fillstored!(m, 0.0) + mul!(v4, m, v, 2.0, 2.0) # y ← 2*(0*x) + 2*[3,...] = [6,...] + @test norm(v4) ≈ 6.0 * sqrt(2*n_global) end end # module diff --git a/test/BlockSparseMatrixAssemblersTests.jl b/test/BlockSparseMatrixAssemblersTests.jl index ecf808e4..13eb3d7e 100644 --- a/test/BlockSparseMatrixAssemblersTests.jl +++ b/test/BlockSparseMatrixAssemblersTests.jl @@ -21,7 +21,7 @@ function is_same_vector(x::BlockPVector,y::PVector,Ub,U) res = map(1:num_fields(Ub)) do i xi = restrict_to_field(Ub,x_fespace,i) yi = restrict_to_field(U,y_fespace,i) - xi ≈ yi + norm(xi - yi) < 1.0e-10 end return all(res) end