Added sequential implementations, removed some commented out old code, added all required SparseMatrixCSR functionality not provided by SparseMatricesCSR package, added Sequential tests

Author: jop611

src/PartitionedArrays.jl

@@ -187,4 +187,26 @@ export nullspace_linear_elasticity!
 export near_nullspace_linear_elasticity
 include("gallery.jl")
 
+export ascsc
+export ascsr
+export copy
+export similar
+export sparsecsr
+export strictly_equivalent
+export approx_equivalent
+export pointer_array
+export index_array
+include("sparse_helpers.jl")
+
+import Base: +,-,*,/
+export RAP
+export RAP!
+export RAPadd!
+export RAP_search!
+export +,-,*,/
+export mul!
+export mul_search!
+include("sequential_implementations.jl")
+
+
 end # module

src/p_sparse_matrix.jl

@@ -1357,177 +1357,6 @@ function psparse_assemble_impl(A,::Type,rows)
     error("Case not implemented yet")
 end
 
-# function psparse_assemble_impl(
-#         A,
-#         ::Type{<:AbstractSplitMatrix},
-#         rows;
-#         reuse=Val(false),
-#         assembly_neighbors_options_cols=(;))
-
-#     function setup_cache_snd(A,parts_snd,rows_sa,cols_sa)
-#         A_ghost_own   = A.blocks.ghost_own
-#         A_ghost_ghost = A.blocks.ghost_ghost
-#         gen = ( owner=>i for (i,owner) in enumerate(parts_snd) )
-#         owner_to_p = Dict(gen)
-#         ptrs = zeros(Int32,length(parts_snd)+1)
-#         ghost_to_owner_row = ghost_to_owner(rows_sa)
-#         ghost_to_global_row = ghost_to_global(rows_sa)
-#         own_to_global_col = own_to_global(cols_sa)
-#         ghost_to_global_col = ghost_to_global(cols_sa)
-#         for (i,_,_) in nziterator(A_ghost_own)
-#             owner = ghost_to_owner_row[i]
-#             ptrs[owner_to_p[owner]+1] += 1
-#         end
-#         for (i,_,_) in nziterator(A_ghost_ghost)
-#             owner = ghost_to_owner_row[i]
-#             ptrs[owner_to_p[owner]+1] += 1
-#         end
-#         length_to_ptrs!(ptrs)
-#         Tv = eltype(A_ghost_own)
-#         ndata = ptrs[end]-1
-#         I_snd_data = zeros(Int,ndata)
-#         J_snd_data = zeros(Int,ndata)
-#         V_snd_data = zeros(Tv,ndata)
-#         k_snd_data = zeros(Int32,ndata)
-#         nnz_ghost_own = 0
-#         for (k,(i,j,v)) in enumerate(nziterator(A_ghost_own))
-#             owner = ghost_to_owner_row[i]
-#             p = ptrs[owner_to_p[owner]]
-#             I_snd_data[p] = ghost_to_global_row[i]
-#             J_snd_data[p] = own_to_global_col[j]
-#             V_snd_data[p] = v
-#             k_snd_data[p] = k
-#             ptrs[owner_to_p[owner]] += 1
-#             nnz_ghost_own += 1
-#         end
-#         for (k,(i,j,v)) in enumerate(nziterator(A_ghost_ghost))
-#             owner = ghost_to_owner_row[i]
-#             p = ptrs[owner_to_p[owner]]
-#             I_snd_data[p] = ghost_to_global_row[i]
-#             J_snd_data[p] = ghost_to_global_col[j]
-#             V_snd_data[p] = v
-#             k_snd_data[p] = k+nnz_ghost_own
-#             ptrs[owner_to_p[owner]] += 1
-#         end
-#         rewind_ptrs!(ptrs)
-#         I_snd = JaggedArray(I_snd_data,ptrs)
-#         J_snd = JaggedArray(J_snd_data,ptrs)
-#         V_snd = JaggedArray(V_snd_data,ptrs)
-#         k_snd = JaggedArray(k_snd_data,ptrs)
-#         (;I_snd,J_snd,V_snd,k_snd,parts_snd)
-#     end
-#     function setup_cache_rcv(I_rcv,J_rcv,V_rcv,parts_rcv)
-#         k_rcv_data = zeros(Int32,length(I_rcv.data))
-#         k_rcv = JaggedArray(k_rcv_data,I_rcv.ptrs)
-#         (;I_rcv,J_rcv,V_rcv,k_rcv,parts_rcv)
-#     end
-#     function setup_own_triplets(A,cache_rcv,rows_sa,cols_sa)
-#         nz_own_own = findnz(A.blocks.own_own)
-#         nz_own_ghost = findnz(A.blocks.own_ghost)
-#         I_rcv_data = cache_rcv.I_rcv.data
-#         J_rcv_data = cache_rcv.J_rcv.data
-#         V_rcv_data = cache_rcv.V_rcv.data
-#         k_rcv_data = cache_rcv.k_rcv.data
-#         global_to_own_col = global_to_own(cols_sa)
-#         is_ghost = findall(j->global_to_own_col[j]==0,J_rcv_data)
-#         is_own = findall(j->global_to_own_col[j]!=0,J_rcv_data)
-#         I_rcv_own = view(I_rcv_data,is_own)
-#         J_rcv_own = view(J_rcv_data,is_own)
-#         V_rcv_own = view(V_rcv_data,is_own)
-#         k_rcv_own = view(k_rcv_data,is_own)
-#         I_rcv_ghost = view(I_rcv_data,is_ghost)
-#         J_rcv_ghost = view(J_rcv_data,is_ghost)
-#         V_rcv_ghost = view(V_rcv_data,is_ghost)
-#         k_rcv_ghost = view(k_rcv_data,is_ghost)
-#         # After this col ids in own_ghost triplet remain global
-#         map_global_to_own!(I_rcv_own,rows_sa)
-#         map_global_to_own!(J_rcv_own,cols_sa)
-#         map_global_to_own!(I_rcv_ghost,rows_sa)
-#         map_ghost_to_global!(nz_own_ghost[2],cols_sa)
-#         own_own_I = vcat(nz_own_own[1],I_rcv_own)
-#         own_own_J = vcat(nz_own_own[2],J_rcv_own)
-#         own_own_V = vcat(nz_own_own[3],V_rcv_own)
-#         own_own_triplet = (own_own_I,own_own_J,own_own_V)
-#         own_ghost_I = vcat(nz_own_ghost[1],I_rcv_ghost)
-#         own_ghost_J = vcat(nz_own_ghost[2],J_rcv_ghost)
-#         own_ghost_V = vcat(nz_own_ghost[3],V_rcv_ghost)
-#         map_global_to_ghost!(nz_own_ghost[2],cols_sa)
-#         own_ghost_triplet = (own_ghost_I,own_ghost_J,own_ghost_V)
-#         triplets = (own_own_triplet,own_ghost_triplet)
-#         aux = (I_rcv_own,J_rcv_own,k_rcv_own,I_rcv_ghost,J_rcv_ghost,k_rcv_ghost,nz_own_own,nz_own_ghost)
-#         triplets, own_ghost_J, aux
-#     end
-#     function finalize_values(A,rows_fa,cols_fa,cache_snd,cache_rcv,triplets,aux)
-#         (own_own_triplet,own_ghost_triplet) = triplets
-#         (I_rcv_own,J_rcv_own,k_rcv_own,I_rcv_ghost,J_rcv_ghost,k_rcv_ghost,nz_own_own,nz_own_ghost) = aux
-#         map_global_to_ghost!(own_ghost_triplet[2],cols_fa)
-#         map_global_to_ghost!(J_rcv_ghost,cols_fa)
-#         TA = typeof(A.blocks.own_own)
-#         n_own_rows = own_length(rows_fa)
-#         n_own_cols = own_length(cols_fa)
-#         n_ghost_rows = ghost_length(rows_fa)
-#         n_ghost_cols = ghost_length(cols_fa)
-#         Ti = indextype(A.blocks.own_own)
-#         Tv = eltype(A.blocks.own_own)
-#         own_own = compresscoo(TA,own_own_triplet...,n_own_rows,n_own_cols)
-#         own_ghost = compresscoo(TA,own_ghost_triplet...,n_own_rows,n_ghost_cols)
-#         ghost_own = compresscoo(TA,Ti[],Ti[],Tv[],n_ghost_rows,n_own_cols)
-#         ghost_ghost = compresscoo(TA,Ti[],Ti[],Tv[],n_ghost_rows,n_ghost_cols)
-#         blocks = split_matrix_blocks(own_own,own_ghost,ghost_own,ghost_ghost)
-#         values = split_matrix(blocks,local_permutation(rows_fa),local_permutation(cols_fa))
-#         nnz_own_own = nnz(own_own)
-#         k_own_sa = precompute_nzindex(own_own,own_own_triplet[1:2]...)
-#         k_ghost_sa = precompute_nzindex(own_ghost,own_ghost_triplet[1:2]...)
-#         for p in 1:length(I_rcv_own)
-#             i = I_rcv_own[p]
-#             j = J_rcv_own[p]
-#             k_rcv_own[p] = nzindex(own_own,i,j)
-#         end
-#         for p in 1:length(I_rcv_ghost)
-#             i = I_rcv_ghost[p]
-#             j = J_rcv_ghost[p]
-#             k_rcv_ghost[p] = nzindex(own_ghost,i,j) + nnz_own_own
-#         end
-#         cache = (;k_own_sa,k_ghost_sa,cache_snd...,cache_rcv...)
-#         values, cache
-#     end
-#     rows_sa = partition(axes(A,1))
-#     cols_sa = partition(axes(A,2))
-#     #rows = map(remove_ghost,rows_sa)
-#     cols = map(remove_ghost,cols_sa)
-#     parts_snd, parts_rcv = assembly_neighbors(rows_sa)
-#     cache_snd = map(setup_cache_snd,partition(A),parts_snd,rows_sa,cols_sa)
-#     I_snd = map(i->i.I_snd,cache_snd)
-#     J_snd = map(i->i.J_snd,cache_snd)
-#     V_snd = map(i->i.V_snd,cache_snd)
-#     graph = ExchangeGraph(parts_snd,parts_rcv)
-#     t_I = exchange(I_snd,graph)
-#     t_J = exchange(J_snd,graph)
-#     t_V = exchange(V_snd,graph)
-#     @fake_async begin
-#         I_rcv = fetch(t_I)
-#         J_rcv = fetch(t_J)
-#         V_rcv = fetch(t_V)
-#         cache_rcv = map(setup_cache_rcv,I_rcv,J_rcv,V_rcv,parts_rcv)
-#         triplets,J,aux = map(setup_own_triplets,partition(A),cache_rcv,rows_sa,cols_sa) |> tuple_of_arrays
-#         J_owner = find_owner(cols_sa,J)
-#         rows_fa = rows
-#         cols_fa = map(union_ghost,cols,J,J_owner)
-#         assembly_neighbors(cols_fa;assembly_neighbors_options_cols...)
-#         vals_fa, cache = map(finalize_values,partition(A),rows_fa,cols_fa,cache_snd,cache_rcv,triplets,aux) |> tuple_of_arrays
-#         assembled = true
-#         B = PSparseMatrix(vals_fa,rows_fa,cols_fa,assembled)
-#         if val_parameter(reuse) == false
-#             B
-#         else
-#             B, cache
-#         end
-#     end
-# end
-
-# New assemble
-####################
-
 function psparse_assemble_impl(A::PSparseMatrix{V,B,C,D,Tv} where {V,B,C,D},
                                ::Type{T},
                                rows;