41 #ifndef __Tpetra_MultiVectorFiller_hpp 42 #define __Tpetra_MultiVectorFiller_hpp 44 #include "Tpetra_MultiVector.hpp" 45 #include "Tpetra_Vector.hpp" 46 #include "Teuchos_CommHelpers.hpp" 74 sortAndMergeIn (Teuchos::Array<T>& allEntries,
75 Teuchos::ArrayView<T> currentEntries,
76 Teuchos::ArrayView<T> newEntries)
78 using Teuchos::ArrayView;
80 typedef typename ArrayView<T>::iterator iter_type;
81 typedef typename ArrayView<T>::size_type size_type;
83 std::sort (newEntries.begin(), newEntries.end());
84 iter_type it = std::unique (newEntries.begin(), newEntries.end());
85 const size_type numNew = as<size_type> (it - newEntries.begin());
87 ArrayView<T> newEntriesView = newEntries.view (0, numNew);
89 const size_type numCur = currentEntries.size();
90 if (allEntries.size() < numCur + numNew) {
91 allEntries.resize (numCur + numNew);
93 ArrayView<T> allView = allEntries.view (0, numCur + numNew);
94 ArrayView<T> newView = allEntries.view (numCur, numNew);
96 std::copy (newEntries.begin(), newEntries.end(), newView.begin());
97 std::inplace_merge (allView.begin(), newView.begin(), allView.end());
98 iter_type it2 = std::unique (allView.begin(), allView.end());
99 const size_type numTotal = as<size_type> (it2 - allView.begin());
101 return allEntries.view (0, numTotal);
112 typedef typename MV::scalar_type scalar_type;
113 typedef typename MV::local_ordinal_type local_ordinal_type;
114 typedef typename MV::global_ordinal_type global_ordinal_type;
115 typedef typename MV::node_type node_type;
117 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>
map_type;
152 const size_t numColumns) :
154 numCols_ (numColumns),
155 sourceIndices_ (numColumns),
156 sourceValues_ (numColumns)
163 const size_t oldNumColumns = getNumColumns();
164 if (newNumColumns >= oldNumColumns) {
165 for (
size_t j = oldNumColumns; j < newNumColumns; ++j) {
166 sourceIndices_.push_back (Teuchos::Array<global_ordinal_type> ());
167 sourceValues_.push_back (Teuchos::Array<scalar_type> ());
172 sourceIndices_.resize (newNumColumns);
173 sourceValues_.resize (newNumColumns);
175 numCols_ = oldNumColumns;
179 sumIntoGlobalValues (
const Teuchos::ArrayView<const global_ordinal_type>& rows,
181 const Teuchos::ArrayView<const scalar_type>& values)
183 if (column >= getNumColumns()) {
184 for (
size_t j = column; j < getNumColumns(); ++j) {
185 sourceIndices_.push_back (Teuchos::Array<global_ordinal_type> ());
186 sourceValues_.push_back (Teuchos::Array<scalar_type> ());
189 std::copy (rows.begin(), rows.end(), std::back_inserter (sourceIndices_[column]));
190 std::copy (values.begin(), values.end(), std::back_inserter (sourceValues_[column]));
202 const Teuchos::ArrayView<const scalar_type>& values)
204 typedef global_ordinal_type GO;
205 typedef scalar_type ST;
206 typedef typename Teuchos::ArrayView<const GO>::const_iterator GoIter;
207 typedef typename Teuchos::ArrayView<const ST>::const_iterator StIter;
209 const size_t numColumns = getNumColumns();
210 GoIter rowIter = rows.begin();
211 StIter valIter = values.begin();
212 for (
size_t j = 0; j < numColumns; ++j) {
213 GoIter rowIterNext = rowIter + numColumns;
214 StIter valIterNext = valIter + numColumns;
215 std::copy (rowIter, rowIterNext, std::back_inserter (sourceIndices_[j]));
216 std::copy (valIter, valIterNext, std::back_inserter (sourceValues_[j]));
217 rowIter = rowIterNext;
218 valIter = valIterNext;
246 template<
class BinaryFunction>
250 using Teuchos::Array;
251 using Teuchos::ArrayRCP;
252 using Teuchos::ArrayView;
254 typedef local_ordinal_type LO;
255 typedef global_ordinal_type GO;
256 typedef scalar_type ST;
257 typedef node_type NT;
259 RCP<const ::Tpetra::Map<LO, GO, NT> > map = X.getMap();
260 Array<LO> localIndices;
261 const size_t numColumns = getNumColumns();
262 for (
size_t j = 0; j < numColumns; ++j) {
263 const typename Array<const GO>::size_type numIndices = sourceIndices_[j].size();
268 if (localIndices.size() < numIndices) {
269 localIndices.resize (numIndices);
271 ArrayView<LO> localIndicesView = localIndices.view (0, numIndices);
272 ArrayView<const GO> globalIndicesView = sourceIndices_[j].view (0, numIndices);
273 for (
typename ArrayView<const GO>::size_type i = 0; i < numIndices; ++i) {
274 localIndices[i] = map->getLocalElement (globalIndicesView[i]);
277 ArrayRCP<ST> X_j = X.getDataNonConst (j);
278 ArrayView<const ST> localValues = sourceValues_[j].view (0, numIndices);
279 for (
typename ArrayView<const GO>::size_type i = 0; i < numIndices; ++i) {
280 const LO localInd = localIndices[i];
281 X_j[localInd] = f (X_j[localInd], localValues[i]);
291 locallyAssemble<std::plus<scalar_type> > (X, f);
296 Teuchos::Array<Teuchos::Array<global_ordinal_type> > newSourceIndices;
297 Teuchos::Array<Teuchos::Array<scalar_type> > newSourceValues;
299 std::swap (sourceIndices_, newSourceIndices);
300 std::swap (sourceValues_, newSourceValues);
304 Teuchos::Array<global_ordinal_type>
307 using Teuchos::Array;
308 using Teuchos::ArrayView;
310 typedef global_ordinal_type GO;
311 typedef typename Array<GO>::size_type size_type;
313 Array<GO> allInds (0);
314 const size_t numCols = getNumColumns();
321 const size_type numNew = sourceIndices_[0].size();
322 allInds.resize (allInds.size() + numNew);
323 std::copy (sourceIndices_[0].begin(), sourceIndices_[0].end(),
325 std::sort (allInds.begin(), allInds.end());
326 typename Array<GO>::iterator it =
327 std::unique (allInds.begin(), allInds.end());
328 const size_type numFinal = as<size_type> (it - allInds.begin());
329 allInds.resize (numFinal);
337 ArrayView<GO> curIndsView = allInds.view (0, 0);
339 for (
size_t j = 0; j < numCols; ++j) {
340 const size_type numNew = sourceIndices_[j].size();
341 if (numNew > newInds.size()) {
342 newInds.resize (numNew);
344 ArrayView<GO> newIndsView = newInds.view (0, numNew);
345 std::copy (sourceIndices_[j].begin(), sourceIndices_[j].end(),
346 newIndsView.begin());
347 curIndsView = sortAndMergeIn<GO> (allInds, curIndsView, newIndsView);
354 Teuchos::RCP<const map_type> map_;
356 Teuchos::Array<Teuchos::Array<global_ordinal_type> > sourceIndices_;
357 Teuchos::Array<Teuchos::Array<scalar_type> > sourceValues_;
359 size_t getNumColumns()
const {
return numCols_; }
371 typedef typename MV::scalar_type scalar_type;
372 typedef typename MV::local_ordinal_type local_ordinal_type;
373 typedef typename MV::global_ordinal_type global_ordinal_type;
374 typedef typename MV::node_type node_type;
376 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>
map_type;
389 const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT,
390 const Teuchos::RCP<Teuchos::FancyOStream>& out=Teuchos::null) :
393 verbLevel_ (verbLevel),
415 const size_t numColumns,
416 const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT,
417 const Teuchos::RCP<Teuchos::FancyOStream>& out=Teuchos::null) :
419 numCols_ (numColumns),
420 localVec_ (new MV (map, numColumns)),
422 nonlocalIndices_ (numColumns),
423 nonlocalValues_ (numColumns),
425 verbLevel_ (verbLevel),
432 std::ostringstream oss;
433 oss <<
"Tpetra::MultiVectorFillerData2<" 434 << Teuchos::TypeNameTraits<MV>::name () <<
">";
439 describe (Teuchos::FancyOStream& out,
440 const Teuchos::EVerbosityLevel verbLevel =
441 Teuchos::Describable::verbLevel_default)
const 444 using Teuchos::Array;
445 using Teuchos::ArrayRCP;
446 using Teuchos::ArrayView;
448 using Teuchos::VERB_DEFAULT;
449 using Teuchos::VERB_NONE;
450 using Teuchos::VERB_LOW;
451 using Teuchos::VERB_MEDIUM;
452 using Teuchos::VERB_HIGH;
453 using Teuchos::VERB_EXTREME;
456 const Teuchos::EVerbosityLevel vl =
457 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
459 RCP<const Teuchos::Comm<int> > comm = map_->getComm();
460 const int myImageID = comm->getRank();
461 const int numImages = comm->getSize();
463 if (vl != VERB_NONE) {
465 Teuchos::OSTab tab (out);
467 if (myImageID == 0) {
468 out << this->description() << endl;
470 for (
int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
471 if (myImageID == imageCtr) {
472 if (vl != VERB_LOW) {
474 out <<
"Process " << myImageID <<
":" << endl;
476 Teuchos::OSTab procTab (out);
478 out <<
"local length=" << localVec_->getLocalLength();
479 if (vl != VERB_MEDIUM) {
481 if (localVec_->isConstantStride()) {
482 out <<
", constant stride=" << localVec_->getStride() << endl;
485 out <<
", not constant stride" << endl;
487 if (vl == VERB_EXTREME) {
489 out <<
"Local values:" << endl;
490 ArrayRCP<ArrayRCP<const scalar_type> > X = localVec_->get2dView();
491 for (
size_t i = 0; i < localVec_->getLocalLength(); ++i) {
492 for (
size_t j = 0; j < localVec_->getNumVectors(); ++j) {
494 if (j < localVec_->getNumVectors() - 1) {
502 out <<
"Nonlocal indices and values:" << endl;
503 for (
size_t j = 0; j < (size_t)nonlocalIndices_.size(); ++j) {
504 ArrayView<const global_ordinal_type> inds = nonlocalIndices_[j]();
505 ArrayView<const scalar_type> vals = nonlocalValues_[j]();
506 for (
typename ArrayView<const global_ordinal_type>::size_type k = 0;
507 k < inds.size(); ++k) {
508 out <<
"X(" << inds[k] <<
"," << j <<
") = " << vals[k] << endl;
524 Teuchos::Array<global_ordinal_type>
526 const Teuchos::RCP<Teuchos::FancyOStream>& out = Teuchos::null,
527 const bool debug =
false)
const 529 using Teuchos::Array;
530 using Teuchos::ArrayView;
534 typedef global_ordinal_type GO;
535 const char prefix[] =
"Tpetra::MultiVectorFiller::getSourceIndices: ";
537 if (debug && ! out.is_null ()) {
538 std::ostringstream os;
539 os <<
"Proc " << comm.getRank () <<
": getSourceIndices" << std::endl;
545 Array<GO> nonlocalIndices = getSortedUniqueNonlocalIndices (comm, out, debug);
548 ArrayView<const GO> localIndices = getLocalIndices ();
554 Array<GO> indices (localIndices.size () + nonlocalIndices.size ());
555 ArrayView<GO> localIndView = indices.view (0, localIndices.size ());
556 TEUCHOS_TEST_FOR_EXCEPTION
557 (localIndView.size () > indices.size (), std::logic_error,
558 prefix <<
"localIndView.size() = " << localIndView.size ()
559 <<
" > indices.size() = " << indices.size () <<
".");
560 TEUCHOS_TEST_FOR_EXCEPTION
561 (localIndView.size () != localIndices.size (), std::logic_error,
562 prefix <<
"localIndView.size() = " << localIndView.size ()
563 <<
" != localIndices.size() = " << localIndices.size () <<
".");
565 std::copy (localIndices.begin (), localIndices.end (), localIndView.begin ());
566 std::sort (localIndView.begin (), localIndView.end ());
568 if (debug && ! out.is_null ()) {
569 std::ostringstream os;
570 os <<
"Proc " << comm.getRank () <<
": Right after copying and sorting" << std::endl;
575 if (nonlocalIndices.size () > 0) {
576 ArrayView<GO> nonlclIndView =
577 indices.view (localIndices.size (), nonlocalIndices.size ());
585 ArrayView<GO> indView = indices.view (0, indices.size ());
586 if (debug && ! out.is_null ()) {
587 std::ostringstream os;
588 os <<
"Right before std::copy" << std::endl;
591 std::copy (nonlocalIndices.begin (),
592 nonlocalIndices.end (),
593 nonlclIndView.begin ());
594 if (debug && ! out.is_null ()) {
595 std::ostringstream os;
596 os <<
"Proc " << comm.getRank () <<
": Right before std::inplace_merge" 600 std::inplace_merge (indView.begin (),
601 indView.begin () + localIndices.size (),
605 if (debug && ! out.is_null ()) {
606 std::ostringstream os;
607 os <<
"Proc " << comm.getRank () <<
": Done with getSourceIndices" 626 using Teuchos::Array;
627 using Teuchos::Range1D;
630 const size_t oldNumColumns = numCols_;
631 if (newNumColumns == oldNumColumns) {
636 if (newNumColumns > oldNumColumns) {
637 newLclVec = rcp (
new MV (map_, newNumColumns));
641 RCP<MV> newLclVecView =
642 newLclVec->subViewNonConst (Range1D (0, oldNumColumns-1));
643 *newLclVecView = *localVec_;
646 if (newNumColumns == 0) {
649 newLclVec = Teuchos::null;
652 newLclVec = localVec_->subViewNonConst (Range1D (0, newNumColumns-1));
657 localVec_ = newLclVec;
658 numCols_ = newNumColumns;
668 const size_t columnIndex,
669 const Teuchos::ArrayView<const scalar_type>& values)
671 using Teuchos::ArrayView;
672 typedef local_ordinal_type LO;
673 typedef global_ordinal_type GO;
674 typedef scalar_type ST;
675 typedef decltype (rows.size ()) size_type;
676 const char prefix[] =
"Tpetra::MultiVectorFiller::sumIntoGlobalValues: ";
678 if (map_.is_null ()) {
681 const size_type numEnt = rows.size ();
682 TEUCHOS_TEST_FOR_EXCEPTION
683 (numEnt != values.size (), std::invalid_argument, prefix
684 <<
"rows.size() = " << numEnt <<
" != values.size() = " 685 << values.size () <<
".");
688 if (columnIndex >= getNumColumns()) {
696 auto X_j = localVec_->getVector (columnIndex);
697 auto X_j_2d = X_j->template getLocalView<typename MV::dual_view_type::t_host::memory_space> ();
698 auto X_j_1d = Kokkos::subview (X_j_2d, Kokkos::ALL (), 0);
700 const map_type& theMap = *map_;
701 for (size_type k = 0; k < numEnt; ++k) {
702 const ST val = values[k];
703 const GO gblRowInd = rows[k];
705 if (lclRowInd == Teuchos::OrdinalTraits<LO>::invalid ()) {
709 auto& innerMap = nonlocals_[columnIndex];
710 auto innerIter = innerMap.find (gblRowInd);
711 if (innerIter == innerMap.end ()) {
712 innerMap.insert (std::make_pair (gblRowInd, values[k]));
714 innerIter->second += val;
721 X_j_1d[lclRowInd] += val;
737 const Teuchos::ArrayView<const scalar_type>& values)
739 using Teuchos::ArrayView;
740 typedef typename ArrayView<const global_ordinal_type>::size_type size_type;
742 const size_t numCols = getNumColumns();
743 for (
size_t j = 0; j < numCols; ++j) {
744 const size_type offset = numCols*j;
745 const size_type len = numCols;
746 this->sumIntoGlobalValues (rows.view (offset, len), j, values.view (offset, len));
774 template<
class BinaryFunction>
778 typedef scalar_type ST;
779 typedef local_ordinal_type LO;
780 typedef global_ordinal_type GO;
781 typedef typename MV::dual_view_type::t_host::memory_space host_memory_space;
783 if (X.getMap ().is_null ()) {
786 const map_type& srcMap = * (X.getMap ());
790 for (
size_t j = 0; j < X.getNumVectors (); ++j) {
791 auto X_j = X.getVector (j);
792 auto X_j_2d = X_j->template getLocalView<host_memory_space> ();
793 auto X_j_1d = Kokkos::subview (X_j_2d, Kokkos::ALL (), 0);
796 if (! localVec_.is_null () && ! localVec_->getMap ().is_null ()
797 && j < localVec_->getNumVectors ()) {
798 auto lcl_j = localVec_->getVector (j);
799 auto lcl_j_2d = lcl_j->template getLocalView<host_memory_space> ();
800 auto lcl_j_1d = Kokkos::subview (lcl_j_2d, Kokkos::ALL (), 0);
805 const map_type& lclMap = * (localVec_->getMap ());
806 const LO lclNumRows =
static_cast<LO
> (lcl_j->getLocalLength ());
807 for (LO i_lcl = 0; i_lcl < lclNumRows; ++i_lcl) {
810 X_j_1d(i_X) = f (X_j_1d(i_X), lcl_j_1d(i_lcl));
815 auto outerIter = nonlocals_.find (j);
816 if (outerIter != nonlocals_.end ()) {
817 auto beg = outerIter->second.begin ();
818 auto end = outerIter->second.end ();
819 for (
auto innerIter = beg; innerIter != end; ++innerIter) {
820 const GO gblRowInd = innerIter->first;
823 if (lclRowInd != Teuchos::OrdinalTraits<LO>::invalid ()) {
824 const ST val = innerIter->second;
825 X_j_1d(lclRowInd) = f (X_j_1d(lclRowInd), val);
841 locallyAssemble<std::plus<scalar_type> > (X, f);
851 std::map<size_t, std::map<global_ordinal_type, scalar_type> > newNonlcls;
855 std::swap (nonlocals_, newNonlcls);
860 if (! localVec_.is_null ()) {
861 localVec_->putScalar (Teuchos::ScalarTraits<scalar_type>::zero());
867 Teuchos::RCP<const map_type> map_;
873 Teuchos::RCP<MV> localVec_;
876 std::map<size_t, std::map<global_ordinal_type, scalar_type> > nonlocals_;
879 Teuchos::EVerbosityLevel verbLevel_;
882 Teuchos::RCP<Teuchos::FancyOStream> out_;
885 size_t getNumColumns()
const {
return numCols_; }
888 Teuchos::ArrayView<const global_ordinal_type>
889 getLocalIndices()
const 891 return map_->getNodeElementList ();
895 Teuchos::Array<global_ordinal_type>
896 getSortedUniqueNonlocalIndices (
const Teuchos::Comm<int>& comm,
897 const Teuchos::RCP<Teuchos::FancyOStream>& out = Teuchos::null,
898 const bool debug =
false)
const 900 using Teuchos::Array;
901 using Teuchos::ArrayView;
903 typedef global_ordinal_type GO;
905 if (debug && ! out.is_null ()) {
906 std::ostringstream os;
907 os <<
"Proc " << comm.getRank () <<
": getSortedUniqueNonlocalIndices" 915 std::set<GO> indsOut;
916 const size_t numCols = getNumColumns ();
917 for (
size_t j = 0; j < numCols; ++j) {
918 auto outerIter = nonlocals_.find (j);
919 if (outerIter != nonlocals_.end ()) {
920 auto beg = outerIter->second.begin ();
921 auto end = outerIter->second.end ();
922 for (
auto innerIter = beg; innerIter != end; ++innerIter) {
924 indsOut.insert (innerIter->first);
929 if (debug && ! out.is_null ()) {
930 std::ostringstream os;
931 os <<
"Proc " << comm.getRank () <<
": Made nonlocals set" << std::endl;
935 Array<GO> allNonlocals (indsOut.begin (), indsOut.end ());
936 if (debug && ! out.is_null ()) {
937 std::ostringstream os;
938 os <<
"Proc " << comm.getRank () <<
": Done with " 939 "getSortedUniqueNonlocalIndices" << std::endl;
977 typedef typename MV::scalar_type scalar_type;
978 typedef typename MV::local_ordinal_type local_ordinal_type;
979 typedef typename MV::global_ordinal_type global_ordinal_type;
980 typedef typename MV::node_type node_type;
981 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>
map_type;
1013 const size_t numCols);
1035 globalAssemble (MV& X_out,
const bool forceReuseMap =
false);
1043 describe (Teuchos::FancyOStream& out,
1044 const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default)
const 1046 data_.describe (out, verbLevel);
1062 Teuchos::ArrayView<const scalar_type> values)
1064 data_.sumIntoGlobalValues (rows, column, values);
1084 Teuchos::ArrayView<const scalar_type> values)
1086 data_.sumIntoGlobalValues (rows, values);
1091 Teuchos::RCP<const map_type> ctorMap_;
1097 Teuchos::RCP<const map_type> sourceMap_;
1104 Teuchos::RCP<const map_type> targetMap_;
1118 Teuchos::RCP<MV> sourceVec_;
1127 Teuchos::RCP<export_type> exporter_;
1139 Teuchos::Array<global_ordinal_type>
1141 const Teuchos::RCP<Teuchos::FancyOStream>& out = Teuchos::null,
1142 const bool debug =
false)
const 1155 Teuchos::RCP<const map_type>
1156 computeSourceMap (
const global_ordinal_type indexBase,
1157 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1158 const Teuchos::RCP<Teuchos::FancyOStream>& out = Teuchos::null,
1159 const bool debug =
false);
1167 data_ (map, numCols),
1172 Teuchos::RCP<const typename MultiVectorFiller<MV>::map_type>
1175 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1176 const Teuchos::RCP<Teuchos::FancyOStream>& out,
1179 using Teuchos::Array;
1180 using Teuchos::ArrayView;
1182 typedef global_ordinal_type GO;
1184 if (debug && ! out.is_null ()) {
1186 std::ostringstream os;
1187 const int myRank = comm->getRank ();
1188 os <<
"Proc " << myRank <<
": computeSourceMap" << std::endl;
1192 Array<GO> indices = getSourceIndices (*comm, out, debug);
1194 if (debug && ! out.is_null ()) {
1195 std::ostringstream os;
1196 const int myRank = comm->getRank ();
1197 os <<
"Proc " << myRank <<
": computeSourceMap: About to create Map" 1207 const GST INV = Teuchos::OrdinalTraits<GST>::invalid ();
1208 return rcp (
new map_type (INV, indices, indexBase, comm));
1216 using Teuchos::ArrayView;
1217 using Teuchos::Array;
1218 using Teuchos::Range1D;
1221 RCP<Teuchos::FancyOStream> outPtr;
1222 const bool debug =
false;
1225 outPtr = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
1229 const size_t numVecs = X_out.getNumVectors ();
1243 RCP<const map_type> targetMap;
1244 bool assumeSameTargetMap =
false;
1245 if (targetMap_.is_null()) {
1246 targetMap_ = X_out.getMap();
1247 targetMap = targetMap_;
1248 assumeSameTargetMap =
false;
1251 if (forceReuseMap) {
1252 targetMap = targetMap_;
1253 assumeSameTargetMap =
true;
1261 if (targetMap_->isSameAs (*(X_out.getMap()))) {
1262 assumeSameTargetMap =
true;
1263 targetMap = targetMap_;
1268 if (debug && ! outPtr.is_null ()) {
1269 std::ostringstream os;
1270 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1271 const int myRank = comm.getRank ();
1272 os <<
"Proc " << myRank <<
": Right before getting source Map" << std::endl;
1273 *outPtr << os.str ();
1281 RCP<const map_type> sourceMap;
1282 bool computedSourceMap =
false;
1284 if (forceReuseMap && ! sourceMap_.is_null()) {
1285 sourceMap = sourceMap_;
1288 sourceMap = computeSourceMap (ctorMap_->getIndexBase (),
1289 ctorMap_->getComm (),
1291 computedSourceMap =
true;
1294 if (computedSourceMap) {
1295 sourceMap_ = sourceMap;
1298 if (debug && ! outPtr.is_null ()) {
1299 std::ostringstream os;
1300 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1301 const int myRank = comm.getRank ();
1302 os <<
"Proc " << myRank <<
": Right before computing Export" << std::endl;
1303 *outPtr << os.str ();
1310 const bool mustComputeExport =
1311 (exporter_.is_null() || (assumeSameTargetMap && ! computedSourceMap));
1312 if (mustComputeExport) {
1313 exporter_ = rcp (
new export_type (sourceMap_, targetMap_));
1316 if (debug && ! outPtr.is_null ()) {
1317 std::ostringstream os;
1318 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1319 const int myRank = comm.getRank ();
1320 os <<
"Proc " << myRank <<
": Right after computing Export" << std::endl;
1321 *outPtr << os.str ();
1326 const bool mustReallocateVec = sourceVec_.is_null() ||
1327 sourceVec_->getNumVectors() < numVecs || ! assumeSameTargetMap;
1329 if (mustReallocateVec) {
1331 sourceVec_ = rcp (
new MV (sourceMap_, numVecs));
1334 if (sourceVec_->getNumVectors() == numVecs) {
1337 X_in = sourceVec_->subViewNonConst (Range1D (0, numVecs-1));
1341 if (debug && ! outPtr.is_null ()) {
1342 std::ostringstream os;
1343 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1344 const int myRank = comm.getRank ();
1345 os <<
"Proc " << myRank <<
": Right before locallyAssemble" << std::endl;
1346 *outPtr << os.str ();
1351 locallyAssemble (*X_in);
1353 if (debug && ! outPtr.is_null ()) {
1354 std::ostringstream os;
1355 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1356 const int myRank = comm.getRank ();
1357 os <<
"Proc " << myRank <<
": Right after locallyAssemble" << std::endl;
1358 *outPtr << os.str ();
1363 X_out.doExport (*X_in, *exporter_, combineMode);
1365 if (debug && ! outPtr.is_null ()) {
1366 std::ostringstream os;
1367 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1368 const int myRank = comm.getRank ();
1369 os <<
"Proc " << myRank <<
": Right after Export" << std::endl;
1370 *outPtr << os.str ();
1374 X_in->putScalar (Teuchos::ScalarTraits<scalar_type>::zero ());
1379 if (debug && ! outPtr.is_null ()) {
1380 std::ostringstream os;
1381 const Teuchos::Comm<int>& comm = * (X_out.getMap ()->getComm ());
1382 const int myRank = comm.getRank ();
1383 os <<
"Proc " << myRank <<
": Done with globalAssemble" << std::endl;
1384 *outPtr << os.str ();
1399 typedef typename MV::scalar_type scalar_type;
1400 typedef typename MV::local_ordinal_type local_ordinal_type;
1401 typedef typename MV::global_ordinal_type global_ordinal_type;
1402 typedef typename MV::node_type node_type;
1403 typedef ::Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type>
map_type;
1415 const global_ordinal_type eltSize,
1416 const size_t numCols,
1417 const Teuchos::RCP<Teuchos::FancyOStream>& outStream=Teuchos::null,
1418 const Teuchos::EVerbosityLevel verbosityLevel=Teuchos::VERB_DEFAULT)
1420 using Teuchos::Array;
1421 using Teuchos::ArrayRCP;
1422 using Teuchos::ArrayView;
1424 using Teuchos::Comm;
1425 using Teuchos::FancyOStream;
1426 using Teuchos::getFancyOStream;
1427 using Teuchos::oblackholestream;
1431 using Teuchos::rcpFromRef;
1432 using Teuchos::REDUCE_SUM;
1433 using Teuchos::reduceAll;
1437 typedef local_ordinal_type LO;
1438 typedef global_ordinal_type GO;
1439 typedef scalar_type ST;
1440 typedef Teuchos::ScalarTraits<ST> STS;
1442 TEUCHOS_TEST_FOR_EXCEPTION(eltSize % 2 == 0, std::invalid_argument,
1443 "Element size (eltSize) argument must be odd.");
1444 TEUCHOS_TEST_FOR_EXCEPTION(numCols == 0, std::invalid_argument,
1445 "Number of columns (numCols) argument must be nonzero.");
1447 RCP<FancyOStream> out = outStream.is_null() ?
1448 getFancyOStream (rcp (
new oblackholestream)) : outStream;
1449 const Teuchos::EVerbosityLevel verbLevel =
1450 (verbosityLevel == Teuchos::VERB_DEFAULT) ?
1451 Teuchos::VERB_NONE : verbosityLevel;
1455 const GO indexBase = targetMap->getIndexBase();
1456 Array<GO> rows (eltSize);
1457 Array<ST> values (eltSize);
1458 std::fill (values.begin(), values.end(), STS::one());
1462 RCP<MultiVectorFiller<MV> > filler =
1465 TEUCHOS_TEST_FOR_EXCEPTION(! targetMap->isContiguous(),
1466 std::invalid_argument,
"MultiVectorFiller test currently only works " 1467 "for contiguous Maps.");
1469 const GO minGlobalIndex = targetMap->getMinGlobalIndex();
1470 const GO maxGlobalIndex = targetMap->getMaxGlobalIndex();
1471 const GO minAllGlobalIndex = targetMap->getMinAllGlobalIndex();
1472 const GO maxAllGlobalIndex = targetMap->getMaxAllGlobalIndex();
1473 for (
size_t j = 0; j < numCols; ++j) {
1474 for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
1476 const GO start = std::max (i - eltSize/2, minAllGlobalIndex);
1477 const GO end = std::min (i + eltSize/2, maxAllGlobalIndex);
1478 const GO len = end - start + 1;
1480 TEUCHOS_TEST_FOR_EXCEPTION(len > eltSize, std::logic_error,
1481 "At start,end = " << start <<
"," << end <<
", len = " << len
1482 <<
" > eltSize = " << eltSize <<
".");
1484 for (GO k = 0; k < len; ++k) {
1485 rows[k] = start + k;
1487 if (verbLevel == Teuchos::VERB_EXTREME) {
1488 *out <<
"Inserting: " 1489 << Teuchos::toString (rows.view(0,len)) <<
", " 1490 << Teuchos::toString (values.view(0, len)) << std::endl;
1492 filler->sumIntoGlobalValues (rows.view(0, len), j, values.view(0, len));
1496 if (verbLevel == Teuchos::VERB_EXTREME) {
1497 *out <<
"Filler:" << std::endl;
1498 filler->describe (*out, verbLevel);
1502 MV X_out (targetMap, numCols);
1503 filler->globalAssemble (X_out);
1504 filler = Teuchos::null;
1506 if (verbLevel == Teuchos::VERB_EXTREME) {
1507 *out <<
"X_out:" << std::endl;
1508 X_out.describe (*out, verbLevel);
1513 MV X_expected (targetMap, numCols);
1514 const scalar_type two = STS::one() + STS::one();
1515 for (
size_t j = 0; j < numCols; ++j) {
1516 ArrayRCP<ST> X_j = X_expected.getDataNonConst (j);
1521 for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
1522 const LO localIndex = targetMap->getLocalElement (i);
1523 TEUCHOS_TEST_FOR_EXCEPTION(i == Teuchos::OrdinalTraits<LO>::invalid(),
1524 std::logic_error,
"Global index " << i <<
" is not in the " 1525 "multivector's Map.");
1527 if (i <= minAllGlobalIndex + eltSize/2) {
1528 X_j[localIndex] = two + as<ST>(i) - as<ST>(indexBase);
1530 else if (i >= maxAllGlobalIndex - eltSize/2) {
1531 X_j[localIndex] = two + as<ST>(maxAllGlobalIndex) - as<ST>(i);
1534 X_j[localIndex] = as<ST>(eltSize);
1539 if (verbLevel == Teuchos::VERB_EXTREME) {
1540 *out <<
"X_expected:" << std::endl;
1541 X_expected.describe (*out, verbLevel);
1545 Array<GO> errorLocations;
1546 for (
size_t j = 0; j < numCols; ++j) {
1547 ArrayRCP<const ST> X_out_j = X_out.getData (j);
1548 ArrayRCP<const ST> X_expected_j = X_expected.getData (j);
1553 for (GO i = minGlobalIndex; i <= maxGlobalIndex; ++i) {
1554 const LO localIndex = targetMap->getLocalElement (i);
1555 TEUCHOS_TEST_FOR_EXCEPTION(i == Teuchos::OrdinalTraits<LO>::invalid(),
1556 std::logic_error,
"Global index " << i <<
" is not in the " 1557 "multivector's Map.");
1561 if (X_out_j[localIndex] != X_expected_j[localIndex]) {
1562 errorLocations.push_back (i);
1566 const typename Array<GO>::size_type localNumErrors = errorLocations.size();
1567 typename Array<GO>::size_type globalNumErrors = 0;
1568 RCP<const Comm<int> > comm = targetMap->getComm();
1569 reduceAll (*comm, REDUCE_SUM, localNumErrors, ptr (&globalNumErrors));
1571 if (globalNumErrors != 0) {
1572 std::ostringstream os;
1573 os <<
"Proc " << comm->getRank() <<
": " << localNumErrors
1574 <<
" incorrect value" << (localNumErrors != 1 ?
"s" :
"")
1575 <<
". Error locations: [ ";
1576 std::copy (errorLocations.begin(), errorLocations.end(),
1577 std::ostream_iterator<GO> (os,
" "));
1581 for (
int p = 0; p < comm->getSize(); ++p) {
1582 if (p == comm->getRank()) {
1583 cerr << os.str() << endl;
1590 TEUCHOS_TEST_FOR_EXCEPTION(globalNumErrors != 0, std::logic_error,
1591 "Over all procs: " << globalNumErrors <<
" total error" 1592 << (globalNumErrors != 1 ?
"s" :
"") <<
".");
1599 template<
class ScalarType,
1600 class LocalOrdinalType,
1601 class GlobalOrdinalType,
1604 testMultiVectorFiller (
const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1605 const size_t unknownsPerNode,
1606 const GlobalOrdinalType unknownsPerElt,
1607 const size_t numCols,
1608 const Teuchos::RCP<Teuchos::FancyOStream>& outStream,
1609 const Teuchos::EVerbosityLevel verbLevel)
1612 using Teuchos::FancyOStream;
1613 using Teuchos::getFancyOStream;
1614 using Teuchos::oblackholestream;
1615 using Teuchos::ParameterList;
1621 typedef ScalarType ST;
1622 typedef LocalOrdinalType LO;
1623 typedef GlobalOrdinalType GO;
1624 typedef NodeType NT;
1625 typedef ::Tpetra::Map<LO, GO, NT>
map_type;
1629 RCP<FancyOStream> out = outStream.is_null() ?
1630 getFancyOStream (rcp (
new oblackholestream)) : outStream;
1631 const GST INV = Teuchos::OrdinalTraits<GST>::invalid ();
1632 const GO indexBase = 0;
1633 RCP<const map_type> targetMap =
1634 rcp (
new map_type (INV, unknownsPerNode, indexBase, comm));
1636 std::ostringstream os;
1639 }
catch (std::exception& e) {
1643 for (
int p = 0; p < comm->getSize(); ++p) {
1644 if (p == comm->getRank()) {
1645 *out <<
"On process " << comm->getRank() <<
": " << os.str() << endl;
1659 testSortAndMergeIn ();
1665 #endif // __Tpetra_MultiVectorFiller_hpp void locallyAssemble(MV &X)
Do the "local" (not MPI) part of assembly.
Namespace Tpetra contains the class and methods constituting the Tpetra library.
static void testSameMap(const Teuchos::RCP< const map_type > &targetMap, const global_ordinal_type eltSize, const size_t numCols, const Teuchos::RCP< Teuchos::FancyOStream > &outStream=Teuchos::null, const Teuchos::EVerbosityLevel verbosityLevel=Teuchos::VERB_DEFAULT)
Test global assembly when constructor Map = target Map.
Second implementation of fill and local assembly for MultiVectorFiller.
void clear()
Clear the contents of the vector, making it implicitly a vector of zeros.
MultiVectorFiller(const Teuchos::RCP< const map_type > &map, const size_t numCols)
Constructor.
void locallyAssemble(MV &X, BinaryFunction &f)
Locally assemble into X, with user-specified combine mode.
Teuchos::Array< global_ordinal_type > getSourceIndices() const
All source indices (local and nonlocal) of the source Map, sorted and unique.
One or more distributed dense vectors.
void locallyAssemble(MV &X)
Overload of locallyAssemble() for the usual ADD combine mode.
Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > createContigMapWithNode(global_size_t numElements, size_t localNumElements, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, const Teuchos::RCP< Node > &node=defaultArgNode< Node >())
Non-member constructor for a (potentially) non-uniformly distributed, contiguous Map with a user-spec...
Tests for MultiVectorFiller.
void sumIntoGlobalValues(const Teuchos::ArrayView< const global_ordinal_type > &rows, const Teuchos::ArrayView< const scalar_type > &values)
void sumIntoGlobalValues(Teuchos::ArrayView< const global_ordinal_type > rows, Teuchos::ArrayView< const scalar_type > values)
Sum data into the multivector.
Adds nonlocal sum-into functionality to Tpetra::MultiVector.
Implementation details of Tpetra.
MultiVectorFillerData(const Teuchos::RCP< const map_type > &map, const size_t numColumns)
Constructor.
void sumIntoGlobalValues(const Teuchos::ArrayView< const global_ordinal_type > &rows, const size_t columnIndex, const Teuchos::ArrayView< const scalar_type > &values)
Set entry (rows[i],columnIndex) to values[i], for all i.
void locallyAssemble(MV &X, BinaryFunction &f)
Locally assemble into X.
size_t global_size_t
Global size_t object.
void sumIntoGlobalValues(const Teuchos::ArrayView< const global_ordinal_type > &rows, const Teuchos::ArrayView< const scalar_type > &values)
Set entry (rows[i],j) to values[i], for all i and j.
Implementation of fill and local assembly for MultiVectorFiller.
MultiVectorFillerData2(const Teuchos::RCP< const map_type > &map, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT, const Teuchos::RCP< Teuchos::FancyOStream > &out=Teuchos::null)
Default constructor (sets number of columns to zero).
Teuchos::Array< global_ordinal_type > getSourceIndices(const Teuchos::Comm< int > &comm, const Teuchos::RCP< Teuchos::FancyOStream > &out=Teuchos::null, const bool debug=false) const
All source indices (local and nonlocal) of the source Map, sorted and unique.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
void globalAssemble(MV &X_out, const bool forceReuseMap=false)
Assemble all the data (local and nonlocal) into X_out.
GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const
The global index corresponding to the given local index.
CombineMode
Rule for combining data in an Import or Export.
Sum new values into existing values.
void clear()
Clear contents, in preparation for another fill cycle.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
Describes a parallel distribution of objects over processes.
MultiVectorFillerData2(const Teuchos::RCP< const map_type > &map, const size_t numColumns, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_DEFAULT, const Teuchos::RCP< Teuchos::FancyOStream > &out=Teuchos::null)
Constructor.
MultiVectorFillerData(const Teuchos::RCP< const map_type > &map)
Default constructor (sets number of columns to zero).
void sumIntoGlobalValues(Teuchos::ArrayView< const global_ordinal_type > rows, size_t column, Teuchos::ArrayView< const scalar_type > values)
Sum data into the multivector.
void setNumColumns(const size_t newNumColumns)
Set the number of columns in the output multivector.
void setNumColumns(const size_t newNumColumns)
Set the number of columns in the output multivector.
void clear()
Clear the contents of the multivector.