From b11c7768035f98e62acc74d6c037ce9ff800eebb Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 7 Aug 2024 07:29:29 +0200 Subject: [PATCH 01/14] Improvements in renumber for PVector --- src/p_vector.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/p_vector.jl b/src/p_vector.jl index ca7a7fac..072d82c8 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -1482,13 +1482,14 @@ function renumber(a::PVector;kwargs...) renumber(a,row_partition_2;kwargs...) end -function renumber(a::PVector,row_partition_2;renumber_local_indices=true) - # TODO storing the vector in split format would help to return a view - # of the input - if renumber_local_indices - values = map(vcat,own_values(a),ghost_values(a)) +function renumber(a::PVector,row_partition_2;renumber_local_indices=Val(true)) + if val_parameter(renumber_local_indices) + perms = map(row_partition_2) do myrows + Int32(1):Int32(local_length(myrows)) + end + values = map(split_vector,own_values(a),ghost_values(a),perms) else - values = map(copy,local_values(a)) + values = local_values(a) end PVector(values,row_partition_2) end From 0b07602a8e9db5dec23cd3d3bc737aafb479e97a Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 7 Aug 2024 07:31:40 +0200 Subject: [PATCH 02/14] Exported local_permutation --- CHANGELOG.md | 1 + src/PartitionedArrays.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fcb6aa41..9a770b65 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Split format support for `PVector`. - Helper functions to build partitioned sparse matrices and vectors in split format, `pvector_from_split_blocks` and `psparse_from_split_blocks`. +- `export local_permutation` ### Deprecated diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index bec2f615..e0d082ab 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -94,6 +94,7 @@ export own_to_local export ghost_to_local export local_to_own export local_to_ghost +export local_permutation export replace_ghost export remove_ghost export union_ghost From bf500f670ee5b483ff109e87937a5ab5117974ef Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 7 Aug 2024 10:59:22 +0200 Subject: [PATCH 03/14] Added array_of_tuples --- CHANGELOG.md | 3 ++- docs/src/reference/arraymethods.md | 1 + src/PartitionedArrays.jl | 1 + src/primitives.jl | 9 +++++++++ test/primitives_tests.jl | 4 ++++ 5 files changed, 17 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9a770b65..6118ccb8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Split format support for `PVector`. - Helper functions to build partitioned sparse matrices and vectors in split format, `pvector_from_split_blocks` and `psparse_from_split_blocks`. -- `export local_permutation` +- Function `array_of_tuples`. +- Export statement for `local_permutation`. ### Deprecated diff --git a/docs/src/reference/arraymethods.md b/docs/src/reference/arraymethods.md index ea830d7d..2c6f02d5 100644 --- a/docs/src/reference/arraymethods.md +++ b/docs/src/reference/arraymethods.md @@ -14,4 +14,5 @@ MAIN i_am_main map_main tuple_of_arrays +array_of_tuples ``` diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index e0d082ab..ed02b884 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -27,6 +27,7 @@ include("sparse_utils.jl") export linear_indices export cartesian_indices export tuple_of_arrays +export array_of_tuples export i_am_main export MAIN export map_main diff --git a/src/primitives.jl b/src/primitives.jl index 0c32f620..37fb99ab 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -96,6 +96,15 @@ function tuple_of_arrays(a) take(a,eltype(a)) end +""" + array_of_tuples(a) +""" +function array_of_tuples(a) + map(a...) do b... + b + end +end + # We don't need a real task since MPI already is able to do # fake_asynchronous (nonblocking) operations. # We want to avoid heap allocations of standard tasks. diff --git a/test/primitives_tests.jl b/test/primitives_tests.jl index 20980bd7..79ffa364 100644 --- a/test/primitives_tests.jl +++ b/test/primitives_tests.jl @@ -15,6 +15,10 @@ function primitives_tests(distribute) end a, b = tuple_of_arrays(a_and_b) + a_and_b_2 = array_of_tuples((a,b)) + map(a_and_b,a_and_b) do v1,v2 + @test v1 == v2 + end map(a,b,rank) do a,b,rank @test a == 2*rank From cb7291c20dbed20b34cda87cfb76cbc91297697f Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 7 Aug 2024 11:15:44 +0200 Subject: [PATCH 04/14] Added BlockedPRange --- Project.toml | 1 + src/PartitionedArrays.jl | 5 ++++ src/p_range.jl | 54 ++++++++++++++++++++++++++++++++++++++++ test/p_range_tests.jl | 15 +++++++++++ 4 files changed, 75 insertions(+) diff --git a/Project.toml b/Project.toml index 98fad2bf..912766c7 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Francesc Verdugo and contributors"] version = "0.5.1" [deps] +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" CircularArrays = "7a955b69-7140-5f4e-a0ed-f168c5e2e749" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index ed02b884..561f9c80 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -8,6 +8,7 @@ using CircularArrays import MPI import IterativeSolvers import Distances +using BlockArrays export length_to_ptrs! export rewind_ptrs! @@ -113,6 +114,10 @@ export map_ghost_to_global! export map_global_to_ghost! export map_own_to_global! export map_global_to_own! +export BlockedPRange +export local_block_ranges +export own_block_ranges +export ghost_block_ranges include("p_range.jl") export local_values diff --git a/src/p_range.jl b/src/p_range.jl index 4b4794a4..f2baa68f 100644 --- a/src/p_range.jl +++ b/src/p_range.jl @@ -1797,3 +1797,57 @@ ghost_to_local(pr::PRange) = map(ghost_to_local,partition(pr)) local_to_own(pr::PRange) = map(local_to_own,partition(pr)) local_to_ghost(pr::PRange) = map(local_to_ghost,partition(pr)) + +struct BlockedPRange{A} <: AbstractUnitRange{Int} + blocks::A + blocklasts::Vector{Int} + function BlockedPRange(blocks) + @assert all(block->isa(block,PRange),blocks) + nblocks = length(blocks) + blocklasts = zeros(Int,nblocks) + prev = 0 + for i in 1:nblocks + prev += length(blocks[i]) + blocklasts[i] = prev + end + A = typeof(blocks) + new{A}(blocks,blocklasts) + end +end + +BlockArrays.blocks(a::BlockedPRange) = a.blocks +BlockArrays.blocklasts(a::BlockedPRange) = a.blocklasts +BlockArrays.blockaxes(a::BlockedPRange) = (Block.(Base.OneTo(length(blocks(a)))),) +Base.getindex(a::BlockedPRange,k::Block{1}) = blocks(a)[first(k.n)] +function BlockArrays.findblock(b::BlockedPRange, k::Integer) + @boundscheck k in b || throw(BoundsError(b,k)) + Block(searchsortedfirst(blocklasts(b), k)) +end +Base.first(a::BlockedPRange) = 1 +Base.last(a::BlockedPRange) = sum(length,blocks(a)) + +function PartitionedArrays.partition(a::BlockedPRange) + local_block_ranges(a) +end + +function local_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(local_length,myids) |> blockedrange + end +end + +function own_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(own_length,myids) |> blockedrange + end +end + +function ghost_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(ghost_length,myids) |> blockedrange + end +end + diff --git a/test/p_range_tests.jl b/test/p_range_tests.jl index 05dd1d32..d59f3635 100644 --- a/test/p_range_tests.jl +++ b/test/p_range_tests.jl @@ -1,6 +1,7 @@ using PartitionedArrays using Test +using BlockArrays using PartitionedArrays: local_range @@ -295,4 +296,18 @@ function p_range_tests(distribute) @test ids1 == ids2 end + ids1 = uniform_partition(parts,(2,2),(5,4)) + ids2 = uniform_partition(parts,(2,2),(6,3)) + ax = BlockedPRange(PRange.([ids1,ids2])) + blockaxes(ax,1) + @test isa(ax[Block(1)],PRange) + @test isa(ax[Block(2)],PRange) + @test blocklasts(ax) == [20, 38] + @test findblock(ax,4) == Block(1) + @test findblock(ax,25) == Block(2) + partition(ax) + local_block_ranges(ax) + own_block_ranges(ax) + ghost_block_ranges(ax) + end From 57bfa5f359736421180d0a99d368c271183f17f8 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 7 Aug 2024 12:07:31 +0200 Subject: [PATCH 05/14] Added initial implementation of BlockPVector --- src/PartitionedArrays.jl | 1 + src/p_range.jl | 6 +-- src/p_vector.jl | 79 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 3 deletions(-) diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index 561f9c80..5542289c 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -145,6 +145,7 @@ export SplitVector export split_vector export split_vector_blocks export pvector_from_split_blocks +export BlockPVector include("p_vector.jl") export SplitMatrix diff --git a/src/p_range.jl b/src/p_range.jl index f2baa68f..5c326045 100644 --- a/src/p_range.jl +++ b/src/p_range.jl @@ -1833,21 +1833,21 @@ end function local_block_ranges(a::BlockedPRange) ids = map(partition,blocks(a)) |> array_of_tuples map(ids) do myids - map(local_length,myids) |> blockedrange + map(local_length,myids) |> collect |> blockedrange end end function own_block_ranges(a::BlockedPRange) ids = map(partition,blocks(a)) |> array_of_tuples map(ids) do myids - map(own_length,myids) |> blockedrange + map(own_length,myids) |> collect |> blockedrange end end function ghost_block_ranges(a::BlockedPRange) ids = map(partition,blocks(a)) |> array_of_tuples map(ids) do myids - map(ghost_length,myids) |> blockedrange + map(ghost_length,myids) |> collect |> blockedrange end end diff --git a/src/p_vector.jl b/src/p_vector.jl index 072d82c8..c2757afe 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -1494,3 +1494,82 @@ function renumber(a::PVector,row_partition_2;renumber_local_indices=Val(true)) PVector(values,row_partition_2) end +struct BlockPVector{A,T} <: BlockArrays.AbstractBlockVector{T} + blocks::A + function BlockPVector(blocks) + T = eltype(first(blocks)) + @assert all(block->isa(block,PVector),blocks) + @assert all(block->eltype(block)==T,blocks) + A = typeof(blocks) + new{A,T}(blocks) + end +end +BlockArrays.blocks(a::BlockPVector) = a.blocks +BlockArrays.viewblock(a::BlockPVector, i::Block) = blocks(a)[i.n...] +Base.axes(a::BlockPVector) = (BlockedPRange(map(block->axes(block,1),blocks(a))),) +Base.length(a::BlockPVector) = sum(length,blocks(a)) +Base.IndexStyle(::Type{<:BlockPVector}) = IndexLinear() +function Base.getindex(a::BlockPVector,gid::Int) + scalar_indexing_action(a) +end +function Base.setindex(a::BlockPVector,v,gid::Int) + scalar_indexing_action(a) +end +function Base.show(io::IO,k::MIME"text/plain",data::BlockPVector) + T = eltype(data) + n = length(data) + ps = partition(axes(data,1)) + np = length(ps) + nb = blocklength(data) + map_main(ps) do _ + println(io,"$nb-blocked $n-element BlockPVector with eltype $T partitioned into $np parts.") + end +end +function Base.show(io::IO,data::BlockPVector) + print(io,"BlockPVector(…)") +end + +function partition(a::BlockPVector) + local_values(a) +end + +function local_values(a::BlockPVector) + r = local_block_ranges(axes(a,1)) + v = map(local_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function own_values(a::BlockPVector) + r = own_block_ranges(axes(a,1)) + v = map(own_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function ghost_values(a::BlockPVector) + r = ghost_block_ranges(axes(a,1)) + v = map(ghost_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function consistent!(a::BlockPVector) + ts = map(consistent!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + +function assemble!(a::BlockPVector) + ts = map(assemble!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + From 9d3863d33a2b611f603b0e60e5db17490ca3accc Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 14 Aug 2024 15:06:39 +0200 Subject: [PATCH 06/14] More code related with BlockPVector --- src/p_vector.jl | 269 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 229 insertions(+), 40 deletions(-) diff --git a/src/p_vector.jl b/src/p_vector.jl index c2757afe..bf3a850c 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -758,6 +758,10 @@ function PVector(::UndefInitializer,index_partition) PVector{Vector{Float64}}(undef,index_partition) end +function PVector(::UndefInitializer,r::PRange) + PVector(undef,partition(r)) +end + """ PVector{V}(undef,index_partition) PVector(undef,index_partition) @@ -772,6 +776,10 @@ function PVector{V}(::UndefInitializer,index_partition) where V PVector(vector_partition,index_partition) end +function PVector{V}(::UndefInitializer,r::PRange) where V + PVector{V}(undef,partition(r)) +end + function Base.copy!(a::PVector,b::PVector) @assert length(a) == length(b) copyto!(a,b) @@ -1180,15 +1188,15 @@ function LinearAlgebra.norm(a::PVector,p::Real=2) reduce(+,contibs;init=zero(eltype(contibs)))^(1/p) end -struct PBroadcasted{A,B,C} +struct BroadcastedPVector{A,B,C} own_values::A ghost_values::B index_partition::C end -own_values(a::PBroadcasted) = a.own_values -ghost_values(a::PBroadcasted) = a.ghost_values +own_values(a::BroadcastedPVector) = a.own_values +ghost_values(a::BroadcastedPVector) = a.ghost_values -function Base.broadcasted(f, args::Union{PVector,PBroadcasted}...) +function Base.broadcasted(f, args::Union{PVector,BroadcastedPVector}...) a1 = first(args) @boundscheck @assert all(ai->matching_own_indices(PRange(ai.index_partition),PRange(a1.index_partition)),args) own_values_in = map(own_values,args) @@ -1199,31 +1207,31 @@ function Base.broadcasted(f, args::Union{PVector,PBroadcasted}...) else ghost_values_out = nothing end - PBroadcasted(own_values_out,ghost_values_out,a1.index_partition) + BroadcastedPVector(own_values_out,ghost_values_out,a1.index_partition) end -function Base.broadcasted( f, a::Number, b::Union{PVector,PBroadcasted}) +function Base.broadcasted( f, a::Number, b::Union{PVector,BroadcastedPVector}) own_values_out = map(b->Base.broadcasted(f,a,b),own_values(b)) if ghost_values(b) !== nothing ghost_values_out = map(b->Base.broadcasted(f,a,b),ghost_values(b)) else ghost_values_out = nothing end - PBroadcasted(own_values_out,ghost_values_out,b.index_partition) + BroadcastedPVector(own_values_out,ghost_values_out,b.index_partition) end -function Base.broadcasted( f, a::Union{PVector,PBroadcasted}, b::Number) +function Base.broadcasted( f, a::Union{PVector,BroadcastedPVector}, b::Number) own_values_out = map(a->Base.broadcasted(f,a,b),own_values(a)) if ghost_values(a) !== nothing ghost_values_out = map(a->Base.broadcasted(f,a,b),ghost_values(a)) else ghost_values_out = nothing end - PBroadcasted(own_values_out,ghost_values_out,a.index_partition) + BroadcastedPVector(own_values_out,ghost_values_out,a.index_partition) end function Base.broadcasted(f, - a::Union{PVector,PBroadcasted}, + a::Union{PVector,BroadcastedPVector}, b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) Base.broadcasted(f,a,Base.materialize(b)) end @@ -1231,11 +1239,11 @@ end function Base.broadcasted( f, a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, - b::Union{PVector,PBroadcasted}) + b::Union{PVector,BroadcastedPVector}) Base.broadcasted(f,Base.materialize(a),b) end -function Base.materialize(b::PBroadcasted) +function Base.materialize(b::BroadcastedPVector) own_values_out = map(Base.materialize,b.own_values) T = eltype(eltype(own_values_out)) a = PVector{Vector{T}}(undef,b.index_partition) @@ -1243,7 +1251,7 @@ function Base.materialize(b::PBroadcasted) a end -function Base.materialize!(a::PVector,b::PBroadcasted) +function Base.materialize!(a::PVector,b::BroadcastedPVector) map(Base.materialize!,own_values(a),own_values(b)) if b.ghost_values !== nothing && a.index_partition === b.index_partition map(Base.materialize!,ghost_values(a),ghost_values(b)) @@ -1254,38 +1262,43 @@ end for M in Distances.metrics @eval begin function (d::$M)(a::PVector,b::PVector) - if Distances.parameters(d) !== nothing - error("Only distances without parameters are implemented at this moment") - end - partials = map(own_values(a),own_values(b)) do a,b - @boundscheck if length(a) != length(b) - throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b)).")) - end - if length(a) == 0 - return zero(Distances.result_type(d, a, b)) + s = distance_eval_body(d,a,b) + Distances.eval_end(d,s) + end + end +end + +function distance_eval_body(d,a::PVector,b::PVector) + if Distances.parameters(d) !== nothing + error("Only distances without parameters are implemented at this moment") + end + partials = map(own_values(a),own_values(b)) do a,b + @boundscheck if length(a) != length(b) + throw(DimensionMismatch("first array has length $(length(a)) which does not match the length of the second, $(length(b)).")) + end + if length(a) == 0 + return zero(Distances.result_type(d, a, b)) + end + @inbounds begin + s = Distances.eval_start(d, a, b) + if (IndexStyle(a, b) === IndexLinear() && eachindex(a) == eachindex(b)) || axes(a) == axes(b) + for I in eachindex(a, b) + ai = a[I] + bi = b[I] + s = Distances.eval_reduce(d, s, Distances.eval_op(d, ai, bi)) end - @inbounds begin - s = Distances.eval_start(d, a, b) - if (IndexStyle(a, b) === IndexLinear() && eachindex(a) == eachindex(b)) || axes(a) == axes(b) - @simd for I in eachindex(a, b) - ai = a[I] - bi = b[I] - s = Distances.eval_reduce(d, s, Distances.eval_op(d, ai, bi)) - end - else - for (ai, bi) in zip(a, b) - s = Distances.eval_reduce(d, s, Distances.eval_op(d, ai, bi)) - end - end - return s + else + for (ai, bi) in zip(a, b) + s = Distances.eval_reduce(d, s, Distances.eval_op(d, ai, bi)) end end - s = reduce((i,j)->Distances.eval_reduce(d,i,j), - partials, - init=Distances.eval_start(d, a, b)) - Distances.eval_end(d,s) + return s end end + s = reduce((i,j)->Distances.eval_reduce(d,i,j), + partials, + init=Distances.eval_start(d, a, b)) + s end # New stuff @@ -1494,6 +1507,8 @@ function renumber(a::PVector,row_partition_2;renumber_local_indices=Val(true)) PVector(values,row_partition_2) end +# BlockPVector + struct BlockPVector{A,T} <: BlockArrays.AbstractBlockVector{T} blocks::A function BlockPVector(blocks) @@ -1573,3 +1588,177 @@ function assemble!(a::BlockPVector) end end +function Base.similar(a::BlockPVector,::Type{T},inds::Tuple{<:BlockedPRange}) where T + r = first(inds) + bs = map((ai,ri)->similar(ai,T,ri),blocks(a),blocks(r)) + BlockPVector(bs) +end + +function Base.similar(::Type{<:BlockPVector{A}},inds::Tuple{<:BlockedPRange}) where A + V = eltype(A) + r = first(inds) + bs = map(ri->similar(V,ri),blocks(r)) + BlockPVector(bs) +end + +function BlockPVector(::UndefInitializer,r::BlockedPRange) + bs = map(ri->PVector(undef,ri),blocks(r)) + BlockPVector(bs) +end + +function BlockPVector{A}(::UndefInitializer,r::BlockedPRange) where A + similar(BlockPVector{A},(r,)) +end + +function Base.copy!(a::BlockPVector,b::BlockPVector) + foreach(copy!,blocks(a),blocks(b)) + a +end + +function Base.copyto!(a::BlockPVector,b::BlockPVector) + foreach(copyto!,blocks(a),blocks(b)) + a +end + +function Base.fill!(a::BlockPVector,v) + foreach(ai->fill!(ai,v),blocks(a)) +end + + +function Base.:(==)(a::BlockPVector,b::BlockPVector) + all(map(==,blocks(a),blocks(b))) +end + +function Base.any(f::Function,x::BlockPVector) + any(xi->any(f,xi),blocks(x)) +end + +function Base.all(f::Function,x::BlockPVector) + all(xi->all(f,xi),blocks(x)) +end + +Base.maximum(x::BlockPVector) = maximum(identity,x) +function Base.maximum(f::Function,x::BlockPVector) + maximum(xi->maximum(f,xi),blocks(x)) +end + +Base.minimum(x::BlockPVector) = minimum(identity,x) +function Base.minimum(f::Function,x::BlockPVector) + minimum(xi->minimum(f,xi),blocks(x)) +end + +function Base.collect(v::BlockPVector) + reduce(vcat,map(collect,blocks(v))) +end + +function Base.:*(a::Number,b::BlockPVector) + bs = map(bi->a*bi,blocks(b)) + BlockPVector(bs) +end + +function Base.:*(b::BlockPVector,a::Number) + a*b +end + +function Base.:/(b::BlockPVector,a::Number) + (1/a)*b +end + +for op in (:+,:-) + @eval begin + function Base.$op(a::BlockPVector) + bs = map(ai->$op(ai),blocks(a)) + BlockPVector(bs) + end + function Base.$op(a::BlockPVector,b::BlockPVector) + $op.(a,b) + end + end +end + +function Base.reduce(op,a::BlockPVector;neutral=neutral_element(op,eltype(a)),kwargs...) + rs = map(ai->reduce(op,ai;neutral,kwargs...),blocks(a)) + reduce(rs;kwargs...) +end + +function Base.sum(a::BlockPVector) + reduce(+,a,init=zero(eltype(a))) +end + +function LinearAlgebra.dot(a::BlockPVector,b::BlockPVector) + c = map(dot,blocks(a),blocks(b)) + sum(c) +end + +function LinearAlgebra.rmul!(a::BlockPVector,v::Number) + map(blocks(a)) do l + rmul!(l,v) + end + a +end + +function LinearAlgebra.norm(a::BlockPVector,p::Real=2) + contibs = map(blocks(a)) do oid_to_value + norm(oid_to_value,p)^p + end + reduce(+,contibs;init=zero(eltype(contibs)))^(1/p) +end + +for M in Distances.metrics + @eval begin + function (d::$M)(a::BlockPVector,b::BlockPVector) + s = distance_eval_body(d,a,b) + Distances.eval_end(d,s) + end + end +end + +function distance_eval_body(d,a::BlockPVector,b::BlockPVector) + partials = map(blocks(a),blocks(b)) do (ai,bi) + distance_eval_body(d,ai,bi) + end + s = reduce((i,j)->Distances.eval_reduce(d,i,j), + partials, + init=Distances.eval_start(d, a, b)) + s +end + +struct BroadcastedBlockPVector{A} + blocks::A +end +BlockArrays.blocks(a::BroadcastedBlockPVector) = a.blocks + +function Base.broadcasted(f, args::Union{BlockPVector,BroadcastedBlockPVector}...) + map( (bs...) -> Base.broadcasted(f,bs...) , map(blocks,args)...) |> BroadcastedBlockPVector +end + +function Base.broadcasted( f, a::Number, b::Union{BlockPVector,BroadcastedBlockPVector}) + map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBlockPVector +end + +function Base.broadcasted( f, a::Union{BlockPVector,BroadcastedBlockPVector}, b::Number) + map( ai -> Base.broadcasted(f,ai,b), blocks(a)) |> BroadcastedBlockPVector +end + +function Base.broadcasted(f, + a::Union{BlockPVector,BroadcastedBlockPVector}, + b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) + Base.broadcasted(f,a,Base.materialize(b)) +end + +function Base.broadcasted( + f, + a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, + b::Union{BlockPVector,BroadcastedBlockPVector}) + Base.broadcasted(f,Base.materialize(a),b) + end + +function Base.materialize(b::BroadcastedBlockPVector) + map(Base.materialize,blocks(b)) |> BlockPVector +end + +function Base.materialize!(a::BlockPVector,b::BroadcastedBlockPVector) + foreach(Base.materialize!,blocks(a),blocks(b)) + a +end + From 82a8f5984b8d60828dd5c9c961dfcda8b437c293 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 14 Aug 2024 15:44:19 +0200 Subject: [PATCH 07/14] Moving code related to block arrays to its own file. --- src/PartitionedArrays.jl | 12 +- src/block_arrays.jl | 316 +++++++++++++++++++++++++ src/p_range.jl | 54 ----- src/p_vector.jl | 255 -------------------- test/block_arrays_tests.jl | 73 ++++++ test/debug_array/block_arrays_tests.jl | 8 + test/debug_array/runtests.jl | 2 + test/mpi_array/block_arrays_tests.jl | 14 ++ test/mpi_array/runtests.jl | 1 + 9 files changed, 421 insertions(+), 314 deletions(-) create mode 100644 src/block_arrays.jl create mode 100644 test/block_arrays_tests.jl create mode 100644 test/debug_array/block_arrays_tests.jl create mode 100644 test/mpi_array/block_arrays_tests.jl diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index b521f8d8..df472960 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -116,10 +116,6 @@ export map_ghost_to_global! export map_global_to_ghost! export map_own_to_global! export map_global_to_own! -export BlockedPRange -export local_block_ranges -export own_block_ranges -export ghost_block_ranges include("p_range.jl") export local_values @@ -147,7 +143,6 @@ export SplitVector export split_vector export split_vector_blocks export pvector_from_split_blocks -export BlockPVector include("p_vector.jl") export SplitMatrix @@ -179,6 +174,13 @@ export spmtm! export centralize include("p_sparse_matrix.jl") +export BlockedPRange +export local_block_ranges +export own_block_ranges +export ghost_block_ranges +export BlockPVector +include("block_arrays.jl") + export PTimer export tic! export toc! diff --git a/src/block_arrays.jl b/src/block_arrays.jl new file mode 100644 index 00000000..e0db22f4 --- /dev/null +++ b/src/block_arrays.jl @@ -0,0 +1,316 @@ + +function BlockArrays.mortar(blocks::Vector{<:PRange}) + BlockedPRange(blocks) +end + +struct BlockedPRange{A} <: AbstractUnitRange{Int} + blocks::A + blocklasts::Vector{Int} + function BlockedPRange(blocks) + @assert all(block->isa(block,PRange),blocks) + nblocks = length(blocks) + blocklasts = zeros(Int,nblocks) + prev = 0 + for i in 1:nblocks + prev += length(blocks[i]) + blocklasts[i] = prev + end + A = typeof(blocks) + new{A}(blocks,blocklasts) + end +end +BlockArrays.blocks(a::BlockedPRange) = a.blocks +BlockArrays.blocklasts(a::BlockedPRange) = a.blocklasts +BlockArrays.blockaxes(a::BlockedPRange) = (Block.(Base.OneTo(length(blocks(a)))),) +Base.getindex(a::BlockedPRange,k::Block{1}) = blocks(a)[first(k.n)] +function BlockArrays.findblock(b::BlockedPRange, k::Integer) + @boundscheck k in b || throw(BoundsError(b,k)) + Block(searchsortedfirst(blocklasts(b), k)) +end +Base.first(a::BlockedPRange) = 1 +Base.last(a::BlockedPRange) = sum(length,blocks(a)) + +function PartitionedArrays.partition(a::BlockedPRange) + local_block_ranges(a) +end + +function local_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(local_length,myids) |> collect |> blockedrange + end +end + +function own_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(own_length,myids) |> collect |> blockedrange + end +end + +function ghost_block_ranges(a::BlockedPRange) + ids = map(partition,blocks(a)) |> array_of_tuples + map(ids) do myids + map(ghost_length,myids) |> collect |> blockedrange + end +end + +# BlockPVector + +function BlockArrays.mortar(blocks::Vector{<:PVector}) + BlockPVector(blocks) +end + +struct BlockPVector{A,T} <: BlockArrays.AbstractBlockVector{T} + blocks::A + function BlockPVector(blocks) + T = eltype(first(blocks)) + @assert all(block->isa(block,PVector),blocks) + @assert all(block->eltype(block)==T,blocks) + A = typeof(blocks) + new{A,T}(blocks) + end +end +BlockArrays.blocks(a::BlockPVector) = a.blocks +BlockArrays.viewblock(a::BlockPVector, i::Block) = blocks(a)[i.n...] +Base.axes(a::BlockPVector) = (BlockedPRange(map(block->axes(block,1),blocks(a))),) +Base.length(a::BlockPVector) = sum(length,blocks(a)) +Base.IndexStyle(::Type{<:BlockPVector}) = IndexLinear() +function Base.getindex(a::BlockPVector,gid::Int) + scalar_indexing_action(a) +end +function Base.setindex(a::BlockPVector,v,gid::Int) + scalar_indexing_action(a) +end +function Base.show(io::IO,k::MIME"text/plain",data::BlockPVector) + T = eltype(data) + n = length(data) + ps = partition(axes(data,1)) + np = length(ps) + nb = blocklength(data) + map_main(ps) do _ + println(io,"$nb-blocked $n-element BlockPVector with eltype $T partitioned into $np parts.") + end +end +function Base.show(io::IO,data::BlockPVector) + print(io,"BlockPVector(…)") +end + +function partition(a::BlockPVector) + local_values(a) +end + +function local_values(a::BlockPVector) + r = local_block_ranges(axes(a,1)) + v = map(local_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function own_values(a::BlockPVector) + r = own_block_ranges(axes(a,1)) + v = map(own_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function ghost_values(a::BlockPVector) + r = ghost_block_ranges(axes(a,1)) + v = map(ghost_values,blocks(a)) |> array_of_tuples + map(v,r) do myv, myr + mortar(collect(myv),(myr,)) + end +end + +function consistent!(a::BlockPVector) + ts = map(consistent!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + +function assemble!(a::BlockPVector) + ts = map(assemble!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + +function Base.similar(a::BlockPVector,::Type{T},inds::Tuple{<:BlockedPRange}) where T + r = first(inds) + bs = map((ai,ri)->similar(ai,T,ri),blocks(a),blocks(r)) + BlockPVector(bs) +end + +function Base.similar(::Type{<:BlockPVector{A}},inds::Tuple{<:BlockedPRange}) where A + V = eltype(A) + r = first(inds) + bs = map(ri->similar(V,ri),blocks(r)) + BlockPVector(bs) +end + +function BlockPVector(::UndefInitializer,r::BlockedPRange) + bs = map(ri->PVector(undef,ri),blocks(r)) + BlockPVector(bs) +end + +function BlockPVector{A}(::UndefInitializer,r::BlockedPRange) where A + similar(BlockPVector{A},(r,)) +end + +function Base.copy!(a::BlockPVector,b::BlockPVector) + foreach(copy!,blocks(a),blocks(b)) + a +end + +function Base.copyto!(a::BlockPVector,b::BlockPVector) + foreach(copyto!,blocks(a),blocks(b)) + a +end + +function Base.fill!(a::BlockPVector,v) + foreach(ai->fill!(ai,v),blocks(a)) +end + + +function Base.:(==)(a::BlockPVector,b::BlockPVector) + all(map(==,blocks(a),blocks(b))) +end + +function Base.any(f::Function,x::BlockPVector) + any(xi->any(f,xi),blocks(x)) +end + +function Base.all(f::Function,x::BlockPVector) + all(xi->all(f,xi),blocks(x)) +end + +Base.maximum(x::BlockPVector) = maximum(identity,x) +function Base.maximum(f::Function,x::BlockPVector) + maximum(xi->maximum(f,xi),blocks(x)) +end + +Base.minimum(x::BlockPVector) = minimum(identity,x) +function Base.minimum(f::Function,x::BlockPVector) + minimum(xi->minimum(f,xi),blocks(x)) +end + +function Base.collect(v::BlockPVector) + reduce(vcat,map(collect,blocks(v))) +end + +function Base.:*(a::Number,b::BlockPVector) + bs = map(bi->a*bi,blocks(b)) + BlockPVector(bs) +end + +function Base.:*(b::BlockPVector,a::Number) + a*b +end + +function Base.:/(b::BlockPVector,a::Number) + (1/a)*b +end + +for op in (:+,:-) + @eval begin + function Base.$op(a::BlockPVector) + bs = map(ai->$op(ai),blocks(a)) + BlockPVector(bs) + end + function Base.$op(a::BlockPVector,b::BlockPVector) + $op.(a,b) + end + end +end + +function Base.reduce(op,a::BlockPVector;neutral=neutral_element(op,eltype(a)),kwargs...) + rs = map(ai->reduce(op,ai;neutral,kwargs...),blocks(a)) + reduce(op,rs;kwargs...) +end + +function Base.sum(a::BlockPVector) + reduce(+,a,init=zero(eltype(a))) +end + +function LinearAlgebra.dot(a::BlockPVector,b::BlockPVector) + c = map(dot,blocks(a),blocks(b)) + sum(c) +end + +function LinearAlgebra.rmul!(a::BlockPVector,v::Number) + map(blocks(a)) do l + rmul!(l,v) + end + a +end + +function LinearAlgebra.norm(a::BlockPVector,p::Real=2) + contibs = map(blocks(a)) do oid_to_value + norm(oid_to_value,p)^p + end + reduce(+,contibs;init=zero(eltype(contibs)))^(1/p) +end + +for M in Distances.metrics + @eval begin + function (d::$M)(a::BlockPVector,b::BlockPVector) + s = distance_eval_body(d,a,b) + Distances.eval_end(d,s) + end + end +end + +function distance_eval_body(d,a::BlockPVector,b::BlockPVector) + partials = map(blocks(a),blocks(b)) do ai,bi + distance_eval_body(d,ai,bi) + end + s = reduce((i,j)->Distances.eval_reduce(d,i,j), + partials, + init=Distances.eval_start(d, a, b)) + s +end + +struct BroadcastedBlockPVector{A} + blocks::A +end +BlockArrays.blocks(a::BroadcastedBlockPVector) = a.blocks + +function Base.broadcasted(f, args::Union{BlockPVector,BroadcastedBlockPVector}...) + map( (bs...) -> Base.broadcasted(f,bs...) , map(blocks,args)...) |> BroadcastedBlockPVector +end + +function Base.broadcasted( f, a::Number, b::Union{BlockPVector,BroadcastedBlockPVector}) + map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBlockPVector +end + +function Base.broadcasted( f, a::Union{BlockPVector,BroadcastedBlockPVector}, b::Number) + map( ai -> Base.broadcasted(f,ai,b), blocks(a)) |> BroadcastedBlockPVector +end + +function Base.broadcasted(f, + a::Union{BlockPVector,BroadcastedBlockPVector}, + b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) + Base.broadcasted(f,a,Base.materialize(b)) +end + +function Base.broadcasted( + f, + a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, + b::Union{BlockPVector,BroadcastedBlockPVector}) + Base.broadcasted(f,Base.materialize(a),b) + end + +function Base.materialize(b::BroadcastedBlockPVector) + map(Base.materialize,blocks(b)) |> BlockPVector +end + +function Base.materialize!(a::BlockPVector,b::BroadcastedBlockPVector) + foreach(Base.materialize!,blocks(a),blocks(b)) + a +end + diff --git a/src/p_range.jl b/src/p_range.jl index 5d626f1c..cfa330b0 100644 --- a/src/p_range.jl +++ b/src/p_range.jl @@ -1849,57 +1849,3 @@ ghost_to_local(pr::PRange) = map(ghost_to_local,partition(pr)) local_to_own(pr::PRange) = map(local_to_own,partition(pr)) local_to_ghost(pr::PRange) = map(local_to_ghost,partition(pr)) - -struct BlockedPRange{A} <: AbstractUnitRange{Int} - blocks::A - blocklasts::Vector{Int} - function BlockedPRange(blocks) - @assert all(block->isa(block,PRange),blocks) - nblocks = length(blocks) - blocklasts = zeros(Int,nblocks) - prev = 0 - for i in 1:nblocks - prev += length(blocks[i]) - blocklasts[i] = prev - end - A = typeof(blocks) - new{A}(blocks,blocklasts) - end -end - -BlockArrays.blocks(a::BlockedPRange) = a.blocks -BlockArrays.blocklasts(a::BlockedPRange) = a.blocklasts -BlockArrays.blockaxes(a::BlockedPRange) = (Block.(Base.OneTo(length(blocks(a)))),) -Base.getindex(a::BlockedPRange,k::Block{1}) = blocks(a)[first(k.n)] -function BlockArrays.findblock(b::BlockedPRange, k::Integer) - @boundscheck k in b || throw(BoundsError(b,k)) - Block(searchsortedfirst(blocklasts(b), k)) -end -Base.first(a::BlockedPRange) = 1 -Base.last(a::BlockedPRange) = sum(length,blocks(a)) - -function PartitionedArrays.partition(a::BlockedPRange) - local_block_ranges(a) -end - -function local_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(local_length,myids) |> collect |> blockedrange - end -end - -function own_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(own_length,myids) |> collect |> blockedrange - end -end - -function ghost_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(ghost_length,myids) |> collect |> blockedrange - end -end - diff --git a/src/p_vector.jl b/src/p_vector.jl index bf3a850c..61fb4b8d 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -1507,258 +1507,3 @@ function renumber(a::PVector,row_partition_2;renumber_local_indices=Val(true)) PVector(values,row_partition_2) end -# BlockPVector - -struct BlockPVector{A,T} <: BlockArrays.AbstractBlockVector{T} - blocks::A - function BlockPVector(blocks) - T = eltype(first(blocks)) - @assert all(block->isa(block,PVector),blocks) - @assert all(block->eltype(block)==T,blocks) - A = typeof(blocks) - new{A,T}(blocks) - end -end -BlockArrays.blocks(a::BlockPVector) = a.blocks -BlockArrays.viewblock(a::BlockPVector, i::Block) = blocks(a)[i.n...] -Base.axes(a::BlockPVector) = (BlockedPRange(map(block->axes(block,1),blocks(a))),) -Base.length(a::BlockPVector) = sum(length,blocks(a)) -Base.IndexStyle(::Type{<:BlockPVector}) = IndexLinear() -function Base.getindex(a::BlockPVector,gid::Int) - scalar_indexing_action(a) -end -function Base.setindex(a::BlockPVector,v,gid::Int) - scalar_indexing_action(a) -end -function Base.show(io::IO,k::MIME"text/plain",data::BlockPVector) - T = eltype(data) - n = length(data) - ps = partition(axes(data,1)) - np = length(ps) - nb = blocklength(data) - map_main(ps) do _ - println(io,"$nb-blocked $n-element BlockPVector with eltype $T partitioned into $np parts.") - end -end -function Base.show(io::IO,data::BlockPVector) - print(io,"BlockPVector(…)") -end - -function partition(a::BlockPVector) - local_values(a) -end - -function local_values(a::BlockPVector) - r = local_block_ranges(axes(a,1)) - v = map(local_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end -end - -function own_values(a::BlockPVector) - r = own_block_ranges(axes(a,1)) - v = map(own_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end -end - -function ghost_values(a::BlockPVector) - r = ghost_block_ranges(axes(a,1)) - v = map(ghost_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end -end - -function consistent!(a::BlockPVector) - ts = map(consistent!,blocks(a)) - @fake_async begin - foreach(wait,ts) - a - end -end - -function assemble!(a::BlockPVector) - ts = map(assemble!,blocks(a)) - @fake_async begin - foreach(wait,ts) - a - end -end - -function Base.similar(a::BlockPVector,::Type{T},inds::Tuple{<:BlockedPRange}) where T - r = first(inds) - bs = map((ai,ri)->similar(ai,T,ri),blocks(a),blocks(r)) - BlockPVector(bs) -end - -function Base.similar(::Type{<:BlockPVector{A}},inds::Tuple{<:BlockedPRange}) where A - V = eltype(A) - r = first(inds) - bs = map(ri->similar(V,ri),blocks(r)) - BlockPVector(bs) -end - -function BlockPVector(::UndefInitializer,r::BlockedPRange) - bs = map(ri->PVector(undef,ri),blocks(r)) - BlockPVector(bs) -end - -function BlockPVector{A}(::UndefInitializer,r::BlockedPRange) where A - similar(BlockPVector{A},(r,)) -end - -function Base.copy!(a::BlockPVector,b::BlockPVector) - foreach(copy!,blocks(a),blocks(b)) - a -end - -function Base.copyto!(a::BlockPVector,b::BlockPVector) - foreach(copyto!,blocks(a),blocks(b)) - a -end - -function Base.fill!(a::BlockPVector,v) - foreach(ai->fill!(ai,v),blocks(a)) -end - - -function Base.:(==)(a::BlockPVector,b::BlockPVector) - all(map(==,blocks(a),blocks(b))) -end - -function Base.any(f::Function,x::BlockPVector) - any(xi->any(f,xi),blocks(x)) -end - -function Base.all(f::Function,x::BlockPVector) - all(xi->all(f,xi),blocks(x)) -end - -Base.maximum(x::BlockPVector) = maximum(identity,x) -function Base.maximum(f::Function,x::BlockPVector) - maximum(xi->maximum(f,xi),blocks(x)) -end - -Base.minimum(x::BlockPVector) = minimum(identity,x) -function Base.minimum(f::Function,x::BlockPVector) - minimum(xi->minimum(f,xi),blocks(x)) -end - -function Base.collect(v::BlockPVector) - reduce(vcat,map(collect,blocks(v))) -end - -function Base.:*(a::Number,b::BlockPVector) - bs = map(bi->a*bi,blocks(b)) - BlockPVector(bs) -end - -function Base.:*(b::BlockPVector,a::Number) - a*b -end - -function Base.:/(b::BlockPVector,a::Number) - (1/a)*b -end - -for op in (:+,:-) - @eval begin - function Base.$op(a::BlockPVector) - bs = map(ai->$op(ai),blocks(a)) - BlockPVector(bs) - end - function Base.$op(a::BlockPVector,b::BlockPVector) - $op.(a,b) - end - end -end - -function Base.reduce(op,a::BlockPVector;neutral=neutral_element(op,eltype(a)),kwargs...) - rs = map(ai->reduce(op,ai;neutral,kwargs...),blocks(a)) - reduce(rs;kwargs...) -end - -function Base.sum(a::BlockPVector) - reduce(+,a,init=zero(eltype(a))) -end - -function LinearAlgebra.dot(a::BlockPVector,b::BlockPVector) - c = map(dot,blocks(a),blocks(b)) - sum(c) -end - -function LinearAlgebra.rmul!(a::BlockPVector,v::Number) - map(blocks(a)) do l - rmul!(l,v) - end - a -end - -function LinearAlgebra.norm(a::BlockPVector,p::Real=2) - contibs = map(blocks(a)) do oid_to_value - norm(oid_to_value,p)^p - end - reduce(+,contibs;init=zero(eltype(contibs)))^(1/p) -end - -for M in Distances.metrics - @eval begin - function (d::$M)(a::BlockPVector,b::BlockPVector) - s = distance_eval_body(d,a,b) - Distances.eval_end(d,s) - end - end -end - -function distance_eval_body(d,a::BlockPVector,b::BlockPVector) - partials = map(blocks(a),blocks(b)) do (ai,bi) - distance_eval_body(d,ai,bi) - end - s = reduce((i,j)->Distances.eval_reduce(d,i,j), - partials, - init=Distances.eval_start(d, a, b)) - s -end - -struct BroadcastedBlockPVector{A} - blocks::A -end -BlockArrays.blocks(a::BroadcastedBlockPVector) = a.blocks - -function Base.broadcasted(f, args::Union{BlockPVector,BroadcastedBlockPVector}...) - map( (bs...) -> Base.broadcasted(f,bs...) , map(blocks,args)...) |> BroadcastedBlockPVector -end - -function Base.broadcasted( f, a::Number, b::Union{BlockPVector,BroadcastedBlockPVector}) - map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBlockPVector -end - -function Base.broadcasted( f, a::Union{BlockPVector,BroadcastedBlockPVector}, b::Number) - map( ai -> Base.broadcasted(f,ai,b), blocks(a)) |> BroadcastedBlockPVector -end - -function Base.broadcasted(f, - a::Union{BlockPVector,BroadcastedBlockPVector}, - b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) - Base.broadcasted(f,a,Base.materialize(b)) -end - -function Base.broadcasted( - f, - a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, - b::Union{BlockPVector,BroadcastedBlockPVector}) - Base.broadcasted(f,Base.materialize(a),b) - end - -function Base.materialize(b::BroadcastedBlockPVector) - map(Base.materialize,blocks(b)) |> BlockPVector -end - -function Base.materialize!(a::BlockPVector,b::BroadcastedBlockPVector) - foreach(Base.materialize!,blocks(a),blocks(b)) - a -end - diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl new file mode 100644 index 00000000..97e02609 --- /dev/null +++ b/test/block_arrays_tests.jl @@ -0,0 +1,73 @@ +using PartitionedArrays +using Test +using LinearAlgebra +using Random +using Distances +using BlockArrays + +function block_arrays_tests(distribute) + + np = 4 + rank = distribute(LinearIndices((np,))) + row_partition = uniform_partition(rank,(2,2),(6,6)) + + a1 = pones(row_partition) + a2 = pzeros(row_partition) + a = mortar([a1,a2]) + display(a) + @test a[Block(1)] == a1 + @test a[Block(2)] == a2 + display(a[Block(1)]) + @test view(a,Block(1)) === a1 + @test view(a,Block(2)) === a2 + + b = similar(a) + b = similar(a,Int) + b = similar(a,Int,axes(a,1)) + b = similar(typeof(a),axes(a,1)) + copy!(b,a) + b = copy(a) + @test typeof(b) == typeof(a) + fill!(b,5.) + rows = axes(a,1) + @test length(a) == length(rows) + assemble!(a) |> wait + consistent!(a) |> wait + + rmul!(a,-1) + fill!(a,0) + @test any(i->i>0,a) == false + @test all(i->i==0,a) == true + @test minimum(a) <= maximum(a) + + b = 2*a + b = a*2 + b = a/2 + c = a .+ a + c = a .+ b .+ a + @test isa(c,BlockPVector) + c = a - b + c = a + b + + fill!(a,1) + r = reduce(+,a) + @test sum(a) == r + @test norm(a) > 0 + @test sqrt(a⋅a) ≈ norm(a) + euclidean(a,a) + @test euclidean(a,a) + 1 ≈ 1 + + u = a + v = b + w = 1 .+ v + @test isa(w,BlockPVector) + w = v .+ 1 + @test isa(w,BlockPVector) + w = v .+ w .- u + @test isa(w,BlockPVector) + w = v .+ 1 .- u + @test isa(w,BlockPVector) + w .= u + +end + diff --git a/test/debug_array/block_arrays_tests.jl b/test/debug_array/block_arrays_tests.jl new file mode 100644 index 00000000..69e1ece5 --- /dev/null +++ b/test/debug_array/block_arrays_tests.jl @@ -0,0 +1,8 @@ +module BlockArrayTests + +include("../block_arrays_tests.jl") + +with_debug(block_arrays_tests) + +end #module + diff --git a/test/debug_array/runtests.jl b/test/debug_array/runtests.jl index b5015a85..2c1a61ab 100644 --- a/test/debug_array/runtests.jl +++ b/test/debug_array/runtests.jl @@ -13,6 +13,8 @@ using PartitionedArrays @testset "p_sparse_matrix" begin include("p_sparse_matrix_tests.jl") end +@testset "block_arrays" begin include("block_arrays_tests.jl") end + @testset "gallery" begin include("gallery_tests.jl") end @testset "p_timer" begin include("p_timer_tests.jl") end diff --git a/test/mpi_array/block_arrays_tests.jl b/test/mpi_array/block_arrays_tests.jl new file mode 100644 index 00000000..28375586 --- /dev/null +++ b/test/mpi_array/block_arrays_tests.jl @@ -0,0 +1,14 @@ +module BlockArraysTests + +using MPI + +project = joinpath(@__DIR__,"..","..") +code = quote + include(joinpath($project,"test","block_arrays_tests.jl")) + with_mpi(block_arrays_tests) +end + +run(`$(mpiexec()) -np 4 julia --project=$project -e $code `) + +end # module + diff --git a/test/mpi_array/runtests.jl b/test/mpi_array/runtests.jl index 87849bea..26a3a5d3 100644 --- a/test/mpi_array/runtests.jl +++ b/test/mpi_array/runtests.jl @@ -9,6 +9,7 @@ using PartitionedArrays @testset "p_vector_tests" begin include("p_vector_tests.jl") end @testset "p_sparse_matrix_tests" begin include("p_sparse_matrix_tests.jl") end @testset "gallery" begin include("gallery_tests.jl") end +@testset "block_arrays" begin include("block_arrays_tests.jl") end @testset "p_timer_tests" begin include("p_timer_tests.jl") end @testset "fdm_example" begin include("fdm_example.jl") end @testset "fem_example" begin include("fem_example.jl") end From 14126e47ca2219be4252c564a391b0d4628e4613 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 14 Aug 2024 15:50:32 +0200 Subject: [PATCH 08/14] Enhancements in BlockPVector --- src/p_vector.jl | 4 ++-- test/block_arrays_tests.jl | 17 +++++++++++++++-- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/p_vector.jl b/src/p_vector.jl index 61fb4b8d..6ed45e45 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -237,8 +237,8 @@ function Base.fill!(a::SplitVector,v) end function Base.:*(a::Number,b::SplitVector) - own_own = a*b.blocks.own - ghost_ghost = a*b.blocks.ghost + own = a*b.blocks.own + ghost = a*b.blocks.ghost blocks = split_vector_blocks(own,ghost) split_vector(blocks,b.permutation) end diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl index 97e02609..841d99e4 100644 --- a/test/block_arrays_tests.jl +++ b/test/block_arrays_tests.jl @@ -11,16 +11,29 @@ function block_arrays_tests(distribute) rank = distribute(LinearIndices((np,))) row_partition = uniform_partition(rank,(2,2),(6,6)) - a1 = pones(row_partition) - a2 = pzeros(row_partition) + a1 = pones(row_partition,split_format=true) + a2 = pzeros(row_partition,split_format=true) a = mortar([a1,a2]) display(a) + rows = axes(a,1) + display(rows) + partition(rows) + local_block_ranges(rows) + own_block_ranges(rows) + ghost_block_ranges(rows) + @test a[Block(1)] == a1 @test a[Block(2)] == a2 display(a[Block(1)]) @test view(a,Block(1)) === a1 @test view(a,Block(2)) === a2 + partition(a) + #map(display,own_values(a)) + local_values(a) + own_values(a) + ghost_values(a) + b = similar(a) b = similar(a,Int) b = similar(a,Int,axes(a,1)) From df88a87bba918b8ccb9fc02f90869a67e2cb1c14 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 14 Aug 2024 15:53:18 +0200 Subject: [PATCH 09/14] Testing block_arrays with and without split format --- test/block_arrays_tests.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl index 841d99e4..c5e4f2f9 100644 --- a/test/block_arrays_tests.jl +++ b/test/block_arrays_tests.jl @@ -6,13 +6,18 @@ using Distances using BlockArrays function block_arrays_tests(distribute) + block_arrays_tests(distribute,false) + block_arrays_tests(distribute,true) +end + +function block_arrays_tests(distribute,split_format) np = 4 rank = distribute(LinearIndices((np,))) row_partition = uniform_partition(rank,(2,2),(6,6)) - a1 = pones(row_partition,split_format=true) - a2 = pzeros(row_partition,split_format=true) + a1 = pones(row_partition;split_format) + a2 = pzeros(row_partition;split_format) a = mortar([a1,a2]) display(a) rows = axes(a,1) From 3d230790fa94f66d151e9e5df3c467cd4f192d41 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 14 Aug 2024 18:52:18 +0200 Subject: [PATCH 10/14] First draft of BlockPSparseMatrix --- src/block_arrays.jl | 198 +++++++++++++++++++++++++++++++++++++ src/p_vector.jl | 33 +++++-- src/primitives.jl | 7 ++ test/block_arrays_tests.jl | 39 ++++++++ 4 files changed, 269 insertions(+), 8 deletions(-) diff --git a/src/block_arrays.jl b/src/block_arrays.jl index e0db22f4..16973208 100644 --- a/src/block_arrays.jl +++ b/src/block_arrays.jl @@ -74,6 +74,7 @@ end BlockArrays.blocks(a::BlockPVector) = a.blocks BlockArrays.viewblock(a::BlockPVector, i::Block) = blocks(a)[i.n...] Base.axes(a::BlockPVector) = (BlockedPRange(map(block->axes(block,1),blocks(a))),) +Base.size(a::BlockPVector) = (length(a),) Base.length(a::BlockPVector) = sum(length,blocks(a)) Base.IndexStyle(::Type{<:BlockPVector}) = IndexLinear() function Base.getindex(a::BlockPVector,gid::Int) @@ -174,6 +175,7 @@ end function Base.fill!(a::BlockPVector,v) foreach(ai->fill!(ai,v),blocks(a)) + a end @@ -314,3 +316,199 @@ function Base.materialize!(a::BlockPVector,b::BroadcastedBlockPVector) a end +# BlockPSparseMatrix + +function BlockArrays.mortar(blocks::Matrix{<:PSparseMatrix}) + BlockPSparseMatrix(blocks) +end + +struct BlockPSparseMatrix{A,T} <: BlockArrays.AbstractBlockMatrix{T} + blocks::A + function BlockPSparseMatrix(blocks) + T = eltype(first(blocks)) + @assert all(block->isa(block,PSparseMatrix),blocks) + @assert all(block->eltype(block)==T,blocks) + A = typeof(blocks) + new{A,T}(blocks) + end +end +BlockArrays.blocks(a::BlockPSparseMatrix) = a.blocks +BlockArrays.viewblock(a::BlockPSparseMatrix, i::Block) = blocks(a)[i.n...] +function Base.axes(a::BlockPSparseMatrix) + rows = BlockedPRange(map(block->axes(block,1),blocks(a)[:,1])) + cols = BlockedPRange(map(block->axes(block,2),blocks(a)[1,:])) + (rows,cols) +end +Base.size(a::BlockPSparseMatrix) = map(length,axes(a)) +function Base.getindex(a::BlockPSparseMatrix,gi::Int,gj::Int) + scalar_indexing_action(a) +end +function Base.setindex!(a::BlockPSparseMatrix,v,gi::Int,gj::Int) + scalar_indexing_action(a) +end +function Base.show(io::IO,k::MIME"text/plain",data::BlockPSparseMatrix) + T = eltype(data) + m,n = size(data) + ps = partition(axes(data,1)) + np = length(ps) + mb,nb = blocksize(data) + map_main(ps) do _ + println(io,"($(mb)×$(nb))-blocked ($(m)×$(n))-element BlockPSparseMatrix with eltype $T partitioned into $np parts.") + end +end +function Base.show(io::IO,data::BlockPSparseMatrix) + print(io,"BlockPSparseMatrix(…)") +end +function partition(a::BlockPSparseMatrix) + local_values(a) +end +function local_values(a::BlockPSparseMatrix) + rows = local_block_ranges(axes(a,1)) + cols = local_block_ranges(axes(a,2)) + v = map(local_values,blocks(a)) |> permute_nesting + map(v,rows,cols) do myv, myrows,mycols + mortar(myv,(myrows,mycols)) + end +end +function own_own_values(a::BlockPSparseMatrix) + rows = own_block_ranges(axes(a,1)) + cols = own_block_ranges(axes(a,2)) + v = map(own_own_values,blocks(a)) |> permute_nesting + map(v,rows,cols) do myv, myrows,mycols + mortar(myv,(myrows,mycols)) + end +end +function own_ghost_values(a::BlockPSparseMatrix) + rows = own_block_ranges(axes(a,1)) + cols = ghost_block_ranges(axes(a,2)) + v = map(own_ghost_values,blocks(a)) |> permute_nesting + map(v,rows,cols) do myv, myrows,mycols + mortar(myv,(myrows,mycols)) + end +end +function ghost_own_values(a::BlockPSparseMatrix) + rows = ghost_block_ranges(axes(a,1)) + cols = own_block_ranges(axes(a,2)) + v = map(ghost_own_values,blocks(a)) |> permute_nesting + map(v,rows,cols) do myv, myrows,mycols + mortar(myv,(myrows,mycols)) + end +end +function ghost_ghost_values(a::BlockPSparseMatrix) + rows = ghost_block_ranges(axes(a,1)) + cols = ghost_block_ranges(axes(a,2)) + v = map(ghost_ghost_values,blocks(a)) |> permute_nesting + map(v,rows,cols) do myv, myrows,mycols + mortar(myv,(myrows,mycols)) + end +end + +function Base.fill!(a::BlockPSparseMatrix,v) + foreach(ai->fill!(ai,v),blocks(a)) + a +end + +function LinearAlgebra.fillstored!(a::BlockPSparseMatrix,v) + foreach(ai->LinearAlgebra.fillstored!(ai,v),blocks(a)) + a +end + +Base.maximum(x::BlockPSparseMatrix) = maximum(identity,x) +function Base.maximum(f::Function,x::BlockPSparseMatrix) + maximum(xi->maximum(f,xi),blocks(x)) +end + +Base.minimum(x::BlockPSparseMatrix) = minimum(identity,x) +function Base.minimum(f::Function,x::BlockPSparseMatrix) + minimum(xi->minimum(f,xi),blocks(x)) +end + +function Base.collect(v::BlockPSparseMatrix) + reduce(vcat,map(collect,blocks(v))) +end + +function centralize(v::BlockPSparseMatrix) + reduce((ai...)->cat(ai...,dims=(1,2)),map(centralize,blocks(v))) +end + +function Base.copy(a::BlockPSparseMatrix) + bs = map(copy,blocks(a)) + BlockPSparseMatrix(bs) +end + +function Base.copy!(a::BlockPSparseMatrix,b::BlockPSparseMatrix) + foreach(copy!,blocks(a),blocks(b)) + a +end + +function Base.copyto!(a::BlockPSparseMatrix,b::BlockPSparseMatrix) + foreach(copyto!,blocks(a),blocks(b)) + a +end + +function SparseArrays.nnz(a::BlockPSparseMatrix) + ns = map(nnz,blocks(a)) + sum(ns) +end + +# This function could be removed if IterativeSolvers was implemented in terms +# of axes(A,d) instead of size(A,d) +function IterativeSolvers.zerox(A::BlockPSparseMatrix,b::BlockPVector) + T = IterativeSolvers.Adivtype(A, b) + x = similar(b, T, axes(A, 2)) + fill!(x, zero(T)) + return x +end + +function Base.:*(a::BlockPSparseMatrix,b::BlockPVector) + Ta = eltype(a) + Tb = eltype(b) + T = typeof(zero(Ta)*zero(Tb)+zero(Ta)*zero(Tb)) + c = similar(b,T,axes(a,1)) + mul!(c,a,b) + c +end + +function Base.:*(a::Number,b::BlockPSparseMatrix) + bs = map(bi->a*bi,blocks(b)) + BlockPSparseMatrix(bs) +end + +function Base.:*(b::BlockPSparseMatrix,a::Number) + a*b +end + +for op in (:+,:-) + @eval begin + function Base.$op(a::BlockPSparseMatrix) + bs = map($op,blocks(a)) + BlockPSparseMatrix(bs) + end + function Base.$op(a::BlockPSparseMatrix,b::BlockPSparseMatrix) + bs = map($op,blocks(a),blocks(b)) + BlockPSparseMatrix(bs) + end + end +end + +function LinearAlgebra.mul!(b::BlockPVector,A::BlockPSparseMatrix,x::BlockPVector) + mul!(b,A,x,1,0) +end + +function LinearAlgebra.mul!(b::BlockPVector,A::BlockPSparseMatrix,x::BlockPVector,α::Number,β::Number) + if β != 1 + β != 0 ? rmul!(b, β) : fill!(b,zero(eltype(b))) + end + bb = blocks(b) + Ab = blocks(A) + xb = blocks(x) + o = one(eltype(b)) + for i in 1:blocksize(A,1) + for j in 1:blocksize(A,2) + mul!(bb[i],Ab[i,j],xb[j],α,o) + end + end + b +end + + diff --git a/src/p_vector.jl b/src/p_vector.jl index 6ed45e45..80da6eac 100644 --- a/src/p_vector.jl +++ b/src/p_vector.jl @@ -527,6 +527,7 @@ struct SplitVectorAssemblyCache{A,B,C,D} buffer_snd::C # NB buffer_rcv::C exchange_setup::D + reversed::Bool end function Base.reverse(a::SplitVectorAssemblyCache) SplitVectorAssemblyCache( @@ -537,6 +538,7 @@ function Base.reverse(a::SplitVectorAssemblyCache) a.buffer_rcv, a.buffer_snd, a.exchange_setup, + !(a.reversed), ) end function copy_cache(a::SplitVectorAssemblyCache) @@ -549,23 +551,29 @@ function copy_cache(a::SplitVectorAssemblyCache) a.own_indices_rcv, buffer_snd, buffer_rcv, - a.exchange_setup + a.exchange_setup, + a.reversed, ) end function p_vector_cache_impl(::Type{<:SplitVector},vector_partition,index_partition) neighbors_snd,neighbors_rcv= assembly_neighbors(index_partition) indices_snd,indices_rcv = assembly_local_indices(index_partition,neighbors_snd,neighbors_rcv) - map(indices_snd,indices_rcv,index_partition) do ids_snd,ids_rcv,myids + ghost_indices_snd = map(indices_snd) do ids + JaggedArray(copy(ids.data),ids.ptrs) + end + own_indices_rcv = map(indices_rcv) do ids + JaggedArray(copy(ids.data),ids.ptrs) + end + foreach(ghost_indices_snd,own_indices_rcv,index_partition) do ids_snd,ids_rcv,myids map_local_to_ghost!(ids_snd.data,myids) map_local_to_own!(ids_rcv.data,myids) end - ghost_indices_snd = indices_snd - own_indices_rcv = indices_rcv - buffers_snd,buffers_rcv = map(assembly_buffers,vector_partition,indices_snd,indices_rcv) |> tuple_of_arrays + buffers_snd,buffers_rcv = map(assembly_buffers,vector_partition,ghost_indices_snd,own_indices_rcv) |> tuple_of_arrays graph = ExchangeGraph(neighbors_snd,neighbors_rcv) exchange_setup = setup_exchange(buffers_rcv,buffers_snd,graph) - SplitVectorAssemblyCache(neighbors_snd,neighbors_rcv,ghost_indices_snd,own_indices_rcv,buffers_snd,buffers_rcv,exchange_setup) + reversed = false + SplitVectorAssemblyCache(neighbors_snd,neighbors_rcv,ghost_indices_snd,own_indices_rcv,buffers_snd,buffers_rcv,exchange_setup,reversed) end function p_vector_cache_impl(::Type{<:SplitVector{<:JaggedArray}},vector_partition,index_partition) @@ -610,6 +618,7 @@ function assemble_impl!(f,vector_partition,cache::JaggedArrayAssemblyCache) end function assemble_impl!(f,vector_partition,cache::SplitVectorAssemblyCache) + reversed = cache.reversed ghost_indices_snd=cache.ghost_indices_snd own_indices_rcv=cache.own_indices_rcv neighbors_snd=cache.neighbors_snd @@ -618,7 +627,11 @@ function assemble_impl!(f,vector_partition,cache::SplitVectorAssemblyCache) buffer_rcv=cache.buffer_rcv exchange_setup=cache.exchange_setup foreach(vector_partition,ghost_indices_snd,buffer_snd) do values,ghost_indices_snd,buffer_snd - ghost_vals = values.blocks.ghost + if reversed + ghost_vals = values.blocks.own + else + ghost_vals = values.blocks.ghost + end for (p,hid) in enumerate(ghost_indices_snd.data) buffer_snd.data[p] = ghost_vals[hid] end @@ -629,7 +642,11 @@ function assemble_impl!(f,vector_partition,cache::SplitVectorAssemblyCache) @fake_async begin wait(t) foreach(vector_partition,own_indices_rcv,buffer_rcv) do values,own_indices_rcv,buffer_rcv - own_vals = values.blocks.own + if reversed + own_vals = values.blocks.ghost + else + own_vals = values.blocks.own + end for (p,oid) in enumerate(own_indices_rcv.data) own_vals[oid] = f(own_vals[oid],buffer_rcv.data[p]) end diff --git a/src/primitives.jl b/src/primitives.jl index 37fb99ab..7b322dfc 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -105,6 +105,13 @@ function array_of_tuples(a) end end +function permute_nesting(a::AbstractArray{<:AbstractArray}) + s = size(a) + map(a...) do b... + reshape(collect(b),s) + end +end + # We don't need a real task since MPI already is able to do # fake_asynchronous (nonblocking) operations. # We want to avoid heap allocations of standard tasks. diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl index c5e4f2f9..404ee988 100644 --- a/test/block_arrays_tests.jl +++ b/test/block_arrays_tests.jl @@ -4,6 +4,7 @@ using LinearAlgebra using Random using Distances using BlockArrays +using SparseArrays function block_arrays_tests(distribute) block_arrays_tests(distribute,false) @@ -20,6 +21,7 @@ function block_arrays_tests(distribute,split_format) a2 = pzeros(row_partition;split_format) a = mortar([a1,a2]) display(a) + collect(a) rows = axes(a,1) display(rows) partition(rows) @@ -87,5 +89,42 @@ function block_arrays_tests(distribute,split_format) @test isa(w,BlockPVector) w .= u + nodes_per_dir = (4,4) + parts_per_dir = (2,2) + args = laplacian_fem(nodes_per_dir,parts_per_dir,rank) + A11 = psparse(args...) |> fetch + x1 = pones(axes(A11,2);split_format) + assemble!(x1) |> wait + consistent!(x1) |> wait + + A = mortar(fill(A11,(2,2))) + display(A) + @show A + #display(centralize(A)) + own_own_values(A) + own_ghost_values(A) + ghost_own_values(A) + ghost_ghost_values(A) + B = copy(A) + copy!(B,A) + copyto!(B,A) + nnz(A) + ax = axes(A,2) + axb = blocks(ax) + x = similar(a,axes(A,2)) + fill!(x,1) + assemble!(x) |> wait + consistent!(x) |> wait + b = similar(x,axes(A,1)) + mul!(b,A,x) + b = A*x + mul!(b,A,x) + B = 2*A + B = A*2 + B = +A + B = -A + C = B+A + D = B-A + end From 7c8e59eba7543402aa00d43dd9c4186ca9fe364d Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 21 Aug 2024 16:44:55 +0200 Subject: [PATCH 11/14] First draft of BRange and BArray --- src/PartitionedArrays.jl | 9 +- src/block_arrays.jl | 470 +++++++++++++++---------------------- src/p_range.jl | 6 +- src/primitives.jl | 4 + test/block_arrays_tests.jl | 57 ++++- 5 files changed, 241 insertions(+), 305 deletions(-) diff --git a/src/PartitionedArrays.jl b/src/PartitionedArrays.jl index df472960..97df913e 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -174,11 +174,10 @@ export spmtm! export centralize include("p_sparse_matrix.jl") -export BlockedPRange -export local_block_ranges -export own_block_ranges -export ghost_block_ranges -export BlockPVector +export BRange +export BArray +export BVector +export BMatrix include("block_arrays.jl") export PTimer diff --git a/src/block_arrays.jl b/src/block_arrays.jl index 16973208..555aeeaa 100644 --- a/src/block_arrays.jl +++ b/src/block_arrays.jl @@ -1,13 +1,8 @@ -function BlockArrays.mortar(blocks::Vector{<:PRange}) - BlockedPRange(blocks) -end - -struct BlockedPRange{A} <: AbstractUnitRange{Int} +struct BRange{A} <: AbstractUnitRange{Int} blocks::A blocklasts::Vector{Int} - function BlockedPRange(blocks) - @assert all(block->isa(block,PRange),blocks) + function BRange(blocks) nblocks = length(blocks) blocklasts = zeros(Int,nblocks) prev = 0 @@ -19,113 +14,150 @@ struct BlockedPRange{A} <: AbstractUnitRange{Int} new{A}(blocks,blocklasts) end end -BlockArrays.blocks(a::BlockedPRange) = a.blocks -BlockArrays.blocklasts(a::BlockedPRange) = a.blocklasts -BlockArrays.blockaxes(a::BlockedPRange) = (Block.(Base.OneTo(length(blocks(a)))),) -Base.getindex(a::BlockedPRange,k::Block{1}) = blocks(a)[first(k.n)] -function BlockArrays.findblock(b::BlockedPRange, k::Integer) +BlockArrays.blocks(a::BRange) = a.blocks +BlockArrays.blocklasts(a::BRange) = a.blocklasts +BlockArrays.blockaxes(a::BRange) = (Block.(Base.OneTo(length(blocks(a)))),) +Base.getindex(a::BRange,k::Block{1}) = blocks(a)[first(k.n)] +function BlockArrays.findblock(b::BRange, k::Integer) @boundscheck k in b || throw(BoundsError(b,k)) Block(searchsortedfirst(blocklasts(b), k)) end -Base.first(a::BlockedPRange) = 1 -Base.last(a::BlockedPRange) = sum(length,blocks(a)) - -function PartitionedArrays.partition(a::BlockedPRange) - local_block_ranges(a) -end +Base.first(a::BRange) = 1 +Base.last(a::BRange) = sum(length,blocks(a)) -function local_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(local_length,myids) |> collect |> blockedrange +function Base.show(io::IO,k::MIME"text/plain",data::BRange) + function tostr(a) + "$a" end -end - -function own_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(own_length,myids) |> collect |> blockedrange + function tostr(a::PRange) + np = length(partition(a)) + "PRange 1:$(length(a)) partitioned into $(np) parts" end -end - -function ghost_block_ranges(a::BlockedPRange) - ids = map(partition,blocks(a)) |> array_of_tuples - map(ids) do myids - map(ghost_length,myids) |> collect |> blockedrange + nb = blocklength(data) + println(io,"BRange $(first(data)):$(last(data)) with $nb blocks:") + for ib in 1:nb + t = ib != nb ? "├──" : "└──" + println(io,"$(t) Block($ib) = $(tostr(blocks(data)[ib]))") end end -# BlockPVector +function Base.show(io::IO,data::BRange) + print(io,"BRange(…)") +end -function BlockArrays.mortar(blocks::Vector{<:PVector}) - BlockPVector(blocks) +function partition(a::BRange) + ps = map(partition,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) end -struct BlockPVector{A,T} <: BlockArrays.AbstractBlockVector{T} +struct BArray{A,T,N} <: BlockArrays.AbstractBlockArray{T,N} blocks::A - function BlockPVector(blocks) + function BArray(blocks) T = eltype(first(blocks)) - @assert all(block->isa(block,PVector),blocks) + N = ndims(first(blocks)) @assert all(block->eltype(block)==T,blocks) + @assert all(block->ndims(block)==N,blocks) A = typeof(blocks) - new{A,T}(blocks) + new{A,T,N}(blocks) end end -BlockArrays.blocks(a::BlockPVector) = a.blocks -BlockArrays.viewblock(a::BlockPVector, i::Block) = blocks(a)[i.n...] -Base.axes(a::BlockPVector) = (BlockedPRange(map(block->axes(block,1),blocks(a))),) -Base.size(a::BlockPVector) = (length(a),) -Base.length(a::BlockPVector) = sum(length,blocks(a)) -Base.IndexStyle(::Type{<:BlockPVector}) = IndexLinear() -function Base.getindex(a::BlockPVector,gid::Int) +const BVector{A,T} = BArray{A,T,1} +const BMatrix{A,T} = BArray{A,T,2} +function BVector(blocks) + N = ndims(first(blocks)) + @assert N==1 + BArray(blocks) +end +function BMatrix(blocks) + N = ndims(first(blocks)) + @assert N==2 + BArray(blocks) +end +BlockArrays.blocks(a::BArray) = a.blocks +BlockArrays.viewblock(a::BArray, i::Block) = blocks(a)[i.n...] +BlockArrays.blockaxes(a::BArray,d::Int) = Block.(Base.OneTo(size(blocks(a),d))) +BlockArrays.blockaxes(a::BArray) = map(i->blockaxes(a,i),ntuple(j->Int(j),ndims(a))) +function Base.axes(a::BArray,d::Int) + N = ndims(a) + I = ntuple(i-> (i == d) ? (:) : 1, Val(N)) + BRange(map(block->axes(block,d),blocks(a)[I...])) +end +Base.axes(a::BArray) = map(i->axes(a,i),ntuple(j->Int(j),ndims(a))) +Base.size(a::BArray) = map(length,axes(a)) +Base.length(a::BArray) = sum(length,blocks(a)) +Base.IndexStyle(::Type{<:BArray}) = IndexCartesian() +function Base.getindex(a::BArray{A,T,N} where {A,T},gid::Vararg{Int,N}) where N + # This could be implemented as it makes sense when the blocks + # are sequential arrays scalar_indexing_action(a) end -function Base.setindex(a::BlockPVector,v,gid::Int) +function Base.setindex(a::BArray{A,T,N} where {A,T},v,gid::Vararg{Int,N}) where N + # This could be implemented as it makes sense when the blocks + # are sequential arrays scalar_indexing_action(a) end -function Base.show(io::IO,k::MIME"text/plain",data::BlockPVector) - T = eltype(data) - n = length(data) - ps = partition(axes(data,1)) - np = length(ps) - nb = blocklength(data) - map_main(ps) do _ - println(io,"$nb-blocked $n-element BlockPVector with eltype $T partitioned into $np parts.") +function Base.show(io::IO,k::MIME"text/plain",data::BArray) + N = ndims(data) + if N == 1 + at = "BVector" + elseif N==2 + at = "BMatrix" + else + at = "BArray" + end + bs = blocksize(data) + bst = map(i->"$(i)×",collect(bs)) + bst[end] = string(bs[end]) + s = size(data) + st = map(i->"$(i)×",collect(s)) + st[end] = string(s[end]) + println(io,"$(join(st))-element $at with $(join(bst)) blocks:") + for cb in CartesianIndices(bs) + ib = LinearIndices(bs)[cb] + nb = prod(bs) + t = ib != nb ? "├──" : "└──" + tcb =Tuple(cb) + b = map(i->"$i, ",collect(tcb)) + b[end] = string(tcb[end]) + println(io,"$(t) Block($(join(b))) = $(blocks(data)[cb])") end end -function Base.show(io::IO,data::BlockPVector) - print(io,"BlockPVector(…)") +function Base.show(io::IO,data::BArray) + print(io,"BArray(…)") +end +function Base.show(io::IO,data::BVector) + print(io,"BVector(…)") +end +function Base.show(io::IO,data::BMatrix) + print(io,"BMatrix(…)") end -function partition(a::BlockPVector) - local_values(a) +function partition(a::BArray) + ps = map(partition,blocks(a)) + ps2 = permute_nesting(ps) + map(BArray,ps2) end -function local_values(a::BlockPVector) - r = local_block_ranges(axes(a,1)) - v = map(local_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end +function local_values(a::BVector) + ps = map(local_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) end -function own_values(a::BlockPVector) - r = own_block_ranges(axes(a,1)) - v = map(own_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end +function own_values(a::BVector) + ps = map(own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) end -function ghost_values(a::BlockPVector) - r = ghost_block_ranges(axes(a,1)) - v = map(ghost_values,blocks(a)) |> array_of_tuples - map(v,r) do myv, myr - mortar(collect(myv),(myr,)) - end +function ghost_values(a::BVector) + ps = map(ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) end -function consistent!(a::BlockPVector) +function consistent!(a::BVector) ts = map(consistent!,blocks(a)) @fake_async begin foreach(wait,ts) @@ -133,7 +165,7 @@ function consistent!(a::BlockPVector) end end -function assemble!(a::BlockPVector) +function assemble!(a::BVector) ts = map(assemble!,blocks(a)) @fake_async begin foreach(wait,ts) @@ -141,117 +173,105 @@ function assemble!(a::BlockPVector) end end -function Base.similar(a::BlockPVector,::Type{T},inds::Tuple{<:BlockedPRange}) where T +function Base.similar(a::BVector,::Type{T},inds::Tuple{<:BRange}) where T r = first(inds) bs = map((ai,ri)->similar(ai,T,ri),blocks(a),blocks(r)) - BlockPVector(bs) + BVector(bs) end -function Base.similar(::Type{<:BlockPVector{A}},inds::Tuple{<:BlockedPRange}) where A - V = eltype(A) - r = first(inds) - bs = map(ri->similar(V,ri),blocks(r)) - BlockPVector(bs) +function Base.copy(a::BArray) + map(copy,blocks(a)) |> BArray end -function BlockPVector(::UndefInitializer,r::BlockedPRange) - bs = map(ri->PVector(undef,ri),blocks(r)) - BlockPVector(bs) -end - -function BlockPVector{A}(::UndefInitializer,r::BlockedPRange) where A - similar(BlockPVector{A},(r,)) -end - -function Base.copy!(a::BlockPVector,b::BlockPVector) +function Base.copy!(a::BArray,b::BArray) foreach(copy!,blocks(a),blocks(b)) a end -function Base.copyto!(a::BlockPVector,b::BlockPVector) +function Base.copyto!(a::BArray,b::BArray) foreach(copyto!,blocks(a),blocks(b)) a end -function Base.fill!(a::BlockPVector,v) +function Base.fill!(a::BArray,v) foreach(ai->fill!(ai,v),blocks(a)) a end - -function Base.:(==)(a::BlockPVector,b::BlockPVector) +function Base.:(==)(a::BArray,b::BArray) all(map(==,blocks(a),blocks(b))) end -function Base.any(f::Function,x::BlockPVector) +function Base.any(f::Function,x::BArray) any(xi->any(f,xi),blocks(x)) end -function Base.all(f::Function,x::BlockPVector) +function Base.all(f::Function,x::BArray) all(xi->all(f,xi),blocks(x)) end -Base.maximum(x::BlockPVector) = maximum(identity,x) -function Base.maximum(f::Function,x::BlockPVector) +Base.maximum(x::BArray) = maximum(identity,x) +function Base.maximum(f::Function,x::BArray) maximum(xi->maximum(f,xi),blocks(x)) end -Base.minimum(x::BlockPVector) = minimum(identity,x) -function Base.minimum(f::Function,x::BlockPVector) +Base.minimum(x::BArray) = minimum(identity,x) +function Base.minimum(f::Function,x::BArray) minimum(xi->minimum(f,xi),blocks(x)) end -function Base.collect(v::BlockPVector) +function Base.collect(v::BVector) reduce(vcat,map(collect,blocks(v))) end -function Base.:*(a::Number,b::BlockPVector) +function Base.:*(a::Number,b::BArray) bs = map(bi->a*bi,blocks(b)) - BlockPVector(bs) + BArray(bs) end -function Base.:*(b::BlockPVector,a::Number) +function Base.:*(b::BArray,a::Number) a*b end -function Base.:/(b::BlockPVector,a::Number) +function Base.:/(b::BArray,a::Number) (1/a)*b end for op in (:+,:-) @eval begin - function Base.$op(a::BlockPVector) - bs = map(ai->$op(ai),blocks(a)) - BlockPVector(bs) + function Base.$op(a::BArray) + bs = map($op,blocks(a)) + BArray(bs) end - function Base.$op(a::BlockPVector,b::BlockPVector) - $op.(a,b) + function Base.$op(a::BArray,b::BArray) + bs = map($op,blocks(a),blocks(b)) + BArray(bs) end end end -function Base.reduce(op,a::BlockPVector;neutral=neutral_element(op,eltype(a)),kwargs...) +function Base.reduce(op,a::BArray;neutral=neutral_element(op,eltype(a)),kwargs...) rs = map(ai->reduce(op,ai;neutral,kwargs...),blocks(a)) reduce(op,rs;kwargs...) end -function Base.sum(a::BlockPVector) +function Base.sum(a::BArray) reduce(+,a,init=zero(eltype(a))) end -function LinearAlgebra.dot(a::BlockPVector,b::BlockPVector) +function LinearAlgebra.dot(a::BVector,b::BVector) c = map(dot,blocks(a),blocks(b)) sum(c) end -function LinearAlgebra.rmul!(a::BlockPVector,v::Number) +function LinearAlgebra.rmul!(a::BArray,v::Number) map(blocks(a)) do l rmul!(l,v) end a end -function LinearAlgebra.norm(a::BlockPVector,p::Real=2) +function LinearAlgebra.norm(a::BVector,p::Real=2) contibs = map(blocks(a)) do oid_to_value norm(oid_to_value,p)^p end @@ -260,14 +280,14 @@ end for M in Distances.metrics @eval begin - function (d::$M)(a::BlockPVector,b::BlockPVector) + function (d::$M)(a::BVector,b::BVector) s = distance_eval_body(d,a,b) Distances.eval_end(d,s) end end end -function distance_eval_body(d,a::BlockPVector,b::BlockPVector) +function distance_eval_body(d,a::BVector,b::BVector) partials = map(blocks(a),blocks(b)) do ai,bi distance_eval_body(d,ai,bi) end @@ -277,25 +297,25 @@ function distance_eval_body(d,a::BlockPVector,b::BlockPVector) s end -struct BroadcastedBlockPVector{A} +struct BroadcastedBArray{A} blocks::A end -BlockArrays.blocks(a::BroadcastedBlockPVector) = a.blocks +BlockArrays.blocks(a::BroadcastedBArray) = a.blocks -function Base.broadcasted(f, args::Union{BlockPVector,BroadcastedBlockPVector}...) - map( (bs...) -> Base.broadcasted(f,bs...) , map(blocks,args)...) |> BroadcastedBlockPVector +function Base.broadcasted(f, args::Union{BArray,BroadcastedBArray}...) + map( (bs...) -> Base.broadcasted(f,bs...) , map(blocks,args)...) |> BroadcastedBArray end -function Base.broadcasted( f, a::Number, b::Union{BlockPVector,BroadcastedBlockPVector}) - map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBlockPVector +function Base.broadcasted( f, a::Number, b::Union{BArray,BroadcastedBArray}) + map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBArray end -function Base.broadcasted( f, a::Union{BlockPVector,BroadcastedBlockPVector}, b::Number) - map( ai -> Base.broadcasted(f,ai,b), blocks(a)) |> BroadcastedBlockPVector +function Base.broadcasted( f, a::Union{BArray,BroadcastedBArray}, b::Number) + map( ai -> Base.broadcasted(f,ai,b), blocks(a)) |> BroadcastedBArray end function Base.broadcasted(f, - a::Union{BlockPVector,BroadcastedBlockPVector}, + a::Union{BArray,BroadcastedBArray}, b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) Base.broadcasted(f,a,Base.materialize(b)) end @@ -303,164 +323,65 @@ end function Base.broadcasted( f, a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, - b::Union{BlockPVector,BroadcastedBlockPVector}) + b::Union{BArray,BroadcastedBArray}) Base.broadcasted(f,Base.materialize(a),b) end -function Base.materialize(b::BroadcastedBlockPVector) - map(Base.materialize,blocks(b)) |> BlockPVector +function Base.materialize(b::BroadcastedBArray) + map(Base.materialize,blocks(b)) |> BArray end -function Base.materialize!(a::BlockPVector,b::BroadcastedBlockPVector) +function Base.materialize!(a::BArray,b::BroadcastedBArray) foreach(Base.materialize!,blocks(a),blocks(b)) a end -# BlockPSparseMatrix - -function BlockArrays.mortar(blocks::Matrix{<:PSparseMatrix}) - BlockPSparseMatrix(blocks) -end - -struct BlockPSparseMatrix{A,T} <: BlockArrays.AbstractBlockMatrix{T} - blocks::A - function BlockPSparseMatrix(blocks) - T = eltype(first(blocks)) - @assert all(block->isa(block,PSparseMatrix),blocks) - @assert all(block->eltype(block)==T,blocks) - A = typeof(blocks) - new{A,T}(blocks) - end -end -BlockArrays.blocks(a::BlockPSparseMatrix) = a.blocks -BlockArrays.viewblock(a::BlockPSparseMatrix, i::Block) = blocks(a)[i.n...] -function Base.axes(a::BlockPSparseMatrix) - rows = BlockedPRange(map(block->axes(block,1),blocks(a)[:,1])) - cols = BlockedPRange(map(block->axes(block,2),blocks(a)[1,:])) - (rows,cols) +function own_own_values(a::BMatrix) + ps = map(own_own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) end -Base.size(a::BlockPSparseMatrix) = map(length,axes(a)) -function Base.getindex(a::BlockPSparseMatrix,gi::Int,gj::Int) - scalar_indexing_action(a) -end -function Base.setindex!(a::BlockPSparseMatrix,v,gi::Int,gj::Int) - scalar_indexing_action(a) -end -function Base.show(io::IO,k::MIME"text/plain",data::BlockPSparseMatrix) - T = eltype(data) - m,n = size(data) - ps = partition(axes(data,1)) - np = length(ps) - mb,nb = blocksize(data) - map_main(ps) do _ - println(io,"($(mb)×$(nb))-blocked ($(m)×$(n))-element BlockPSparseMatrix with eltype $T partitioned into $np parts.") - end -end -function Base.show(io::IO,data::BlockPSparseMatrix) - print(io,"BlockPSparseMatrix(…)") +function own_ghost_values(a::BMatrix) + ps = map(own_ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) end -function partition(a::BlockPSparseMatrix) - local_values(a) +function ghost_own_values(a::BMatrix) + ps = map(ghost_own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) end -function local_values(a::BlockPSparseMatrix) - rows = local_block_ranges(axes(a,1)) - cols = local_block_ranges(axes(a,2)) - v = map(local_values,blocks(a)) |> permute_nesting - map(v,rows,cols) do myv, myrows,mycols - mortar(myv,(myrows,mycols)) - end -end -function own_own_values(a::BlockPSparseMatrix) - rows = own_block_ranges(axes(a,1)) - cols = own_block_ranges(axes(a,2)) - v = map(own_own_values,blocks(a)) |> permute_nesting - map(v,rows,cols) do myv, myrows,mycols - mortar(myv,(myrows,mycols)) - end -end -function own_ghost_values(a::BlockPSparseMatrix) - rows = own_block_ranges(axes(a,1)) - cols = ghost_block_ranges(axes(a,2)) - v = map(own_ghost_values,blocks(a)) |> permute_nesting - map(v,rows,cols) do myv, myrows,mycols - mortar(myv,(myrows,mycols)) - end -end -function ghost_own_values(a::BlockPSparseMatrix) - rows = ghost_block_ranges(axes(a,1)) - cols = own_block_ranges(axes(a,2)) - v = map(ghost_own_values,blocks(a)) |> permute_nesting - map(v,rows,cols) do myv, myrows,mycols - mortar(myv,(myrows,mycols)) - end -end -function ghost_ghost_values(a::BlockPSparseMatrix) - rows = ghost_block_ranges(axes(a,1)) - cols = ghost_block_ranges(axes(a,2)) - v = map(ghost_ghost_values,blocks(a)) |> permute_nesting - map(v,rows,cols) do myv, myrows,mycols - mortar(myv,(myrows,mycols)) - end +function ghost_ghost_values(a::BMatrix) + ps = map(ghost_ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) end -function Base.fill!(a::BlockPSparseMatrix,v) - foreach(ai->fill!(ai,v),blocks(a)) - a -end - -function LinearAlgebra.fillstored!(a::BlockPSparseMatrix,v) +function LinearAlgebra.fillstored!(a::BMatrix,v) foreach(ai->LinearAlgebra.fillstored!(ai,v),blocks(a)) a end -Base.maximum(x::BlockPSparseMatrix) = maximum(identity,x) -function Base.maximum(f::Function,x::BlockPSparseMatrix) - maximum(xi->maximum(f,xi),blocks(x)) -end - -Base.minimum(x::BlockPSparseMatrix) = minimum(identity,x) -function Base.minimum(f::Function,x::BlockPSparseMatrix) - minimum(xi->minimum(f,xi),blocks(x)) -end - -function Base.collect(v::BlockPSparseMatrix) - reduce(vcat,map(collect,blocks(v))) -end - -function centralize(v::BlockPSparseMatrix) - reduce((ai...)->cat(ai...,dims=(1,2)),map(centralize,blocks(v))) -end - -function Base.copy(a::BlockPSparseMatrix) - bs = map(copy,blocks(a)) - BlockPSparseMatrix(bs) -end - -function Base.copy!(a::BlockPSparseMatrix,b::BlockPSparseMatrix) - foreach(copy!,blocks(a),blocks(b)) - a -end - -function Base.copyto!(a::BlockPSparseMatrix,b::BlockPSparseMatrix) - foreach(copyto!,blocks(a),blocks(b)) - a -end +#function centralize(v::BMatrix) +# # TODO not correct +# reduce((ai...)->cat(ai...,dims=(1,2)),map(centralize,blocks(v))) +#end -function SparseArrays.nnz(a::BlockPSparseMatrix) +function SparseArrays.nnz(a::BMatrix) ns = map(nnz,blocks(a)) sum(ns) end # This function could be removed if IterativeSolvers was implemented in terms # of axes(A,d) instead of size(A,d) -function IterativeSolvers.zerox(A::BlockPSparseMatrix,b::BlockPVector) +function IterativeSolvers.zerox(A::BMatrix,b::BVector) T = IterativeSolvers.Adivtype(A, b) x = similar(b, T, axes(A, 2)) fill!(x, zero(T)) return x end -function Base.:*(a::BlockPSparseMatrix,b::BlockPVector) +function Base.:*(a::BMatrix,b::BVector) Ta = eltype(a) Tb = eltype(b) T = typeof(zero(Ta)*zero(Tb)+zero(Ta)*zero(Tb)) @@ -469,33 +390,11 @@ function Base.:*(a::BlockPSparseMatrix,b::BlockPVector) c end -function Base.:*(a::Number,b::BlockPSparseMatrix) - bs = map(bi->a*bi,blocks(b)) - BlockPSparseMatrix(bs) -end - -function Base.:*(b::BlockPSparseMatrix,a::Number) - a*b -end - -for op in (:+,:-) - @eval begin - function Base.$op(a::BlockPSparseMatrix) - bs = map($op,blocks(a)) - BlockPSparseMatrix(bs) - end - function Base.$op(a::BlockPSparseMatrix,b::BlockPSparseMatrix) - bs = map($op,blocks(a),blocks(b)) - BlockPSparseMatrix(bs) - end - end -end - -function LinearAlgebra.mul!(b::BlockPVector,A::BlockPSparseMatrix,x::BlockPVector) +function LinearAlgebra.mul!(b::BVector,A::BMatrix,x::BVector) mul!(b,A,x,1,0) end -function LinearAlgebra.mul!(b::BlockPVector,A::BlockPSparseMatrix,x::BlockPVector,α::Number,β::Number) +function LinearAlgebra.mul!(b::BVector,A::BMatrix,x::BVector,α::Number,β::Number) if β != 1 β != 0 ? rmul!(b, β) : fill!(b,zero(eltype(b))) end @@ -511,4 +410,3 @@ function LinearAlgebra.mul!(b::BlockPVector,A::BlockPSparseMatrix,x::BlockPVecto b end - diff --git a/src/p_range.jl b/src/p_range.jl index cfa330b0..aa7e7b7f 100644 --- a/src/p_range.jl +++ b/src/p_range.jl @@ -1808,10 +1808,14 @@ Base.last(a::PRange) = getany(map(global_length,partition(a))) function Base.show(io::IO,k::MIME"text/plain",data::PRange) np = length(partition(data)) map_main(partition(data)) do indices - println(io,"1:$(global_length(indices)) partitioned into $(np) parts") + println(io,"PRange 1:$(global_length(indices)) partitioned into $(np) parts") end end +function Base.show(io::IO,data::PRange) + print(io,"PRange(…)") +end + function matching_local_indices(a::PRange,b::PRange) partition(a) === partition(b) && return true c = map(matching_local_indices,partition(a),partition(b)) diff --git a/src/primitives.jl b/src/primitives.jl index 7b322dfc..8fa98f0b 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -112,6 +112,10 @@ function permute_nesting(a::AbstractArray{<:AbstractArray}) end end +function permute_nesting(a::Tuple{Vararg{<:AbstractArray}}) + array_of_tuples(a) +end + # We don't need a real task since MPI already is able to do # fake_asynchronous (nonblocking) operations. # We want to avoid heap allocations of standard tasks. diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl index 404ee988..b1567aaa 100644 --- a/test/block_arrays_tests.jl +++ b/test/block_arrays_tests.jl @@ -17,17 +17,43 @@ function block_arrays_tests(distribute,split_format) rank = distribute(LinearIndices((np,))) row_partition = uniform_partition(rank,(2,2),(6,6)) + r1 = PRange(row_partition) + r2 = r1 + + display(r1) + @show r1 + + r = BRange([r1,r2]) + + display(r) + @show r + + partition(r) |> display + + r = BRange([r1,r2]) + + display(r) + @show r + + partition(r) |> display + + + a1 = pones(row_partition;split_format) a2 = pzeros(row_partition;split_format) - a = mortar([a1,a2]) + a = BVector([a1,a2]) display(a) + + b = BVector([[1,2,3],[4,5,6,7]]) + display(b) + @test size(b) == (7,) + @test blocksize(b) == (2,) + @test blocklength(b) == 2 + collect(a) rows = axes(a,1) - display(rows) + @test isa(rows,BRange) partition(rows) - local_block_ranges(rows) - own_block_ranges(rows) - ghost_block_ranges(rows) @test a[Block(1)] == a1 @test a[Block(2)] == a2 @@ -44,7 +70,6 @@ function block_arrays_tests(distribute,split_format) b = similar(a) b = similar(a,Int) b = similar(a,Int,axes(a,1)) - b = similar(typeof(a),axes(a,1)) copy!(b,a) b = copy(a) @test typeof(b) == typeof(a) @@ -65,7 +90,7 @@ function block_arrays_tests(distribute,split_format) b = a/2 c = a .+ a c = a .+ b .+ a - @test isa(c,BlockPVector) + @test isa(c,BVector) c = a - b c = a + b @@ -80,13 +105,13 @@ function block_arrays_tests(distribute,split_format) u = a v = b w = 1 .+ v - @test isa(w,BlockPVector) + @test isa(w,BVector) w = v .+ 1 - @test isa(w,BlockPVector) + @test isa(w,BVector) w = v .+ w .- u - @test isa(w,BlockPVector) + @test isa(w,BVector) w = v .+ 1 .- u - @test isa(w,BlockPVector) + @test isa(w,BVector) w .= u nodes_per_dir = (4,4) @@ -97,10 +122,16 @@ function block_arrays_tests(distribute,split_format) assemble!(x1) |> wait consistent!(x1) |> wait - A = mortar(fill(A11,(2,2))) + @test size(A11) == (16,16) + A = BMatrix(fill(A11,(2,2))) display(A) @show A - #display(centralize(A)) + + display(axes(A,1)) + + @test blocksize(A) == (2,2) + @test size(A) == (32,32) + own_own_values(A) own_ghost_values(A) ghost_own_values(A) From ce9cf8aae415e9f8edd4823f5c4c4021b242f057 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 21 Aug 2024 16:55:44 +0200 Subject: [PATCH 12/14] More tests --- test/block_arrays_tests.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/test/block_arrays_tests.jl b/test/block_arrays_tests.jl index b1567aaa..5909edc9 100644 --- a/test/block_arrays_tests.jl +++ b/test/block_arrays_tests.jl @@ -5,6 +5,7 @@ using Random using Distances using BlockArrays using SparseArrays +using IterativeSolvers function block_arrays_tests(distribute) block_arrays_tests(distribute,false) @@ -143,13 +144,14 @@ function block_arrays_tests(distribute,split_format) ax = axes(A,2) axb = blocks(ax) x = similar(a,axes(A,2)) + @test isa(x,BVector) fill!(x,1) assemble!(x) |> wait consistent!(x) |> wait b = similar(x,axes(A,1)) mul!(b,A,x) b = A*x - mul!(b,A,x) + @test isa(b,BVector) B = 2*A B = A*2 B = +A @@ -157,5 +159,11 @@ function block_arrays_tests(distribute,split_format) C = B+A D = B-A + y = copy(x) + fill!(y,0) + IterativeSolvers.cg!(y,A,b,verbose=i_am_main(rank)) + y = IterativeSolvers.cg(A,b,verbose=i_am_main(rank)) + @test isa(y,BVector) + end From 28793603aab15254e421e855e3fec4787ac0106a Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 21 Aug 2024 17:02:51 +0200 Subject: [PATCH 13/14] Adding compat entry for BlockArrays --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index e8f39040..f3ef6ef2 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] +BlockArrays = "0.16" CircularArrays = "1" Distances = "0.10" FillArrays = "0.10, 0.11, 0.12, 0.13, 1" From 37ad378e7cf6650ce4ce67c2962599c53c4519a9 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 21 Aug 2024 17:10:11 +0200 Subject: [PATCH 14/14] Fixing tests --- test/p_range_tests.jl | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/test/p_range_tests.jl b/test/p_range_tests.jl index d59f3635..05dd1d32 100644 --- a/test/p_range_tests.jl +++ b/test/p_range_tests.jl @@ -1,7 +1,6 @@ using PartitionedArrays using Test -using BlockArrays using PartitionedArrays: local_range @@ -296,18 +295,4 @@ function p_range_tests(distribute) @test ids1 == ids2 end - ids1 = uniform_partition(parts,(2,2),(5,4)) - ids2 = uniform_partition(parts,(2,2),(6,3)) - ax = BlockedPRange(PRange.([ids1,ids2])) - blockaxes(ax,1) - @test isa(ax[Block(1)],PRange) - @test isa(ax[Block(2)],PRange) - @test blocklasts(ax) == [20, 38] - @test findblock(ax,4) == Block(1) - @test findblock(ax,25) == Block(2) - partition(ax) - local_block_ranges(ax) - own_block_ranges(ax) - ghost_block_ranges(ax) - end