fverdugo / PartitionedArrays.jl

Large-scale, distributed, sparse linear algebra in Julia.
MIT License
119 stars 14 forks source link

Bug: PermutedLocalIndices #124

Open JordiManyer opened 10 months ago

JordiManyer commented 10 months ago

@fverdugo I think I found some issues for PermutedLocalIndices. See the MWE:

using PartitionedArrays

ranks = with_debug() do distribute
  distribute(LinearIndices((2,)))
end

indices = map(ranks) do r
  n_global = 4
  if r == 1
    local_to_global = [1,2,3]
    local_to_owner  = [1,1,2]
  else
    local_to_global = [2,3,4]
    local_to_owner  = [1,2,2]
  end
  LocalIndices(n_global,r,local_to_global,local_to_owner)
end

perm = map(ranks) do r
  if r == 1
    [2,1,3]
  else
    [1,3,2]
  end
end

perm_indices = map(permute_indices,indices,perm) # Error 1
v = pfill(0.0,perm_indices) # Error 2

Trace for error 1:

ERROR: MethodError: no method matching PartitionedArrays.LocalToGlobal(::SubArray{Int64, 1, Vector{Int64}, Tuple{Vector{Int32}}, false}, ::SubArray{Int64, 1, Vector{Int64}, Tuple{Vector{Int32}}, false}, ::Vector{Int32})

Closest candidates are:
  PartitionedArrays.LocalToGlobal(::A, ::Vector{Int64}, ::C) where {A, C}
   @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/p_range.jl:955
  PartitionedArrays.LocalToGlobal(::PartitionedArrays.BlockPartitionOwnToGlobal, ::Vector{Int64}, ::Any)
   @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/p_range.jl:1428

Trace for error 2:

ERROR: MethodError: no method matching PartitionedArrays.LocalToOwner(::PartitionedArrays.OwnToOwner, ::SubArray{Int32, 1, Vector{Int32}, Tuple{Vector{Int32}}, false}, ::Vector{Int32})

Closest candidates are:
  PartitionedArrays.LocalToOwner(::PartitionedArrays.OwnToOwner, ::Vector{Int32}, ::C) where C
   @ PartitionedArrays ~/.julia/packages/PartitionedArrays/py6uo/src/p_range.jl:972

It seems the LocalToGlobal and LocalToOwner structs are enforcing datatypes which PermutedLocalIndices does not support, when created from a LocalIndices object.


I would also think that this is related to the fact that LocalToGlobal assumes we have our indices ordered as [own...,ghost...] like for OwnAndGhostIndices (which is the only type for which permutation is supported). Indeed, if we solve the issue with something like

function local_to_global(a::PermutedLocalIndices)
  LocalToGlobal(own_to_global(a),collect(ghost_to_global(a)),a.perm)
end

we can create the object, but the result is wrong:

2-element DebugArray{PartitionedArrays.LocalToGlobal{SubArray{Int64, 1, Vector{Int64}, Tuple{Vector{Int32}}, false}, Vector{Int32}}, 1}:
[1] = [2, 1, 3]
[2] = [3, 4, 2] # This should be still [2,3,4], but it gets reordered as owned then ghost.

All in all, I think something along these lines should be considered to fix this issue while preserving performance and the current API:

Alternatively, we could avoid a new structure by just dispatching to

function permute_indices(indices::LocalIndices,perm)
  id = part_id(indices)
  n_glob = global_length(indices)
  l2g = view(local_to_global(indices),perm)
  l2o = view(local_to_owner(indices),perm)
  return LocalIndices(n_glob,id,l2g,l2o)
end

Also, and just as a comment, I think it would be nice to add to the documentation in which direction the permutation is going, i.e from the old local indices to the new or vice-versa. I believe it;s meant to be a map new -> old.

fverdugo commented 10 months ago

Hi @JordiManyer I will take a look hopefully next week.

If you have a patch feel free to open a PR.