Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Tpetra_Details_CooMatrix.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_COOMATRIX_HPP
11 #define TPETRA_DETAILS_COOMATRIX_HPP
12 
17 
18 #include "TpetraCore_config.h"
19 #include "Tpetra_Details_PackTriples.hpp"
20 #include "Tpetra_DistObject.hpp"
23 #include "Teuchos_TypeNameTraits.hpp"
24 
25 #include <initializer_list>
26 #include <map>
27 #include <memory>
28 #include <string>
29 #include <type_traits>
30 #include <vector>
31 
32 namespace Tpetra {
33 namespace Details {
34 
35 // Implementation details of Tpetra::Details.
36 // So, users REALLY should not use anything in here.
37 namespace Impl {
38 
41 template <class IndexType>
42 struct CooGraphEntry {
43  IndexType row;
44  IndexType col;
45 };
46 
51 template <class IndexType>
53  bool
54  operator()(const CooGraphEntry<IndexType>& x,
55  const CooGraphEntry<IndexType>& y) const {
56  return (x.row < y.row) || (x.row == y.row && x.col < y.col);
57  }
58 };
59 
62 template <class SC, class GO>
64  private:
72  typedef std::map<CooGraphEntry<GO>, SC,
74  entries_type;
75 
78  entries_type entries_;
79 
80  public:
82  typedef char packet_type;
83 
85  CooMatrixImpl() = default;
86 
93  void
94  sumIntoGlobalValue(const GO gblRowInd,
95  const GO gblColInd,
96  const SC& val) {
97  // There's no sense in worrying about the insertion hint, since
98  // indices may be all over the place. Users make no promises of
99  // sorting or locality of input.
100  CooGraphEntry<GO> ij{gblRowInd, gblColInd};
101  auto result = this->entries_.insert({ij, val});
102  if (!result.second) { // already in the map
103  result.first->second += val; // sum in the new value
104  }
105  }
106 
115  void
116  sumIntoGlobalValues(const GO gblRowInds[],
117  const GO gblColInds[],
118  const SC vals[],
119  const std::size_t numEnt) {
120  for (std::size_t k = 0; k < numEnt; ++k) {
121  // There's no sense in worrying about the insertion hint, since
122  // indices may be all over the place. Users make no promises of
123  // sorting or locality of input.
124  CooGraphEntry<GO> ij{gblRowInds[k], gblColInds[k]};
125  const SC val = vals[k];
126  auto result = this->entries_.insert({ij, val});
127  if (!result.second) { // already in the map
128  result.first->second += val; // sum in the new value
129  }
130  }
131  }
132 
134  std::size_t
136  return this->entries_.size();
137  }
138 
142  void
143  forAllEntries(std::function<void(const GO, const GO, const SC&)> f) const {
144  for (auto iter = this->entries_.begin();
145  iter != this->entries_.end(); ++iter) {
146  f(iter->first.row, iter->first.col, iter->second);
147  }
148  }
149 
155  void
156  mergeIntoRow(const GO tgtGblRow,
157  const CooMatrixImpl<SC, GO>& src,
158  const GO srcGblRow) {
159  // const char prefix[] =
160  // "Tpetra::Details::Impl::CooMatrixImpl::mergeIntoRow: ";
161 
162  entries_type& tgtEntries = this->entries_;
163  const entries_type& srcEntries = src.entries_;
164 
165  // Find all entries with the given global row index. The min GO
166  // value is guaranteed to be the least possible column index, so
167  // beg1 is a lower bound for all (row index, column index) pairs.
168  // Lower bound is inclusive; upper bound is exclusive.
169  auto srcBeg = srcEntries.lower_bound({srcGblRow, std::numeric_limits<GO>::min()});
170  auto srcEnd = srcEntries.upper_bound({srcGblRow + 1, std::numeric_limits<GO>::min()});
171  auto tgtBeg = tgtEntries.lower_bound({tgtGblRow, std::numeric_limits<GO>::min()});
172  auto tgtEnd = tgtEntries.upper_bound({tgtGblRow + 1, std::numeric_limits<GO>::min()});
173 
174  // Don't waste time iterating over the current row of *this, if
175  // the current row of src is empty.
176  if (srcBeg != srcEnd) {
177  for (auto tgtCur = tgtBeg; tgtCur != tgtEnd; ++tgtCur) {
178  auto srcCur = srcBeg;
179  while (srcCur != srcEnd && srcCur->first.col < tgtCur->first.col) {
180  ++srcCur;
181  }
182  // At this point, one of the following is true:
183  //
184  // 1. srcCur == srcEnd, thus nothing more to insert.
185  // 2. srcCur->first.col > tgtCur->first.col, thus the current
186  // row of src has no entry matching tgtCur->first, so we
187  // have to insert it. Insertion does not invalidate
188  // tgtEntries iterators, and neither srcEntries nor
189  // tgtEntries have duplicates, so this is safe.
190  // 3. srcCur->first.col == tgtCur->first.col, so add the two entries.
191  if (srcCur != srcEnd) {
192  if (srcCur->first.col == tgtCur->first.col) {
193  tgtCur->second += srcCur->second;
194  } else {
195  // tgtCur is (optimally) right before where we want to put
196  // the new entry, since srcCur points to the first entry
197  // in srcEntries whose column index is greater than that
198  // of the entry to which tgtCur points.
199  (void)tgtEntries.insert(tgtCur, *srcCur);
200  }
201  } // if srcCur != srcEnd
202  } // for each entry in the current row (lclRow) of *this
203  } // if srcBeg != srcEnd
204  }
205 
215  int countPackRow(int& numPackets,
216  const GO gblRow,
217  const ::Teuchos::Comm<int>& comm,
218  std::ostream* errStrm = NULL) const {
219  using std::endl;
220  using ::Tpetra::Details::countPackTriples;
221  using ::Tpetra::Details::countPackTriplesCount;
222  const char prefix[] = "Tpetra::Details::Impl::CooMatrixImpl::countPackRow: ";
223 #ifdef HAVE_TPETRACORE_MPI
224  int errCode = MPI_SUCCESS;
225 #else
226  int errCode = 0;
227 #endif // HAVE_TPETRACORE_MPI
228 
229  // Count the number of entries in the given row.
230  const GO minGO = std::numeric_limits<GO>::min();
231  auto beg = this->entries_.lower_bound({gblRow, minGO});
232  auto end = this->entries_.upper_bound({gblRow + 1, minGO});
233  int numEnt = 0;
234  for (auto iter = beg; iter != end; ++iter) {
235  ++numEnt;
236  if (numEnt == std::numeric_limits<int>::max()) {
237  if (errStrm != NULL) {
238  *errStrm << prefix << "In (global) row " << gblRow << ", the number "
239  "of entries thus far, numEnt = "
240  << numEnt << ", has reached the "
241  "maximum int value. We need to store this count as int for MPI's "
242  "sake."
243  << endl;
244  }
245  return -1;
246  }
247  }
248 
249  int numPacketsForCount = 0; // output argument of countPackTriplesCount
250  {
251  errCode =
252  countPackTriplesCount(comm, numPacketsForCount, errStrm);
253  if (errCode != 0) {
254  if (errStrm != NULL) {
255  *errStrm << prefix << "countPackTriplesCount "
256  "returned errCode = "
257  << errCode << " != 0." << endl;
258  }
259  return errCode;
260  }
261  if (numPacketsForCount < 0) {
262  if (errStrm != NULL) {
263  *errStrm << prefix << "countPackTriplesCount returned "
264  "numPacketsForCount = "
265  << numPacketsForCount << " < 0." << endl;
266  }
267  return -1;
268  }
269  }
270 
271  int numPacketsForTriples = 0; // output argument of countPackTriples
272  {
273  errCode = countPackTriples<SC, GO>(numEnt, comm, numPacketsForTriples);
274  TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix << "countPackTriples "
275  "returned errCode = "
276  << errCode << " != 0.");
277  TEUCHOS_TEST_FOR_EXCEPTION(numPacketsForTriples < 0, std::logic_error, prefix << "countPackTriples "
278  "returned numPacketsForTriples = "
279  << numPacketsForTriples << " < 0.");
280  }
281 
282  numPackets = numPacketsForCount + numPacketsForTriples;
283  return errCode;
284  }
285 
300  void
302  const int outBufSize,
303  int& outBufCurPos, // in/out argument
304  const ::Teuchos::Comm<int>& comm,
305  std::vector<GO>& gblRowInds,
306  std::vector<GO>& gblColInds,
307  std::vector<SC>& vals,
308  const GO gblRow) const {
309  using ::Tpetra::Details::packTriples;
310  using ::Tpetra::Details::packTriplesCount;
311  const char prefix[] = "Tpetra::Details::Impl::CooMatrixImpl::packRow: ";
312 
313  const GO minGO = std::numeric_limits<GO>::min();
314  auto beg = this->entries_.lower_bound({gblRow, minGO});
315  auto end = this->entries_.upper_bound({gblRow + 1, minGO});
316 
317  // This doesn't actually deallocate. Only swapping with an empty
318  // std::vector does that.
319  gblRowInds.resize(0);
320  gblColInds.resize(0);
321  vals.resize(0);
322 
323  int numEnt = 0;
324  for (auto iter = beg; iter != end; ++iter) {
325  gblRowInds.push_back(iter->first.row);
326  gblColInds.push_back(iter->first.col);
327  vals.push_back(iter->second);
328  ++numEnt;
329  TEUCHOS_TEST_FOR_EXCEPTION(numEnt == std::numeric_limits<int>::max(), std::runtime_error, prefix << "In (global) row " << gblRow << ", the number of entries thus far, "
330  "numEnt = "
331  << numEnt << ", has reached the maximum int value. "
332  "We need to store this count as int for MPI's sake.");
333  }
334 
335  {
336  const int errCode =
337  packTriplesCount(numEnt, outBuf, outBufSize, outBufCurPos, comm);
338  TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix << "In (global) row " << gblRow << ", packTriplesCount returned "
339  "errCode = "
340  << errCode << " != 0.");
341  }
342  {
343  const int errCode =
344  packTriples(gblRowInds.data(),
345  gblColInds.data(),
346  vals.data(),
347  numEnt,
348  outBuf,
349  outBufSize,
350  outBufCurPos, // in/out argument
351  comm);
352  TEUCHOS_TEST_FOR_EXCEPTION(errCode != 0, std::runtime_error, prefix << "In (global) row " << gblRow << ", packTriples returned errCode = " << errCode << " != 0.");
353  }
354  }
355 
363  GO getMyGlobalRowIndices(std::vector<GO>& rowInds) const {
364  rowInds.clear();
365 
366  GO lclMinRowInd = std::numeric_limits<GO>::max(); // compute local min
367  for (typename entries_type::const_iterator iter = this->entries_.begin();
368  iter != this->entries_.end(); ++iter) {
369  const GO rowInd = iter->first.row;
370  if (rowInd < lclMinRowInd) {
371  lclMinRowInd = rowInd;
372  }
373  if (rowInds.empty() || rowInds.back() != rowInd) {
374  rowInds.push_back(rowInd); // don't insert duplicates
375  }
376  }
377  return lclMinRowInd;
378  }
379 
380  template <class OffsetType>
381  void
382  buildCrs(std::vector<OffsetType>& rowOffsets,
383  GO gblColInds[],
384  SC vals[]) const {
385  static_assert(std::is_integral<OffsetType>::value,
386  "OffsetType must be a built-in integer type.");
387 
388  // clear() doesn't generally free storage; it just resets the
389  // length. Thus, we reuse any existing storage here.
390  rowOffsets.clear();
391 
392  const std::size_t numEnt = this->getLclNumEntries();
393  if (numEnt == 0) {
394  rowOffsets.push_back(0);
395  } else {
396  typename entries_type::const_iterator iter = this->entries_.begin();
397  GO prevGblRowInd = iter->first.row;
398 
399  OffsetType k = 0;
400  for (; iter != this->entries_.end(); ++iter, ++k) {
401  const GO gblRowInd = iter->first.row;
402  if (k == 0 || gblRowInd != prevGblRowInd) {
403  // The row offsets array always has at least one entry. The
404  // first entry is always zero, and the last entry is always
405  // the number of matrix entries.
406  rowOffsets.push_back(k);
407  prevGblRowInd = gblRowInd;
408  }
409  gblColInds[k] = iter->first.col;
410 
411  static_assert(std::is_same<typename std::decay<decltype(iter->second)>::type, SC>::value,
412  "The type of iter->second != SC.");
413  vals[k] = iter->second;
414  }
415  rowOffsets.push_back(static_cast<OffsetType>(numEnt));
416  }
417  }
418 
431  template <class OffsetType, class LO>
432  void
433  buildLocallyIndexedCrs(std::vector<OffsetType>& rowOffsets,
434  LO lclColInds[],
435  SC vals[],
436  std::function<LO(const GO)> gblToLcl) const {
437  static_assert(std::is_integral<OffsetType>::value,
438  "OffsetType must be a built-in integer type.");
439  static_assert(std::is_integral<LO>::value,
440  "LO must be a built-in integer type.");
441 
442  // clear() doesn't generally free storage; it just resets the
443  // length. Thus, we reuse any existing storage here.
444  rowOffsets.clear();
445 
446  const std::size_t numEnt = this->getLclNumEntries();
447  if (numEnt == 0) {
448  rowOffsets.push_back(0);
449  } else {
450  typename entries_type::const_iterator iter = this->entries_.begin();
451  GO prevGblRowInd = iter->first.row;
452 
453  OffsetType k = 0;
454  for (; iter != this->entries_.end(); ++iter, ++k) {
455  const GO gblRowInd = iter->first.row;
456  if (k == 0 || gblRowInd != prevGblRowInd) {
457  // The row offsets array always has at least one entry. The
458  // first entry is always zero, and the last entry is always
459  // the number of matrix entries.
460  rowOffsets.push_back(k);
461  prevGblRowInd = gblRowInd;
462  }
463  lclColInds[k] = gblToLcl(iter->first.col);
464  vals[k] = iter->second;
465  }
466  rowOffsets.push_back(static_cast<OffsetType>(numEnt));
467  }
468  }
469 };
470 
471 } // namespace Impl
472 
521 template <class SC,
525 class CooMatrix : public ::Tpetra::DistObject<char, LO, GO, NT> {
526  public:
528  typedef char packet_type;
530  typedef SC scalar_type;
531  typedef LO local_ordinal_type;
532  typedef GO global_ordinal_type;
533  typedef NT node_type;
534  typedef typename NT::device_type device_type;
536  typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
537 
538  private:
540  typedef ::Tpetra::DistObject<packet_type, local_ordinal_type,
541  global_ordinal_type, node_type>
542  base_type;
543 
546 
547  public:
554  : base_type(::Teuchos::null)
555  , localError_(new bool(false))
556  , errs_(new std::shared_ptr<std::ostringstream>()) // ptr to a null ptr
557  {}
558 
566  CooMatrix(const ::Teuchos::RCP<const map_type>& map)
567  : base_type(map)
568  , localError_(new bool(false))
569  , errs_(new std::shared_ptr<std::ostringstream>()) // ptr to a null ptr
570  {}
571 
573  virtual ~CooMatrix() {}
574 
581  void
582  sumIntoGlobalValue(const GO gblRowInd,
583  const GO gblColInd,
584  const SC& val) {
585  this->impl_.sumIntoGlobalValue(gblRowInd, gblColInd, val);
586  }
587 
596  void
597  sumIntoGlobalValues(const GO gblRowInds[],
598  const GO gblColInds[],
599  const SC vals[],
600  const std::size_t numEnt) {
601  this->impl_.sumIntoGlobalValues(gblRowInds, gblColInds, vals, numEnt);
602  }
603 
605  void
606  sumIntoGlobalValues(std::initializer_list<GO> gblRowInds,
607  std::initializer_list<GO> gblColInds,
608  std::initializer_list<SC> vals,
609  const std::size_t numEnt) {
610  this->impl_.sumIntoGlobalValues(gblRowInds.begin(), gblColInds.begin(),
611  vals.begin(), numEnt);
612  }
613 
615  std::size_t
617  return this->impl_.getLclNumEntries();
618  }
619 
620  template <class OffsetType>
621  void
622  buildCrs(::Kokkos::View<OffsetType*, device_type>& rowOffsets,
623  ::Kokkos::View<GO*, device_type>& gblColInds,
624 #if KOKKOS_VERSION >= 40799
625  ::Kokkos::View<typename ::KokkosKernels::ArithTraits<SC>::val_type*, device_type>& vals) {
626 #else
627  ::Kokkos::View<typename ::Kokkos::ArithTraits<SC>::val_type*, device_type>& vals) {
628 #endif
629  static_assert(std::is_integral<OffsetType>::value,
630  "OffsetType must be a built-in integer type.");
631  using ::Kokkos::create_mirror_view;
633  using ::Kokkos::View;
634 #if KOKKOS_VERSION >= 40799
635  typedef typename ::KokkosKernels::ArithTraits<SC>::val_type ISC;
636 #else
637  typedef typename ::Kokkos::ArithTraits<SC>::val_type ISC;
638 #endif
639 
640  const std::size_t numEnt = this->getLclNumEntries();
641 
642  gblColInds = View<GO*, device_type>("gblColInds", numEnt);
643  vals = View<ISC*, device_type>("vals", numEnt);
644  auto gblColInds_h = create_mirror_view(gblColInds);
645  auto vals_h = create_mirror_view(vals);
646 
647  std::vector<std::size_t> rowOffsetsSV;
648  this->impl_.buildCrs(rowOffsetsSV,
649  gblColInds_h.data(),
650  vals_h.data());
651  rowOffsets =
652  View<OffsetType*, device_type>("rowOffsets", rowOffsetsSV.size());
653  typename View<OffsetType*, device_type>::host_mirror_type
654  rowOffsets_h(rowOffsetsSV.data(), rowOffsetsSV.size());
655  deep_copy(rowOffsets, rowOffsets_h);
656 
657  deep_copy(gblColInds, gblColInds_h);
658  deep_copy(vals, vals_h);
659  }
660 
671  void
672  fillComplete(const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm) {
673  if (comm.is_null()) {
674  this->map_ = ::Teuchos::null;
675  } else {
676  this->map_ = this->buildMap(comm);
677  }
678  }
679 
688  void
690  TEUCHOS_TEST_FOR_EXCEPTION(this->getMap().is_null(), std::runtime_error,
691  "Tpetra::Details::"
692  "CooMatrix::fillComplete: This object does not yet have a Map. "
693  "You must call the version of fillComplete "
694  "that takes an input communicator.");
695  this->fillComplete(this->getMap()->getComm());
696  }
697 
712  bool localError() const {
713  return *localError_;
714  }
715 
733  std::string errorMessages() const {
734  return ((*errs_).get() == NULL) ? std::string("") : (*errs_)->str();
735  }
736 
737  private:
750  std::shared_ptr<bool> localError_;
751 
759  std::shared_ptr<std::shared_ptr<std::ostringstream> > errs_;
760 
762  std::ostream&
763  markLocalErrorAndGetStream() {
764  *(this->localError_) = true;
765  if ((*errs_).get() == NULL) {
766  *errs_ = std::shared_ptr<std::ostringstream>(new std::ostringstream());
767  }
768  return **errs_;
769  }
770 
771  public:
774  virtual std::string description() const {
775  using Teuchos::TypeNameTraits;
776 
777  std::ostringstream os;
778  os << "\"Tpetra::Details::CooMatrix\": {"
779  << "SC: " << TypeNameTraits<SC>::name()
780  << ", LO: " << TypeNameTraits<LO>::name()
781  << ", GO: " << TypeNameTraits<GO>::name()
782  << ", NT: " << TypeNameTraits<NT>::name();
783  if (this->getObjectLabel() != "") {
784  os << ", Label: \"" << this->getObjectLabel() << "\"";
785  }
786  os << ", Has Map: " << (this->map_.is_null() ? "false" : "true")
787  << "}";
788  return os.str();
789  }
790 
793  virtual void
794  describe(Teuchos::FancyOStream& out,
795  const Teuchos::EVerbosityLevel verbLevel =
796  Teuchos::Describable::verbLevel_default) const {
797  using std::endl;
798  using ::Teuchos::EVerbosityLevel;
799  using ::Teuchos::OSTab;
800  using ::Teuchos::TypeNameTraits;
801  using ::Teuchos::VERB_DEFAULT;
802  using ::Teuchos::VERB_LOW;
803  using ::Teuchos::VERB_MEDIUM;
804  using ::Tpetra::Details::gathervPrint;
805 
806  const auto vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
807 
808  auto comm = this->getMap().is_null() ? Teuchos::null : this->getMap()->getComm();
809  const int myRank = comm.is_null() ? 0 : comm->getRank();
810  // const int numProcs = comm.is_null () ? 1 : comm->getSize ();
811 
812  if (vl != Teuchos::VERB_NONE) {
813  // Convention is to start describe() implementations with a tab.
814  OSTab tab0(out);
815  if (myRank == 0) {
816  out << "\"Tpetra::Details::CooMatrix\":" << endl;
817  }
818  OSTab tab1(out);
819  if (myRank == 0) {
820  out << "Template parameters:" << endl;
821  {
822  OSTab tab2(out);
823  out << "SC: " << TypeNameTraits<SC>::name() << endl
824  << "LO: " << TypeNameTraits<LO>::name() << endl
825  << "GO: " << TypeNameTraits<GO>::name() << endl
826  << "NT: " << TypeNameTraits<NT>::name() << endl;
827  }
828  if (this->getObjectLabel() != "") {
829  out << "Label: \"" << this->getObjectLabel() << "\"" << endl;
830  }
831  out << "Has Map: " << (this->map_.is_null() ? "false" : "true") << endl;
832  } // if myRank == 0
833 
834  // Describe the Map, if it exists.
835  if (!this->map_.is_null()) {
836  if (myRank == 0) {
837  out << "Map:" << endl;
838  }
839  OSTab tab2(out);
840  this->map_->describe(out, vl);
841  }
842 
843  // At verbosity > VERB_LOW, each process prints something.
844  if (vl > VERB_LOW) {
845  std::ostringstream os;
846  os << "Process " << myRank << ":" << endl;
847  // OSTab tab3 (os);
848 
849  const std::size_t numEnt = this->impl_.getLclNumEntries();
850  os << " Local number of entries: " << numEnt << endl;
851  if (vl > VERB_MEDIUM) {
852  os << " Local entries:" << endl;
853  // OSTab tab4 (os);
854  this->impl_.forAllEntries([&os](const GO row, const GO col, const SC& val) {
855  os << " {"
856  << "row: " << row
857  << ", col: " << col
858  << ", val: " << val
859  << "}" << endl;
860  });
861  }
862  gathervPrint(out, os.str(), *comm);
863  }
864  } // vl != Teuchos::VERB_NONE
865  }
866 
867  private:
876  Teuchos::RCP<const map_type>
877  buildMap(const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm) {
878  using ::Teuchos::outArg;
879  using ::Teuchos::rcp;
880  using ::Teuchos::REDUCE_MIN;
881  using ::Teuchos::reduceAll;
882  typedef ::Tpetra::global_size_t GST;
883  // const char prefix[] = "Tpetra::Details::CooMatrix::buildMap: ";
884 
885  // Processes where comm is null don't participate in the Map.
886  if (comm.is_null()) {
887  return ::Teuchos::null;
888  }
889 
890  // mfh 17 Jan 2017: We just happen to use row indices, because
891  // that's what Tpetra::CrsMatrix currently uses. That's probably
892  // not the best thing to use, but it's not bad for commonly
893  // encountered matrices. A better more general approach might be
894  // to hash (row index, column index) pairs to a global index. One
895  // could make that unique by doing a parallel scan at map
896  // construction time.
897 
898  std::vector<GO> rowInds;
899  const GO lclMinGblRowInd = this->impl_.getMyGlobalRowIndices(rowInds);
900 
901  // Compute the globally min row index for the "index base."
902  GO gblMinGblRowInd = 0; // output argument
903  reduceAll<int, GO>(*comm, REDUCE_MIN, lclMinGblRowInd,
904  outArg(gblMinGblRowInd));
905  const GO indexBase = gblMinGblRowInd;
906  const GST INV = Tpetra::Details::OrdinalTraits<GST>::invalid();
907  return rcp(new map_type(INV, rowInds.data(), rowInds.size(),
908  indexBase, comm));
909  }
910 
911  protected:
915  virtual size_t constantNumberOfPackets() const {
916  return static_cast<size_t>(0);
917  }
918 
922  virtual bool
923  checkSizes(const ::Tpetra::SrcDistObject& source) {
924  using std::endl;
925  typedef CooMatrix<SC, LO, GO, NT> this_COO_type;
926  const char prefix[] = "Tpetra::Details::CooMatrix::checkSizes: ";
927 
928  const this_COO_type* src = dynamic_cast<const this_COO_type*>(&source);
929 
930  if (src == NULL) {
931  std::ostream& err = markLocalErrorAndGetStream();
932  err << prefix << "The source object of the Import or Export "
933  "must be a CooMatrix with the same template parameters as the "
934  "target object."
935  << endl;
936  } else if (this->map_.is_null()) {
937  std::ostream& err = markLocalErrorAndGetStream();
938  err << prefix << "The target object of the Import or Export "
939  "must be a CooMatrix with a nonnull Map."
940  << endl;
941  }
942  return !(*(this->localError_));
943  }
944 
946  using buffer_device_type =
947  typename ::Tpetra::DistObject<char, LO, GO, NT>::buffer_device_type;
948 
949  virtual void
950  copyAndPermute(const ::Tpetra::SrcDistObject& sourceObject,
951  const size_t numSameIDs,
952  const Kokkos::DualView<const LO*,
953  buffer_device_type>& permuteToLIDs,
954  const Kokkos::DualView<const LO*,
955  buffer_device_type>& permuteFromLIDs,
956  const CombineMode /* CM */) {
957  using std::endl;
958  using this_COO_type = CooMatrix<SC, LO, GO, NT>;
959  const char prefix[] = "Tpetra::Details::CooMatrix::copyAndPermute: ";
960 
961  // There's no communication in this method, so it's safe just to
962  // return on error.
963 
964  if (*(this->localError_)) {
965  std::ostream& err = this->markLocalErrorAndGetStream();
966  err << prefix << "The target object of the Import or Export is "
967  "already in an error state."
968  << endl;
969  return;
970  }
971 
972  const this_COO_type* src = dynamic_cast<const this_COO_type*>(&sourceObject);
973  if (src == nullptr) {
974  std::ostream& err = this->markLocalErrorAndGetStream();
975  err << prefix << "Input argument 'sourceObject' is not a CooMatrix."
976  << endl;
977  return;
978  }
979 
980  const size_t numPermuteIDs =
981  static_cast<size_t>(permuteToLIDs.extent(0));
982  if (numPermuteIDs != static_cast<size_t>(permuteFromLIDs.extent(0))) {
983  std::ostream& err = this->markLocalErrorAndGetStream();
984  err << prefix << "permuteToLIDs.extent(0) = "
985  << numPermuteIDs << " != permuteFromLIDs.extent(0) = "
986  << permuteFromLIDs.extent(0) << "." << endl;
987  return;
988  }
989  if (sizeof(int) <= sizeof(size_t) &&
990  numPermuteIDs > static_cast<size_t>(std::numeric_limits<int>::max())) {
991  std::ostream& err = this->markLocalErrorAndGetStream();
992  err << prefix << "numPermuteIDs = " << numPermuteIDs
993  << ", a size_t, overflows int." << endl;
994  return;
995  }
996 
997  // Even though this is an std::set, we can start with numSameIDs
998  // just by iterating through the first entries of the set.
999 
1000  if (sizeof(int) <= sizeof(size_t) &&
1001  numSameIDs > static_cast<size_t>(std::numeric_limits<int>::max())) {
1002  std::ostream& err = this->markLocalErrorAndGetStream();
1003  err << prefix << "numSameIDs = " << numSameIDs
1004  << ", a size_t, overflows int." << endl;
1005  return;
1006  }
1007 
1008  //
1009  // Copy in entries from any initial rows with the same global row indices.
1010  //
1011  const LO numSame = static_cast<int>(numSameIDs);
1012  // Count of local row indices encountered here with invalid global
1013  // row indices. If nonzero, something went wrong. If something
1014  // did go wrong, we'll defer responding until the end of this
1015  // method, so we can print as much useful info as possible.
1016  LO numInvalidSameRows = 0;
1017  for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1018  // All numSame initial rows have the same global index in both
1019  // source and target Maps, so we only need to convert to global
1020  // once.
1021  const GO gblRow = this->map_->getGlobalElement(lclRow);
1022  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1023  ++numInvalidSameRows;
1024  continue;
1025  } else {
1026  this->impl_.mergeIntoRow(gblRow, src->impl_, gblRow);
1027  }
1028  }
1029 
1030  //
1031  // Copy in entries from remaining rows that are permutations, that
1032  // is, that live in both the source and target Maps, but aren't
1033  // included in the "same" list (see above).
1034  //
1035  const LO numPermute = static_cast<int>(numPermuteIDs);
1036  // Count of local "from" row indices encountered here with invalid
1037  // global row indices. If nonzero, something went wrong. If
1038  // something did go wrong, we'll defer responding until the end of
1039  // this method, so we can print as much useful info as possible.
1040  LO numInvalidRowsFrom = 0;
1041  // Count of local "to" row indices encountered here with invalid
1042  // global row indices. If nonzero, something went wrong. If
1043  // something did go wrong, we'll defer responding until the end of
1044  // this method, so we can print as much useful info as possible.
1045  LO numInvalidRowsTo = 0;
1046 
1047  TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_host());
1048  TEUCHOS_ASSERT(!permuteToLIDs.need_sync_host());
1049  auto permuteFromLIDs_h = permuteFromLIDs.view_host();
1050  auto permuteToLIDs_h = permuteToLIDs.view_host();
1051 
1052  for (LO k = 0; k < numPermute; ++k) {
1053  const LO lclRowFrom = permuteFromLIDs_h[k];
1054  const LO lclRowTo = permuteToLIDs_h[k];
1055  const GO gblRowFrom = src->map_->getGlobalElement(lclRowFrom);
1056  const GO gblRowTo = this->map_->getGlobalElement(lclRowTo);
1057 
1058  bool bothConversionsValid = true;
1059  if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1060  ++numInvalidRowsFrom;
1061  bothConversionsValid = false;
1062  }
1063  if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1064  ++numInvalidRowsTo;
1065  bothConversionsValid = false;
1066  }
1067  if (bothConversionsValid) {
1068  this->impl_.mergeIntoRow(gblRowTo, src->impl_, gblRowFrom);
1069  }
1070  }
1071 
1072  // Print info if any errors occurred.
1073  if (numInvalidSameRows != 0 || numInvalidRowsFrom != 0 ||
1074  numInvalidRowsTo != 0) {
1075  // Collect and print all the invalid input row indices, for the
1076  // "same," "from," and "to" lists.
1077  std::vector<std::pair<LO, GO> > invalidSameRows;
1078  invalidSameRows.reserve(numInvalidSameRows);
1079  std::vector<std::pair<LO, GO> > invalidRowsFrom;
1080  invalidRowsFrom.reserve(numInvalidRowsFrom);
1081  std::vector<std::pair<LO, GO> > invalidRowsTo;
1082  invalidRowsTo.reserve(numInvalidRowsTo);
1083 
1084  for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1085  // All numSame initial rows have the same global index in both
1086  // source and target Maps, so we only need to convert to global
1087  // once.
1088  const GO gblRow = this->map_->getGlobalElement(lclRow);
1089  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1090  invalidSameRows.push_back({lclRow, gblRow});
1091  }
1092  }
1093 
1094  for (LO k = 0; k < numPermute; ++k) {
1095  const LO lclRowFrom = permuteFromLIDs_h[k];
1096  const LO lclRowTo = permuteToLIDs_h[k];
1097  const GO gblRowFrom = src->map_->getGlobalElement(lclRowFrom);
1098  const GO gblRowTo = this->map_->getGlobalElement(lclRowTo);
1099 
1100  if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1101  invalidRowsFrom.push_back({lclRowFrom, gblRowFrom});
1102  }
1103  if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1104  invalidRowsTo.push_back({lclRowTo, gblRowTo});
1105  }
1106  }
1107 
1108  std::ostringstream os;
1109  if (numInvalidSameRows != 0) {
1110  os << "Invalid permute \"same\" (local, global) index pairs: [";
1111  for (std::size_t k = 0; k < invalidSameRows.size(); ++k) {
1112  const auto& p = invalidSameRows[k];
1113  os << "(" << p.first << "," << p.second << ")";
1114  if (k + 1 < invalidSameRows.size()) {
1115  os << ", ";
1116  }
1117  }
1118  os << "]" << std::endl;
1119  }
1120  if (numInvalidRowsFrom != 0) {
1121  os << "Invalid permute \"from\" (local, global) index pairs: [";
1122  for (std::size_t k = 0; k < invalidRowsFrom.size(); ++k) {
1123  const auto& p = invalidRowsFrom[k];
1124  os << "(" << p.first << "," << p.second << ")";
1125  if (k + 1 < invalidRowsFrom.size()) {
1126  os << ", ";
1127  }
1128  }
1129  os << "]" << std::endl;
1130  }
1131  if (numInvalidRowsTo != 0) {
1132  os << "Invalid permute \"to\" (local, global) index pairs: [";
1133  for (std::size_t k = 0; k < invalidRowsTo.size(); ++k) {
1134  const auto& p = invalidRowsTo[k];
1135  os << "(" << p.first << "," << p.second << ")";
1136  if (k + 1 < invalidRowsTo.size()) {
1137  os << ", ";
1138  }
1139  }
1140  os << "]" << std::endl;
1141  }
1142 
1143  std::ostream& err = this->markLocalErrorAndGetStream();
1144  err << prefix << os.str();
1145  return;
1146  }
1147  }
1148 
1149  virtual void
1150  packAndPrepare(const ::Tpetra::SrcDistObject& sourceObject,
1151  const Kokkos::DualView<const local_ordinal_type*,
1152  buffer_device_type>& exportLIDs,
1153  Kokkos::DualView<packet_type*,
1154  buffer_device_type>& exports,
1155  Kokkos::DualView<size_t*,
1157  numPacketsPerLID,
1158  size_t& constantNumPackets) {
1159  using std::endl;
1160  using Teuchos::Comm;
1161  using Teuchos::RCP;
1162  using this_COO_type = CooMatrix<SC, LO, GO, NT>;
1163  const char prefix[] = "Tpetra::Details::CooMatrix::packAndPrepare: ";
1164  const char suffix[] =
1165  " This should never happen. "
1166  "Please report this bug to the Tpetra developers.";
1167 
1168  // Tell the caller that different rows may have different numbers
1169  // of matrix entries.
1170  constantNumPackets = 0;
1171 
1172  const this_COO_type* src = dynamic_cast<const this_COO_type*>(&sourceObject);
1173  if (src == nullptr) {
1174  std::ostream& err = this->markLocalErrorAndGetStream();
1175  err << prefix << "Input argument 'sourceObject' is not a CooMatrix."
1176  << endl;
1177  } else if (*(src->localError_)) {
1178  std::ostream& err = this->markLocalErrorAndGetStream();
1179  err << prefix << "The source (input) object of the Import or Export "
1180  "is already in an error state on this process."
1181  << endl;
1182  } else if (*(this->localError_)) {
1183  std::ostream& err = this->markLocalErrorAndGetStream();
1184  err << prefix << "The target (output, \"this\") object of the Import "
1185  "or Export is already in an error state on this process."
1186  << endl;
1187  }
1188  // Respond to detected error(s) by resizing 'exports' to zero (so
1189  // we won't be tempted to read it later), and filling
1190  // numPacketsPerLID with zeros.
1191  if (*(this->localError_)) {
1192  // Resize 'exports' to zero, so we won't be tempted to read it.
1193  Details::reallocDualViewIfNeeded(exports, 0, "CooMatrix exports");
1194  // Trick to get around const DualView& being const.
1195  {
1196  auto numPacketsPerLID_tmp = numPacketsPerLID;
1197  numPacketsPerLID_tmp.sync_host();
1198  numPacketsPerLID_tmp.modify_host();
1199  }
1200  // Fill numPacketsPerLID with zeros.
1201  Kokkos::deep_copy(numPacketsPerLID.view_host(), static_cast<size_t>(0));
1202  return;
1203  }
1204 
1205  const size_t numExports = exportLIDs.extent(0);
1206  if (numExports == 0) {
1207  Details::reallocDualViewIfNeeded(exports, 0, exports.view_host().label());
1208  return; // nothing to send
1209  }
1210  RCP<const Comm<int> > comm = src->getMap().is_null() ? Teuchos::null : src->getMap()->getComm();
1211  if (comm.is_null() || comm->getSize() == 1) {
1212  if (numExports != static_cast<size_t>(0)) {
1213  std::ostream& err = this->markLocalErrorAndGetStream();
1214  err << prefix << "The input communicator is either null or "
1215  "has only one process, but numExports = "
1216  << numExports << " != 0."
1217  << suffix << endl;
1218  return;
1219  }
1220  // Don't go into the rest of this method unless there are
1221  // actually processes other than the calling process. This is
1222  // because the pack and unpack functions only have nonstub
1223  // implementations if building with MPI.
1224  return;
1225  }
1226 
1227  numPacketsPerLID.sync_host();
1228  numPacketsPerLID.modify_host();
1229 
1230  TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
1231  auto exportLIDs_h = exportLIDs.view_host();
1232 
1233  int totalNumPackets = 0;
1234  size_t numInvalidRowInds = 0;
1235  std::ostringstream errStrm; // current loop iteration's error messages
1236  for (size_t k = 0; k < numExports; ++k) {
1237  const LO lclRow = exportLIDs_h[k];
1238  // We're packing the source object's data, so we need to use the
1239  // source object's Map to convert from local to global indices.
1240  const GO gblRow = src->map_->getGlobalElement(lclRow);
1241  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1242  // Mark the error later; just count for now.
1243  ++numInvalidRowInds;
1244  numPacketsPerLID.view_host()[k] = 0;
1245  continue;
1246  }
1247 
1248  // Count the number of bytes needed to pack the current row of
1249  // the source object.
1250  int numPackets = 0;
1251  const int errCode =
1252  src->impl_.countPackRow(numPackets, gblRow, *comm, &errStrm);
1253  if (errCode != 0) {
1254  std::ostream& err = this->markLocalErrorAndGetStream();
1255  err << prefix << errStrm.str() << endl;
1256  numPacketsPerLID.view_host()[k] = 0;
1257  continue;
1258  }
1259 
1260  // Make sure that the total number of packets fits in int.
1261  // MPI requires this.
1262  const long long newTotalNumPackets =
1263  static_cast<long long>(totalNumPackets) +
1264  static_cast<long long>(numPackets);
1265  if (newTotalNumPackets >
1266  static_cast<long long>(std::numeric_limits<int>::max())) {
1267  std::ostream& err = this->markLocalErrorAndGetStream();
1268  err << prefix << "The new total number of packets "
1269  << newTotalNumPackets << " does not fit in int." << endl;
1270  // At this point, we definitely cannot continue. In order to
1271  // leave the output arguments in a rational state, we zero out
1272  // all remaining entries of numPacketsPerLID before returning.
1273  for (size_t k2 = k; k2 < numExports; ++k2) {
1274  numPacketsPerLID.view_host()[k2] = 0;
1275  }
1276  return;
1277  }
1278  numPacketsPerLID.view_host()[k] = static_cast<size_t>(numPackets);
1279  totalNumPackets = static_cast<int>(newTotalNumPackets);
1280  }
1281 
1282  // If we found invalid row indices in exportLIDs, go back,
1283  // collect, and report them.
1284  if (numInvalidRowInds != 0) {
1285  std::vector<std::pair<LO, GO> > invalidRowInds;
1286  for (size_t k = 0; k < numExports; ++k) {
1287  const LO lclRow = exportLIDs_h[k];
1288  // We're packing the source object's data, so we need to use
1289  // the source object's Map to convert from local to global
1290  // indices.
1291  const GO gblRow = src->map_->getGlobalElement(lclRow);
1292  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1293  invalidRowInds.push_back({lclRow, gblRow});
1294  }
1295  }
1296  std::ostringstream os;
1297  os << prefix << "We found " << numInvalidRowInds << " invalid row ind"
1298  << (numInvalidRowInds != static_cast<size_t>(1) ? "ices" : "ex")
1299  << " out of " << numExports << " in exportLIDs. Here is the list "
1300  << " of invalid row indices: [";
1301  for (size_t k = 0; k < invalidRowInds.size(); ++k) {
1302  os << "(LID: " << invalidRowInds[k].first << ", GID: "
1303  << invalidRowInds[k].second << ")";
1304  if (k + 1 < invalidRowInds.size()) {
1305  os << ", ";
1306  }
1307  }
1308  os << "].";
1309 
1310  std::ostream& err = this->markLocalErrorAndGetStream();
1311  err << prefix << os.str() << std::endl;
1312  return;
1313  }
1314 
1315  {
1316  const bool reallocated =
1317  Details::reallocDualViewIfNeeded(exports, totalNumPackets,
1318  "CooMatrix exports");
1319  if (reallocated) {
1320  exports.sync_host(); // make sure alloc'd on host
1321  }
1322  }
1323  exports.modify_host();
1324 
1325  // FIXME (mfh 17 Jan 2017) packTriples wants three arrays, not a
1326  // single array of structs. For now, we just copy.
1327  std::vector<GO> gblRowInds;
1328  std::vector<GO> gblColInds;
1329  std::vector<SC> vals;
1330 
1331  int outBufCurPos = 0;
1332  packet_type* outBuf = exports.view_host().data();
1333  for (size_t k = 0; k < numExports; ++k) {
1334  const LO lclRow = exportLIDs.view_host()[k];
1335  // We're packing the source object's data, so we need to use the
1336  // source object's Map to convert from local to global indices.
1337  const GO gblRow = src->map_->getGlobalElement(lclRow);
1338  // Pack the current row of the source object.
1339  src->impl_.packRow(outBuf, totalNumPackets, outBufCurPos, *comm,
1340  gblRowInds, gblColInds, vals, gblRow);
1341  }
1342  }
1343 
1344  virtual void
1345  unpackAndCombine(const Kokkos::DualView<const local_ordinal_type*,
1346  buffer_device_type>& importLIDs,
1347  Kokkos::DualView<packet_type*,
1349  imports,
1350  Kokkos::DualView<size_t*,
1352  numPacketsPerLID,
1353  const size_t /* constantNumPackets */,
1354  const ::Tpetra::CombineMode /* combineMode */) {
1355  using std::endl;
1356  using Teuchos::Comm;
1357  using Teuchos::RCP;
1358  const char prefix[] = "Tpetra::Details::CooMatrix::unpackAndCombine: ";
1359  const char suffix[] =
1360  " This should never happen. "
1361  "Please report this bug to the Tpetra developers.";
1362 
1363  TEUCHOS_ASSERT(!importLIDs.need_sync_host());
1364  auto importLIDs_h = importLIDs.view_host();
1365 
1366  const std::size_t numImports = importLIDs.extent(0);
1367  if (numImports == 0) {
1368  return; // nothing to receive
1369  } else if (imports.extent(0) == 0) {
1370  std::ostream& err = this->markLocalErrorAndGetStream();
1371  err << prefix << "importLIDs.extent(0) = " << numImports << " != 0, "
1372  << "but imports.extent(0) = 0. This doesn't make sense, because "
1373  << "for every incoming LID, CooMatrix packs at least the count of "
1374  << "triples associated with that LID, even if the count is zero. "
1375  << "importLIDs = [";
1376  for (std::size_t k = 0; k < numImports; ++k) {
1377  err << importLIDs_h[k];
1378  if (k + 1 < numImports) {
1379  err << ", ";
1380  }
1381  }
1382  err << "]. " << suffix << endl;
1383  return;
1384  }
1385 
1386  RCP<const Comm<int> > comm = this->getMap().is_null() ? Teuchos::null : this->getMap()->getComm();
1387  if (comm.is_null() || comm->getSize() == 1) {
1388  if (numImports != static_cast<size_t>(0)) {
1389  std::ostream& err = this->markLocalErrorAndGetStream();
1390  err << prefix << "The input communicator is either null or "
1391  "has only one process, but numImports = "
1392  << numImports << " != 0."
1393  << suffix << endl;
1394  return;
1395  }
1396  // Don't go into the rest of this method unless there are
1397  // actually processes other than the calling process. This is
1398  // because the pack and unpack functions only have nonstub
1399  // implementations if building with MPI.
1400  return;
1401  }
1402 
1403  // Make sure that the length of 'imports' fits in int.
1404  // This is ultimately an MPI restriction.
1405  if (static_cast<size_t>(imports.extent(0)) >
1406  static_cast<size_t>(std::numeric_limits<int>::max())) {
1407  std::ostream& err = this->markLocalErrorAndGetStream();
1408  err << prefix << "imports.extent(0) = "
1409  << imports.extent(0) << " does not fit in int." << endl;
1410  return;
1411  }
1412  const int inBufSize = static_cast<int>(imports.extent(0));
1413 
1414  if (imports.need_sync_host()) {
1415  imports.sync_host();
1416  }
1417  if (numPacketsPerLID.need_sync_host()) {
1418  numPacketsPerLID.sync_host();
1419  }
1420  auto imports_h = imports.view_host();
1421  auto numPacketsPerLID_h = numPacketsPerLID.view_host();
1422 
1423  // FIXME (mfh 17,24 Jan 2017) packTriples wants three arrays, not a
1424  // single array of structs. For now, we just copy.
1425  std::vector<GO> gblRowInds;
1426  std::vector<GO> gblColInds;
1427  std::vector<SC> vals;
1428 
1429  const packet_type* inBuf = imports_h.data();
1430  int inBufCurPos = 0;
1431  size_t numInvalidRowInds = 0;
1432  int errCode = 0;
1433  std::ostringstream errStrm; // for unpack* error output.
1434  for (size_t k = 0; k < numImports; ++k) {
1435  const LO lclRow = importLIDs_h(k);
1436  const GO gblRow = this->map_->getGlobalElement(lclRow);
1437  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1438  ++numInvalidRowInds;
1439  continue;
1440  }
1441 
1442  // Remember where we were, so we don't overrun the buffer
1443  // length. inBufCurPos is an in/out argument of unpackTriples*.
1444  const int origInBufCurPos = inBufCurPos;
1445 
1446  int numEnt = 0; // output argument of unpackTriplesCount
1447  errCode = unpackTriplesCount(inBuf, inBufSize, inBufCurPos,
1448  numEnt, *comm, &errStrm);
1449  if (errCode != 0 || numEnt < 0 || inBufCurPos < origInBufCurPos) {
1450  std::ostream& err = this->markLocalErrorAndGetStream();
1451 
1452  err << prefix << "In unpack loop, k=" << k << ": ";
1453  if (errCode != 0) {
1454  err << " unpackTriplesCount returned errCode = " << errCode
1455  << " != 0." << endl;
1456  }
1457  if (numEnt < 0) {
1458  err << " unpackTriplesCount returned errCode = 0, but numEnt = "
1459  << numEnt << " < 0." << endl;
1460  }
1461  if (inBufCurPos < origInBufCurPos) {
1462  err << " After unpackTriplesCount, inBufCurPos = " << inBufCurPos
1463  << " < origInBufCurPos = " << origInBufCurPos << "." << endl;
1464  }
1465  err << " unpackTriplesCount report: " << errStrm.str() << endl;
1466  err << suffix << endl;
1467 
1468  // We only continue in a debug build, because the above error
1469  // messages could consume too much memory and cause an
1470  // out-of-memory error, without actually printing. Printing
1471  // everything is useful in a debug build, but not so much in a
1472  // release build.
1473 #ifdef HAVE_TPETRA_DEBUG
1474  // Clear out the current error stream, so we don't accumulate
1475  // over loop iterations.
1476  errStrm.str("");
1477  continue;
1478 #else
1479  return;
1480 #endif // HAVE_TPETRA_DEBUG
1481  }
1482 
1483  // FIXME (mfh 17,24 Jan 2017) packTriples wants three arrays,
1484  // not a single array of structs. For now, we just copy.
1485  gblRowInds.resize(numEnt);
1486  gblColInds.resize(numEnt);
1487  vals.resize(numEnt);
1488 
1489  errCode = unpackTriples(inBuf, inBufSize, inBufCurPos,
1490  gblRowInds.data(), gblColInds.data(),
1491  vals.data(), numEnt, *comm, &errStrm);
1492  if (errCode != 0) {
1493  std::ostream& err = this->markLocalErrorAndGetStream();
1494  err << prefix << "unpackTriples returned errCode = "
1495  << errCode << " != 0. It reports: " << errStrm.str()
1496  << endl;
1497  // We only continue in a debug build, because the above error
1498  // messages could consume too much memory and cause an
1499  // out-of-memory error, without actually printing. Printing
1500  // everything is useful in a debug build, but not so much in a
1501  // release build.
1502 #ifdef HAVE_TPETRA_DEBUG
1503  // Clear out the current error stream, so we don't accumulate
1504  // over loop iterations.
1505  errStrm.str("");
1506  continue;
1507 #else
1508  return;
1509 #endif // HAVE_TPETRA_DEBUG
1510  }
1511  this->sumIntoGlobalValues(gblRowInds.data(), gblColInds.data(),
1512  vals.data(), numEnt);
1513  }
1514 
1515  // If we found invalid row indices in exportLIDs, go back,
1516  // collect, and report them.
1517  if (numInvalidRowInds != 0) {
1518  // Mark the error now, before we possibly run out of memory.
1519  // The latter could raise an exception (e.g., std::bad_alloc),
1520  // but at least we would get the error state right.
1521  std::ostream& err = this->markLocalErrorAndGetStream();
1522 
1523  std::vector<std::pair<LO, GO> > invalidRowInds;
1524  for (size_t k = 0; k < numImports; ++k) {
1525  const LO lclRow = importLIDs_h(k);
1526  const GO gblRow = this->map_->getGlobalElement(lclRow);
1527  if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid()) {
1528  invalidRowInds.push_back({lclRow, gblRow});
1529  }
1530  }
1531 
1532  err << prefix << "We found " << numInvalidRowInds << " invalid row ind"
1533  << (numInvalidRowInds != static_cast<size_t>(1) ? "ices" : "ex")
1534  << " out of " << numImports << " in importLIDs. Here is the list "
1535  << " of invalid row indices: [";
1536  for (size_t k = 0; k < invalidRowInds.size(); ++k) {
1537  err << "(LID: " << invalidRowInds[k].first << ", GID: "
1538  << invalidRowInds[k].second << ")";
1539  if (k + 1 < invalidRowInds.size()) {
1540  err << ", ";
1541  }
1542  }
1543  err << "].";
1544  return;
1545  }
1546  }
1547 };
1548 
1549 } // namespace Details
1550 } // namespace Tpetra
1551 
1552 #endif // TPETRA_DETAILS_COOMATRIX_HPP
char packet_type
Type for packing and unpacking data.
int packTriplesCount(const int, char[], const int, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm)
Pack the count (number) of matrix triples.
void sumIntoGlobalValues(std::initializer_list< GO > gblRowInds, std::initializer_list< GO > gblColInds, std::initializer_list< SC > vals, const std::size_t numEnt)
Initializer-list overload of the above method (which see).
Node node_type
The Node type. If you don&#39;t know what this is, don&#39;t use it.
CooMatrixImpl()=default
Default constructor.
std::string errorMessages() const
The current stream of error messages.
Declaration of a function that prints strings from each process.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a descriptiion of this object to the given output stream; overrides Teuchos::Describable method...
void packRow(packet_type outBuf[], const int outBufSize, int &outBufCurPos, const ::Teuchos::Comm< int > &comm, std::vector< GO > &gblRowInds, std::vector< GO > &gblColInds, std::vector< SC > &vals, const GO gblRow) const
Pack the given row of the matrix.
Function comparing two CooGraphEntry structs, lexicographically, first by row index, then by column index.
void sumIntoGlobalValue(const GO gblRowInd, const GO gblColInd, const SC &val)
Insert one entry locally into the sparse matrix, if it does not exist there yet. If it does exist...
Sparse matrix used only for file input / output.
void buildLocallyIndexedCrs(std::vector< OffsetType > &rowOffsets, LO lclColInds[], SC vals[], std::function< LO(const GO)> gblToLcl) const
Build a locally indexed version of CRS storage.
void fillComplete()
Special version of fillComplete that assumes that the matrix already has a Map, and reuses its commun...
GlobalOrdinal global_ordinal_type
The type of global indices.
virtual ~CooMatrix()
Destructor (virtual for memory safety of derived classes).
void gathervPrint(std::ostream &out, const std::string &s, const Teuchos::Comm< int > &comm)
On Process 0 in the given communicator, print strings from each process in that communicator, in rank order.
void fillComplete(const ::Teuchos::RCP< const ::Teuchos::Comm< int > > &comm)
Tell the matrix that you are done inserting entries locally, and that the matrix should build its Map...
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.
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling process.
void sumIntoGlobalValue(const GO gblRowInd, const GO gblColInd, const SC &val)
Insert one entry locally into the sparse matrix, if it does not exist there yet. If it does exist...
Type of each (row index, column index) pair in the Tpetra::Details::CooMatrix (see below)...
bool localError() const
Whether this object had an error on the calling process.
int unpackTriplesCount(const char[], const int, int &, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm)
Unpack just the count of triples from the given input buffer.
SC scalar_type
Type of each entry (value) in the sparse matrix.
virtual size_t constantNumberOfPackets() const
By returning 0, tell DistObject that this class may not necessarily have a constant number of &quot;packet...
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
CooMatrix(const ::Teuchos::RCP< const map_type > &map)
Constructor that takes a Map.
virtual bool checkSizes(const ::Tpetra::SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
CombineMode
Rule for combining data in an Import or Export.
Teuchos::RCP< const map_type > map_
The Map over which this object is distributed.
bool reallocDualViewIfNeeded(Kokkos::DualView< ValueType *, DeviceType > &dv, const size_t newSize, const char newLabel[], const size_t tooBigFactor=2, const bool needFenceBeforeRealloc=true)
Reallocate the DualView in/out argument, if needed.
int countPackTriplesCount(const ::Teuchos::Comm< int > &, int &size, std::ostream *errStrm)
Compute the buffer size required by packTriples for packing the number of matrix entries (&quot;triples&quot;)...
::Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Type of the Map specialization to give to the constructor.
Declaration and definition of Tpetra::Details::reallocDualViewIfNeeded, an implementation detail of T...
void forAllEntries(std::function< void(const GO, const GO, const SC &)> f) const
Execute the given function for all entries of the sparse matrix, sequentially (no thread parallelism)...
LocalOrdinal local_ordinal_type
The type of local indices.
int countPackRow(int &numPackets, const GO gblRow, const ::Teuchos::Comm< int > &comm, std::ostream *errStrm=NULL) const
Count the number of packets (bytes, in this case) needed to pack the given row of the matrix...
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling process.
GO getMyGlobalRowIndices(std::vector< GO > &rowInds) const
Get the global row indices on this process, sorted and made unique, and return the minimum global row...
typename::Tpetra::DistObject< char, LO, GO, NT >::buffer_device_type buffer_device_type
Kokkos::Device specialization for DistObject communication buffers.
void mergeIntoRow(const GO tgtGblRow, const CooMatrixImpl< SC, GO > &src, const GO srcGblRow)
Into global row tgtGblRow of *this, merge global row srcGblRow of src.
char packet_type
This class transfers data as bytes (MPI_BYTE).
virtual std::string description() const
One-line descriptiion of this object; overrides Teuchos::Describable method.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
int unpackTriples(const char[], const int, int &, OrdinalType[], OrdinalType[], ScalarType[], const int, const ::Teuchos::Comm< int > &, std::ostream *errStrm=NULL)
Unpack matrix entries (&quot;triples&quot; (i, j, A(i,j))) from the given input buffer.
Base class for distributed Tpetra objects that support data redistribution.
Implementation detail of Tpetra::Details::CooMatrix (which see below).
int packTriples(const OrdinalType[], const OrdinalType[], const ScalarType[], const int, char[], const int, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm=NULL)
Pack matrix entries (&quot;triples&quot; (i, j, A(i,j))) into the given output buffer.