Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Details_packCrsMatrix_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Tpetra: Templated Linear Algebra Services Package
4 //
5 // Copyright 2008 NTESS and the Tpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef TPETRA_DETAILS_PACKCRSMATRIX_DEF_HPP
11 #define TPETRA_DETAILS_PACKCRSMATRIX_DEF_HPP
12 
13 #include "TpetraCore_config.h"
14 #include "Teuchos_Array.hpp"
15 #include "Teuchos_ArrayView.hpp"
24 #include <memory>
25 #include <sstream>
26 #include <stdexcept>
27 #include <string>
28 
51 
52 namespace Tpetra {
53 
54 //
55 // Users must never rely on anything in the Details namespace.
56 //
57 namespace Details {
58 
59 namespace PackCrsMatrixImpl {
67 template<class OutputOffsetsViewType,
68  class CountsViewType,
69  class InputOffsetsViewType,
70  class InputLocalRowIndicesViewType,
71  class InputLocalRowPidsViewType,
72  const bool debug =
73 #ifdef HAVE_TPETRA_DEBUG
74  true
75 #else
76  false
77 #endif // HAVE_TPETRA_DEBUG
78  >
80 public:
81  typedef typename OutputOffsetsViewType::non_const_value_type output_offset_type;
82  typedef typename CountsViewType::non_const_value_type count_type;
83  typedef typename InputOffsetsViewType::non_const_value_type input_offset_type;
84  typedef typename InputLocalRowIndicesViewType::non_const_value_type local_row_index_type;
85  typedef typename InputLocalRowPidsViewType::non_const_value_type local_row_pid_type;
86  // output Views drive where execution happens.
87  typedef typename OutputOffsetsViewType::device_type device_type;
88  static_assert (std::is_same<typename CountsViewType::device_type::execution_space,
89  typename device_type::execution_space>::value,
90  "OutputOffsetsViewType and CountsViewType must have the same execution space.");
91  static_assert (Kokkos::is_view<OutputOffsetsViewType>::value,
92  "OutputOffsetsViewType must be a Kokkos::View.");
93  static_assert (std::is_same<typename OutputOffsetsViewType::value_type, output_offset_type>::value,
94  "OutputOffsetsViewType must be a nonconst Kokkos::View.");
95  static_assert (std::is_integral<output_offset_type>::value,
96  "The type of each entry of OutputOffsetsViewType must be a built-in integer type.");
97  static_assert (Kokkos::is_view<CountsViewType>::value,
98  "CountsViewType must be a Kokkos::View.");
99  static_assert (std::is_same<typename CountsViewType::value_type, output_offset_type>::value,
100  "CountsViewType must be a nonconst Kokkos::View.");
101  static_assert (std::is_integral<count_type>::value,
102  "The type of each entry of CountsViewType must be a built-in integer type.");
103  static_assert (Kokkos::is_view<InputOffsetsViewType>::value,
104  "InputOffsetsViewType must be a Kokkos::View.");
105  static_assert (std::is_integral<input_offset_type>::value,
106  "The type of each entry of InputOffsetsViewType must be a built-in integer type.");
107  static_assert (Kokkos::is_view<InputLocalRowIndicesViewType>::value,
108  "InputLocalRowIndicesViewType must be a Kokkos::View.");
109  static_assert (std::is_integral<local_row_index_type>::value,
110  "The type of each entry of InputLocalRowIndicesViewType must be a built-in integer type.");
111 
112  NumPacketsAndOffsetsFunctor (const OutputOffsetsViewType& outputOffsets,
113  const CountsViewType& counts,
114  const InputOffsetsViewType& rowOffsets,
115  const InputLocalRowIndicesViewType& lclRowInds,
116  const InputLocalRowPidsViewType& lclRowPids,
117  const count_type sizeOfLclCount,
118  const count_type sizeOfGblColInd,
119  const count_type sizeOfPid,
120  const count_type sizeOfValue) :
121  outputOffsets_ (outputOffsets),
122  counts_ (counts),
123  rowOffsets_ (rowOffsets),
124  lclRowInds_ (lclRowInds),
125  lclRowPids_ (lclRowPids),
126  sizeOfLclCount_ (sizeOfLclCount),
127  sizeOfGblColInd_ (sizeOfGblColInd),
128  sizeOfPid_ (sizeOfPid),
129  sizeOfValue_ (sizeOfValue),
130  error_ ("error") // don't forget this, or you'll get segfaults!
131  {
132  if (debug) {
133  const size_t numRowsToPack = static_cast<size_t> (lclRowInds_.extent (0));
134 
135  if (numRowsToPack != static_cast<size_t> (counts_.extent (0))) {
136  std::ostringstream os;
137  os << "lclRowInds.extent(0) = " << numRowsToPack
138  << " != counts.extent(0) = " << counts_.extent (0)
139  << ".";
140  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str ());
141  }
142  if (static_cast<size_t> (numRowsToPack + 1) !=
143  static_cast<size_t> (outputOffsets_.extent (0))) {
144  std::ostringstream os;
145  os << "lclRowInds.extent(0) + 1 = " << (numRowsToPack + 1)
146  << " != outputOffsets.extent(0) = " << outputOffsets_.extent (0)
147  << ".";
148  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str ());
149  }
150  }
151  }
152 
153  KOKKOS_INLINE_FUNCTION void
154  operator() (const local_row_index_type& curInd,
155  output_offset_type& update,
156  const bool final) const
157  {
158  if (debug) {
159  if (curInd < static_cast<local_row_index_type> (0)) {
160  error_ () = 1;
161  return;
162  }
163  }
164 
165  if (final) {
166  if (debug) {
167  if (curInd >= static_cast<local_row_index_type> (outputOffsets_.extent (0))) {
168  error_ () = 2;
169  return;
170  }
171  }
172  outputOffsets_(curInd) = update;
173  }
174 
175  if (curInd < static_cast<local_row_index_type> (counts_.extent (0))) {
176  const auto lclRow = lclRowInds_(curInd);
177  if (static_cast<size_t> (lclRow + 1) >= static_cast<size_t> (rowOffsets_.extent (0)) ||
178  static_cast<local_row_index_type> (lclRow) < static_cast<local_row_index_type> (0)) {
179  error_ () = 3;
180  return;
181  }
182  // count_type could differ from the type of each row offset.
183  // For example, row offsets might each be 64 bits, but if their
184  // difference always fits in 32 bits, we may then safely use a
185  // 32-bit count_type.
186  const count_type count =
187  static_cast<count_type> (rowOffsets_(lclRow+1) - rowOffsets_(lclRow));
188 
189  // We pack first the number of entries in the row, then that
190  // many global column indices, then that many pids (if any),
191  // then that many values. However, if the number of entries in
192  // the row is zero, we pack nothing.
193  const count_type numBytes = (count == 0) ?
194  static_cast<count_type> (0) :
195  sizeOfLclCount_ + count * (sizeOfGblColInd_ +
196  (lclRowPids_.size() > 0 ? sizeOfPid_ : 0) +
197  sizeOfValue_);
198 
199  if (final) {
200  counts_(curInd) = numBytes;
201  }
202  update += numBytes;
203  }
204  }
205 
206  // mfh 31 May 2017: Don't need init or join. If you have join, MUST
207  // have join both with and without volatile! Otherwise intrawarp
208  // joins are really slow on GPUs.
209 
211  int getError () const {
212  typedef typename device_type::execution_space execution_space;
213  auto error_h = Kokkos::create_mirror_view (error_);
214  // DEEP_COPY REVIEW - DEVICE-TO-HOSTMIRROR
215  Kokkos::deep_copy (execution_space(), error_h, error_);
216  return error_h ();
217  }
218 
219 private:
220  OutputOffsetsViewType outputOffsets_;
221  CountsViewType counts_;
222  typename InputOffsetsViewType::const_type rowOffsets_;
223  typename InputLocalRowIndicesViewType::const_type lclRowInds_;
224  typename InputLocalRowPidsViewType::const_type lclRowPids_;
225  count_type sizeOfLclCount_;
226  count_type sizeOfGblColInd_;
227  count_type sizeOfPid_;
228  count_type sizeOfValue_;
229  Kokkos::View<int, device_type> error_;
230 };
231 
241 template<class OutputOffsetsViewType,
242  class CountsViewType,
243  class InputOffsetsViewType,
244  class InputLocalRowIndicesViewType,
245  class InputLocalRowPidsViewType>
246 typename CountsViewType::non_const_value_type
247 computeNumPacketsAndOffsets (const OutputOffsetsViewType& outputOffsets,
248  const CountsViewType& counts,
249  const InputOffsetsViewType& rowOffsets,
250  const InputLocalRowIndicesViewType& lclRowInds,
251  const InputLocalRowPidsViewType& lclRowPids,
252  const typename CountsViewType::non_const_value_type sizeOfLclCount,
253  const typename CountsViewType::non_const_value_type sizeOfGblColInd,
254  const typename CountsViewType::non_const_value_type sizeOfPid,
255  const typename CountsViewType::non_const_value_type sizeOfValue)
256 {
257  typedef NumPacketsAndOffsetsFunctor<OutputOffsetsViewType,
258  CountsViewType, typename InputOffsetsViewType::const_type,
259  typename InputLocalRowIndicesViewType::const_type,
260  typename InputLocalRowPidsViewType::const_type> functor_type;
261  typedef typename CountsViewType::non_const_value_type count_type;
262  typedef typename OutputOffsetsViewType::size_type size_type;
263  typedef typename OutputOffsetsViewType::execution_space execution_space;
264  typedef typename functor_type::local_row_index_type LO;
265  typedef Kokkos::RangePolicy<execution_space, LO> range_type;
266  const char prefix[] = "computeNumPacketsAndOffsets: ";
267 
268  count_type count = 0;
269  const count_type numRowsToPack = lclRowInds.extent (0);
270 
271  if (numRowsToPack == 0) {
272  return count;
273  }
274  else {
275  TEUCHOS_TEST_FOR_EXCEPTION
276  (rowOffsets.extent (0) <= static_cast<size_type> (1),
277  std::invalid_argument, prefix << "There is at least one row to pack, "
278  "but the matrix has no rows. lclRowInds.extent(0) = " <<
279  numRowsToPack << ", but rowOffsets.extent(0) = " <<
280  rowOffsets.extent (0) << " <= 1.");
281  TEUCHOS_TEST_FOR_EXCEPTION
282  (outputOffsets.extent (0) !=
283  static_cast<size_type> (numRowsToPack + 1), std::invalid_argument,
284  prefix << "Output dimension does not match number of rows to pack. "
285  << "outputOffsets.extent(0) = " << outputOffsets.extent (0)
286  << " != lclRowInds.extent(0) + 1 = "
287  << static_cast<size_type> (numRowsToPack + 1) << ".");
288  TEUCHOS_TEST_FOR_EXCEPTION
289  (counts.extent (0) != numRowsToPack, std::invalid_argument,
290  prefix << "counts.extent(0) = " << counts.extent (0)
291  << " != numRowsToPack = " << numRowsToPack << ".");
292 
293  functor_type f (outputOffsets, counts, rowOffsets,
294  lclRowInds, lclRowPids, sizeOfLclCount,
295  sizeOfGblColInd, sizeOfPid, sizeOfValue);
296  Kokkos::parallel_scan (range_type (0, numRowsToPack + 1), f);
297 
298  // At least in debug mode, this functor checks for errors.
299  const int errCode = f.getError ();
300  TEUCHOS_TEST_FOR_EXCEPTION
301  (errCode != 0, std::runtime_error, prefix << "parallel_scan error code "
302  << errCode << " != 0.");
303 
304 #if 0
305  size_t total = 0;
306  for (LO k = 0; k < numRowsToPack; ++k) {
307  total += counts[k];
308  }
309  if (outputOffsets(numRowsToPack) != total) {
310  if (errStr.get () == NULL) {
311  errStr = std::unique_ptr<std::ostringstream> (new std::ostringstream ());
312  }
313  std::ostringstream& os = *errStr;
314  os << prefix
315  << "outputOffsets(numRowsToPack=" << numRowsToPack << ") "
316  << outputOffsets(numRowsToPack) << " != sum of counts = "
317  << total << "." << std::endl;
318  if (numRowsToPack != 0) {
319  // Only print the array if it's not too long.
320  if (numRowsToPack < static_cast<LO> (10)) {
321  os << "outputOffsets: [";
322  for (LO i = 0; i <= numRowsToPack; ++i) {
323  os << outputOffsets(i);
324  if (static_cast<LO> (i + 1) <= numRowsToPack) {
325  os << ",";
326  }
327  }
328  os << "]" << std::endl;
329  os << "counts: [";
330  for (LO i = 0; i < numRowsToPack; ++i) {
331  os << counts(i);
332  if (static_cast<LO> (i + 1) < numRowsToPack) {
333  os << ",";
334  }
335  }
336  os << "]" << std::endl;
337  }
338  else {
339  os << "outputOffsets(" << (numRowsToPack-1) << ") = "
340  << outputOffsets(numRowsToPack-1) << "." << std::endl;
341  }
342  }
343  count = outputOffsets(numRowsToPack);
344  return {false, errStr};
345  }
346 #endif // HAVE_TPETRA_DEBUG
347 
348  // Get last entry of outputOffsets, which is the sum of the entries
349  // of counts. Don't assume UVM.
350  using Tpetra::Details::getEntryOnHost;
351  return static_cast<count_type> (getEntryOnHost (outputOffsets,
352  numRowsToPack));
353  }
354 }
355 
371 template<class ST, class ColumnMap, class BufferDeviceType>
372 KOKKOS_FUNCTION
373 Kokkos::pair<int, size_t>
374 packCrsMatrixRow (const ColumnMap& col_map,
375  const Kokkos::View<char*, BufferDeviceType>& exports,
376  const typename PackTraits<typename ColumnMap::local_ordinal_type>::input_array_type& lids_in,
377  const typename PackTraits<int>::input_array_type& pids_in,
378  const typename PackTraits<ST>::input_array_type& vals_in,
379  const size_t offset,
380  const size_t num_ent,
381  const size_t num_bytes_per_value,
382  const bool pack_pids)
383 {
384  using Kokkos::subview;
385  using LO = typename ColumnMap::local_ordinal_type;
386  using GO = typename ColumnMap::global_ordinal_type;
387  using return_type = Kokkos::pair<int, size_t>;
388 
389  if (num_ent == 0) {
390  // Empty rows always take zero bytes, to ensure sparsity.
391  return return_type (0, 0);
392  }
393 
394  const LO num_ent_LO = static_cast<LO> (num_ent); // packValueCount wants this
395  const size_t num_ent_beg = offset;
396  const size_t num_ent_len = PackTraits<LO>::packValueCount (num_ent_LO);
397 
398  const size_t gids_beg = num_ent_beg + num_ent_len;
399  const size_t gids_len = num_ent * PackTraits<GO>::packValueCount (GO (0));
400 
401  const size_t pids_beg = gids_beg + gids_len;
402  const size_t pids_len = pack_pids ?
403  num_ent * PackTraits<int>::packValueCount (int (0)) :
404  static_cast<size_t> (0);
405 
406  const size_t vals_beg = gids_beg + gids_len + pids_len;
407  const size_t vals_len = num_ent * num_bytes_per_value;
408 
409  char* const num_ent_out = exports.data () + num_ent_beg;
410  char* const gids_out = exports.data () + gids_beg;
411  char* const pids_out = pack_pids ? exports.data () + pids_beg : NULL;
412  char* const vals_out = exports.data () + vals_beg;
413 
414  size_t num_bytes_out = 0;
415  int error_code = 0;
416  num_bytes_out += PackTraits<LO>::packValue (num_ent_out, num_ent_LO);
417 
418  {
419  // Copy column indices one at a time, so that we don't need
420  // temporary storage.
421  for (size_t k = 0; k < num_ent; ++k) {
422  const LO lid = lids_in[k];
423  const GO gid = col_map.getGlobalElement (lid);
424  num_bytes_out += PackTraits<GO>::packValue (gids_out, k, gid);
425  }
426  // Copy PIDs one at a time, so that we don't need temporary storage.
427  if (pack_pids) {
428  for (size_t k = 0; k < num_ent; ++k) {
429  const LO lid = lids_in[k];
430  const int pid = pids_in[lid];
431  num_bytes_out += PackTraits<int>::packValue (pids_out, k, pid);
432  }
433  }
434  const auto p =
435  PackTraits<ST>::packArray (vals_out, vals_in.data (), num_ent);
436  error_code += p.first;
437  num_bytes_out += p.second;
438  }
439 
440  if (error_code != 0) {
441  return return_type (10, num_bytes_out);
442  }
443 
444  const size_t expected_num_bytes =
445  num_ent_len + gids_len + pids_len + vals_len;
446  if (num_bytes_out != expected_num_bytes) {
447  return return_type (11, num_bytes_out);
448  }
449  return return_type (0, num_bytes_out);
450 }
451 
452 template<class LocalMatrix, class LocalMap, class BufferDeviceType>
453 struct PackCrsMatrixFunctor {
454  typedef LocalMatrix local_matrix_device_type;
455  typedef LocalMap local_map_type;
456  typedef typename local_matrix_device_type::value_type ST;
457  typedef typename local_map_type::local_ordinal_type LO;
458  typedef typename local_map_type::global_ordinal_type GO;
459  typedef typename local_matrix_device_type::device_type DT;
460 
461  typedef Kokkos::View<const size_t*, BufferDeviceType>
462  num_packets_per_lid_view_type;
463  typedef Kokkos::View<const size_t*, BufferDeviceType> offsets_view_type;
464  typedef Kokkos::View<char*, BufferDeviceType> exports_view_type;
465  using export_lids_view_type = typename PackTraits<LO>::input_array_type;
466  using source_pids_view_type = typename PackTraits<int>::input_array_type;
467 
468  typedef typename num_packets_per_lid_view_type::non_const_value_type
469  count_type;
470  typedef typename offsets_view_type::non_const_value_type
471  offset_type;
472  typedef Kokkos::pair<int, LO> value_type;
473 
474  static_assert (std::is_same<LO, typename local_matrix_device_type::ordinal_type>::value,
475  "local_map_type::local_ordinal_type and "
476  "local_matrix_device_type::ordinal_type must be the same.");
477 
478  local_matrix_device_type local_matrix;
479  local_map_type local_col_map;
480  exports_view_type exports;
481  num_packets_per_lid_view_type num_packets_per_lid;
482  export_lids_view_type export_lids;
483  source_pids_view_type source_pids;
484  offsets_view_type offsets;
485  size_t num_bytes_per_value;
486  bool pack_pids;
487 
488  PackCrsMatrixFunctor (const local_matrix_device_type& local_matrix_in,
489  const local_map_type& local_col_map_in,
490  const exports_view_type& exports_in,
491  const num_packets_per_lid_view_type& num_packets_per_lid_in,
492  const export_lids_view_type& export_lids_in,
493  const source_pids_view_type& source_pids_in,
494  const offsets_view_type& offsets_in,
495  const size_t num_bytes_per_value_in,
496  const bool pack_pids_in) :
497  local_matrix (local_matrix_in),
498  local_col_map (local_col_map_in),
499  exports (exports_in),
500  num_packets_per_lid (num_packets_per_lid_in),
501  export_lids (export_lids_in),
502  source_pids (source_pids_in),
503  offsets (offsets_in),
504  num_bytes_per_value (num_bytes_per_value_in),
505  pack_pids (pack_pids_in)
506  {
507  const LO numRows = local_matrix_in.numRows ();
508  const LO rowMapDim =
509  static_cast<LO> (local_matrix.graph.row_map.extent (0));
510  TEUCHOS_TEST_FOR_EXCEPTION
511  (numRows != 0 && rowMapDim != numRows + static_cast<LO> (1),
512  std::logic_error, "local_matrix.graph.row_map.extent(0) = "
513  << rowMapDim << " != numRows (= " << numRows << " ) + 1.");
514  }
515 
516  KOKKOS_INLINE_FUNCTION void init (value_type& dst) const
517  {
518  using ::Tpetra::Details::OrdinalTraits;
519  dst = Kokkos::make_pair (0, OrdinalTraits<LO>::invalid ());
520  }
521 
522  KOKKOS_INLINE_FUNCTION void
523  join (value_type& dst, const value_type& src) const
524  {
525  // `dst` should reflect the first (least) bad index and all other
526  // associated error codes and data, so prefer keeping it.
527  if (src.first != 0 && dst.first == 0) {
528  dst = src;
529  }
530  }
531 
532  KOKKOS_INLINE_FUNCTION
533  void operator() (const LO i, value_type& dst) const
534  {
535  const size_t offset = offsets[i];
536  const LO export_lid = export_lids[i];
537  const size_t buf_size = exports.size();
538  const size_t num_bytes = num_packets_per_lid(i);
539  const size_t num_ent =
540  static_cast<size_t> (local_matrix.graph.row_map[export_lid+1]
541  - local_matrix.graph.row_map[export_lid]);
542 
543  // Only pack this row's data if it has a nonzero number of
544  // entries. We can do this because receiving processes get the
545  // number of packets, and will know that zero packets means zero
546  // entries.
547  if (num_ent == 0) {
548  return;
549  }
550 
551  if (export_lid >= local_matrix.numRows ()) {
552  if (dst.first != 0) { // keep only the first error
553  dst = Kokkos::make_pair (1, i); // invalid row
554  }
555  return;
556  }
557  else if ((offset > buf_size || offset + num_bytes > buf_size)) {
558  if (dst.first != 0) { // keep only the first error
559  dst = Kokkos::make_pair (2, i); // out of bounds
560  }
561  return;
562  }
563 
564  // We can now pack this row
565 
566  // Since the matrix is locally indexed on the calling process, we
567  // have to use its column Map (which it _must_ have in this case)
568  // to convert to global indices.
569  const auto row_beg = local_matrix.graph.row_map[export_lid];
570  const auto row_end = local_matrix.graph.row_map[export_lid + 1];
571  auto vals_in = subview (local_matrix.values,
572  Kokkos::make_pair (row_beg, row_end));
573  auto lids_in = subview (local_matrix.graph.entries,
574  Kokkos::make_pair (row_beg, row_end));
575  typedef local_map_type LMT;
576  typedef BufferDeviceType BDT;
577  auto p = packCrsMatrixRow<ST, LMT, BDT> (local_col_map, exports, lids_in,
578  source_pids, vals_in, offset,
579  num_ent, num_bytes_per_value,
580  pack_pids);
581  int error_code_this_row = p.first;
582  size_t num_bytes_packed_this_row = p.second;
583  if (error_code_this_row != 0) {
584  if (dst.first != 0) { // keep only the first error
585  dst = Kokkos::make_pair (error_code_this_row, i); // bad pack
586  }
587  }
588  else if (num_bytes_packed_this_row != num_bytes) {
589  if (dst.first != 0) { // keep only the first error
590  dst = Kokkos::make_pair (3, i);
591  }
592  }
593  }
594 };
595 
603 template<class LocalMatrix, class LocalMap, class BufferDeviceType>
604 void
605 do_pack (const LocalMatrix& local_matrix,
606  const LocalMap& local_map,
607  const Kokkos::View<char*, BufferDeviceType>& exports,
608  const typename PackTraits<size_t>::input_array_type& num_packets_per_lid,
609  const typename PackTraits<typename LocalMap::local_ordinal_type>::input_array_type& export_lids,
610  const typename PackTraits<int>::input_array_type& source_pids,
611  const Kokkos::View<const size_t*, BufferDeviceType>& offsets,
612  const size_t num_bytes_per_value,
613  const bool pack_pids)
614 {
615  using LO = typename LocalMap::local_ordinal_type;
616  using DT = typename LocalMatrix::device_type;
617  using range_type = Kokkos::RangePolicy<typename DT::execution_space, LO>;
618  const char prefix[] = "Tpetra::Details::do_pack: ";
619 
620  if (export_lids.extent (0) != 0) {
621  TEUCHOS_TEST_FOR_EXCEPTION
622  (static_cast<size_t> (offsets.extent (0)) !=
623  static_cast<size_t> (export_lids.extent (0) + 1),
624  std::invalid_argument, prefix << "offsets.extent(0) = "
625  << offsets.extent (0) << " != export_lids.extent(0) (= "
626  << export_lids.extent (0) << ") + 1.");
627  TEUCHOS_TEST_FOR_EXCEPTION
628  (export_lids.extent (0) != num_packets_per_lid.extent (0),
629  std::invalid_argument, prefix << "export_lids.extent(0) = " <<
630  export_lids.extent (0) << " != num_packets_per_lid.extent(0) = "
631  << num_packets_per_lid.extent (0) << ".");
632  // If exports has nonzero length at this point, then the matrix
633  // has at least one entry to pack. Thus, if packing process
634  // ranks, we had better have at least one process rank to pack.
635  TEUCHOS_TEST_FOR_EXCEPTION
636  (pack_pids && exports.extent (0) != 0 &&
637  source_pids.extent (0) == 0, std::invalid_argument, prefix <<
638  "pack_pids is true, and exports.extent(0) = " <<
639  exports.extent (0) << " != 0, meaning that we need to pack at "
640  "least one matrix entry, but source_pids.extent(0) = 0.");
641  }
642 
643  using pack_functor_type =
644  PackCrsMatrixFunctor<LocalMatrix, LocalMap, BufferDeviceType>;
645  pack_functor_type f (local_matrix, local_map, exports,
646  num_packets_per_lid, export_lids,
647  source_pids, offsets, num_bytes_per_value,
648  pack_pids);
649 
650  typename pack_functor_type::value_type result;
651  range_type range (0, num_packets_per_lid.extent (0));
652  Kokkos::parallel_reduce (range, f, result);
653 
654  if (result.first != 0) {
655  // We can't deep_copy from AnonymousSpace Views, so we can't print
656  // out any information from them in case of error.
657  TEUCHOS_TEST_FOR_EXCEPTION
658  (true, std::runtime_error, prefix << "PackCrsMatrixFunctor "
659  "reported error code " << result.first << " for the first "
660  "bad row " << result.second << ".");
661  }
662 }
663 
693 template<typename ST, typename LO, typename GO, typename NT, typename BufferDeviceType>
694 void
695 packCrsMatrix (const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
696  Kokkos::DualView<char*, BufferDeviceType>& exports,
697  const Kokkos::View<size_t*, BufferDeviceType>& num_packets_per_lid,
698  const Kokkos::View<const LO*, BufferDeviceType>& export_lids,
699  const Kokkos::View<const int*, typename NT::device_type>& export_pids,
700  size_t& constant_num_packets,
701  const bool pack_pids)
702 {
703  ::Tpetra::Details::ProfilingRegion region_pack_crs_matrix(
704  "Tpetra::Details::PackCrsMatrixImpl::packCrsMatrix",
705  "Import/Export"
706  );
707  using Kokkos::View;
708  typedef BufferDeviceType DT;
709  typedef Kokkos::DualView<char*, BufferDeviceType> exports_view_type;
710  const char prefix[] = "Tpetra::Details::PackCrsMatrixImpl::packCrsMatrix: ";
711  constexpr bool debug = false;
712 
713  auto local_matrix = sourceMatrix.getLocalMatrixDevice ();
714  auto local_col_map = sourceMatrix.getColMap ()->getLocalMap ();
715 
716  // Setting this to zero tells the caller to expect a possibly
717  // different ("nonconstant") number of packets per local index
718  // (i.e., a possibly different number of entries per row).
719  constant_num_packets = 0;
720 
721  const size_t num_export_lids =
722  static_cast<size_t> (export_lids.extent (0));
723  TEUCHOS_TEST_FOR_EXCEPTION
724  (num_export_lids !=
725  static_cast<size_t> (num_packets_per_lid.extent (0)),
726  std::invalid_argument, prefix << "num_export_lids.extent(0) = "
727  << num_export_lids << " != num_packets_per_lid.extent(0) = "
728  << num_packets_per_lid.extent (0) << ".");
729  if (num_export_lids != 0) {
730  TEUCHOS_TEST_FOR_EXCEPTION
731  (num_packets_per_lid.data () == NULL, std::invalid_argument,
732  prefix << "num_export_lids = "<< num_export_lids << " != 0, but "
733  "num_packets_per_lid.data() = "
734  << num_packets_per_lid.data () << " == NULL.");
735  }
736 
737  const size_t num_bytes_per_lid = PackTraits<LO>::packValueCount (LO (0));
738  const size_t num_bytes_per_gid = PackTraits<GO>::packValueCount (GO (0));
739  const size_t num_bytes_per_pid = PackTraits<int>::packValueCount (int (0));
740 
741  size_t num_bytes_per_value = 0;
742  if (PackTraits<ST>::compileTimeSize) {
743  // Assume ST is default constructible; packValueCount wants an instance.
744  num_bytes_per_value = PackTraits<ST>::packValueCount (ST ());
745  }
746  else {
747  // Since the packed data come from the source matrix, we can use
748  // the source matrix to get the number of bytes per Scalar value
749  // stored in the matrix. This assumes that all Scalar values in
750  // the source matrix require the same number of bytes. If the
751  // source matrix has no entries on the calling process, then we
752  // hope that some process does have some idea how big a Scalar
753  // value is. Of course, if no processes have any entries, then no
754  // values should be packed (though this does assume that in our
755  // packing scheme, rows with zero entries take zero bytes).
756  size_t num_bytes_per_value_l = 0;
757  if (local_matrix.values.extent(0) > 0) {
758  const ST& val = local_matrix.values(0);
759  num_bytes_per_value_l = PackTraits<ST>::packValueCount (val);
760  }
761  using Teuchos::reduceAll;
762  reduceAll<int, size_t> (* (sourceMatrix.getComm ()),
763  Teuchos::REDUCE_MAX,
764  num_bytes_per_value_l,
765  Teuchos::outArg (num_bytes_per_value));
766  }
767 
768  if (num_export_lids == 0) {
769  exports = exports_view_type ("exports", 0);
770  return;
771  }
772 
773  // Array of offsets into the pack buffer.
774  Kokkos::View<size_t*, DT> offsets ("offsets", num_export_lids + 1);
775 
776  // Compute number of packets per LID (row to send), as well as
777  // corresponding offsets (the prefix sum of the packet counts).
778  const size_t count =
779  computeNumPacketsAndOffsets (offsets, num_packets_per_lid,
780  local_matrix.graph.row_map, export_lids,
781  export_pids,
782  num_bytes_per_lid, num_bytes_per_gid,
783  num_bytes_per_pid, num_bytes_per_value);
784 
785  // Resize the output pack buffer if needed.
786  if (count > static_cast<size_t> (exports.extent (0))) {
787  exports = exports_view_type ("exports", count);
788  if (debug) {
789  std::ostringstream os;
790  os << "*** exports resized to " << count << std::endl;
791  std::cerr << os.str ();
792  }
793  }
794  if (debug) {
795  std::ostringstream os;
796  os << "*** count: " << count << ", exports.extent(0): "
797  << exports.extent (0) << std::endl;
798  std::cerr << os.str ();
799  }
800 
801  // If exports has nonzero length at this point, then the matrix has
802  // at least one entry to pack. Thus, if packing process ranks, we
803  // had better have at least one process rank to pack.
804  TEUCHOS_TEST_FOR_EXCEPTION
805  (pack_pids && exports.extent (0) != 0 &&
806  export_pids.extent (0) == 0, std::invalid_argument, prefix <<
807  "pack_pids is true, and exports.extent(0) = " <<
808  exports.extent (0) << " != 0, meaning that we need to pack at least "
809  "one matrix entry, but export_pids.extent(0) = 0.");
810 
811  typedef typename std::decay<decltype (local_matrix)>::type
812  local_matrix_device_type;
813  typedef typename std::decay<decltype (local_col_map)>::type
814  local_map_type;
815 
816  exports.modify_device ();
817  auto exports_d = exports.view_device ();
818  do_pack<local_matrix_device_type, local_map_type, DT>
819  (local_matrix, local_col_map, exports_d, num_packets_per_lid,
820  export_lids, export_pids, offsets, num_bytes_per_value,
821  pack_pids);
822  // If we got this far, we succeeded.
823 }
824 
825 } // namespace PackCrsMatrixImpl
826 
827 template<typename ST, typename LO, typename GO, typename NT>
828 void
830  Teuchos::Array<char>& exports,
831  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
832  const Teuchos::ArrayView<const LO>& exportLIDs,
833  size_t& constantNumPackets)
834 {
835 
837  using buffer_device_type = typename DistObject<char, LO, GO, NT>::buffer_device_type;
838  using host_exec_space = typename Kokkos::View<size_t*, device_type>::HostMirror::execution_space;
839  using device_exec_space = typename device_type::execution_space;
840  using host_dev_type = Kokkos::Device<host_exec_space, Kokkos::HostSpace>;
841 
842  // Convert all Teuchos::Array to Kokkos::View
843 
844  // This is an output array, so we don't have to copy to device here.
845  // However, we'll have to remember to copy back to host when done.
846  Kokkos::View<size_t*, buffer_device_type> num_packets_per_lid_d =
847  create_mirror_view_from_raw_host_array (buffer_device_type (),
848  numPacketsPerLID.getRawPtr (),
849  numPacketsPerLID.size (), false,
850  "num_packets_per_lid");
851  // FIXME (mfh 05 Feb 2019) We should just pass the exportLIDs
852  // DualView through here, instead of recreating a device View from a
853  // host ArrayView that itself came from a DualView.
854  //
855  // This is an input array, so we have to copy to device here.
856  // However, we never need to copy it back to host.
857  Kokkos::View<const LO*, buffer_device_type> export_lids_d =
858  create_mirror_view_from_raw_host_array (buffer_device_type (),
859  exportLIDs.getRawPtr (),
860  exportLIDs.size (), true,
861  "export_lids");
862 
863  Kokkos::View<int*, device_type> export_pids_d; // output arg
864  Kokkos::DualView<char*, buffer_device_type> exports_dv; // output arg
865  constexpr bool pack_pids = false;
866  PackCrsMatrixImpl::packCrsMatrix<ST, LO, GO, NT, buffer_device_type> (
867  sourceMatrix, exports_dv, num_packets_per_lid_d, export_lids_d,
868  export_pids_d, constantNumPackets, pack_pids);
869 
870  // The counts are an output of PackCrsMatrixImpl::packCrsMatrix, so we have to
871  // copy them back to host.
872  Kokkos::View<size_t*, host_dev_type> num_packets_per_lid_h
873  (numPacketsPerLID.getRawPtr (),
874  numPacketsPerLID.size ());
875  // DEEP_COPY REVIEW - DEVICE-TO-HOST
876  Kokkos::deep_copy (device_exec_space(), num_packets_per_lid_h, num_packets_per_lid_d);
877 
878  // FIXME (mfh 23 Aug 2017) If we're forced to use a DualView for
879  // exports_dv above, then we have two host copies for exports_h.
880 
881  // The exports are an output of PackCrsMatrixImpl::packCrsMatrix, so we have
882  // to copy them back to host.
883  if (static_cast<size_t> (exports.size ()) !=
884  static_cast<size_t> (exports_dv.extent (0))) {
885  exports.resize (exports_dv.extent (0));
886  }
887  Kokkos::View<char*, host_dev_type> exports_h (exports.getRawPtr (),
888  exports.size ());
889  // DEEP_COPY REVIEW - DEVICE-TO-HOST
890  Kokkos::deep_copy (device_exec_space(), exports_h, exports_dv.d_view);
891 }
892 
893 template<typename ST, typename LO, typename GO, typename NT>
894 void
896  const CrsMatrix<ST, LO, GO, NT>& sourceMatrix,
897  Kokkos::DualView<char*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exports,
898  const Kokkos::DualView<size_t*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& numPacketsPerLID,
899  const Kokkos::DualView<const LO*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exportLIDs,
900  size_t& constantNumPackets)
901 {
902  using device_type = typename CrsMatrix<ST, LO, GO, NT>::device_type;
903  using buffer_device_type = typename DistObject<char, LO, GO, NT>::buffer_device_type;
904 
905  // Create an empty array of PIDs, since the interface needs it.
906  Kokkos::View<int*, device_type> exportPIDs_d ("exportPIDs", 0);
907  constexpr bool pack_pids = false;
908 
909  // Write-only device access
910  auto numPacketsPerLID_nc = numPacketsPerLID; // const DV& -> DV
911  numPacketsPerLID_nc.clear_sync_state ();
912  numPacketsPerLID_nc.modify_device ();
913  auto numPacketsPerLID_d = numPacketsPerLID.view_device ();
914 
915  // Read-only device access
916  TEUCHOS_ASSERT( ! exportLIDs.need_sync_device () );
917  auto exportLIDs_d = exportLIDs.view_device ();
918 
919  ::Tpetra::Details::ProfilingRegion region_pack_crs_matrix_new(
920  "Tpetra::Details::packCrsMatrixNew",
921  "Import/Export"
922  );
923  PackCrsMatrixImpl::packCrsMatrix<ST,LO,GO,NT,buffer_device_type> (
924  sourceMatrix, exports, numPacketsPerLID_d, exportLIDs_d,
925  exportPIDs_d, constantNumPackets, pack_pids);
926 }
927 
928 template<typename ST, typename LO, typename GO, typename NT>
929 void
931  Kokkos::DualView<char*, typename DistObject<char, LO, GO, NT>::buffer_device_type>& exports_dv,
932  const Teuchos::ArrayView<size_t>& numPacketsPerLID,
933  const Teuchos::ArrayView<const LO>& exportLIDs,
934  const Teuchos::ArrayView<const int>& sourcePIDs,
935  size_t& constantNumPackets)
936 {
937  typedef typename CrsMatrix<ST,LO,GO,NT>::local_matrix_device_type local_matrix_device_type;
938  typedef typename DistObject<char, LO, GO, NT>::buffer_device_type buffer_device_type;
939  typedef typename Kokkos::DualView<char*, buffer_device_type>::t_host::execution_space host_exec_space;
940  typedef Kokkos::Device<host_exec_space, Kokkos::HostSpace> host_dev_type;
941 
942  typename local_matrix_device_type::device_type outputDevice;
943  typedef typename NT::execution_space execution_space;
944 
945 
946  const bool verbose = ::Tpetra::Details::Behavior::verbose ();
947  std::unique_ptr<std::string> prefix;
948  if (verbose) {
949  const int myRank = [&] () {
950  auto map = sourceMatrix.getMap ();
951  if (map.get () == nullptr) {
952  return -1;
953  }
954  auto comm = map->getComm ();
955  if (comm.get () == nullptr) {
956  return -2;
957  }
958  return comm->getRank ();
959  } ();
960  std::ostringstream os;
961  os << "Proc " << myRank << ": packCrsMatrixWithOwningPIDs: ";
962  prefix = std::unique_ptr<std::string> (new std::string (os.str ()));
963 
964  std::ostringstream os2;
965  os2 << *prefix << "start" << std::endl;
966  std::cerr << os2.str ();
967  }
968 
969  // Convert all Teuchos::Array to Kokkos::View
970 
971  // This is an output array, so we don't have to copy to device here.
972  // However, we'll have to remember to copy back to host when done.
973  auto num_packets_per_lid_d =
974  create_mirror_view_from_raw_host_array (buffer_device_type (),
975  numPacketsPerLID.getRawPtr (),
976  numPacketsPerLID.size (), false,
977  "num_packets_per_lid");
978 
979  // This is an input array, so we have to copy to device here.
980  // However, we never need to copy it back to host.
981  auto export_lids_d =
982  create_mirror_view_from_raw_host_array (buffer_device_type (),
983  exportLIDs.getRawPtr (),
984  exportLIDs.size (), true,
985  "export_lids");
986  // This is an input array, so we have to copy to device here.
987  // However, we never need to copy it back to host.
988  auto export_pids_d =
990  sourcePIDs.getRawPtr (),
991  sourcePIDs.size (), true,
992  "export_pids");
993  constexpr bool pack_pids = true;
994  try {
996  (sourceMatrix, exports_dv, num_packets_per_lid_d, export_lids_d,
997  export_pids_d, constantNumPackets, pack_pids);
998  }
999  catch (std::exception& e) {
1000  if (verbose) {
1001  std::ostringstream os;
1002  os << *prefix << "PackCrsMatrixImpl::packCrsMatrix threw: "
1003  << e.what () << std::endl;
1004  std::cerr << os.str ();
1005  }
1006  throw;
1007  }
1008  catch (...) {
1009  if (verbose) {
1010  std::ostringstream os;
1011  os << *prefix << "PackCrsMatrixImpl::packCrsMatrix threw an exception "
1012  "not a subclass of std::exception" << std::endl;
1013  std::cerr << os.str ();
1014  }
1015  throw;
1016  }
1017 
1018  if (numPacketsPerLID.size () != 0) {
1019  try {
1020  // The counts are an output of PackCrsMatrixImpl::packCrsMatrix,
1021  // so we have to copy them back to host.
1022  Kokkos::View<size_t*, host_dev_type> num_packets_per_lid_h
1023  (numPacketsPerLID.getRawPtr (), numPacketsPerLID.size ());
1024  // DEEP_COPY REVIEW - DEVICE-TO-HOST
1025  Kokkos::deep_copy (execution_space(), num_packets_per_lid_h, num_packets_per_lid_d);
1026  }
1027  catch (std::exception& e) {
1028  if (verbose) {
1029  std::ostringstream os;
1030  os << *prefix << "Kokkos::deep_copy threw: " << e.what () << std::endl;
1031  std::cerr << os.str ();
1032  }
1033  throw;
1034  }
1035  catch (...) {
1036  if (verbose) {
1037  std::ostringstream os;
1038  os << *prefix << "Kokkos::deep_copy threw an exception not a subclass "
1039  "of std::exception" << std::endl;
1040  std::cerr << os.str ();
1041  }
1042  throw;
1043  }
1044  }
1045 
1046  if (verbose) {
1047  std::ostringstream os;
1048  os << *prefix << "done" << std::endl;
1049  std::cerr << os.str ();
1050  }
1051 }
1052 
1053 } // namespace Details
1054 } // namespace Tpetra
1055 
1056 #define TPETRA_DETAILS_PACKCRSMATRIX_INSTANT( ST, LO, GO, NT ) \
1057  template void \
1058  Details::packCrsMatrix<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
1059  Teuchos::Array<char>&, \
1060  const Teuchos::ArrayView<size_t>&, \
1061  const Teuchos::ArrayView<const LO>&, \
1062  size_t&); \
1063  template void \
1064  Details::packCrsMatrixNew<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
1065  Kokkos::DualView<char*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
1066  const Kokkos::DualView<size_t*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
1067  const Kokkos::DualView<const LO*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
1068  size_t&); \
1069  template void \
1070  Details::packCrsMatrixWithOwningPIDs<ST, LO, GO, NT> (const CrsMatrix<ST, LO, GO, NT>&, \
1071  Kokkos::DualView<char*, DistObject<char, LO, GO, NT>::buffer_device_type>&, \
1072  const Teuchos::ArrayView<size_t>&, \
1073  const Teuchos::ArrayView<const LO>&, \
1074  const Teuchos::ArrayView<const int>&, \
1075  size_t&);
1076 
1077 #endif // TPETRA_DETAILS_PACKCRSMATRIX_DEF_HPP
GlobalOrdinal global_ordinal_type
The type of global indices.
Impl::CreateMirrorViewFromUnmanagedHostArray< ValueType, OutputDeviceType >::output_view_type create_mirror_view_from_raw_host_array(const OutputDeviceType &, ValueType *inPtr, const size_t inSize, const bool copy=true, const char label[]="")
Variant of Kokkos::create_mirror_view that takes a raw host 1-d array as input.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Import KokkosSparse::OrdinalTraits, a traits class for &quot;invalid&quot; (flag) values of integer types...
Declaration of Tpetra::Details::Profiling, a scope guard for Kokkos Profiling.
static KOKKOS_INLINE_FUNCTION Kokkos::pair< int, size_t > packArray(char outBuf[], const value_type inBuf[], const size_t numEnt)
Pack the first numEnt entries of the given input buffer of value_type, into the output buffer of byte...
Declaration of the Tpetra::CrsMatrix class.
Declaration and generic definition of traits class that tells Tpetra::CrsMatrix how to pack and unpac...
KokkosSparse::CrsMatrix< impl_scalar_type, local_ordinal_type, device_type, void, typename local_graph_device_type::size_type > local_matrix_device_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
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.
static KOKKOS_INLINE_FUNCTION size_t packValue(char outBuf[], const T &inVal)
Pack the given value of type value_type into the given output buffer of bytes (char).
Compute the number of packets and offsets for the pack procedure.
void packCrsMatrixNew(const CrsMatrix< ST, LO, GO, NT > &sourceMatrix, Kokkos::DualView< char *, typename DistObject< char, LO, GO, NT >::buffer_device_type > &exports, const Kokkos::DualView< size_t *, typename DistObject< char, LO, GO, NT >::buffer_device_type > &numPacketsPerLID, const Kokkos::DualView< const LO *, typename DistObject< char, LO, GO, NT >::buffer_device_type > &exportLIDs, size_t &constantNumPackets)
Pack specified entries of the given local sparse matrix for communication, for &quot;new&quot; DistObject inter...
static bool verbose()
Whether Tpetra is in verbose mode.
typename Node::device_type device_type
The Kokkos device type.
void packCrsMatrix(const CrsMatrix< ST, LO, GO, NT > &sourceMatrix, Teuchos::Array< char > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, const Teuchos::ArrayView< const LO > &exportLIDs, size_t &constantNumPackets)
Pack specified entries of the given local sparse matrix for communication.
Kokkos::View< const value_type *, Kokkos::AnonymousSpace > input_array_type
The type of an input array of value_type.
static KOKKOS_INLINE_FUNCTION size_t packValueCount(const T &)
Number of bytes required to pack or unpack the given value of type value_type.
Declaration and definition of Tpetra::Details::castAwayConstDualView, an implementation detail of Tpe...
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
Declaration and definition of Tpetra::Details::getEntryOnHost.
void packCrsMatrixWithOwningPIDs(const CrsMatrix< ST, LO, GO, NT > &sourceMatrix, Kokkos::DualView< char *, typename DistObject< char, LO, GO, NT >::buffer_device_type > &exports_dv, const Teuchos::ArrayView< size_t > &numPacketsPerLID, const Teuchos::ArrayView< const LO > &exportLIDs, const Teuchos::ArrayView< const int > &sourcePIDs, size_t &constantNumPackets)
Pack specified entries of the given local sparse matrix for communication.
Base class for distributed Tpetra objects that support data redistribution.
LocalOrdinal local_ordinal_type
The type of local indices.
Functions that wrap Kokkos::create_mirror_view, in order to avoid deep copies when not necessary...
Declaration of Tpetra::Details::Behavior, a class that describes Tpetra&#39;s behavior.