42 #ifndef TPETRA_DETAILS_CRSUTILS_HPP
43 #define TPETRA_DETAILS_CRSUTILS_HPP
46 #include <type_traits>
48 #include "TpetraCore_config.h"
49 #include "Kokkos_Core.hpp"
50 #include "Kokkos_UnorderedMap.hpp"
63 template<
class ViewType>
65 uninitialized_view(
const std::string& name,
const size_t& size)
67 return ViewType (Kokkos::view_alloc(name, Kokkos::WithoutInitializing), size);
71 template<
class RowPtr,
class Indices,
class Values,
class Padding>
74 const RowPtr& row_ptr_beg,
75 const RowPtr& row_ptr_end,
78 const Padding& padding)
81 if (padding.size() == 0 || row_ptr_beg.size() == 0) {
86 auto pad_values = values.extent(0) == indices.extent(0);
89 auto num_row = row_ptr_beg.size() - 1;
90 auto policy = Kokkos::RangePolicy<typename Padding::execution_space>(0, num_row);
91 RowPtr entries_this_row(
"entries_this_row", num_row);
93 size_t additional_size_needed = 0;
94 Kokkos::parallel_reduce(
"Determine additional size needed", policy,
95 KOKKOS_LAMBDA(
const int i,
size_t& ladditional_size_needed) {
97 auto allocated_this_row = row_ptr_beg(i+1) - row_ptr_beg(i);
98 auto used_this_row = row_ptr_end(i) - row_ptr_beg(i);
99 auto free_this_row = allocated_this_row - used_this_row;
100 entries_this_row(i) = allocated_this_row;
102 auto k = padding.find(static_cast<typename Padding::key_type>(i));
103 if (padding.valid_at(k)) {
105 auto num_ent = padding.value_at(k);
106 auto n = (num_ent > free_this_row) ? num_ent - free_this_row : 0;
107 entries_this_row(i) += n;
108 ladditional_size_needed += n;
110 }, additional_size_needed);
112 if (additional_size_needed == 0)
115 using ptrs_value_type =
typename RowPtr::non_const_value_type;
116 using inds_value_type =
typename Indices::non_const_value_type;
117 using vals_value_type =
typename Values::non_const_value_type;
120 auto indices_new = uninitialized_view<Indices>(
"ind new", indices.size()+additional_size_needed);
121 auto values_new = uninitialized_view<Values>(
"val new", pad_values ? values.size()+additional_size_needed : 0);
128 auto this_row_beg = row_ptr_beg(0);
129 auto this_row_end = row_ptr_end(0);
130 using range = Kokkos::pair<ptrs_value_type, ptrs_value_type>;
131 for (
typename RowPtr::size_type i=0; i<num_row-1; i++) {
133 auto used_this_row = this_row_end - this_row_beg;
137 auto indices_old_subview = subview(indices, range(this_row_beg, this_row_beg+used_this_row));
138 auto indices_new_subview = subview(indices_new, range(row_ptr_beg(i), row_ptr_beg(i)+used_this_row));
140 memcpy(indices_new_subview.data(), indices_old_subview.data(), used_this_row *
sizeof(inds_value_type));
145 auto values_old_subview = subview(values, range(this_row_beg, this_row_beg+used_this_row));
146 auto values_new_subview = subview(values_new, range(row_ptr_beg(i), row_ptr_beg(i)+used_this_row));
147 memcpy(values_new_subview.data(), values_old_subview.data(), used_this_row *
sizeof(vals_value_type));
151 this_row_beg = row_ptr_beg(i+1);
152 this_row_end = row_ptr_end(i+1);
154 auto used_next_row = row_ptr_end(i+1) - row_ptr_beg(i+1);
155 row_ptr_beg(i+1) = row_ptr_beg(i) + entries_this_row(i);
156 row_ptr_end(i+1) = row_ptr_beg(i+1) + used_next_row;
161 row_ptr_beg(num_row) = indices_new.size();
162 auto n = num_row - 1;
163 auto used_this_row = row_ptr_end(n) - row_ptr_beg(n);
166 auto indices_old_subview = subview(indices, range(this_row_beg, this_row_beg+used_this_row));
167 auto indices_new_subview = subview(indices_new, range(row_ptr_beg(n), row_ptr_beg(n)+used_this_row));
168 memcpy(indices_new_subview.data(), indices_old_subview.data(), used_this_row *
sizeof(inds_value_type));
172 auto values_old_subview = subview(values, range(this_row_beg, this_row_beg+used_this_row));
173 auto values_new_subview = subview(values_new, range(row_ptr_beg(n), row_ptr_beg(n)+used_this_row));
174 memcpy(values_new_subview.data(), values_old_subview.data(), used_this_row *
sizeof(vals_value_type));
178 indices = indices_new;
183 template <
class Po
inters,
class InOutIndices,
class InIndices,
class IndexMap>
186 typename Pointers::value_type
const row,
187 Pointers
const& row_ptrs,
188 InOutIndices& cur_indices,
189 size_t& num_assigned,
190 InIndices
const& new_indices,
192 std::function<
void(
size_t const,
size_t const,
size_t const)> cb)
194 if (new_indices.size() == 0) {
198 if (cur_indices.size() == 0) {
200 return Teuchos::OrdinalTraits<size_t>::invalid();
203 using offset_type =
typename std::decay<decltype (row_ptrs[0])>::type;
204 using ordinal_type =
typename std::decay<decltype (cur_indices[0])>::type;
206 const offset_type start = row_ptrs[row];
207 offset_type end = start +
static_cast<offset_type
> (num_assigned);
208 const size_t num_avail = (row_ptrs[row + 1] < end) ?
size_t (0) :
209 row_ptrs[row + 1] - end;
210 const size_t num_new_indices =
static_cast<size_t> (new_indices.size ());
211 size_t num_inserted = 0;
213 for (
size_t k = 0; k < num_new_indices; ++k) {
214 const ordinal_type idx = std::forward<IndexMap>(map)(new_indices[k]);
215 offset_type row_offset = start;
216 for (; row_offset < end; ++row_offset) {
217 if (idx == cur_indices[row_offset]) {
222 if (row_offset == end) {
223 if (num_inserted >= num_avail) {
224 return Teuchos::OrdinalTraits<size_t>::invalid();
227 cur_indices[end++] = idx;
231 cb(k, start, row_offset - start);
234 num_assigned += num_inserted;
239 template <
class Po
inters,
class Indices1,
class Indices2,
class IndexMap,
class Callback>
242 typename Pointers::value_type
const row,
243 Pointers
const& row_ptrs,
244 const size_t curNumEntries,
245 Indices1
const& cur_indices,
246 Indices2
const& new_indices,
250 if (new_indices.size() == 0)
253 using ordinal =
typename Indices1::value_type;
254 auto invalid_ordinal = Teuchos::OrdinalTraits<ordinal>::invalid();
256 const size_t start =
static_cast<size_t> (row_ptrs[row]);
257 const size_t end = start + curNumEntries;
258 size_t num_found = 0;
259 for (
size_t k = 0; k < new_indices.size(); k++)
261 auto row_offset = start;
262 auto idx = std::forward<IndexMap>(map)(new_indices[k]);
263 if (idx == invalid_ordinal)
265 for (; row_offset < end; row_offset++)
267 if (idx == cur_indices[row_offset])
269 std::forward<Callback>(cb)(k, start, row_offset - start);
301 template<
class RowPtr,
class Indices,
class Padding>
304 const RowPtr& rowPtrBeg,
305 const RowPtr& rowPtrEnd,
307 const Padding& padding)
309 using impl::pad_crs_arrays;
312 pad_crs_arrays<RowPtr, Indices, Indices, Padding>(rowPtrBeg, rowPtrEnd, indices, values, padding);
315 template<
class RowPtr,
class Indices,
class Values,
class Padding>
318 const RowPtr& rowPtrBeg,
319 const RowPtr& rowPtrEnd,
322 const Padding& padding)
324 using impl::pad_crs_arrays;
325 pad_crs_arrays<RowPtr, Indices, Values, Padding>(rowPtrBeg, rowPtrEnd, indices, values, padding);
373 template <
class Po
inters,
class InOutIndices,
class InIndices>
376 typename Pointers::value_type
const row,
377 Pointers
const& rowPtrs,
378 InOutIndices& curIndices,
380 InIndices
const& newIndices,
381 std::function<
void(
const size_t,
const size_t,
const size_t)> cb =
382 std::function<
void(
const size_t,
const size_t,
const size_t)>())
384 static_assert(std::is_same<
typename std::remove_const<typename InOutIndices::value_type>::type,
385 typename std::remove_const<typename InIndices::value_type>::type>::value,
386 "Expected views to have same value type");
389 using ordinal =
typename InOutIndices::value_type;
390 auto numInserted = impl::insert_crs_indices(row, rowPtrs, curIndices,
391 numAssigned, newIndices, [](ordinal
const idx) {
return idx; }, cb);
395 template <
class Po
inters,
class InOutIndices,
class InIndices>
398 typename Pointers::value_type
const row,
399 Pointers
const& rowPtrs,
400 InOutIndices& curIndices,
402 InIndices
const& newIndices,
403 std::function<
typename InOutIndices::value_type(
const typename InIndices::value_type)> map,
404 std::function<
void(
const size_t,
const size_t,
const size_t)> cb =
405 std::function<
void(
const size_t,
const size_t,
const size_t)>())
407 auto numInserted = impl::insert_crs_indices(row, rowPtrs, curIndices,
408 numAssigned, newIndices, map, cb);
442 template <
class Po
inters,
class Indices1,
class Indices2,
class Callback>
445 typename Pointers::value_type
const row,
446 Pointers
const& rowPtrs,
447 const size_t curNumEntries,
448 Indices1
const& curIndices,
449 Indices2
const& newIndices,
452 static_assert(std::is_same<
typename std::remove_const<typename Indices1::value_type>::type,
453 typename std::remove_const<typename Indices2::value_type>::type>::value,
454 "Expected views to have same value type");
456 using ordinal =
typename Indices2::value_type;
457 auto numFound = impl::find_crs_indices(row, rowPtrs, curNumEntries, curIndices, newIndices,
458 [=](ordinal ind){
return ind; }, cb);
462 template <
class Po
inters,
class Indices1,
class Indices2,
class IndexMap,
class Callback>
465 typename Pointers::value_type
const row,
466 Pointers
const& rowPtrs,
467 const size_t curNumEntries,
468 Indices1
const& curIndices,
469 Indices2
const& newIndices,
473 return impl::find_crs_indices(row, rowPtrs, curNumEntries, curIndices, newIndices, map, cb);
479 #endif // TPETRA_DETAILS_CRSUTILS_HPP
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
Copy the contents of the MultiVector src into dst.
size_t insertCrsIndices(typename Pointers::value_type const row, Pointers const &rowPtrs, InOutIndices &curIndices, size_t &numAssigned, InIndices const &newIndices, std::function< void(const size_t, const size_t, const size_t)> cb=std::function< void(const size_t, const size_t, const size_t)>())
Insert new indices in to current list of indices.
size_t findCrsIndices(typename Pointers::value_type const row, Pointers const &rowPtrs, const size_t curNumEntries, Indices1 const &curIndices, Indices2 const &newIndices, Callback &&cb)
Finds offsets in to current list of indices.
void padCrsArrays(const RowPtr &rowPtrBeg, const RowPtr &rowPtrEnd, Indices &indices, const Padding &padding)
Determine if the row pointers and indices arrays need to be resized to accommodate new entries...