 
         alloc.deallocate((char *)thrust::raw_pointer_cast(out_idx_ptr),0);
         alloc.deallocate((char *)thrust::raw_pointer_cast(n_ptr),0);
+
+        /*
+         * sort by neighbor and compute start and end indices
+         */
+        thrust::sort_by_key(thrust::cuda::par(alloc),
+            out_neighbor_ptr,
+            out_neighbor_ptr + n_out,
+            ghost_tag_ptr);
+
+
         #else
         gpu_expand_neighbors(n_out,
             thrust::raw_pointer_cast(offsets_ptr),

        unique_neighbors_ptr,
        unique_neighbors_ptr + n_unique_neigh,
        ghost_end_ptr);

    // deallocate temporary arrays
    alloc.deallocate((char *)thrust::raw_pointer_cast(out_neighbor_ptr),0);
    }

        alloc.deallocate((char *)thrust::raw_pointer_cast(plan_out_ptr),0);
        }

    #if 0
    thrust::lower_bound(thrust::cuda::par(alloc),
        out_neighbor_ptr,
        out_neighbor_ptr + n_out,
