Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_Details_CooMatrix.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Tpetra: Templated Linear Algebra Services Package
5// Copyright (2008) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40// @HEADER
41
42#ifndef TPETRA_DETAILS_COOMATRIX_HPP
43#define TPETRA_DETAILS_COOMATRIX_HPP
44
49
50#include "TpetraCore_config.h"
51#include "Tpetra_Details_PackTriples.hpp"
52#include "Tpetra_DistObject.hpp"
55#include "Teuchos_TypeNameTraits.hpp"
56
57#include <initializer_list>
58#include <map>
59#include <memory>
60#include <string>
61#include <type_traits>
62#include <vector>
63
64
65namespace Tpetra {
66namespace Details {
67
68// Implementation details of Tpetra::Details.
69// So, users REALLY should not use anything in here.
70namespace Impl {
71
74template<class IndexType>
76 IndexType row;
77 IndexType col;
78};
79
84template<class IndexType>
86 bool
88 const CooGraphEntry<IndexType>& y) const
89 {
90 return (x.row < y.row) || (x.row == y.row && x.col < y.col);
91 }
92};
93
96template<class SC, class GO>
98private:
106 typedef std::map<CooGraphEntry<GO>, SC,
107 CompareCooGraphEntries<GO> > entries_type;
108
111 entries_type entries_;
112
113public:
115 typedef char packet_type;
116
118 CooMatrixImpl () = default;
119
126 void
128 const GO gblColInd,
129 const SC& val)
130 {
131 // There's no sense in worrying about the insertion hint, since
132 // indices may be all over the place. Users make no promises of
133 // sorting or locality of input.
135 auto result = this->entries_.insert ({ij, val});
136 if (! result.second) { // already in the map
137 result.first->second += val; // sum in the new value
138 }
139 }
140
149 void
151 const GO gblColInds[],
152 const SC vals[],
153 const std::size_t numEnt)
154 {
155 for (std::size_t k = 0; k < numEnt; ++k) {
156 // There's no sense in worrying about the insertion hint, since
157 // indices may be all over the place. Users make no promises of
158 // sorting or locality of input.
160 const SC val = vals[k];
161 auto result = this->entries_.insert ({ij, val});
162 if (! result.second) { // already in the map
163 result.first->second += val; // sum in the new value
164 }
165 }
166 }
167
169 std::size_t
171 {
172 return this->entries_.size ();
173 }
174
178 void
179 forAllEntries (std::function<void (const GO, const GO, const SC&)> f) const
180 {
181 for (auto iter = this->entries_.begin ();
182 iter != this->entries_.end (); ++iter) {
183 f (iter->first.row, iter->first.col, iter->second);
184 }
185 }
186
192 void
194 const CooMatrixImpl<SC, GO>& src,
195 const GO srcGblRow)
196 {
197 // const char prefix[] =
198 // "Tpetra::Details::Impl::CooMatrixImpl::mergeIntoRow: ";
199
200 entries_type& tgtEntries = this->entries_;
201 const entries_type& srcEntries = src.entries_;
202
203 // Find all entries with the given global row index. The min GO
204 // value is guaranteed to be the least possible column index, so
205 // beg1 is a lower bound for all (row index, column index) pairs.
206 // Lower bound is inclusive; upper bound is exclusive.
207 auto srcBeg = srcEntries.lower_bound ({srcGblRow, std::numeric_limits<GO>::min ()});
208 auto srcEnd = srcEntries.upper_bound ({srcGblRow+1, std::numeric_limits<GO>::min ()});
209 auto tgtBeg = tgtEntries.lower_bound ({tgtGblRow, std::numeric_limits<GO>::min ()});
210 auto tgtEnd = tgtEntries.upper_bound ({tgtGblRow+1, std::numeric_limits<GO>::min ()});
211
212 // Don't waste time iterating over the current row of *this, if
213 // the current row of src is empty.
214 if (srcBeg != srcEnd) {
215 for (auto tgtCur = tgtBeg; tgtCur != tgtEnd; ++tgtCur) {
216 auto srcCur = srcBeg;
217 while (srcCur != srcEnd && srcCur->first.col < tgtCur->first.col) {
218 ++srcCur;
219 }
220 // At this point, one of the following is true:
221 //
222 // 1. srcCur == srcEnd, thus nothing more to insert.
223 // 2. srcCur->first.col > tgtCur->first.col, thus the current
224 // row of src has no entry matching tgtCur->first, so we
225 // have to insert it. Insertion does not invalidate
226 // tgtEntries iterators, and neither srcEntries nor
227 // tgtEntries have duplicates, so this is safe.
228 // 3. srcCur->first.col == tgtCur->first.col, so add the two entries.
229 if (srcCur != srcEnd) {
230 if (srcCur->first.col == tgtCur->first.col) {
231 tgtCur->second += srcCur->second;
232 }
233 else {
234 // tgtCur is (optimally) right before where we want to put
235 // the new entry, since srcCur points to the first entry
236 // in srcEntries whose column index is greater than that
237 // of the entry to which tgtCur points.
238 (void) tgtEntries.insert (tgtCur, *srcCur);
239 }
240 } // if srcCur != srcEnd
241 } // for each entry in the current row (lclRow) of *this
242 } // if srcBeg != srcEnd
243 }
244
254 int
256 const GO gblRow,
257 const ::Teuchos::Comm<int>& comm,
258 std::ostream* errStrm = NULL) const
259 {
260 using ::Tpetra::Details::countPackTriples;
261 using ::Tpetra::Details::countPackTriplesCount;
262 using std::endl;
263 const char prefix[] = "Tpetra::Details::Impl::CooMatrixImpl::countPackRow: ";
264#ifdef HAVE_TPETRACORE_MPI
265 int errCode = MPI_SUCCESS;
266#else
267 int errCode = 0;
268#endif // HAVE_TPETRACORE_MPI
269
270 // Count the number of entries in the given row.
271 const GO minGO = std::numeric_limits<GO>::min ();
272 auto beg = this->entries_.lower_bound ({gblRow, minGO});
273 auto end = this->entries_.upper_bound ({gblRow+1, minGO});
274 int numEnt = 0;
275 for (auto iter = beg; iter != end; ++iter) {
276 ++numEnt;
277 if (numEnt == std::numeric_limits<int>::max ()) {
278 if (errStrm != NULL) {
279 *errStrm << prefix << "In (global) row " << gblRow << ", the number "
280 "of entries thus far, numEnt = " << numEnt << ", has reached the "
281 "maximum int value. We need to store this count as int for MPI's "
282 "sake." << endl;
283 }
284 return -1;
285 }
286 }
287
288 int numPacketsForCount = 0; // output argument of countPackTriplesCount
289 {
290 errCode =
292 if (errCode != 0) {
293 if (errStrm != NULL) {
294 *errStrm << prefix << "countPackTriplesCount "
295 "returned errCode = " << errCode << " != 0." << endl;
296 }
297 return errCode;
298 }
299 if (numPacketsForCount < 0) {
300 if (errStrm != NULL) {
301 *errStrm << prefix << "countPackTriplesCount returned "
302 "numPacketsForCount = " << numPacketsForCount << " < 0." << endl;
303 }
304 return -1;
305 }
306 }
307
308 int numPacketsForTriples = 0; // output argument of countPackTriples
309 {
312 (errCode != 0, std::runtime_error, prefix << "countPackTriples "
313 "returned errCode = " << errCode << " != 0.");
315 (numPacketsForTriples < 0, std::logic_error, prefix << "countPackTriples "
316 "returned numPacketsForTriples = " << numPacketsForTriples << " < 0.");
317 }
318
320 return errCode;
321 }
322
337 void
339 const int outBufSize,
340 int& outBufCurPos, // in/out argument
341 const ::Teuchos::Comm<int>& comm,
342 std::vector<GO>& gblRowInds,
343 std::vector<GO>& gblColInds,
344 std::vector<SC>& vals,
345 const GO gblRow) const
346 {
347 using ::Tpetra::Details::packTriples;
348 using ::Tpetra::Details::packTriplesCount;
349 const char prefix[] = "Tpetra::Details::Impl::CooMatrixImpl::packRow: ";
350
351 const GO minGO = std::numeric_limits<GO>::min ();
352 auto beg = this->entries_.lower_bound ({gblRow, minGO});
353 auto end = this->entries_.upper_bound ({gblRow+1, minGO});
354
355 // This doesn't actually deallocate. Only swapping with an empty
356 // std::vector does that.
357 gblRowInds.resize (0);
358 gblColInds.resize (0);
359 vals.resize (0);
360
361 int numEnt = 0;
362 for (auto iter = beg; iter != end; ++iter) {
363 gblRowInds.push_back (iter->first.row);
364 gblColInds.push_back (iter->first.col);
365 vals.push_back (iter->second);
366 ++numEnt;
368 (numEnt == std::numeric_limits<int>::max (), std::runtime_error, prefix
369 << "In (global) row " << gblRow << ", the number of entries thus far, "
370 "numEnt = " << numEnt << ", has reached the maximum int value. "
371 "We need to store this count as int for MPI's sake.");
372 }
373
374 {
375 const int errCode =
378 (errCode != 0, std::runtime_error, prefix
379 << "In (global) row " << gblRow << ", packTriplesCount returned "
380 "errCode = " << errCode << " != 0.");
381 }
382 {
383 const int errCode =
384 packTriples (gblRowInds.data (),
385 gblColInds.data (),
386 vals.data (),
387 numEnt,
388 outBuf,
390 outBufCurPos, // in/out argument
391 comm);
393 (errCode != 0, std::runtime_error, prefix << "In (global) row "
394 << gblRow << ", packTriples returned errCode = " << errCode
395 << " != 0.");
396 }
397 }
398
406 GO
407 getMyGlobalRowIndices (std::vector<GO>& rowInds) const
408 {
409 rowInds.clear ();
410
411 GO lclMinRowInd = std::numeric_limits<GO>::max (); // compute local min
412 for (typename entries_type::const_iterator iter = this->entries_.begin ();
413 iter != this->entries_.end (); ++iter) {
414 const GO rowInd = iter->first.row;
415 if (rowInd < lclMinRowInd) {
417 }
418 if (rowInds.empty () || rowInds.back () != rowInd) {
419 rowInds.push_back (rowInd); // don't insert duplicates
420 }
421 }
422 return lclMinRowInd;
423 }
424
425 template<class OffsetType>
426 void
427 buildCrs (std::vector<OffsetType>& rowOffsets,
428 GO gblColInds[],
429 SC vals[]) const
430 {
431 static_assert (std::is_integral<OffsetType>::value,
432 "OffsetType must be a built-in integer type.");
433
434 // clear() doesn't generally free storage; it just resets the
435 // length. Thus, we reuse any existing storage here.
436 rowOffsets.clear ();
437
438 const std::size_t numEnt = this->getLclNumEntries ();
439 if (numEnt == 0) {
440 rowOffsets.push_back (0);
441 }
442 else {
443 typename entries_type::const_iterator iter = this->entries_.begin ();
444 GO prevGblRowInd = iter->first.row;
445
446 OffsetType k = 0;
447 for ( ; iter != this->entries_.end (); ++iter, ++k) {
448 const GO gblRowInd = iter->first.row;
449 if (k == 0 || gblRowInd != prevGblRowInd) {
450 // The row offsets array always has at least one entry. The
451 // first entry is always zero, and the last entry is always
452 // the number of matrix entries.
453 rowOffsets.push_back (k);
455 }
456 gblColInds[k] = iter->first.col;
457
458 static_assert (std::is_same<typename std::decay<decltype (iter->second)>::type, SC>::value,
459 "The type of iter->second != SC.");
460 vals[k] = iter->second;
461 }
462 rowOffsets.push_back (static_cast<OffsetType> (numEnt));
463 }
464 }
465
478 template<class OffsetType, class LO>
479 void
480 buildLocallyIndexedCrs (std::vector<OffsetType>& rowOffsets,
481 LO lclColInds[],
482 SC vals[],
483 std::function<LO (const GO)> gblToLcl) const
484 {
485 static_assert (std::is_integral<OffsetType>::value,
486 "OffsetType must be a built-in integer type.");
487 static_assert (std::is_integral<LO>::value,
488 "LO must be a built-in integer type.");
489
490 // clear() doesn't generally free storage; it just resets the
491 // length. Thus, we reuse any existing storage here.
492 rowOffsets.clear ();
493
494 const std::size_t numEnt = this->getLclNumEntries ();
495 if (numEnt == 0) {
496 rowOffsets.push_back (0);
497 }
498 else {
499 typename entries_type::const_iterator iter = this->entries_.begin ();
500 GO prevGblRowInd = iter->first.row;
501
502 OffsetType k = 0;
503 for ( ; iter != this->entries_.end (); ++iter, ++k) {
504 const GO gblRowInd = iter->first.row;
505 if (k == 0 || gblRowInd != prevGblRowInd) {
506 // The row offsets array always has at least one entry. The
507 // first entry is always zero, and the last entry is always
508 // the number of matrix entries.
509 rowOffsets.push_back (k);
511 }
512 lclColInds[k] = gblToLcl (iter->first.col);
513 vals[k] = iter->second;
514 }
515 rowOffsets.push_back (static_cast<OffsetType> (numEnt));
516 }
517 }
518};
519
520} // namespace Impl
521
570template<class SC,
574class CooMatrix : public ::Tpetra::DistObject<char, LO, GO, NT> {
575public:
577 typedef char packet_type;
580 typedef LO local_ordinal_type;
581 typedef GO global_ordinal_type;
582 typedef NT node_type;
583 typedef typename NT::device_type device_type;
585 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
586
587private:
589 typedef ::Tpetra::DistObject<packet_type, local_ordinal_type,
590 global_ordinal_type, node_type> base_type;
591
594
595public:
602 base_type (::Teuchos::null),
603 localError_ (new bool (false)),
604 errs_ (new std::shared_ptr<std::ostringstream> ()) // ptr to a null ptr
605 {}
606
614 CooMatrix (const ::Teuchos::RCP<const map_type>& map) :
615 base_type (map),
616 localError_ (new bool (false)),
617 errs_ (new std::shared_ptr<std::ostringstream> ()) // ptr to a null ptr
618 {}
619
621 virtual ~CooMatrix () {}
622
629 void
631 const GO gblColInd,
632 const SC& val)
633 {
635 }
636
645 void
647 const GO gblColInds[],
648 const SC vals[],
649 const std::size_t numEnt)
650 {
652 }
653
655 void
656 sumIntoGlobalValues (std::initializer_list<GO> gblRowInds,
657 std::initializer_list<GO> gblColInds,
658 std::initializer_list<SC> vals,
659 const std::size_t numEnt)
660 {
661 this->impl_.sumIntoGlobalValues (gblRowInds.begin (), gblColInds.begin (),
662 vals.begin (), numEnt);
663 }
664
666 std::size_t
668 {
669 return this->impl_.getLclNumEntries ();
670 }
671
672 template<class OffsetType>
673 void
674 buildCrs (::Kokkos::View<OffsetType*, device_type>& rowOffsets,
675 ::Kokkos::View<GO*, device_type>& gblColInds,
676 ::Kokkos::View<typename ::Kokkos::Details::ArithTraits<SC>::val_type*, device_type>& vals)
677 {
678 static_assert (std::is_integral<OffsetType>::value,
679 "OffsetType must be a built-in integer type.");
680 using ::Kokkos::create_mirror_view;
681 using ::Kokkos::deep_copy;
682 using ::Kokkos::View;
683 typedef typename ::Kokkos::Details::ArithTraits<SC>::val_type ISC;
684
685 const std::size_t numEnt = this->getLclNumEntries ();
686
687 gblColInds = View<GO*, device_type> ("gblColInds", numEnt);
689 auto gblColInds_h = create_mirror_view (gblColInds);
690 auto vals_h = create_mirror_view (vals);
691
692 std::vector<std::size_t> rowOffsetsSV;
693 this->impl_.buildCrs (rowOffsetsSV,
694 gblColInds_h.data (),
695 vals_h.data ());
696 rowOffsets =
697 View<OffsetType*, device_type> ("rowOffsets", rowOffsetsSV.size ());
699 rowOffsets_h (rowOffsetsSV.data (), rowOffsetsSV.size ());
701
704 }
705
716 void
717 fillComplete (const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm)
718 {
719 if (comm.is_null ()) {
720 this->map_ = ::Teuchos::null;
721 }
722 else {
723 this->map_ = this->buildMap (comm);
724 }
725 }
726
735 void
737 {
739 (this->getMap ().is_null (), std::runtime_error, "Tpetra::Details::"
740 "CooMatrix::fillComplete: This object does not yet have a Map. "
741 "You must call the version of fillComplete "
742 "that takes an input communicator.");
743 this->fillComplete (this->getMap ()->getComm ());
744 }
745
760 bool localError () const {
761 return *localError_;
762 }
763
781 std::string errorMessages () const {
782 return ((*errs_).get () == NULL) ? std::string ("") : (*errs_)->str ();
783 }
784
785private:
798 std::shared_ptr<bool> localError_;
799
807 std::shared_ptr<std::shared_ptr<std::ostringstream> > errs_;
808
810 std::ostream&
811 markLocalErrorAndGetStream ()
812 {
813 * (this->localError_) = true;
814 if ((*errs_).get () == NULL) {
815 *errs_ = std::shared_ptr<std::ostringstream> (new std::ostringstream ());
816 }
817 return **errs_;
818 }
819
820public:
823 virtual std::string description () const {
824 using Teuchos::TypeNameTraits;
825
826 std::ostringstream os;
827 os << "\"Tpetra::Details::CooMatrix\": {"
828 << "SC: " << TypeNameTraits<SC>::name ()
829 << ", LO: " << TypeNameTraits<LO>::name ()
830 << ", GO: " << TypeNameTraits<GO>::name ()
831 << ", NT: " << TypeNameTraits<NT>::name ();
832 if (this->getObjectLabel () != "") {
833 os << ", Label: \"" << this->getObjectLabel () << "\"";
834 }
835 os << ", Has Map: " << (this->map_.is_null () ? "false" : "true")
836 << "}";
837 return os.str ();
838 }
839
842 virtual void
843 describe (Teuchos::FancyOStream& out,
844 const Teuchos::EVerbosityLevel verbLevel =
845 Teuchos::Describable::verbLevel_default) const
846 {
847 using ::Tpetra::Details::gathervPrint;
848 using ::Teuchos::EVerbosityLevel;
849 using ::Teuchos::OSTab;
850 using ::Teuchos::TypeNameTraits;
851 using ::Teuchos::VERB_DEFAULT;
852 using ::Teuchos::VERB_LOW;
853 using ::Teuchos::VERB_MEDIUM;
854 using std::endl;
855
856 const auto vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
857
858 auto comm = this->getMap ().is_null () ?
859 Teuchos::null : this->getMap ()->getComm ();
860 const int myRank = comm.is_null () ? 0 : comm->getRank ();
861 // const int numProcs = comm.is_null () ? 1 : comm->getSize ();
862
863 if (vl != Teuchos::VERB_NONE) {
864 // Convention is to start describe() implementations with a tab.
865 OSTab tab0 (out);
866 if (myRank == 0) {
867 out << "\"Tpetra::Details::CooMatrix\":" << endl;
868 }
869 OSTab tab1 (out);
870 if (myRank == 0) {
871 out << "Template parameters:" << endl;
872 {
873 OSTab tab2 (out);
874 out << "SC: " << TypeNameTraits<SC>::name () << endl
875 << "LO: " << TypeNameTraits<LO>::name () << endl
876 << "GO: " << TypeNameTraits<GO>::name () << endl
877 << "NT: " << TypeNameTraits<NT>::name () << endl;
878 }
879 if (this->getObjectLabel () != "") {
880 out << "Label: \"" << this->getObjectLabel () << "\"" << endl;
881 }
882 out << "Has Map: " << (this->map_.is_null () ? "false" : "true") << endl;
883 } // if myRank == 0
884
885 // Describe the Map, if it exists.
886 if (! this->map_.is_null ()) {
887 if (myRank == 0) {
888 out << "Map:" << endl;
889 }
890 OSTab tab2 (out);
891 this->map_->describe (out, vl);
892 }
893
894 // At verbosity > VERB_LOW, each process prints something.
895 if (vl > VERB_LOW) {
896 std::ostringstream os;
897 os << "Process " << myRank << ":" << endl;
898 //OSTab tab3 (os);
899
900 const std::size_t numEnt = this->impl_.getLclNumEntries ();
901 os << " Local number of entries: " << numEnt << endl;
902 if (vl > VERB_MEDIUM) {
903 os << " Local entries:" << endl;
904 //OSTab tab4 (os);
905 this->impl_.forAllEntries ([&os] (const GO row, const GO col, const SC& val) {
906 os << " {"
907 << "row: " << row
908 << ", col: " << col
909 << ", val: " << val
910 << "}" << endl;
911 });
912 }
913 gathervPrint (out, os.str (), *comm);
914 }
915 } // vl != Teuchos::VERB_NONE
916 }
917
918private:
927 Teuchos::RCP<const map_type>
928 buildMap (const ::Teuchos::RCP<const ::Teuchos::Comm<int> >& comm)
929 {
930 using ::Teuchos::outArg;
931 using ::Teuchos::rcp;
932 using ::Teuchos::REDUCE_MIN;
933 using ::Teuchos::reduceAll;
935 //const char prefix[] = "Tpetra::Details::CooMatrix::buildMap: ";
936
937 // Processes where comm is null don't participate in the Map.
938 if (comm.is_null ()) {
939 return ::Teuchos::null;
940 }
941
942 // mfh 17 Jan 2017: We just happen to use row indices, because
943 // that's what Tpetra::CrsMatrix currently uses. That's probably
944 // not the best thing to use, but it's not bad for commonly
945 // encountered matrices. A better more general approach might be
946 // to hash (row index, column index) pairs to a global index. One
947 // could make that unique by doing a parallel scan at map
948 // construction time.
949
950 std::vector<GO> rowInds;
951 const GO lclMinGblRowInd = this->impl_.getMyGlobalRowIndices (rowInds);
952
953 // Compute the globally min row index for the "index base."
954 GO gblMinGblRowInd = 0; // output argument
957 const GO indexBase = gblMinGblRowInd;
958 const GST INV = Tpetra::Details::OrdinalTraits<GST>::invalid ();
959 return rcp (new map_type (INV, rowInds.data (), rowInds.size (),
960 indexBase, comm));
961 }
962
963protected:
967 virtual size_t constantNumberOfPackets () const {
968 return static_cast<size_t> (0);
969 }
970
974 virtual bool
975 checkSizes (const ::Tpetra::SrcDistObject& source)
976 {
977 using std::endl;
979 const char prefix[] = "Tpetra::Details::CooMatrix::checkSizes: ";
980
981 const this_COO_type* src = dynamic_cast<const this_COO_type* > (&source);
982
983 if (src == NULL) {
984 std::ostream& err = markLocalErrorAndGetStream ();
985 err << prefix << "The source object of the Import or Export "
986 "must be a CooMatrix with the same template parameters as the "
987 "target object." << endl;
988 }
989 else if (this->map_.is_null ()) {
990 std::ostream& err = markLocalErrorAndGetStream ();
991 err << prefix << "The target object of the Import or Export "
992 "must be a CooMatrix with a nonnull Map." << endl;
993 }
994 return ! (* (this->localError_));
995 }
996
999 typename ::Tpetra::DistObject<char, LO, GO, NT>::buffer_device_type;
1000
1001 virtual void
1002 copyAndPermute
1003 (const ::Tpetra::SrcDistObject& sourceObject,
1004 const size_t numSameIDs,
1005 const Kokkos::DualView<const LO*,
1007 const Kokkos::DualView<const LO*,
1009 const CombineMode /* CM */)
1010 {
1011 using std::endl;
1013 const char prefix[] = "Tpetra::Details::CooMatrix::copyAndPermute: ";
1014
1015 // There's no communication in this method, so it's safe just to
1016 // return on error.
1017
1018 if (* (this->localError_)) {
1019 std::ostream& err = this->markLocalErrorAndGetStream ();
1020 err << prefix << "The target object of the Import or Export is "
1021 "already in an error state." << endl;
1022 return;
1023 }
1024
1025 const this_COO_type* src = dynamic_cast<const this_COO_type*> (&sourceObject);
1026 if (src == nullptr) {
1027 std::ostream& err = this->markLocalErrorAndGetStream ();
1028 err << prefix << "Input argument 'sourceObject' is not a CooMatrix."
1029 << endl;
1030 return;
1031 }
1032
1033 const size_t numPermuteIDs =
1034 static_cast<size_t> (permuteToLIDs.extent (0));
1035 if (numPermuteIDs != static_cast<size_t> (permuteFromLIDs.extent (0))) {
1036 std::ostream& err = this->markLocalErrorAndGetStream ();
1037 err << prefix << "permuteToLIDs.extent(0) = "
1038 << numPermuteIDs << " != permuteFromLIDs.extent(0) = "
1039 << permuteFromLIDs.extent (0) << "." << endl;
1040 return;
1041 }
1042 if (sizeof (int) <= sizeof (size_t) &&
1043 numPermuteIDs > static_cast<size_t> (std::numeric_limits<int>::max ())) {
1044 std::ostream& err = this->markLocalErrorAndGetStream ();
1045 err << prefix << "numPermuteIDs = " << numPermuteIDs
1046 << ", a size_t, overflows int." << endl;
1047 return;
1048 }
1049
1050 // Even though this is an std::set, we can start with numSameIDs
1051 // just by iterating through the first entries of the set.
1052
1053 if (sizeof (int) <= sizeof (size_t) &&
1054 numSameIDs > static_cast<size_t> (std::numeric_limits<int>::max ())) {
1055 std::ostream& err = this->markLocalErrorAndGetStream ();
1056 err << prefix << "numSameIDs = " << numSameIDs
1057 << ", a size_t, overflows int." << endl;
1058 return;
1059 }
1060
1061 //
1062 // Copy in entries from any initial rows with the same global row indices.
1063 //
1064 const LO numSame = static_cast<int> (numSameIDs);
1065 // Count of local row indices encountered here with invalid global
1066 // row indices. If nonzero, something went wrong. If something
1067 // did go wrong, we'll defer responding until the end of this
1068 // method, so we can print as much useful info as possible.
1069 LO numInvalidSameRows = 0;
1070 for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1071 // All numSame initial rows have the same global index in both
1072 // source and target Maps, so we only need to convert to global
1073 // once.
1074 const GO gblRow = this->map_->getGlobalElement (lclRow);
1075 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1077 continue;
1078 }
1079 else {
1080 this->impl_.mergeIntoRow (gblRow, src->impl_, gblRow);
1081 }
1082 }
1083
1084 //
1085 // Copy in entries from remaining rows that are permutations, that
1086 // is, that live in both the source and target Maps, but aren't
1087 // included in the "same" list (see above).
1088 //
1089 const LO numPermute = static_cast<int> (numPermuteIDs);
1090 // Count of local "from" row indices encountered here with invalid
1091 // global row indices. If nonzero, something went wrong. If
1092 // something did go wrong, we'll defer responding until the end of
1093 // this method, so we can print as much useful info as possible.
1094 LO numInvalidRowsFrom = 0;
1095 // Count of local "to" row indices encountered here with invalid
1096 // global row indices. If nonzero, something went wrong. If
1097 // something did go wrong, we'll defer responding until the end of
1098 // this method, so we can print as much useful info as possible.
1099 LO numInvalidRowsTo = 0;
1100
1101 TEUCHOS_ASSERT( ! permuteFromLIDs.need_sync_host () );
1102 TEUCHOS_ASSERT( ! permuteToLIDs.need_sync_host () );
1103 auto permuteFromLIDs_h = permuteFromLIDs.view_host ();
1104 auto permuteToLIDs_h = permuteToLIDs.view_host ();
1105
1106 for (LO k = 0; k < numPermute; ++k) {
1107 const LO lclRowFrom = permuteFromLIDs_h[k];
1108 const LO lclRowTo = permuteToLIDs_h[k];
1109 const GO gblRowFrom = src->map_->getGlobalElement (lclRowFrom);
1110 const GO gblRowTo = this->map_->getGlobalElement (lclRowTo);
1111
1112 bool bothConversionsValid = true;
1113 if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1115 bothConversionsValid = false;
1116 }
1117 if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1119 bothConversionsValid = false;
1120 }
1122 this->impl_.mergeIntoRow (gblRowTo, src->impl_, gblRowFrom);
1123 }
1124 }
1125
1126 // Print info if any errors occurred.
1127 if (numInvalidSameRows != 0 || numInvalidRowsFrom != 0 ||
1128 numInvalidRowsTo != 0) {
1129 // Collect and print all the invalid input row indices, for the
1130 // "same," "from," and "to" lists.
1131 std::vector<std::pair<LO, GO> > invalidSameRows;
1133 std::vector<std::pair<LO, GO> > invalidRowsFrom;
1135 std::vector<std::pair<LO, GO> > invalidRowsTo;
1137
1138 for (LO lclRow = 0; lclRow < numSame; ++lclRow) {
1139 // All numSame initial rows have the same global index in both
1140 // source and target Maps, so we only need to convert to global
1141 // once.
1142 const GO gblRow = this->map_->getGlobalElement (lclRow);
1143 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1144 invalidSameRows.push_back ({lclRow, gblRow});
1145 }
1146 }
1147
1148 for (LO k = 0; k < numPermute; ++k) {
1149 const LO lclRowFrom = permuteFromLIDs_h[k];
1150 const LO lclRowTo = permuteToLIDs_h[k];
1151 const GO gblRowFrom = src->map_->getGlobalElement (lclRowFrom);
1152 const GO gblRowTo = this->map_->getGlobalElement (lclRowTo);
1153
1154 if (gblRowFrom == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1155 invalidRowsFrom.push_back ({lclRowFrom, gblRowFrom});
1156 }
1157 if (gblRowTo == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1158 invalidRowsTo.push_back ({lclRowTo, gblRowTo});
1159 }
1160 }
1161
1162 std::ostringstream os;
1163 if (numInvalidSameRows != 0) {
1164 os << "Invalid permute \"same\" (local, global) index pairs: [";
1165 for (std::size_t k = 0; k < invalidSameRows.size (); ++k) {
1166 const auto& p = invalidSameRows[k];
1167 os << "(" << p.first << "," << p.second << ")";
1168 if (k + 1 < invalidSameRows.size ()) {
1169 os << ", ";
1170 }
1171 }
1172 os << "]" << std::endl;
1173 }
1174 if (numInvalidRowsFrom != 0) {
1175 os << "Invalid permute \"from\" (local, global) index pairs: [";
1176 for (std::size_t k = 0; k < invalidRowsFrom.size (); ++k) {
1177 const auto& p = invalidRowsFrom[k];
1178 os << "(" << p.first << "," << p.second << ")";
1179 if (k + 1 < invalidRowsFrom.size ()) {
1180 os << ", ";
1181 }
1182 }
1183 os << "]" << std::endl;
1184 }
1185 if (numInvalidRowsTo != 0) {
1186 os << "Invalid permute \"to\" (local, global) index pairs: [";
1187 for (std::size_t k = 0; k < invalidRowsTo.size (); ++k) {
1188 const auto& p = invalidRowsTo[k];
1189 os << "(" << p.first << "," << p.second << ")";
1190 if (k + 1 < invalidRowsTo.size ()) {
1191 os << ", ";
1192 }
1193 }
1194 os << "]" << std::endl;
1195 }
1196
1197 std::ostream& err = this->markLocalErrorAndGetStream ();
1198 err << prefix << os.str ();
1199 return;
1200 }
1201 }
1202
1203 virtual void
1204 packAndPrepare
1205 (const ::Tpetra::SrcDistObject& sourceObject,
1206 const Kokkos::DualView<const local_ordinal_type*,
1208 Kokkos::DualView<packet_type*,
1209 buffer_device_type>& exports,
1210 Kokkos::DualView<size_t*,
1212 size_t& constantNumPackets)
1213 {
1214 using Teuchos::Comm;
1215 using Teuchos::RCP;
1216 using std::endl;
1218 const char prefix[] = "Tpetra::Details::CooMatrix::packAndPrepare: ";
1219 const char suffix[] = " This should never happen. "
1220 "Please report this bug to the Tpetra developers.";
1221
1222 // Tell the caller that different rows may have different numbers
1223 // of matrix entries.
1225
1226 const this_COO_type* src = dynamic_cast<const this_COO_type*> (&sourceObject);
1227 if (src == nullptr) {
1228 std::ostream& err = this->markLocalErrorAndGetStream ();
1229 err << prefix << "Input argument 'sourceObject' is not a CooMatrix."
1230 << endl;
1231 }
1232 else if (* (src->localError_)) {
1233 std::ostream& err = this->markLocalErrorAndGetStream ();
1234 err << prefix << "The source (input) object of the Import or Export "
1235 "is already in an error state on this process."
1236 << endl;
1237 }
1238 else if (* (this->localError_)) {
1239 std::ostream& err = this->markLocalErrorAndGetStream ();
1240 err << prefix << "The target (output, \"this\") object of the Import "
1241 "or Export is already in an error state on this process." << endl;
1242 }
1243 // Respond to detected error(s) by resizing 'exports' to zero (so
1244 // we won't be tempted to read it later), and filling
1245 // numPacketsPerLID with zeros.
1246 if (* (this->localError_)) {
1247 // Resize 'exports' to zero, so we won't be tempted to read it.
1248 Details::reallocDualViewIfNeeded (exports, 0, "CooMatrix exports");
1249 // Trick to get around const DualView& being const.
1250 {
1251 auto numPacketsPerLID_tmp = numPacketsPerLID;
1252 numPacketsPerLID_tmp.sync_host ();
1253 numPacketsPerLID_tmp.modify_host ();
1254 }
1255 // Fill numPacketsPerLID with zeros.
1256 Kokkos::deep_copy (numPacketsPerLID.h_view, static_cast<size_t> (0));
1257 return;
1258 }
1259
1260 const size_t numExports = exportLIDs.extent (0);
1261 if (numExports == 0) {
1262 Details::reallocDualViewIfNeeded (exports, 0, exports.h_view.label ());
1263 return; // nothing to send
1264 }
1265 RCP<const Comm<int> > comm = src->getMap ().is_null () ?
1266 Teuchos::null : src->getMap ()->getComm ();
1267 if (comm.is_null () || comm->getSize () == 1) {
1268 if (numExports != static_cast<size_t> (0)) {
1269 std::ostream& err = this->markLocalErrorAndGetStream ();
1270 err << prefix << "The input communicator is either null or "
1271 "has only one process, but numExports = " << numExports << " != 0."
1272 << suffix << endl;
1273 return;
1274 }
1275 // Don't go into the rest of this method unless there are
1276 // actually processes other than the calling process. This is
1277 // because the pack and unpack functions only have nonstub
1278 // implementations if building with MPI.
1279 return;
1280 }
1281
1282 numPacketsPerLID.sync_host ();
1283 numPacketsPerLID.modify_host ();
1284
1285 TEUCHOS_ASSERT( ! exportLIDs.need_sync_host () );
1286 auto exportLIDs_h = exportLIDs.view_host ();
1287
1288 int totalNumPackets = 0;
1289 size_t numInvalidRowInds = 0;
1290 std::ostringstream errStrm; // current loop iteration's error messages
1291 for (size_t k = 0; k < numExports; ++k) {
1292 const LO lclRow = exportLIDs_h[k];
1293 // We're packing the source object's data, so we need to use the
1294 // source object's Map to convert from local to global indices.
1295 const GO gblRow = src->map_->getGlobalElement (lclRow);
1296 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1297 // Mark the error later; just count for now.
1298 ++numInvalidRowInds;
1299 numPacketsPerLID.h_view[k] = 0;
1300 continue;
1301 }
1302
1303 // Count the number of bytes needed to pack the current row of
1304 // the source object.
1305 int numPackets = 0;
1306 const int errCode =
1307 src->impl_.countPackRow (numPackets, gblRow, *comm, &errStrm);
1308 if (errCode != 0) {
1309 std::ostream& err = this->markLocalErrorAndGetStream ();
1310 err << prefix << errStrm.str () << endl;
1311 numPacketsPerLID.h_view[k] = 0;
1312 continue;
1313 }
1314
1315 // Make sure that the total number of packets fits in int.
1316 // MPI requires this.
1317 const long long newTotalNumPackets =
1318 static_cast<long long> (totalNumPackets) +
1319 static_cast<long long> (numPackets);
1320 if (newTotalNumPackets >
1321 static_cast<long long> (std::numeric_limits<int>::max ())) {
1322 std::ostream& err = this->markLocalErrorAndGetStream ();
1323 err << prefix << "The new total number of packets "
1324 << newTotalNumPackets << " does not fit in int." << endl;
1325 // At this point, we definitely cannot continue. In order to
1326 // leave the output arguments in a rational state, we zero out
1327 // all remaining entries of numPacketsPerLID before returning.
1328 for (size_t k2 = k; k2 < numExports; ++k2) {
1329 numPacketsPerLID.h_view[k2] = 0;
1330 }
1331 return;
1332 }
1333 numPacketsPerLID.h_view[k] = static_cast<size_t> (numPackets);
1334 totalNumPackets = static_cast<int> (newTotalNumPackets);
1335 }
1336
1337 // If we found invalid row indices in exportLIDs, go back,
1338 // collect, and report them.
1339 if (numInvalidRowInds != 0) {
1340 std::vector<std::pair<LO, GO> > invalidRowInds;
1341 for (size_t k = 0; k < numExports; ++k) {
1342 const LO lclRow = exportLIDs_h[k];
1343 // We're packing the source object's data, so we need to use
1344 // the source object's Map to convert from local to global
1345 // indices.
1346 const GO gblRow = src->map_->getGlobalElement (lclRow);
1347 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1348 invalidRowInds.push_back ({lclRow, gblRow});
1349 }
1350 }
1351 std::ostringstream os;
1352 os << prefix << "We found " << numInvalidRowInds << " invalid row ind"
1353 << (numInvalidRowInds != static_cast<size_t> (1) ? "ices" : "ex")
1354 << " out of " << numExports << " in exportLIDs. Here is the list "
1355 << " of invalid row indices: [";
1356 for (size_t k = 0; k < invalidRowInds.size (); ++k) {
1357 os << "(LID: " << invalidRowInds[k].first << ", GID: "
1358 << invalidRowInds[k].second << ")";
1359 if (k + 1 < invalidRowInds.size ()) {
1360 os << ", ";
1361 }
1362 }
1363 os << "].";
1364
1365 std::ostream& err = this->markLocalErrorAndGetStream ();
1366 err << prefix << os.str () << std::endl;
1367 return;
1368 }
1369
1370 {
1371 const bool reallocated =
1372 Details::reallocDualViewIfNeeded (exports, totalNumPackets,
1373 "CooMatrix exports");
1374 if (reallocated) {
1375 exports.sync_host (); // make sure alloc'd on host
1376 }
1377 }
1378 exports.modify_host ();
1379
1380 // FIXME (mfh 17 Jan 2017) packTriples wants three arrays, not a
1381 // single array of structs. For now, we just copy.
1382 std::vector<GO> gblRowInds;
1383 std::vector<GO> gblColInds;
1384 std::vector<SC> vals;
1385
1386 int outBufCurPos = 0;
1387 packet_type* outBuf = exports.h_view.data ();
1388 for (size_t k = 0; k < numExports; ++k) {
1389 const LO lclRow = exportLIDs.h_view[k];
1390 // We're packing the source object's data, so we need to use the
1391 // source object's Map to convert from local to global indices.
1392 const GO gblRow = src->map_->getGlobalElement (lclRow);
1393 // Pack the current row of the source object.
1394 src->impl_.packRow (outBuf, totalNumPackets, outBufCurPos, *comm,
1395 gblRowInds, gblColInds, vals, gblRow);
1396 }
1397 }
1398
1399 virtual void
1400 unpackAndCombine
1401 (const Kokkos::DualView<const local_ordinal_type*,
1402 buffer_device_type>& importLIDs,
1403 Kokkos::DualView<packet_type*,
1404 buffer_device_type> imports,
1405 Kokkos::DualView<size_t*,
1406 buffer_device_type> numPacketsPerLID,
1407 const size_t /* constantNumPackets */,
1408 const ::Tpetra::CombineMode /* combineMode */)
1409 {
1410 using Teuchos::Comm;
1411 using Teuchos::RCP;
1412 using std::endl;
1413 const char prefix[] = "Tpetra::Details::CooMatrix::unpackAndCombine: ";
1414 const char suffix[] = " This should never happen. "
1415 "Please report this bug to the Tpetra developers.";
1416
1417 TEUCHOS_ASSERT( ! importLIDs.need_sync_host () );
1418 auto importLIDs_h = importLIDs.view_host ();
1419
1420 const std::size_t numImports = importLIDs.extent (0);
1421 if (numImports == 0) {
1422 return; // nothing to receive
1423 }
1424 else if (imports.extent (0) == 0) {
1425 std::ostream& err = this->markLocalErrorAndGetStream ();
1426 err << prefix << "importLIDs.extent(0) = " << numImports << " != 0, "
1427 << "but imports.extent(0) = 0. This doesn't make sense, because "
1428 << "for every incoming LID, CooMatrix packs at least the count of "
1429 << "triples associated with that LID, even if the count is zero. "
1430 << "importLIDs = [";
1431 for (std::size_t k = 0; k < numImports; ++k) {
1432 err << importLIDs_h[k];
1433 if (k + 1 < numImports) {
1434 err << ", ";
1435 }
1436 }
1437 err << "]. " << suffix << endl;
1438 return;
1439 }
1440
1441 RCP<const Comm<int> > comm = this->getMap ().is_null () ?
1442 Teuchos::null : this->getMap ()->getComm ();
1443 if (comm.is_null () || comm->getSize () == 1) {
1444 if (numImports != static_cast<size_t> (0)) {
1445 std::ostream& err = this->markLocalErrorAndGetStream ();
1446 err << prefix << "The input communicator is either null or "
1447 "has only one process, but numImports = " << numImports << " != 0."
1448 << suffix << endl;
1449 return;
1450 }
1451 // Don't go into the rest of this method unless there are
1452 // actually processes other than the calling process. This is
1453 // because the pack and unpack functions only have nonstub
1454 // implementations if building with MPI.
1455 return;
1456 }
1457
1458 // Make sure that the length of 'imports' fits in int.
1459 // This is ultimately an MPI restriction.
1460 if (static_cast<size_t> (imports.extent (0)) >
1461 static_cast<size_t> (std::numeric_limits<int>::max ())) {
1462 std::ostream& err = this->markLocalErrorAndGetStream ();
1463 err << prefix << "imports.extent(0) = "
1464 << imports.extent (0) << " does not fit in int." << endl;
1465 return;
1466 }
1467 const int inBufSize = static_cast<int> (imports.extent (0));
1468
1469 if (imports.need_sync_host ()) {
1470 imports.sync_host ();
1471 }
1472 if (numPacketsPerLID.need_sync_host ()) {
1473 numPacketsPerLID.sync_host ();
1474 }
1475 auto imports_h = imports.view_host ();
1476 auto numPacketsPerLID_h = numPacketsPerLID.view_host ();
1477
1478 // FIXME (mfh 17,24 Jan 2017) packTriples wants three arrays, not a
1479 // single array of structs. For now, we just copy.
1480 std::vector<GO> gblRowInds;
1481 std::vector<GO> gblColInds;
1482 std::vector<SC> vals;
1483
1484 const packet_type* inBuf = imports_h.data ();
1485 int inBufCurPos = 0;
1486 size_t numInvalidRowInds = 0;
1487 int errCode = 0;
1488 std::ostringstream errStrm; // for unpack* error output.
1489 for (size_t k = 0; k < numImports; ++k) {
1490 const LO lclRow = importLIDs_h(k);
1491 const GO gblRow = this->map_->getGlobalElement (lclRow);
1492 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1493 ++numInvalidRowInds;
1494 continue;
1495 }
1496
1497 // Remember where we were, so we don't overrun the buffer
1498 // length. inBufCurPos is an in/out argument of unpackTriples*.
1499 const int origInBufCurPos = inBufCurPos;
1500
1501 int numEnt = 0; // output argument of unpackTriplesCount
1502 errCode = unpackTriplesCount (inBuf, inBufSize, inBufCurPos,
1503 numEnt, *comm, &errStrm);
1504 if (errCode != 0 || numEnt < 0 || inBufCurPos < origInBufCurPos) {
1505 std::ostream& err = this->markLocalErrorAndGetStream ();
1506
1507 err << prefix << "In unpack loop, k=" << k << ": ";
1508 if (errCode != 0) {
1509 err << " unpackTriplesCount returned errCode = " << errCode
1510 << " != 0." << endl;
1511 }
1512 if (numEnt < 0) {
1513 err << " unpackTriplesCount returned errCode = 0, but numEnt = "
1514 << numEnt << " < 0." << endl;
1515 }
1516 if (inBufCurPos < origInBufCurPos) {
1517 err << " After unpackTriplesCount, inBufCurPos = " << inBufCurPos
1518 << " < origInBufCurPos = " << origInBufCurPos << "." << endl;
1519 }
1520 err << " unpackTriplesCount report: " << errStrm.str () << endl;
1521 err << suffix << endl;
1522
1523 // We only continue in a debug build, because the above error
1524 // messages could consume too much memory and cause an
1525 // out-of-memory error, without actually printing. Printing
1526 // everything is useful in a debug build, but not so much in a
1527 // release build.
1528#ifdef HAVE_TPETRA_DEBUG
1529 // Clear out the current error stream, so we don't accumulate
1530 // over loop iterations.
1531 errStrm.str ("");
1532 continue;
1533#else
1534 return;
1535#endif // HAVE_TPETRA_DEBUG
1536 }
1537
1538 // FIXME (mfh 17,24 Jan 2017) packTriples wants three arrays,
1539 // not a single array of structs. For now, we just copy.
1540 gblRowInds.resize (numEnt);
1541 gblColInds.resize (numEnt);
1542 vals.resize (numEnt);
1543
1544 errCode = unpackTriples (inBuf, inBufSize, inBufCurPos,
1545 gblRowInds.data (), gblColInds.data (),
1546 vals.data (), numEnt, *comm, &errStrm);
1547 if (errCode != 0) {
1548 std::ostream& err = this->markLocalErrorAndGetStream ();
1549 err << prefix << "unpackTriples returned errCode = "
1550 << errCode << " != 0. It reports: " << errStrm.str ()
1551 << endl;
1552 // We only continue in a debug build, because the above error
1553 // messages could consume too much memory and cause an
1554 // out-of-memory error, without actually printing. Printing
1555 // everything is useful in a debug build, but not so much in a
1556 // release build.
1557#ifdef HAVE_TPETRA_DEBUG
1558 // Clear out the current error stream, so we don't accumulate
1559 // over loop iterations.
1560 errStrm.str ("");
1561 continue;
1562#else
1563 return;
1564#endif // HAVE_TPETRA_DEBUG
1565 }
1566 this->sumIntoGlobalValues (gblRowInds.data (), gblColInds.data (),
1567 vals.data (), numEnt);
1568 }
1569
1570 // If we found invalid row indices in exportLIDs, go back,
1571 // collect, and report them.
1572 if (numInvalidRowInds != 0) {
1573 // Mark the error now, before we possibly run out of memory.
1574 // The latter could raise an exception (e.g., std::bad_alloc),
1575 // but at least we would get the error state right.
1576 std::ostream& err = this->markLocalErrorAndGetStream ();
1577
1578 std::vector<std::pair<LO, GO> > invalidRowInds;
1579 for (size_t k = 0; k < numImports; ++k) {
1580 const LO lclRow = importLIDs_h(k);
1581 const GO gblRow = this->map_->getGlobalElement (lclRow);
1582 if (gblRow == ::Tpetra::Details::OrdinalTraits<GO>::invalid ()) {
1583 invalidRowInds.push_back ({lclRow, gblRow});
1584 }
1585 }
1586
1587 err << prefix << "We found " << numInvalidRowInds << " invalid row ind"
1588 << (numInvalidRowInds != static_cast<size_t> (1) ? "ices" : "ex")
1589 << " out of " << numImports << " in importLIDs. Here is the list "
1590 << " of invalid row indices: [";
1591 for (size_t k = 0; k < invalidRowInds.size (); ++k) {
1592 err << "(LID: " << invalidRowInds[k].first << ", GID: "
1593 << invalidRowInds[k].second << ")";
1594 if (k + 1 < invalidRowInds.size ()) {
1595 err << ", ";
1596 }
1597 }
1598 err << "].";
1599 return;
1600 }
1601 }
1602};
1603
1604} // namespace Details
1605} // namespace Tpetra
1606
1607#endif // TPETRA_DETAILS_COOMATRIX_HPP
Declaration of a function that prints strings from each process.
Declaration and definition of Tpetra::Details::reallocDualViewIfNeeded, an implementation detail of T...
Struct that holds views of the contents of a CrsMatrix.
Sparse matrix used only for file input / output.
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
virtual std::string description() const
One-line descriptiion of this object; overrides Teuchos::Describable method.
virtual size_t constantNumberOfPackets() const
By returning 0, tell DistObject that this class may not necessarily have a constant number of "packet...
virtual ~CooMatrix()
Destructor (virtual for memory safety of derived classes).
void fillComplete()
Special version of fillComplete that assumes that the matrix already has a Map, and reuses its commun...
SC scalar_type
Type of each entry (value) in the sparse matrix.
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).
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 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,...
typename ::Tpetra::DistObject< char, LO, GO, NT >::buffer_device_type buffer_device_type
Kokkos::Device specialization for DistObject communication buffers.
CooMatrix(const ::Teuchos::RCP< const map_type > &map)
Constructor that takes a Map.
bool localError() const
Whether this object had an error on the calling process.
virtual bool checkSizes(const ::Tpetra::SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling 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...
std::string errorMessages() const
The current stream of error messages.
char packet_type
This class transfers data as bytes (MPI_BYTE).
::Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > map_type
Type of the Map specialization to give to the constructor.
Implementation detail of Tpetra::Details::CooMatrix (which see below).
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,...
std::size_t getLclNumEntries() const
Number of entries in the sparse matrix on the calling process.
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.
void sumIntoGlobalValues(const GO gblRowInds[], const GO gblColInds[], const SC vals[], const std::size_t numEnt)
Insert multiple entries locally into the sparse matrix.
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.
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)...
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...
char packet_type
Type for packing and unpacking data.
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.
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.
CooMatrixImpl()=default
Default constructor.
Base class for distributed Tpetra objects that support data redistribution.
Teuchos::RCP< const map_type > map_
The Map over which this object is distributed.
Node node_type
The Node type. If you don't know what this is, don't use it.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
Implementation details of Tpetra.
int packTriplesCount(const int, char[], const int, int &, const ::Teuchos::Comm< int > &, std::ostream *errStrm)
Pack the count (number) of matrix triples.
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 ("triples").
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 ("triples" (i, j, A(i,j))) into the given output buffer.
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.
int unpackTriples(const char[], const int, int &, OrdinalType[], OrdinalType[], ScalarType[], const int, const ::Teuchos::Comm< int > &, std::ostream *errStrm=NULL)
Unpack matrix entries ("triples" (i, j, A(i,j))) from the given input buffer.
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.
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,...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
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.
CombineMode
Rule for combining data in an Import or Export.
Function comparing two CooGraphEntry structs, lexicographically, first by row index,...
Type of each (row index, column index) pair in the Tpetra::Details::CooMatrix (see below).