diff --git a/CHANGELOG.md b/CHANGELOG.md index aa8b5d39..b22c5b24 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,9 +5,18 @@ 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). +## "master" branch + +### Added + +- Function `array_of_tuples`. +- Export statement for `local_permutation`. +- Experimental support for block arrays via types `BRange`, `BVector`, `BMatrix`, and `BArray`. + ## [0.5.3] - 2024-08-16 ### Fixed + - Typo: `node_coorinates_unit_cube` -> `node_coordinates_unit_cube`. - Bug in `nullspace_linear_elasticity`. - Bug in `PVector` when working in split format. diff --git a/Project.toml b/Project.toml index 73c8b770..f3ef6ef2 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Francesc Verdugo and contributors"] version = "0.5.3" [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" @@ -17,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" 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 eef27803..a402de0a 100644 --- a/src/PartitionedArrays.jl +++ b/src/PartitionedArrays.jl @@ -9,6 +9,7 @@ using StaticArrays import MPI import IterativeSolvers import Distances +using BlockArrays export length_to_ptrs! export rewind_ptrs! @@ -28,6 +29,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 @@ -95,6 +97,7 @@ export own_to_local export ghost_to_local export local_to_own export local_to_ghost +export local_permutation export global_to_owner export replace_ghost export remove_ghost @@ -171,6 +174,12 @@ export spmtm! export centralize include("p_sparse_matrix.jl") +export BRange +export BArray +export BVector +export BMatrix +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..555aeeaa --- /dev/null +++ b/src/block_arrays.jl @@ -0,0 +1,412 @@ + +struct BRange{A} <: AbstractUnitRange{Int} + blocks::A + blocklasts::Vector{Int} + function BRange(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::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::BRange) = 1 +Base.last(a::BRange) = sum(length,blocks(a)) + +function Base.show(io::IO,k::MIME"text/plain",data::BRange) + function tostr(a) + "$a" + end + function tostr(a::PRange) + np = length(partition(a)) + "PRange 1:$(length(a)) partitioned into $(np) parts" + end + 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 + +function Base.show(io::IO,data::BRange) + print(io,"BRange(…)") +end + +function partition(a::BRange) + ps = map(partition,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) +end + +struct BArray{A,T,N} <: BlockArrays.AbstractBlockArray{T,N} + blocks::A + function BArray(blocks) + T = eltype(first(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,N}(blocks) + end +end +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::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::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::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::BArray) + ps = map(partition,blocks(a)) + ps2 = permute_nesting(ps) + map(BArray,ps2) +end + +function local_values(a::BVector) + ps = map(local_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) +end + +function own_values(a::BVector) + ps = map(own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) +end + +function ghost_values(a::BVector) + ps = map(ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BVector,ps2) +end + +function consistent!(a::BVector) + ts = map(consistent!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + +function assemble!(a::BVector) + ts = map(assemble!,blocks(a)) + @fake_async begin + foreach(wait,ts) + a + end +end + +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)) + BVector(bs) +end + +function Base.copy(a::BArray) + map(copy,blocks(a)) |> BArray +end + +function Base.copy!(a::BArray,b::BArray) + foreach(copy!,blocks(a),blocks(b)) + a +end + +function Base.copyto!(a::BArray,b::BArray) + foreach(copyto!,blocks(a),blocks(b)) + a +end + +function Base.fill!(a::BArray,v) + foreach(ai->fill!(ai,v),blocks(a)) + a +end + +function Base.:(==)(a::BArray,b::BArray) + all(map(==,blocks(a),blocks(b))) +end + +function Base.any(f::Function,x::BArray) + any(xi->any(f,xi),blocks(x)) +end + +function Base.all(f::Function,x::BArray) + all(xi->all(f,xi),blocks(x)) +end + +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::BArray) = minimum(identity,x) +function Base.minimum(f::Function,x::BArray) + minimum(xi->minimum(f,xi),blocks(x)) +end + +function Base.collect(v::BVector) + reduce(vcat,map(collect,blocks(v))) +end + +function Base.:*(a::Number,b::BArray) + bs = map(bi->a*bi,blocks(b)) + BArray(bs) +end + +function Base.:*(b::BArray,a::Number) + a*b +end + +function Base.:/(b::BArray,a::Number) + (1/a)*b +end + +for op in (:+,:-) + @eval begin + function Base.$op(a::BArray) + bs = map($op,blocks(a)) + BArray(bs) + end + function Base.$op(a::BArray,b::BArray) + bs = map($op,blocks(a),blocks(b)) + BArray(bs) + end + end +end + +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::BArray) + reduce(+,a,init=zero(eltype(a))) +end + +function LinearAlgebra.dot(a::BVector,b::BVector) + c = map(dot,blocks(a),blocks(b)) + sum(c) +end + +function LinearAlgebra.rmul!(a::BArray,v::Number) + map(blocks(a)) do l + rmul!(l,v) + end + a +end + +function LinearAlgebra.norm(a::BVector,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::BVector,b::BVector) + s = distance_eval_body(d,a,b) + Distances.eval_end(d,s) + end + end +end + +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 + s = reduce((i,j)->Distances.eval_reduce(d,i,j), + partials, + init=Distances.eval_start(d, a, b)) + s +end + +struct BroadcastedBArray{A} + blocks::A +end +BlockArrays.blocks(a::BroadcastedBArray) = a.blocks + +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{BArray,BroadcastedBArray}) + map( bi -> Base.broadcasted(f,a,bi), blocks(b)) |> BroadcastedBArray +end + +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{BArray,BroadcastedBArray}, + 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{BArray,BroadcastedBArray}) + Base.broadcasted(f,Base.materialize(a),b) + end + +function Base.materialize(b::BroadcastedBArray) + map(Base.materialize,blocks(b)) |> BArray +end + +function Base.materialize!(a::BArray,b::BroadcastedBArray) + foreach(Base.materialize!,blocks(a),blocks(b)) + a +end + +function own_own_values(a::BMatrix) + ps = map(own_own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) +end +function own_ghost_values(a::BMatrix) + ps = map(own_ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) +end +function ghost_own_values(a::BMatrix) + ps = map(ghost_own_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) +end +function ghost_ghost_values(a::BMatrix) + ps = map(ghost_ghost_values,blocks(a)) + ps2 = permute_nesting(ps) + map(BMatrix,ps2) +end + +function LinearAlgebra.fillstored!(a::BMatrix,v) + foreach(ai->LinearAlgebra.fillstored!(ai,v),blocks(a)) + 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::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::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::BMatrix,b::BVector) + 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 LinearAlgebra.mul!(b::BVector,A::BMatrix,x::BVector) + mul!(b,A,x,1,0) +end + +function LinearAlgebra.mul!(b::BVector,A::BMatrix,x::BVector,α::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_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/p_vector.jl b/src/p_vector.jl index ea75e03c..80da6eac 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 @@ -775,6 +775,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) @@ -789,6 +793,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) @@ -1197,15 +1205,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) @@ -1216,31 +1224,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 @@ -1248,11 +1256,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) @@ -1260,7 +1268,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)) @@ -1271,38 +1279,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 @@ -1499,13 +1512,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 diff --git a/src/primitives.jl b/src/primitives.jl index 0c32f620..8fa98f0b 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -96,6 +96,26 @@ 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 + +function permute_nesting(a::AbstractArray{<:AbstractArray}) + s = size(a) + map(a...) do b... + reshape(collect(b),s) + 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 new file mode 100644 index 00000000..5909edc9 --- /dev/null +++ b/test/block_arrays_tests.jl @@ -0,0 +1,169 @@ +using PartitionedArrays +using Test +using LinearAlgebra +using Random +using Distances +using BlockArrays +using SparseArrays +using IterativeSolvers + +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)) + + 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 = 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) + @test isa(rows,BRange) + partition(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)) + 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,BVector) + 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,BVector) + w = v .+ 1 + @test isa(w,BVector) + w = v .+ w .- u + @test isa(w,BVector) + w = v .+ 1 .- u + @test isa(w,BVector) + 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 + + @test size(A11) == (16,16) + A = BMatrix(fill(A11,(2,2))) + display(A) + @show 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) + 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)) + @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 + @test isa(b,BVector) + B = 2*A + B = A*2 + B = +A + B = -A + 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 + 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 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