40#ifndef TPETRA_DETAILS_CRSUTILS_HPP
41#define TPETRA_DETAILS_CRSUTILS_HPP
45#include "TpetraCore_config.h"
46#include "Kokkos_Core.hpp"
48#include "Tpetra_Details_CrsPadding.hpp"
49#include "Tpetra_Details_WrappedDualView.hpp"
52#include <unordered_map>
65template<
class ViewType>
67make_uninitialized_view(
68 const std::string& name,
71 const std::string*
const prefix)
74 std::ostringstream os;
75 os << *prefix <<
"Allocate Kokkos::View " << name
76 <<
": " << size << std::endl;
77 std::cerr << os.str();
79 using Kokkos::view_alloc;
80 using Kokkos::WithoutInitializing;
81 return ViewType(view_alloc(name, WithoutInitializing), size);
84template<
class ViewType>
87 const std::string& name,
90 const std::string*
const prefix)
93 std::ostringstream os;
94 os << *prefix <<
"Allocate & initialize Kokkos::View "
95 << name <<
": " << size << std::endl;
96 std::cerr << os.str();
98 return ViewType(name, size);
101template<
class OutViewType,
class InViewType>
103assign_to_view(OutViewType& out,
104 const InViewType& in,
105 const char viewName[],
107 const std::string*
const prefix)
110 std::ostringstream os;
111 os << *prefix <<
"Assign to Kokkos::View " << viewName
112 <<
": Old size: " << out.extent(0)
113 <<
", New size: " << in.extent(0) << std::endl;
114 std::cerr << os.str();
119template<
class MemorySpace,
class ViewType>
120auto create_mirror_view(
121 const MemorySpace& memSpace,
122 const ViewType& view,
124 const std::string*
const prefix) ->
125 decltype(Kokkos::create_mirror_view(memSpace, view))
128 std::ostringstream os;
129 os << *prefix <<
"Create mirror view: "
130 <<
"view.extent(0): " << view.extent(0) << std::endl;
131 std::cerr << os.str();
133 return Kokkos::create_mirror_view(memSpace, view);
136enum class PadCrsAction {
149template<
class RowPtr,
class Indices,
class Values,
class Padding>
152 const PadCrsAction action,
153 const RowPtr& row_ptr_beg,
154 const RowPtr& row_ptr_end,
155 Indices& indices_wdv,
157 const Padding& padding,
161 using execution_space =
typename Indices::t_dev::execution_space;
162 using Kokkos::view_alloc;
163 using Kokkos::WithoutInitializing;
165 std::unique_ptr<std::string> prefix;
167 const size_t maxNumToPrint = verbose ?
170 std::ostringstream os;
171 os <<
"Proc " << my_rank <<
": Tpetra::...::pad_crs_arrays: ";
172 prefix = std::unique_ptr<std::string>(
new std::string(os.str()));
173 os <<
"Start" << endl;
174 std::cerr << os.str();
176 Kokkos::HostSpace hostSpace;
179 std::ostringstream os;
180 os << *prefix <<
"On input: ";
182 Kokkos::create_mirror_view(hostSpace, row_ptr_beg);
184 Kokkos::deep_copy(row_ptr_beg_h, row_ptr_beg);
189 Kokkos::create_mirror_view(hostSpace, row_ptr_end);
191 Kokkos::deep_copy(row_ptr_end_h, row_ptr_end);
194 os <<
", indices.extent(0): " << indices_wdv.extent(0)
195 <<
", values.extent(0): " << values_wdv.extent(0)
199 std::cerr << os.str();
202 if (row_ptr_beg.size() == 0) {
204 std::ostringstream os;
205 os << *prefix <<
"Done; local matrix has no rows" << endl;
206 std::cerr << os.str();
211 const size_t lclNumRows(row_ptr_beg.size() - 1);
212 RowPtr newAllocPerRow =
213 make_uninitialized_view<RowPtr>(
"newAllocPerRow", lclNumRows,
214 verbose, prefix.get());
216 std::ostringstream os;
217 os << *prefix <<
"Fill newAllocPerRow & compute increase" << endl;
218 std::cerr << os.str();
223 auto row_ptr_end_h = create_mirror_view(
224 hostSpace, row_ptr_end, verbose, prefix.get());
226 Kokkos::deep_copy(execution_space(), row_ptr_end_h, row_ptr_end);
227 auto row_ptr_beg_h = create_mirror_view(
228 hostSpace, row_ptr_beg, verbose, prefix.get());
230 Kokkos::deep_copy(execution_space(), row_ptr_beg_h, row_ptr_beg);
232 auto newAllocPerRow_h = create_mirror_view(
233 hostSpace, newAllocPerRow, verbose, prefix.get());
234 using host_range_type = Kokkos::RangePolicy<
235 Kokkos::DefaultHostExecutionSpace,
size_t>;
236 Kokkos::parallel_reduce
237 (
"Tpetra::CrsGraph: Compute new allocation size per row",
238 host_range_type(0, lclNumRows),
239 [&] (
const size_t lclRowInd,
size_t& lclIncrease) {
240 const size_t start = row_ptr_beg_h[lclRowInd];
241 const size_t end = row_ptr_beg_h[lclRowInd+1];
242 TEUCHOS_ASSERT( end >= start );
243 const size_t oldAllocSize = end - start;
244 const size_t oldNumEnt = row_ptr_end_h[lclRowInd] - start;
245 TEUCHOS_ASSERT( oldNumEnt <= oldAllocSize );
252 auto result = padding.get_result(lclRowInd);
253 const size_t newNumEnt = oldNumEnt + result.numInSrcNotInTgt;
254 if (newNumEnt > oldAllocSize) {
255 lclIncrease += (newNumEnt - oldAllocSize);
256 newAllocPerRow_h[lclRowInd] = newNumEnt;
259 newAllocPerRow_h[lclRowInd] = oldAllocSize;
264 std::ostringstream os;
265 os << *prefix <<
"increase: " << increase <<
", ";
269 std::cerr << os.str();
276 Kokkos::deep_copy(execution_space(), newAllocPerRow, newAllocPerRow_h);
279 using inds_value_type =
280 typename Indices::t_dev::non_const_value_type;
281 using vals_value_type =
typename Values::t_dev::non_const_value_type;
284 auto indices_old = indices_wdv.getDeviceView(Access::ReadOnly);
285 const size_t newIndsSize = size_t(indices_old.size()) + increase;
286 auto indices_new = make_uninitialized_view<typename Indices::t_dev>(
287 "Tpetra::CrsGraph column indices", newIndsSize, verbose,
290 typename Values::t_dev values_new;
291 auto values_old = values_wdv.getDeviceView(Access::ReadOnly);
292 if (action == PadCrsAction::INDICES_AND_VALUES) {
293 const size_t newValsSize = newIndsSize;
296 values_new = make_initialized_view<typename Values::t_dev>(
297 "Tpetra::CrsMatrix values", newValsSize, verbose, prefix.get());
301 std::ostringstream os;
302 os << *prefix <<
"Repack" << endl;
303 std::cerr << os.str();
306 using range_type = Kokkos::RangePolicy<execution_space, size_t>;
307 Kokkos::parallel_scan(
308 "Tpetra::CrsGraph or CrsMatrix repack",
309 range_type(
size_t(0),
size_t(lclNumRows+1)),
310 KOKKOS_LAMBDA (
const size_t lclRow,
size_t& newRowBeg,
311 const bool finalPass)
316 const size_t row_beg = row_ptr_beg[lclRow];
317 const size_t row_end =
318 lclRow < lclNumRows ? row_ptr_end[lclRow] : row_beg;
319 const size_t numEnt = row_end - row_beg;
320 const size_t newRowAllocSize =
321 lclRow < lclNumRows ? newAllocPerRow[lclRow] : size_t(0);
323 if (lclRow < lclNumRows) {
324 const Kokkos::pair<size_t, size_t> oldRange(
325 row_beg, row_beg + numEnt);
326 const Kokkos::pair<size_t, size_t> newRange(
327 newRowBeg, newRowBeg + numEnt);
328 auto oldColInds = Kokkos::subview(indices_old, oldRange);
329 auto newColInds = Kokkos::subview(indices_new, newRange);
332 memcpy(newColInds.data(), oldColInds.data(),
333 numEnt *
sizeof(inds_value_type));
334 if (action == PadCrsAction::INDICES_AND_VALUES) {
336 Kokkos::subview(values_old, oldRange);
337 auto newVals = Kokkos::subview(values_new, newRange);
338 memcpy(newVals.data(), oldVals.data(),
339 numEnt *
sizeof(vals_value_type));
343 row_ptr_beg[lclRow] = newRowBeg;
344 if (lclRow < lclNumRows) {
345 row_ptr_end[lclRow] = newRowBeg + numEnt;
348 newRowBeg += newRowAllocSize;
353 std::ostringstream os;
357 Kokkos::create_mirror_view(hostSpace, row_ptr_beg);
359 Kokkos::deep_copy(row_ptr_beg_h, row_ptr_beg);
366 Kokkos::create_mirror_view(hostSpace, row_ptr_end);
368 Kokkos::deep_copy(row_ptr_end_h, row_ptr_end);
373 std::cout << os.str();
376 indices_wdv = Indices(indices_new);
377 values_wdv = Values(values_new);
381 auto indices_h = indices_wdv.getHostView(Access::ReadOnly);
382 auto values_h = values_wdv.getHostView(Access::ReadOnly);
383 std::ostringstream os;
394 std::ostringstream os;
395 os << *prefix <<
"Done" << endl;
396 std::cerr << os.str();
401template <
class Po
inters,
class InOutIndices,
class InIndices,
class IndexMap>
404 typename Pointers::value_type
const row,
405 Pointers
const& row_ptrs,
406 InOutIndices& cur_indices,
407 size_t& num_assigned,
408 InIndices
const& new_indices,
410 std::function<
void(
size_t const,
size_t const,
size_t const)> cb)
412 if (new_indices.size() == 0) {
416 if (cur_indices.size() == 0) {
418 return Teuchos::OrdinalTraits<size_t>::invalid();
421 using offset_type =
typename std::decay<
decltype (row_ptrs[0])>::type;
422 using ordinal_type =
typename std::decay<
decltype (cur_indices[0])>::type;
424 const offset_type start = row_ptrs[row];
425 offset_type end = start +
static_cast<offset_type
> (num_assigned);
426 const size_t num_avail = (row_ptrs[row + 1] < end) ?
size_t (0) :
427 row_ptrs[row + 1] - end;
428 const size_t num_new_indices =
static_cast<size_t> (new_indices.size ());
429 size_t num_inserted = 0;
431 size_t numIndicesLookup = num_assigned + num_new_indices;
434 const size_t useLookUpTableThreshold = 400;
436 if (numIndicesLookup <= useLookUpTableThreshold || num_new_indices == 1) {
439 for (
size_t k = 0; k < num_new_indices; ++k) {
440 const ordinal_type idx = std::forward<IndexMap>(map)(new_indices[k]);
441 offset_type row_offset = start;
442 for (; row_offset < end; ++row_offset) {
443 if (idx == cur_indices[row_offset]) {
448 if (row_offset == end) {
449 if (num_inserted >= num_avail) {
450 return Teuchos::OrdinalTraits<size_t>::invalid();
453 cur_indices[end++] = idx;
457 cb(k, start, row_offset - start);
463 std::unordered_map<ordinal_type, offset_type> idxLookup(numIndicesLookup);
466 for (
size_t k = 0; k < num_assigned; k++) {
467 idxLookup[cur_indices[start+k]] = start+k;
471 for (
size_t k = 0; k < num_new_indices; k++) {
472 const ordinal_type idx = std::forward<IndexMap>(map)(new_indices[k]);
473 offset_type row_offset;
475 auto it = idxLookup.find(idx);
476 if (it == idxLookup.end()) {
477 if (num_inserted >= num_avail) {
478 return Teuchos::OrdinalTraits<size_t>::invalid();
482 cur_indices[end++] = idx;
483 idxLookup[idx] = row_offset;
488 row_offset = it->second;
491 cb(k, start, row_offset - start);
495 num_assigned += num_inserted;
500template <
class Po
inters,
class Indices1,
class Indices2,
class IndexMap,
class Callback>
503 typename Pointers::value_type
const row,
504 Pointers
const& row_ptrs,
505 const size_t curNumEntries,
506 Indices1
const& cur_indices,
507 Indices2
const& new_indices,
511 if (new_indices.size() == 0)
515 typename std::remove_const<typename Indices1::value_type>::type;
516 auto invalid_ordinal = Teuchos::OrdinalTraits<ordinal>::invalid();
518 const size_t start =
static_cast<size_t> (row_ptrs[row]);
519 const size_t end = start + curNumEntries;
520 size_t num_found = 0;
521 for (
size_t k = 0; k < new_indices.size(); k++)
523 auto row_offset = start;
524 auto idx = std::forward<IndexMap>(map)(new_indices[k]);
525 if (idx == invalid_ordinal)
527 for (; row_offset < end; row_offset++)
529 if (idx == cur_indices[row_offset])
531 std::forward<Callback>(cb)(k, start, row_offset - start);
559template<
class RowPtr,
class Indices,
class Padding>
562 const RowPtr& rowPtrBeg,
563 const RowPtr& rowPtrEnd,
564 Indices& indices_wdv,
565 const Padding& padding,
569 using impl::pad_crs_arrays;
572 pad_crs_arrays<RowPtr, Indices, Indices, Padding>(
573 impl::PadCrsAction::INDICES_ONLY, rowPtrBeg, rowPtrEnd,
574 indices_wdv, values_null, padding, my_rank, verbose);
577template<
class RowPtr,
class Indices,
class Values,
class Padding>
580 const RowPtr& rowPtrBeg,
581 const RowPtr& rowPtrEnd,
582 Indices& indices_wdv,
584 const Padding& padding,
588 using impl::pad_crs_arrays;
589 pad_crs_arrays<RowPtr, Indices, Values, Padding>(
590 impl::PadCrsAction::INDICES_AND_VALUES, rowPtrBeg, rowPtrEnd,
591 indices_wdv, values_wdv, padding, my_rank, verbose);
639template <
class Po
inters,
class InOutIndices,
class InIndices>
642 typename Pointers::value_type
const row,
643 Pointers
const& rowPtrs,
644 InOutIndices& curIndices,
646 InIndices
const& newIndices,
647 std::function<
void(
const size_t,
const size_t,
const size_t)> cb =
648 std::function<
void(
const size_t,
const size_t,
const size_t)>())
650 static_assert(std::is_same<typename std::remove_const<typename InOutIndices::value_type>::type,
651 typename std::remove_const<typename InIndices::value_type>::type>::value,
652 "Expected views to have same value type");
655 using ordinal =
typename InOutIndices::value_type;
656 auto numInserted = impl::insert_crs_indices(row, rowPtrs, curIndices,
657 numAssigned, newIndices, [](ordinal
const idx) {
return idx; }, cb);
661template <
class Po
inters,
class InOutIndices,
class InIndices>
664 typename Pointers::value_type
const row,
665 Pointers
const& rowPtrs,
666 InOutIndices& curIndices,
668 InIndices
const& newIndices,
669 std::function<
typename InOutIndices::value_type(
const typename InIndices::value_type)> map,
670 std::function<
void(
const size_t,
const size_t,
const size_t)> cb =
671 std::function<
void(
const size_t,
const size_t,
const size_t)>())
673 auto numInserted = impl::insert_crs_indices(row, rowPtrs, curIndices,
674 numAssigned, newIndices, map, cb);
708template <
class Po
inters,
class Indices1,
class Indices2,
class Callback>
711 typename Pointers::value_type
const row,
712 Pointers
const& rowPtrs,
713 const size_t curNumEntries,
714 Indices1
const& curIndices,
715 Indices2
const& newIndices,
718 static_assert(std::is_same<typename std::remove_const<typename Indices1::value_type>::type,
719 typename std::remove_const<typename Indices2::value_type>::type>::value,
720 "Expected views to have same value type");
722 using ordinal =
typename Indices2::value_type;
723 auto numFound = impl::find_crs_indices(row, rowPtrs, curNumEntries, curIndices, newIndices,
724 [=](ordinal ind){
return ind; }, cb);
728template <
class Po
inters,
class Indices1,
class Indices2,
class IndexMap,
class Callback>
731 typename Pointers::value_type
const row,
732 Pointers
const& rowPtrs,
733 const size_t curNumEntries,
734 Indices1
const& curIndices,
735 Indices2
const& newIndices,
739 return impl::find_crs_indices(row, rowPtrs, curNumEntries, curIndices, newIndices, map, cb);
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra's behavior.
void pad_crs_arrays(const PadCrsAction action, const RowPtr &row_ptr_beg, const RowPtr &row_ptr_end, Indices &indices_wdv, Values &values_wdv, const Padding &padding, const int my_rank, const bool verbose)
Implementation of padCrsArrays.
size_t find_crs_indices(typename Pointers::value_type const row, Pointers const &row_ptrs, const size_t curNumEntries, Indices1 const &cur_indices, Indices2 const &new_indices, IndexMap &&map, Callback &&cb)
Implementation of findCrsIndices.
size_t insert_crs_indices(typename Pointers::value_type const row, Pointers const &row_ptrs, InOutIndices &cur_indices, size_t &num_assigned, InIndices const &new_indices, IndexMap &&map, std::function< void(size_t const, size_t const, size_t const)> cb)
Implementation of insertCrsIndices.
static size_t verbosePrintCountThreshold()
Number of entries below which arrays, lists, etc. will be printed in debug mode.
Implementation details of Tpetra.
void padCrsArrays(const RowPtr &rowPtrBeg, const RowPtr &rowPtrEnd, Indices &indices_wdv, const Padding &padding, const int my_rank, const bool verbose)
Determine if the row pointers and indices arrays need to be resized to accommodate new entries....
void verbosePrintArray(std::ostream &out, const ArrayType &x, const char name[], const size_t maxNumToPrint)
Print min(x.size(), maxNumToPrint) entries of x.
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.
Namespace Tpetra contains the class and methods constituting the Tpetra library.