3969 RCP<const map_type> nonlocalRowMap;
3970 Teuchos::Array<size_t> numEntPerNonlocalRow(myNumNonlocalRows);
3972 Teuchos::Array<GO> myNonlocalGblRows(myNumNonlocalRows);
3973 size_type curPos = 0;
3975 ++mapIter, ++curPos) {
3976 myNonlocalGblRows[curPos] = mapIter->first;
3979 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
3980 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
3987 sort2(gblCols.begin(), gblCols.end(), vals.begin());
3988 typename Teuchos::Array<GO>::iterator gblCols_newEnd;
3989 typename Teuchos::Array<Scalar>::iterator vals_newEnd;
3990 merge2(gblCols_newEnd, vals_newEnd,
3991 gblCols.begin(), gblCols.end(),
3992 vals.begin(), vals.end());
3993 gblCols.erase(gblCols_newEnd, gblCols.end());
3994 vals.erase(vals_newEnd, vals.end());
3995 numEntPerNonlocalRow[curPos] = gblCols.size();
4006 GO myMinNonlocalGblRow = std::numeric_limits<GO>::max();
4008 auto iter = std::min_element(myNonlocalGblRows.begin(),
4009 myNonlocalGblRows.end());
4010 if (iter != myNonlocalGblRows.end()) {
4011 myMinNonlocalGblRow = *iter;
4014 GO gblMinNonlocalGblRow = 0;
4015 reduceAll<int, GO>(*comm, REDUCE_MIN, myMinNonlocalGblRow,
4016 outArg(gblMinNonlocalGblRow));
4017 const GO indexBase = gblMinNonlocalGblRow;
4018 const global_size_t INV = Teuchos::OrdinalTraits<global_size_t>::invalid();
4019 nonlocalRowMap = rcp(
new map_type(INV, myNonlocalGblRows(), indexBase, comm));
4028 std::ostringstream os;
4029 os << *prefix <<
"Create nonlocal matrix" << endl;
4030 std::cerr << os.str();
4032 RCP<crs_matrix_type> nonlocalMatrix =
4033 rcp(
new crs_matrix_type(nonlocalRowMap, numEntPerNonlocalRow()));
4035 size_type curPos = 0;
4036 for (
auto mapIter = nonlocals_.begin(); mapIter != nonlocals_.end();
4037 ++mapIter, ++curPos) {
4038 const GO gblRow = mapIter->first;
4040 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4041 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4043 nonlocalMatrix->insertGlobalValues(gblRow, gblCols(), vals());
4055 auto origRowMap = this->getRowMap();
4056 const bool origRowMapIsOneToOne = origRowMap->isOneToOne();
4058 int isLocallyComplete = 1;
4060 if (origRowMapIsOneToOne) {
4062 std::ostringstream os;
4063 os << *prefix <<
"Original row Map is 1-to-1" << endl;
4064 std::cerr << os.str();
4066 export_type exportToOrig(nonlocalRowMap, origRowMap);
4067 if (!exportToOrig.isLocallyComplete()) {
4068 isLocallyComplete = 0;
4071 std::ostringstream os;
4072 os << *prefix <<
"doExport from nonlocalMatrix" << endl;
4073 std::cerr << os.str();
4075 this->doExport(*nonlocalMatrix, exportToOrig,
Tpetra::ADD);
4079 std::ostringstream os;
4080 os << *prefix <<
"Original row Map is NOT 1-to-1" << endl;
4081 std::cerr << os.str();
4088 export_type exportToOneToOne(nonlocalRowMap, oneToOneRowMap);
4089 if (!exportToOneToOne.isLocallyComplete()) {
4090 isLocallyComplete = 0;
4098 std::ostringstream os;
4099 os << *prefix <<
"Create & doExport into 1-to-1 matrix"
4101 std::cerr << os.str();
4103 crs_matrix_type oneToOneMatrix(oneToOneRowMap, 0);
4105 oneToOneMatrix.doExport(*nonlocalMatrix, exportToOneToOne,
4111 std::ostringstream os;
4112 os << *prefix <<
"Free nonlocalMatrix" << endl;
4113 std::cerr << os.str();
4115 nonlocalMatrix = Teuchos::null;
4119 std::ostringstream os;
4120 os << *prefix <<
"doImport from 1-to-1 matrix" << endl;
4121 std::cerr << os.str();
4123 import_type importToOrig(oneToOneRowMap, origRowMap);
4124 this->doImport(oneToOneMatrix, importToOrig,
Tpetra::ADD);
4132 std::ostringstream os;
4133 os << *prefix <<
"Free nonlocals_ (std::map)" << endl;
4134 std::cerr << os.str();
4136 decltype(nonlocals_) newNonlocals;
4137 std::swap(nonlocals_, newNonlocals);
4146 int isGloballyComplete = 0;
4147 reduceAll<int, int>(*comm, REDUCE_MIN, isLocallyComplete,
4148 outArg(isGloballyComplete));
4149 TEUCHOS_TEST_FOR_EXCEPTION(isGloballyComplete != 1, std::runtime_error,
4150 "On at least one process, "
4151 "you called insertGlobalValues with a global row index which is not in "
4152 "the matrix's row Map on any process in its communicator.");
4155template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4157 resumeFill(
const Teuchos::RCP<Teuchos::ParameterList>& params) {
4159 myGraph_->resumeFill(params);
4162 applyHelper.reset();
4166template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4169 return getCrsGraphRef().haveGlobalConstants();
4172template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4174 fillComplete(
const Teuchos::RCP<Teuchos::ParameterList>& params) {
4175 const char tfecfFuncName[] =
"fillComplete(params): ";
4177 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->
getCrsGraph().is_null(), std::logic_error,
4178 "getCrsGraph() returns null. This should not happen at this point. "
4179 "Please report this bug to the Tpetra developers.");
4188 Teuchos::RCP<const map_type> rangeMap = graph.
getRowMap();
4189 Teuchos::RCP<const map_type> domainMap = rangeMap;
4194template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4196 fillComplete(
const Teuchos::RCP<const map_type>& domainMap,
4197 const Teuchos::RCP<const map_type>& rangeMap,
4198 const Teuchos::RCP<Teuchos::ParameterList>& params) {
4202 using Teuchos::ArrayRCP;
4205 const char tfecfFuncName[] =
"fillComplete: ";
4206 ProfilingRegion regionFillComplete(
"Tpetra::CrsMatrix::fillComplete");
4207 const bool verbose = Behavior::verbose(
"CrsMatrix");
4208 std::unique_ptr<std::string> prefix;
4210 prefix = this->createPrefix(
"CrsMatrix",
"fillComplete(dom,ran,p)");
4211 std::ostringstream os;
4212 os << *prefix << endl;
4213 std::cerr << os.str();
4216 "Tpetra::CrsMatrix::fillCompete",
4220 "Matrix fill state must be active (isFillActive() "
4221 "must be true) before you may call fillComplete().");
4222 const int numProcs = this->
getComm()->getSize();
4232 bool assertNoNonlocalInserts =
false;
4235 bool sortGhosts =
true;
4237 if (!params.is_null()) {
4238 assertNoNonlocalInserts = params->get(
"No Nonlocal Changes",
4239 assertNoNonlocalInserts);
4240 if (params->isParameter(
"sort column map ghost gids")) {
4241 sortGhosts = params->get(
"sort column map ghost gids", sortGhosts);
4242 }
else if (params->isParameter(
"Sort column Map ghost GIDs")) {
4243 sortGhosts = params->get(
"Sort column Map ghost GIDs", sortGhosts);
4248 const bool needGlobalAssemble = !assertNoNonlocalInserts && numProcs > 1;
4250 if (!this->myGraph_.is_null()) {
4251 this->myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
4254 if (!this->getCrsGraphRef().indicesAreAllocated()) {
4263 if (needGlobalAssemble) {
4266 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numProcs == 1 &&
nonlocals_.size() > 0,
4268 "Cannot have nonlocal entries on a serial run. "
4269 "An invalid entry (i.e., with row index not in the row Map) must have "
4270 "been submitted to the CrsMatrix.");
4281#ifdef HAVE_TPETRA_DEBUG
4299 const bool domainMapsMatch =
4300 this->staticGraph_->getDomainMap()->isSameAs(*domainMap);
4301 const bool rangeMapsMatch =
4302 this->staticGraph_->getRangeMap()->isSameAs(*rangeMap);
4304 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!domainMapsMatch, std::runtime_error,
4305 "The CrsMatrix's domain Map does not match the graph's domain Map. "
4306 "The graph cannot be changed because it was given to the CrsMatrix "
4307 "constructor as const. You can fix this by passing in the graph's "
4308 "domain Map and range Map to the matrix's fillComplete call.");
4310 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!rangeMapsMatch, std::runtime_error,
4311 "The CrsMatrix's range Map does not match the graph's range Map. "
4312 "The graph cannot be changed because it was given to the CrsMatrix "
4313 "constructor as const. You can fix this by passing in the graph's "
4314 "domain Map and range Map to the matrix's fillComplete call.");
4326 this->myGraph_->setDomainRangeMaps(domainMap, rangeMap);
4329 Teuchos::Array<int> remotePIDs(0);
4330 const bool mustBuildColMap = !this->
hasColMap();
4331 if (mustBuildColMap) {
4332 this->myGraph_->makeColMap(remotePIDs);
4337 const std::pair<size_t, std::string> makeIndicesLocalResult =
4338 this->myGraph_->makeIndicesLocal(verbose);
4343 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(makeIndicesLocalResult.first != 0, std::runtime_error,
4344 makeIndicesLocalResult.second);
4346 const bool sorted = this->myGraph_->isSorted();
4347 const bool merged = this->myGraph_->isMerged();
4353 this->myGraph_->makeImportExport(remotePIDs, mustBuildColMap);
4359 const bool callGraphComputeGlobalConstants = params.get() ==
nullptr ||
4360 params->get(
"compute global constants",
true);
4361 if (callGraphComputeGlobalConstants) {
4362 this->myGraph_->computeGlobalConstants();
4364 this->myGraph_->computeLocalConstants();
4366 this->myGraph_->fillComplete_ =
true;
4367 this->myGraph_->checkInternalState();
4375 "Tpetra::CrsMatrix::fillCompete",
"checkInternalState");
4380template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4383 const Teuchos::RCP<const map_type>& rangeMap,
4384 const Teuchos::RCP<const import_type>& importer,
4385 const Teuchos::RCP<const export_type>& exporter,
4386 const Teuchos::RCP<Teuchos::ParameterList>& params) {
4387#ifdef HAVE_TPETRA_MMM_TIMINGS
4389 if (!params.is_null())
4390 label = params->get(
"Timer Label", label);
4391 std::string prefix = std::string(
"Tpetra ") + label + std::string(
": ");
4392 using Teuchos::TimeMonitor;
4394 Teuchos::TimeMonitor all(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-all")));
4397 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
4400 "Matrix fill state must be active (isFillActive() "
4401 "must be true) before calling fillComplete().");
4402 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4403 myGraph_.is_null(), std::logic_error,
"myGraph_ is null. This is not allowed.");
4406#ifdef HAVE_TPETRA_MMM_TIMINGS
4407 Teuchos::TimeMonitor graph(*TimeMonitor::getNewTimer(prefix + std::string(
"eSFC-M-Graph")));
4410 myGraph_->expertStaticFillComplete(domainMap, rangeMap, importer, exporter, params);
4414#ifdef HAVE_TPETRA_MMM_TIMINGS
4415 TimeMonitor fLGAM(*TimeMonitor::getNewTimer(prefix + std::string(
"eSFC-M-fLGAM")));
4426#ifdef HAVE_TPETRA_DEBUG
4427 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillActive(), std::logic_error,
4428 ": We're at the end of fillComplete(), but isFillActive() is true. "
4429 "Please report this bug to the Tpetra developers.");
4430 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete(), std::logic_error,
4431 ": We're at the end of fillComplete(), but isFillActive() is true. "
4432 "Please report this bug to the Tpetra developers.");
4435#ifdef HAVE_TPETRA_MMM_TIMINGS
4436 Teuchos::TimeMonitor cIS(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cIS")));
4443template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4448 LocalOrdinal* beg = cols;
4449 LocalOrdinal* end = cols + rowLen;
4450 LocalOrdinal* newend = beg;
4452 LocalOrdinal* cur = beg + 1;
4456 while (cur != end) {
4457 if (*cur != *newend) {
4473 return newend - beg;
4476template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4480 typedef LocalOrdinal LO;
4481 typedef typename Kokkos::View<LO*, device_type>::host_mirror_type::execution_space
4482 host_execution_space;
4483 typedef Kokkos::RangePolicy<host_execution_space, LO> range_type;
4484 const char tfecfFuncName[] =
"sortAndMergeIndicesAndValues: ";
4485 ProfilingRegion regionSAM(
"Tpetra::CrsMatrix::sortAndMergeIndicesAndValues");
4487 if (!sorted || !merged) {
4488 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->
isStaticGraph(), std::runtime_error,
4489 "Cannot sort or merge with "
4490 "\"static\" (const) graph, since the matrix does not own the graph.");
4491 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->myGraph_.is_null(), std::logic_error,
4492 "myGraph_ is null, but "
4493 "this matrix claims ! isStaticGraph(). "
4494 "Please report this bug to the Tpetra developers.");
4496 "It is invalid to call "
4497 "this method if the graph's storage has already been optimized. "
4498 "Please report this bug to the Tpetra developers.");
4502 size_t totalNumDups = 0;
4507 auto vals_ = this->valuesUnpacked_wdv.getHostView(Access::ReadWrite);
4509 Kokkos::parallel_reduce(
4510 "sortAndMergeIndicesAndValues", range_type(0, lclNumRows),
4511 [=](
const LO lclRow,
size_t& numDups) {
4512 size_t rowBegin = rowBegins_(lclRow);
4513 size_t rowLen = rowLengths_(lclRow);
4514 LO* cols = cols_.data() + rowBegin;
4517 sort2(cols, cols + rowLen, vals);
4521 rowLengths_(lclRow) = newRowLength;
4522 numDups += rowLen - newRowLength;
4536template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4541 Scalar beta)
const {
4544 using Teuchos::rcp_const_cast;
4545 using Teuchos::rcpFromRef;
4547 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
4548 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
4554 if (alpha ==
ZERO) {
4557 }
else if (beta != ONE) {
4571 RCP<const import_type> importer = this->
getGraph()->getImporter();
4572 RCP<const export_type> exporter = this->
getGraph()->getExporter();
4578 const bool Y_is_overwritten = (beta ==
ZERO);
4581 const bool Y_is_replicated =
4590 if (Y_is_replicated && this->
getComm()->getRank() > 0) {
4597 RCP<const MV> X_colMap;
4598 if (importer.is_null()) {
4608 X_colMap = rcp_const_cast<const MV>(X_colMapNonConst);
4612 X_colMap = rcpFromRef(X_in);
4615 ProfilingRegion regionImport(
"Tpetra::CrsMatrix::apply: Import");
4624 X_colMapNonConst->doImport(X_in, *importer,
INSERT);
4625 X_colMap = rcp_const_cast<const MV>(X_colMapNonConst);
4639 if (!exporter.is_null()) {
4640 this->
localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha,
ZERO);
4642 ProfilingRegion regionExport(
"Tpetra::CrsMatrix::apply: Export");
4648 if (Y_is_overwritten) {
4679 this->
localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
4682 this->
localApply(*X_colMap, Y_in, Teuchos::NO_TRANS, alpha, beta);
4690 if (Y_is_replicated) {
4691 ProfilingRegion regionReduce(
"Tpetra::CrsMatrix::apply: Reduce Y");
4696template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4700 const Teuchos::ETransp mode,
4702 Scalar beta)
const {
4703 using Teuchos::null;
4706 using Teuchos::rcp_const_cast;
4707 using Teuchos::rcpFromRef;
4709 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
4712 if (alpha ==
ZERO) {
4734 RCP<const import_type> importer = this->
getGraph()->getImporter();
4735 RCP<const export_type> exporter = this->
getGraph()->getExporter();
4741 const bool Y_is_overwritten = (beta ==
ZERO);
4742 if (Y_is_replicated && this->
getComm()->getRank() > 0) {
4748 X = rcp(
new MV(X_in, Teuchos::Copy));
4750 X = rcpFromRef(X_in);
4754 if (importer != Teuchos::null) {
4762 if (exporter != Teuchos::null) {
4773 if (!exporter.is_null()) {
4774 ProfilingRegion regionImport(
"Tpetra::CrsMatrix::apply (transpose): Import");
4782 if (importer != Teuchos::null) {
4783 ProfilingRegion regionExport(
"Tpetra::CrsMatrix::apply (transpose): Export");
4794 if (Y_is_overwritten) {
4811 MV Y(Y_in, Teuchos::Copy);
4815 this->
localApply(*X, Y_in, mode, alpha, beta);
4822 if (Y_is_replicated) {
4823 ProfilingRegion regionReduce(
"Tpetra::CrsMatrix::apply (transpose): Reduce Y");
4828template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4832 const Teuchos::ETransp mode,
4833 const Scalar& alpha,
4834 const Scalar& beta)
const {
4835 using Teuchos::NO_TRANS;
4837 ProfilingRegion regionLocalApply(
"Tpetra::CrsMatrix::localApply");
4844 const char tfecfFuncName[] =
"localApply: ";
4847 "Y.getNumVectors() = "
4849 const bool transpose = (mode != Teuchos::NO_TRANS);
4850 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!transpose && X.
getLocalLength() !=
4853 "NO_TRANS case: X has the wrong number of local rows. "
4854 "X.getLocalLength() = "
4856 "getColMap()->getLocalNumElements() = "
4857 <<
getColMap()->getLocalNumElements() <<
".");
4858 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!transpose && Y.
getLocalLength() !=
4861 "NO_TRANS case: Y has the wrong number of local rows. "
4862 "Y.getLocalLength() = "
4864 "getRowMap()->getLocalNumElements() = "
4865 <<
getRowMap()->getLocalNumElements() <<
".");
4866 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(transpose && X.
getLocalLength() !=
4869 "TRANS or CONJ_TRANS case: X has the wrong number of local "
4870 "rows. X.getLocalLength() = "
4872 <<
" != getRowMap()->getLocalNumElements() = "
4873 <<
getRowMap()->getLocalNumElements() <<
".");
4874 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(transpose && Y.
getLocalLength() !=
4877 "TRANS or CONJ_TRANS case: X has the wrong number of local "
4878 "rows. Y.getLocalLength() = "
4880 <<
" != getColMap()->getLocalNumElements() = "
4881 <<
getColMap()->getLocalNumElements() <<
".");
4882 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete(), std::runtime_error,
4883 "The matrix is not "
4884 "fill complete. You must call fillComplete() (possibly with "
4885 "domain and range Map arguments) without an intervening "
4886 "resumeFill() call before you may call this method.");
4888 std::runtime_error,
"X and Y must be constant stride.");
4893 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(X_lcl.data() == Y_lcl.data() && X_lcl.data() !=
nullptr && X_lcl.extent(0) != 0,
4894 std::runtime_error,
"X and Y may not alias one another.");
4899 if (!applyHelper.get()) {
4902 bool useMergePath =
false;
4903#ifdef KOKKOSKERNELS_ENABLE_TPL_CUSPARSE
4909 if constexpr (std::is_same_v<execution_space, Kokkos::Cuda>) {
4911 LocalOrdinal maxRowImbalance = 0;
4916 useMergePath =
true;
4919 applyHelper = std::make_shared<ApplyHelper>(A_lcl.nnz(), A_lcl.graph.row_map,
4920 useMergePath ? KokkosSparse::SPMV_MERGE_PATH : KokkosSparse::SPMV_DEFAULT);
4924 const char* modeKK =
nullptr;
4926 case Teuchos::NO_TRANS:
4927 modeKK = KokkosSparse::NoTranspose;
4929 case Teuchos::TRANS:
4930 modeKK = KokkosSparse::Transpose;
4932 case Teuchos::CONJ_TRANS:
4933 modeKK = KokkosSparse::ConjugateTranspose;
4936 throw std::invalid_argument(
"Tpetra::CrsMatrix::localApply: invalid mode");
4939 if (applyHelper->shouldUseIntRowptrs()) {
4940 auto A_lcl_int_rowptrs = applyHelper->getIntRowptrMatrix(A_lcl);
4942 &applyHelper->handle_int, modeKK,
4946 &applyHelper->handle, modeKK,
4951template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4955 Teuchos::ETransp mode,
4957 Scalar beta)
const {
4959 const char fnName[] =
"Tpetra::CrsMatrix::apply";
4961 TEUCHOS_TEST_FOR_EXCEPTION(!
isFillComplete(), std::runtime_error,
4962 fnName <<
": Cannot call apply() until fillComplete() "
4963 "has been called.");
4965 if (mode == Teuchos::NO_TRANS) {
4966 ProfilingRegion regionNonTranspose(fnName);
4969 ProfilingRegion regionTranspose(
"Tpetra::CrsMatrix::apply (transpose)");
4974template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4976Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node>>
4980 typedef CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node> output_matrix_type;
4981 const char tfecfFuncName[] =
"convert: ";
4983 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!this->
isFillComplete(), std::runtime_error,
4984 "This matrix (the source "
4985 "of the conversion) is not fill complete. You must first call "
4986 "fillComplete() (possibly with the domain and range Map) without an "
4987 "intervening call to resumeFill(), before you may call this method.");
4989 RCP<output_matrix_type> newMatrix(
new output_matrix_type(this->
getCrsGraph()));
4993 copyConvert(newMatrix->getLocalMatrixDevice().values,
4994 this->getLocalMatrixDevice().values);
5003template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5008 const char tfecfFuncName[] =
"checkInternalState: ";
5010 "Internal state is not consistent. "
5011 "Please report this bug to the Tpetra developers.";
5015 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(staticGraph_.is_null(), std::logic_error, err);
5019 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!myGraph_.is_null() && myGraph_ != staticGraph_,
5020 std::logic_error, err);
5022 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillComplete() && !staticGraph_->isFillComplete(),
5023 std::logic_error, err <<
" Specifically, the matrix is fill complete, "
5024 "but its graph is NOT fill complete.");
5027 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(staticGraph_->indicesAreAllocated() &&
5028 staticGraph_->getLocalAllocationSize() > 0 &&
5029 staticGraph_->getLocalNumRows() > 0 &&
5030 valuesUnpacked_wdv.extent(0) == 0,
5031 std::logic_error, err);
5035template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5039 std::ostringstream os;
5041 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
5042 if (this->getObjectLabel() !=
"") {
5043 os <<
"Label: \"" << this->getObjectLabel() <<
"\", ";
5046 os <<
"isFillComplete: true"
5052 os <<
"isFillComplete: false"
5059template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5061 describe(Teuchos::FancyOStream& out,
5062 const Teuchos::EVerbosityLevel verbLevel)
const {
5065 using Teuchos::ArrayView;
5066 using Teuchos::Comm;
5068 using Teuchos::TypeNameTraits;
5069 using Teuchos::VERB_DEFAULT;
5070 using Teuchos::VERB_EXTREME;
5071 using Teuchos::VERB_HIGH;
5072 using Teuchos::VERB_LOW;
5073 using Teuchos::VERB_MEDIUM;
5074 using Teuchos::VERB_NONE;
5076 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
5078 if (vl == VERB_NONE) {
5083 Teuchos::OSTab tab0(out);
5085 RCP<const Comm<int>> comm = this->
getComm();
5086 const int myRank = comm->getRank();
5087 const int numProcs = comm->getSize();
5092 width = std::max<size_t>(width,
static_cast<size_t>(11)) + 2;
5102 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
5104 Teuchos::OSTab tab1(out);
5107 if (this->getObjectLabel() !=
"") {
5108 out <<
"Label: \"" << this->getObjectLabel() <<
"\", ";
5111 out <<
"Template parameters:" << endl;
5112 Teuchos::OSTab tab2(out);
5113 out <<
"Scalar: " << TypeNameTraits<Scalar>::name() << endl
5114 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name() << endl
5115 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name() << endl
5116 <<
"Node: " << TypeNameTraits<Node>::name() << endl;
5119 out <<
"isFillComplete: true" << endl
5124 <<
"Global max number of entries in a row: "
5127 out <<
"isFillComplete: false" << endl
5133 if (vl < VERB_MEDIUM) {
5140 <<
"Row Map:" << endl;
5144 out <<
"null" << endl;
5155 out <<
"Column Map: ";
5159 out <<
"null" << endl;
5163 out <<
"same as row Map" << endl;
5174 out <<
"Domain Map: ";
5178 out <<
"null" << endl;
5182 out <<
"same as row Map" << endl;
5186 out <<
"same as column Map" << endl;
5197 out <<
"Range Map: ";
5201 out <<
"null" << endl;
5205 out <<
"same as domain Map" << endl;
5209 out <<
"same as row Map" << endl;
5219 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5220 if (myRank == curRank) {
5221 out <<
"Process rank: " << curRank << endl;
5222 Teuchos::OSTab tab2(out);
5223 if (!staticGraph_->indicesAreAllocated()) {
5224 out <<
"Graph indices not allocated" << endl;
5226 out <<
"Number of allocated entries: "
5227 << staticGraph_->getLocalAllocationSize() << endl;
5239 if (vl < VERB_HIGH) {
5244 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5245 if (myRank == curRank) {
5246 out << std::setw(width) <<
"Proc Rank"
5247 << std::setw(width) <<
"Global Row"
5248 << std::setw(width) <<
"Num Entries";
5249 if (vl == VERB_EXTREME) {
5250 out << std::setw(width) <<
"(Index,Value)";
5255 GlobalOrdinal gid =
getRowMap()->getGlobalElement(r);
5256 out << std::setw(width) << myRank
5257 << std::setw(width) << gid
5258 << std::setw(width) << nE;
5259 if (vl == VERB_EXTREME) {
5261 global_inds_host_view_type rowinds;
5262 values_host_view_type rowvals;
5264 for (
size_t j = 0; j < nE; ++j) {
5265 out <<
" (" << rowinds[j]
5266 <<
", " << rowvals[j]
5270 local_inds_host_view_type rowinds;
5271 values_host_view_type rowvals;
5273 for (
size_t j = 0; j < nE; ++j) {
5274 out <<
" (" <<
getColMap()->getGlobalElement(rowinds[j])
5275 <<
", " << rowvals[j]
5291template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5303 return (srcRowMat !=
nullptr);
5306template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5309 const typename crs_graph_type::padding_type& padding,
5310 const bool verbose) {
5315 using row_ptrs_type =
5316 typename local_graph_device_type::row_map_type::non_const_type;
5317 using range_policy =
5318 Kokkos::RangePolicy<execution_space, Kokkos::IndexType<LO>>;
5319 const char tfecfFuncName[] =
"applyCrsPadding";
5320 const char suffix[] =
5321 ". Please report this bug to the Tpetra developers.";
5322 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::applyCrsPadding");
5324 std::unique_ptr<std::string> prefix;
5326 prefix = this->createPrefix(
"CrsMatrix", tfecfFuncName);
5327 std::ostringstream os;
5328 os << *prefix <<
"padding: ";
5331 std::cerr << os.str();
5333 const int myRank = !verbose ? -1 : [&]() {
5334 auto map = this->getMap();
5335 if (map.is_null()) {
5338 auto comm = map->getComm();
5339 if (comm.is_null()) {
5342 return comm->getRank();
5346 if (!myGraph_->indicesAreAllocated()) {
5348 std::ostringstream os;
5349 os << *prefix <<
"Call allocateIndices" << endl;
5350 std::cerr << os.str();
5352 allocateValues(GlobalIndices, GraphNotYetAllocated, verbose);
5364 std::ostringstream os;
5365 os << *prefix <<
"Allocate row_ptrs_beg: "
5366 << myGraph_->getRowPtrsUnpackedHost().extent(0) << endl;
5367 std::cerr << os.str();
5369 using Kokkos::view_alloc;
5370 using Kokkos::WithoutInitializing;
5371 row_ptrs_type row_ptr_beg(view_alloc(
"row_ptr_beg", WithoutInitializing),
5372 myGraph_->rowPtrsUnpacked_dev_.extent(0));
5374 Kokkos::deep_copy(
execution_space(), row_ptr_beg, myGraph_->rowPtrsUnpacked_dev_);
5376 const size_t N = row_ptr_beg.extent(0) == 0 ? size_t(0) : size_t(row_ptr_beg.extent(0) - 1);
5378 std::ostringstream os;
5379 os << *prefix <<
"Allocate row_ptrs_end: " << N << endl;
5380 std::cerr << os.str();
5382 row_ptrs_type row_ptr_end(
5383 view_alloc(
"row_ptr_end", WithoutInitializing), N);
5385 row_ptrs_type num_row_entries_d;
5387 const bool refill_num_row_entries =
5388 myGraph_->k_numRowEntries_.extent(0) != 0;
5390 if (refill_num_row_entries) {
5393 num_row_entries_d = create_mirror_view_and_copy(memory_space(),
5394 myGraph_->k_numRowEntries_);
5395 Kokkos::parallel_for(
5396 "Fill end row pointers", range_policy(0, N),
5397 KOKKOS_LAMBDA(
const size_t i) {
5398 row_ptr_end(i) = row_ptr_beg(i) + num_row_entries_d(i);
5404 Kokkos::parallel_for(
5405 "Fill end row pointers", range_policy(0, N),
5406 KOKKOS_LAMBDA(
const size_t i) {
5407 row_ptr_end(i) = row_ptr_beg(i + 1);
5411 if (myGraph_->isGloballyIndexed()) {
5413 myGraph_->gblInds_wdv,
5414 valuesUnpacked_wdv, padding, myRank, verbose);
5415 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5416 const auto newColIndsLen = myGraph_->gblInds_wdv.extent(0);
5417 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(newValuesLen != newColIndsLen, std::logic_error,
5418 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5419 <<
" != myGraph_->gblInds_wdv.extent(0)=" << newColIndsLen
5423 myGraph_->lclIndsUnpacked_wdv,
5424 valuesUnpacked_wdv, padding, myRank, verbose);
5425 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5426 const auto newColIndsLen = myGraph_->lclIndsUnpacked_wdv.extent(0);
5427 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(newValuesLen != newColIndsLen, std::logic_error,
5428 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5429 <<
" != myGraph_->lclIndsUnpacked_wdv.extent(0)=" << newColIndsLen
5433 if (refill_num_row_entries) {
5434 Kokkos::parallel_for(
5435 "Fill num entries", range_policy(0, N),
5436 KOKKOS_LAMBDA(
const size_t i) {
5437 num_row_entries_d(i) = row_ptr_end(i) - row_ptr_beg(i);
5439 Kokkos::deep_copy(myGraph_->k_numRowEntries_, num_row_entries_d);
5443 std::ostringstream os;
5444 os << *prefix <<
"Assign myGraph_->rowPtrsUnpacked_; "
5445 <<
"old size: " << myGraph_->rowPtrsUnpacked_host_.extent(0)
5446 <<
", new size: " << row_ptr_beg.extent(0) << endl;
5447 std::cerr << os.str();
5448 TEUCHOS_ASSERT(myGraph_->getRowPtrsUnpackedHost().extent(0) ==
5449 row_ptr_beg.extent(0));
5451 myGraph_->setRowPtrsUnpacked(row_ptr_beg);
5454template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5455void copyAndPermuteStaticGraphNew(
5458 const size_t numSameIDs,
5459 const LocalOrdinal permuteToLIDs[],
5460 const LocalOrdinal permuteFromLIDs[],
5461 const size_t numPermutes);
5463template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5467 const size_t numSameIDs,
5468 const LocalOrdinal permuteToLIDs[],
5469 const LocalOrdinal permuteFromLIDs[],
5470 const size_t numPermutes) {
5473 using Teuchos::Array;
5474 using Teuchos::ArrayView;
5475 using LO = LocalOrdinal;
5476 using GO = GlobalOrdinal;
5477 const char tfecfFuncName[] =
"copyAndPermuteStaticGraph";
5478 const char suffix[] =
5479 " Please report this bug to the Tpetra developers.";
5480 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteStaticGraph");
5484 std::unique_ptr<std::string> prefix;
5486 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5487 std::ostringstream os;
5488 os << *prefix <<
"Start" << endl;
5490 const char*
const prefix_raw =
5491 verbose ? prefix.get()->c_str() :
nullptr;
5493 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed();
5498 const auto& srcRowMap = *(srcMat.getRowMap());
5499 nonconst_global_inds_host_view_type rowInds;
5500 nonconst_values_host_view_type rowVals;
5501 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
5502 if (sourceIsLocallyIndexed) {
5503 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5507 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5508 const GO targetGID = sourceGID;
5510 ArrayView<const GO> rowIndsConstView;
5511 ArrayView<const Scalar> rowValsConstView;
5513 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5514 if (rowLength >
static_cast<size_t>(rowInds.size())) {
5515 Kokkos::resize(rowInds, rowLength);
5516 Kokkos::resize(rowVals, rowLength);
5520 nonconst_global_inds_host_view_type rowIndsView;
5521 nonconst_values_host_view_type rowValsView;
5526 size_t checkRowLength = 0;
5529 const crs_matrix_type* srcMatCrsPtr =
dynamic_cast<const crs_matrix_type*
>(&srcMat);
5530 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(srcMatCrsPtr ==
nullptr, std::runtime_error,
"bad srcMatCrsPtr");
5531 const crs_matrix_type& srcMatCrs = *srcMatCrsPtr;
5533 auto globalRow = sourceGID;
5534 auto StaticGraphRCP = srcMatCrs.getGraph();
5536 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(StaticGraphPtr ==
nullptr, std::runtime_error,
"bad StaticGraphPtr");
5538 const RowInfo rowinfo = StaticGraph.getRowInfoFromGlobalRowIndex(globalRow);
5539 const size_t theNumEntries = rowinfo.numEntries;
5540 checkRowLength = theNumEntries;
5541 auto numEntries = theNumEntries;
5543 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid()) {
5544 if (StaticGraph.isLocallyIndexed()) {
5545 const map_type& colMap = *(StaticGraph.getColMap());
5546 auto curLclInds = StaticGraph.getLocalIndsViewHost(rowinfo);
5547 auto rowValsViewLocal = srcMatCrs.getValuesViewHost(rowinfo);
5548 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5549 reinterpret_cast<const Scalar*
>(rowValsViewLocal.data()),
5550 rowValsViewLocal.extent(0),
5551 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5552 auto rowIndsViewLocal = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5553 rowIndsConstView = Teuchos::ArrayView<const GO>(
5554 rowIndsViewLocal.data(), rowIndsViewLocal.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
5555 bool err = colMap.getGlobalElements(curLclInds.data(), numEntries, rowIndsViewLocal.data());
5556 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(err, std::runtime_error,
"getGlobalElements error");
5557 }
else if (StaticGraph.isGloballyIndexed()) {
5558 auto rowIndsViewLocal = StaticGraph.getGlobalIndsViewHost(rowinfo);
5559 rowIndsConstView = Teuchos::ArrayView<const GO>(
5560 rowIndsViewLocal.data(), rowIndsViewLocal.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
5561 auto rowValsViewLocal = srcMatCrs.getValuesViewHost(rowinfo);
5562 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5563 reinterpret_cast<const Scalar*
>(rowValsViewLocal.data()),
5564 rowValsViewLocal.extent(0),
5565 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5570 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5571 rowLength != checkRowLength,
5573 "For global row index " << sourceGID <<
", the source matrix's getNumEntriesInGlobalRow returns a row length of " << rowLength <<
", but getGlobalRowCopy reports a row length of " << checkRowLength <<
"." << suffix);
5576 combineGlobalValues(
5577 targetGID, rowIndsConstView, rowValsConstView,
REPLACE, prefix_raw, debug, verbose);
5580 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5584 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5585 const GO targetGID = sourceGID;
5587 ArrayView<const GO> rowIndsConstView;
5588 ArrayView<const Scalar> rowValsConstView;
5590 global_inds_host_view_type rowIndsView;
5591 values_host_view_type rowValsView;
5592 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5597 rowIndsConstView = Teuchos::ArrayView<const GO>(
5598 rowIndsView.data(), rowIndsView.extent(0),
5599 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5600 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5601 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5602 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5607 combineGlobalValues(targetGID, rowIndsConstView,
5609 prefix_raw, debug, verbose);
5614 std::ostringstream os;
5615 os << *prefix <<
"Do permutes" << endl;
5622 const map_type& tgtRowMap = *(this->getRowMap());
5623 for (
size_t p = 0; p < numPermutes; ++p) {
5624 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
5625 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
5627 ArrayView<const GO> rowIndsConstView;
5628 ArrayView<const Scalar> rowValsConstView;
5630 if (sourceIsLocallyIndexed) {
5631 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5632 if (rowLength >
static_cast<size_t>(rowInds.size())) {
5633 Kokkos::resize(rowInds, rowLength);
5634 Kokkos::resize(rowVals, rowLength);
5638 nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5639 nonconst_values_host_view_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5644 size_t checkRowLength = 0;
5645 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
5646 rowValsView, checkRowLength);
5648 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5650 "source matrix global row index "
5651 << sourceGID <<
", "
5652 "getNumEntriesInGlobalRow returns a row length of "
5653 << rowLength <<
", but getGlobalRowCopy a row length of "
5654 << checkRowLength <<
"." << suffix);
5661 rowIndsConstView = Teuchos::ArrayView<const GO>(
5662 rowIndsView.data(), rowIndsView.extent(0),
5663 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5664 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5665 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5666 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5670 global_inds_host_view_type rowIndsView;
5671 values_host_view_type rowValsView;
5672 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5677 rowIndsConstView = Teuchos::ArrayView<const GO>(
5678 rowIndsView.data(), rowIndsView.extent(0),
5679 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5680 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5681 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5682 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5687 combineGlobalValues(targetGID, rowIndsConstView,
5689 prefix_raw, debug, verbose);
5693 std::ostringstream os;
5694 os << *prefix <<
"Done" << endl;
5698template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5702 const size_t numSameIDs,
5703 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs_dv,
5704 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs_dv,
5705 const size_t numPermutes) {
5708 using Teuchos::Array;
5709 using Teuchos::ArrayView;
5710 using LO = LocalOrdinal;
5711 using GO = GlobalOrdinal;
5712 const char tfecfFuncName[] =
"copyAndPermuteNonStaticGraph";
5713 const char suffix[] =
5714 " Please report this bug to the Tpetra developers.";
5715 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteNonStaticGraph");
5719 std::unique_ptr<std::string> prefix;
5721 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5722 std::ostringstream os;
5723 os << *prefix <<
"Start" << endl;
5725 const char*
const prefix_raw =
5726 verbose ? prefix.get()->c_str() :
nullptr;
5730 const row_graph_type& srcGraph = *(srcMat.getGraph());
5732 myGraph_->computeCrsPadding(srcGraph, numSameIDs,
5733 permuteToLIDs_dv, permuteFromLIDs_dv, verbose);
5734 applyCrsPadding(*padding, verbose);
5736 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed();
5741 const map_type& srcRowMap = *(srcMat.getRowMap());
5742 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
5743 using gids_type = nonconst_global_inds_host_view_type;
5744 using vals_type = nonconst_values_host_view_type;
5747 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5751 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5752 const GO targetGID = sourceGID;
5754 ArrayView<const GO> rowIndsConstView;
5755 ArrayView<const Scalar> rowValsConstView;
5757 if (sourceIsLocallyIndexed) {
5758 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5759 if (rowLength >
static_cast<size_t>(rowInds.extent(0))) {
5760 Kokkos::resize(rowInds, rowLength);
5761 Kokkos::resize(rowVals, rowLength);
5765 gids_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5766 vals_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5771 size_t checkRowLength = 0;
5772 srcMat.getGlobalRowCopy(sourceGID, rowIndsView, rowValsView,
5775 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5778 << sourceGID <<
", the source "
5779 "matrix's getNumEntriesInGlobalRow returns a row length "
5781 << rowLength <<
", but getGlobalRowCopy reports "
5783 << checkRowLength <<
"." << suffix);
5785 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
5786 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar*
>(rowValsView.data()), rowLength);
5788 global_inds_host_view_type rowIndsView;
5789 values_host_view_type rowValsView;
5790 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5796 rowIndsConstView = Teuchos::ArrayView<const GO>(
5797 rowIndsView.data(), rowIndsView.extent(0),
5798 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5799 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5800 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5801 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5807 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
5808 rowValsConstView, prefix_raw, debug, verbose);
5812 std::ostringstream os;
5813 os << *prefix <<
"Do permutes" << endl;
5815 const LO*
const permuteFromLIDs = permuteFromLIDs_dv.view_host().data();
5816 const LO*
const permuteToLIDs = permuteToLIDs_dv.view_host().data();
5818 const map_type& tgtRowMap = *(this->getRowMap());
5819 for (
size_t p = 0; p < numPermutes; ++p) {
5820 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
5821 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
5823 ArrayView<const GO> rowIndsConstView;
5824 ArrayView<const Scalar> rowValsConstView;
5826 if (sourceIsLocallyIndexed) {
5827 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5828 if (rowLength >
static_cast<size_t>(rowInds.extent(0))) {
5829 Kokkos::resize(rowInds, rowLength);
5830 Kokkos::resize(rowVals, rowLength);
5834 gids_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5835 vals_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5840 size_t checkRowLength = 0;
5841 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
5842 rowValsView, checkRowLength);
5844 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5846 "source matrix global row index "
5847 << sourceGID <<
", "
5848 "getNumEntriesInGlobalRow returns a row length of "
5849 << rowLength <<
", but getGlobalRowCopy a row length of "
5850 << checkRowLength <<
"." << suffix);
5852 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
5853 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar*
>(rowValsView.data()), rowLength);
5855 global_inds_host_view_type rowIndsView;
5856 values_host_view_type rowValsView;
5857 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5863 rowIndsConstView = Teuchos::ArrayView<const GO>(
5864 rowIndsView.data(), rowIndsView.extent(0),
5865 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5866 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5867 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5868 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5874 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
5875 rowValsConstView, prefix_raw, debug, verbose);
5879 std::ostringstream os;
5880 os << *prefix <<
"Done" << endl;
5884template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5888 const size_t numSameIDs,
5889 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs,
5890 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs,
5898 const char tfecfFuncName[] =
"copyAndPermute: ";
5899 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermute");
5901 const bool verbose = Behavior::verbose(
"CrsMatrix");
5902 std::unique_ptr<std::string> prefix;
5904 prefix = this->createPrefix(
"CrsMatrix",
"copyAndPermute");
5905 std::ostringstream os;
5906 os << *prefix << endl
5907 << *prefix <<
" numSameIDs: " << numSameIDs << endl
5908 << *prefix <<
" numPermute: " << permuteToLIDs.extent(0)
5911 << dualViewStatusToString(permuteToLIDs,
"permuteToLIDs")
5914 << dualViewStatusToString(permuteFromLIDs,
"permuteFromLIDs")
5917 <<
"isStaticGraph: " << (
isStaticGraph() ?
"true" :
"false")
5919 std::cerr << os.str();
5922 const auto numPermute = permuteToLIDs.extent(0);
5923 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numPermute != permuteFromLIDs.extent(0),
5924 std::invalid_argument,
"permuteToLIDs.extent(0) = " << numPermute <<
"!= permuteFromLIDs.extent(0) = " << permuteFromLIDs.extent(0) <<
".");
5929 const RMT& srcMat =
dynamic_cast<const RMT&
>(srcObj);
5932 TEUCHOS_ASSERT(!permuteToLIDs.need_sync_device());
5933 auto permuteToLIDs_d = permuteToLIDs.view_device();
5934 TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_device());
5935 auto permuteFromLIDs_d = permuteFromLIDs.view_device();
5936 copyAndPermuteStaticGraphNew(
5937 srcMat, *
this, numSameIDs, permuteToLIDs_d.data(), permuteFromLIDs_d.data(), numPermute);
5940 TEUCHOS_ASSERT(!permuteToLIDs.need_sync_host());
5941 auto permuteToLIDs_h = permuteToLIDs.view_host();
5942 TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_host());
5943 auto permuteFromLIDs_h = permuteFromLIDs.view_host();
5945 copyAndPermuteStaticGraph(srcMat, numSameIDs,
5946 permuteToLIDs_h.data(),
5947 permuteFromLIDs_h.data(),
5951 copyAndPermuteNonStaticGraph(srcMat, numSameIDs, permuteToLIDs,
5952 permuteFromLIDs, numPermute);
5956 std::ostringstream os;
5957 os << *prefix <<
"Done" << endl;
5958 std::cerr << os.str();
5962template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5965 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
5966 Kokkos::DualView<char*, buffer_device_type>& exports,
5967 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
5968 size_t& constantNumPackets) {
5973 using Teuchos::outArg;
5974 using Teuchos::REDUCE_MAX;
5975 using Teuchos::reduceAll;
5976 typedef LocalOrdinal LO;
5977 typedef GlobalOrdinal GO;
5978 const char tfecfFuncName[] =
"packAndPrepare: ";
5979 ProfilingRegion regionPAP(
"Tpetra::CrsMatrix::packAndPrepare");
5981 const bool debug = Behavior::debug(
"CrsMatrix");
5982 const bool verbose = Behavior::verbose(
"CrsMatrix");
5985 Teuchos::RCP<const Teuchos::Comm<int>> pComm = this->getComm();
5986 if (pComm.is_null()) {
5989 const Teuchos::Comm<int>& comm = *pComm;
5990 const int myRank = comm.getSize();
5992 std::unique_ptr<std::string> prefix;
5994 prefix = this->createPrefix(
"CrsMatrix",
"packAndPrepare");
5995 std::ostringstream os;
5996 os << *prefix <<
"Start" << endl
5998 << dualViewStatusToString(exportLIDs,
"exportLIDs")
6001 << dualViewStatusToString(exports,
"exports")
6004 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6006 std::cerr << os.str();
6029 std::ostringstream msg;
6033 const crs_matrix_type* srcCrsMat =
6034 dynamic_cast<const crs_matrix_type*
>(&source);
6035 if (srcCrsMat !=
nullptr) {
6037 std::ostringstream os;
6038 os << *prefix <<
"Source matrix same (CrsMatrix) type as target; "
6041 std::cerr << os.str();
6044 srcCrsMat->packNew(exportLIDs, exports, numPacketsPerLID,
6045 constantNumPackets);
6046 }
catch (std::exception& e) {
6048 msg <<
"Proc " << myRank <<
": " << e.what() << std::endl;
6051 using Kokkos::HostSpace;
6052 using Kokkos::subview;
6053 using exports_type = Kokkos::DualView<char*, buffer_device_type>;
6054 using range_type = Kokkos::pair<size_t, size_t>;
6057 std::ostringstream os;
6058 os << *prefix <<
"Source matrix NOT same (CrsMatrix) type as target"
6060 std::cerr << os.str();
6063 const row_matrix_type* srcRowMat =
6064 dynamic_cast<const row_matrix_type*
>(&source);
6065 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(srcRowMat ==
nullptr, std::invalid_argument,
6066 "The source object of the Import or Export operation is neither a "
6067 "CrsMatrix (with the same template parameters as the target object), "
6068 "nor a RowMatrix (with the same first four template parameters as the "
6079 TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
6080 auto exportLIDs_h = exportLIDs.view_host();
6081 Teuchos::ArrayView<const LO> exportLIDs_av(exportLIDs_h.data(),
6082 exportLIDs_h.size());
6086 Teuchos::Array<char> exports_a;
6092 numPacketsPerLID.clear_sync_state();
6093 numPacketsPerLID.modify_host();
6094 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
6095 Teuchos::ArrayView<size_t> numPacketsPerLID_av(numPacketsPerLID_h.data(),
6096 numPacketsPerLID_h.size());
6101 srcRowMat->pack(exportLIDs_av, exports_a, numPacketsPerLID_av,
6102 constantNumPackets);
6103 }
catch (std::exception& e) {
6105 msg <<
"Proc " << myRank <<
": " << e.what() << std::endl;
6109 const size_t newAllocSize =
static_cast<size_t>(exports_a.size());
6110 if (
static_cast<size_t>(exports.extent(0)) < newAllocSize) {
6111 const std::string oldLabel = exports.view_device().label();
6112 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6113 exports = exports_type(newLabel, newAllocSize);
6118 exports.modify_host();
6120 auto exports_h = exports.view_host();
6121 auto exports_h_sub = subview(exports_h, range_type(0, newAllocSize));
6125 typedef typename exports_type::t_host::execution_space HES;
6126 typedef Kokkos::Device<HES, HostSpace> host_device_type;
6127 Kokkos::View<const char*, host_device_type>
6128 exports_a_kv(exports_a.getRawPtr(), newAllocSize);
6130 Kokkos::deep_copy(exports_h_sub, exports_a_kv);
6135 reduceAll<int, int>(comm, REDUCE_MAX, lclBad, outArg(gblBad));
6138 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error,
6139 "packNew() or pack() threw an exception on "
6140 "one or more participating processes.");
6143 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(lclBad != 0, std::logic_error,
6144 "packNew threw an exception on one "
6145 "or more participating processes. Here is this process' error "
6151 std::ostringstream os;
6152 os << *prefix <<
"packAndPrepare: Done!" << endl
6162 std::cerr << os.str();
6166template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6170 const size_t offset,
6171 const size_t numEnt,
6172 const GlobalOrdinal gidsIn[],
6174 const size_t numBytesPerValue)
const {
6175 using Kokkos::subview;
6177 using Tpetra::Details::PackTraits;
6178 typedef LocalOrdinal LO;
6179 typedef GlobalOrdinal GO;
6188 const LO numEntLO =
static_cast<size_t>(numEnt);
6190 const size_t numEntBeg = offset;
6191 const size_t numEntLen = PackTraits<LO>::packValueCount(numEntLO);
6192 const size_t gidsBeg = numEntBeg + numEntLen;
6193 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
6194 const size_t valsBeg = gidsBeg + gidsLen;
6195 const size_t valsLen = numEnt * numBytesPerValue;
6197 char*
const numEntOut = exports + numEntBeg;
6198 char*
const gidsOut = exports + gidsBeg;
6199 char*
const valsOut = exports + valsBeg;
6201 size_t numBytesOut = 0;
6203 numBytesOut += PackTraits<LO>::packValue(numEntOut, numEntLO);
6206 Kokkos::pair<int, size_t> p;
6207 p = PackTraits<GO>::packArray(gidsOut, gidsIn, numEnt);
6208 errorCode += p.first;
6209 numBytesOut += p.second;
6211 p = PackTraits<ST>::packArray(valsOut, valsIn, numEnt);
6212 errorCode += p.first;
6213 numBytesOut += p.second;
6216 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6217 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
6220 << numBytesOut <<
" != expectedNumBytes = "
6221 << expectedNumBytes <<
".");
6222 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
6224 "PackTraits::packArray returned a nonzero error code");
6229template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6234 const char imports[],
6235 const size_t offset,
6236 const size_t numBytes,
6237 const size_t numEnt,
6238 const size_t numBytesPerValue) {
6239 using Kokkos::subview;
6241 using Tpetra::Details::PackTraits;
6242 typedef LocalOrdinal LO;
6243 typedef GlobalOrdinal GO;
6247 "Tpetra::CrsMatrix::unpackRow",
6250 if (numBytes == 0) {
6253 const int myRank = this->getMap()->getComm()->getRank();
6254 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6255 "unpackRow: The number of bytes to unpack numBytes=0, but the "
6256 "number of entries to unpack (as reported by numPacketsPerLID) "
6257 "for this row numEnt="
6258 << numEnt <<
" != 0.");
6263 if (numEnt == 0 && numBytes != 0) {
6264 const int myRank = this->getMap()->getComm()->getRank();
6265 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6266 "unpackRow: The number of entries to unpack (as reported by "
6267 "numPacketsPerLID) numEnt=0, but the number of bytes to unpack "
6269 << numBytes <<
" != 0.");
6275 const size_t numEntBeg = offset;
6276 const size_t numEntLen = PackTraits<LO>::packValueCount(lid);
6277 const size_t gidsBeg = numEntBeg + numEntLen;
6278 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
6279 const size_t valsBeg = gidsBeg + gidsLen;
6280 const size_t valsLen = numEnt * numBytesPerValue;
6282 const char*
const numEntIn = imports + numEntBeg;
6283 const char*
const gidsIn = imports + gidsBeg;
6284 const char*
const valsIn = imports + valsBeg;
6286 size_t numBytesOut = 0;
6289 numBytesOut += PackTraits<LO>::unpackValue(numEntOut, numEntIn);
6290 if (
static_cast<size_t>(numEntOut) != numEnt ||
6291 numEntOut ==
static_cast<LO
>(0)) {
6292 const int myRank = this->getMap()->getComm()->getRank();
6293 std::ostringstream os;
6294 os <<
"(Proc " << myRank <<
") CrsMatrix::unpackRow: ";
6295 bool firstErrorCondition =
false;
6296 if (
static_cast<size_t>(numEntOut) != numEnt) {
6297 os <<
"Number of entries from numPacketsPerLID numEnt=" << numEnt
6298 <<
" does not equal number of entries unpacked from imports "
6300 << numEntOut <<
".";
6301 firstErrorCondition =
true;
6303 if (numEntOut ==
static_cast<LO
>(0)) {
6304 if (firstErrorCondition) {
6307 os <<
"Number of entries unpacked from imports buffer numEntOut=0, "
6308 "but number of bytes to unpack for this row numBytes="
6310 <<
" != 0. This should never happen, since packRow should only "
6311 "ever pack rows with a nonzero number of entries. In this case, "
6312 "the number of entries from numPacketsPerLID is numEnt="
6316 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error, os.str());
6320 Kokkos::pair<int, size_t> p;
6321 p = PackTraits<GO>::unpackArray(gidsOut, gidsIn, numEnt);
6322 errorCode += p.first;
6323 numBytesOut += p.second;
6325 p = PackTraits<ST>::unpackArray(valsOut, valsIn, numEnt);
6326 errorCode += p.first;
6327 numBytesOut += p.second;
6330 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != numBytes, std::logic_error,
"unpackRow: numBytesOut = " << numBytesOut <<
" != numBytes = " << numBytes <<
".");
6332 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6333 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
6336 << numBytesOut <<
" != expectedNumBytes = "
6337 << expectedNumBytes <<
".");
6339 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
6341 "PackTraits::unpackArray returned a nonzero error code");
6346template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6349 size_t& totalNumEntries,
6350 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs)
const {
6355 typedef LocalOrdinal LO;
6356 typedef GlobalOrdinal GO;
6362 const bool verbose = Behavior::verbose(
"CrsMatrix");
6363 std::unique_ptr<std::string> prefix;
6365 prefix = this->
createPrefix(
"CrsMatrix",
"allocatePackSpaceNew");
6366 std::ostringstream os;
6367 os << *prefix <<
"Before:"
6375 std::cerr << os.str();
6380 const LO numExportLIDs =
static_cast<LO
>(exportLIDs.extent(0));
6382 TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
6383 auto exportLIDs_h = exportLIDs.view_host();
6386 totalNumEntries = 0;
6387 for (LO i = 0; i < numExportLIDs; ++i) {
6388 const LO lclRow = exportLIDs_h[i];
6389 size_t curNumEntries = this->getNumEntriesInLocalRow(lclRow);
6392 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid()) {
6395 totalNumEntries += curNumEntries;
6406 const size_t allocSize =
6407 static_cast<size_t>(numExportLIDs) *
sizeof(LO) +
6408 totalNumEntries * (
sizeof(IST) +
sizeof(GO));
6409 if (
static_cast<size_t>(exports.extent(0)) < allocSize) {
6410 using exports_type = Kokkos::DualView<char*, buffer_device_type>;
6412 const std::string oldLabel = exports.view_device().label();
6413 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6414 exports = exports_type(newLabel, allocSize);
6418 std::ostringstream os;
6419 os << *prefix <<
"After:"
6427 std::cerr << os.str();
6431template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6433 packNew(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6434 Kokkos::DualView<char*, buffer_device_type>& exports,
6435 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
6436 size_t& constantNumPackets)
const {
6441 packCrsMatrixNew(*
this, exports, numPacketsPerLID, exportLIDs,
6442 constantNumPackets);
6444 this->packNonStaticNew(exportLIDs, exports, numPacketsPerLID,
6445 constantNumPackets);
6449template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6451 packNonStaticNew(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6452 Kokkos::DualView<char*, buffer_device_type>& exports,
6453 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
6454 size_t& constantNumPackets)
const {
6461 using LO = LocalOrdinal;
6462 using GO = GlobalOrdinal;
6464 const char tfecfFuncName[] =
"packNonStaticNew: ";
6466 const bool verbose = Behavior::verbose(
"CrsMatrix");
6467 std::unique_ptr<std::string> prefix;
6469 prefix = this->createPrefix(
"CrsMatrix",
"packNonStaticNew");
6470 std::ostringstream os;
6471 os << *prefix <<
"Start" << endl;
6472 std::cerr << os.str();
6475 const size_t numExportLIDs =
static_cast<size_t>(exportLIDs.extent(0));
6476 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numExportLIDs !=
static_cast<size_t>(numPacketsPerLID.extent(0)),
6477 std::invalid_argument,
"exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size() = " << numPacketsPerLID.extent(0) <<
".");
6482 constantNumPackets = 0;
6487 size_t totalNumEntries = 0;
6488 this->allocatePackSpaceNew(exports, totalNumEntries, exportLIDs);
6489 const size_t bufSize =
static_cast<size_t>(exports.extent(0));
6492 exports.clear_sync_state();
6493 exports.modify_host();
6494 auto exports_h = exports.view_host();
6496 std::ostringstream os;
6497 os << *prefix <<
"After marking exports as modified on host, "
6498 << dualViewStatusToString(exports,
"exports") << endl;
6499 std::cerr << os.str();
6503 auto exportLIDs_h = exportLIDs.view_host();
6506 const_cast<Kokkos::DualView<size_t*, buffer_device_type>*
>(&numPacketsPerLID)->clear_sync_state();
6507 const_cast<Kokkos::DualView<size_t*, buffer_device_type>*
>(&numPacketsPerLID)->modify_host();
6508 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
6513 auto maxRowNumEnt = this->getLocalMaxNumRowEntries();
6516 typename global_inds_host_view_type::non_const_type gidsIn_k;
6517 if (this->isLocallyIndexed()) {
6519 typename global_inds_host_view_type::non_const_type(
"packGids",
6524 for (
size_t i = 0; i < numExportLIDs; ++i) {
6525 const LO lclRow = exportLIDs_h[i];
6527 size_t numBytes = 0;
6528 size_t numEnt = this->getNumEntriesInLocalRow(lclRow);
6535 numPacketsPerLID_h[i] = 0;
6539 if (this->isLocallyIndexed()) {
6540 typename global_inds_host_view_type::non_const_type gidsIn;
6541 values_host_view_type valsIn;
6545 local_inds_host_view_type lidsIn;
6546 this->getLocalRowView(lclRow, lidsIn, valsIn);
6547 const map_type& colMap = *(this->getColMap());
6548 for (
size_t k = 0; k < numEnt; ++k) {
6549 gidsIn_k[k] = colMap.getGlobalElement(lidsIn[k]);
6551 gidsIn = Kokkos::subview(gidsIn_k, Kokkos::make_pair(GO(0), GO(numEnt)));
6553 const size_t numBytesPerValue =
6554 PackTraits<ST>::packValueCount(valsIn[0]);
6555 numBytes = this->
packRow(exports_h.data(), offset, numEnt,
6556 gidsIn.data(), valsIn.data(),
6558 }
else if (this->isGloballyIndexed()) {
6559 global_inds_host_view_type gidsIn;
6560 values_host_view_type valsIn;
6566 const map_type& rowMap = *(this->getRowMap());
6568 this->getGlobalRowView(gblRow, gidsIn, valsIn);
6570 const size_t numBytesPerValue =
6571 PackTraits<ST>::packValueCount(valsIn[0]);
6572 numBytes = this->
packRow(exports_h.data(), offset, numEnt,
6573 gidsIn.data(), valsIn.data(),
6580 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(offset > bufSize || offset + numBytes > bufSize, std::logic_error,
6581 "First invalid offset into 'exports' pack buffer at index i = " << i
6582 <<
". exportLIDs_h[i]: " << exportLIDs_h[i] <<
", bufSize: " << bufSize <<
", offset: " << offset <<
", numBytes: " << numBytes <<
".");
6586 numPacketsPerLID_h[i] = numBytes;
6591 std::ostringstream os;
6592 os << *prefix <<
"Tpetra::CrsMatrix::packNonStaticNew: After:" << endl
6599 std::cerr << os.str();
6603template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6607 const LocalOrdinal numEnt,
6609 const GlobalOrdinal cols[],
6611 const char*
const prefix,
6613 const bool verbose) {
6614 using GO = GlobalOrdinal;
6618 const GO gblRow = myGraph_->rowMap_->getGlobalElement(lclRow);
6619 Teuchos::ArrayView<const GO> cols_av(numEnt == 0 ?
nullptr : cols, numEnt);
6620 Teuchos::ArrayView<const Scalar> vals_av(numEnt == 0 ?
nullptr :
reinterpret_cast<const Scalar*
>(vals), numEnt);
6625 combineGlobalValues(gblRow, cols_av, vals_av, combMode,
6626 prefix, debug, verbose);
6630template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6633 const GlobalOrdinal globalRowIndex,
6634 const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
6635 const Teuchos::ArrayView<const Scalar>& values,
6637 const char*
const prefix,
6639 const bool verbose) {
6640 const char tfecfFuncName[] =
"combineGlobalValues: ";
6642 if (isStaticGraph()) {
6646 if (combineMode ==
ADD) {
6647 sumIntoGlobalValues(globalRowIndex, columnIndices, values);
6648 }
else if (combineMode ==
REPLACE) {
6649 replaceGlobalValues(globalRowIndex, columnIndices, values);
6650 }
else if (combineMode ==
ABSMAX) {
6651 using ::Tpetra::Details::AbsMax;
6653 this->
template transformGlobalValues<AbsMax<Scalar>>(globalRowIndex,
6656 }
else if (combineMode ==
INSERT) {
6657 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isStaticGraph() && combineMode ==
INSERT,
6658 std::invalid_argument,
6659 "INSERT combine mode is forbidden "
6660 "if the matrix has a static (const) graph (i.e., was "
6661 "constructed with the CrsMatrix constructor that takes a "
6662 "const CrsGraph pointer).");
6664 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error,
6665 "Invalid combine mode; should "
6667 "Please report this bug to the Tpetra developers.");
6670 if (combineMode ==
ADD || combineMode ==
INSERT) {
6677 insertGlobalValuesFilteredChecked(globalRowIndex,
6678 columnIndices, values, prefix, debug, verbose);
6689 else if (combineMode ==
ABSMAX) {
6690 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6691 !isStaticGraph() && combineMode ==
ABSMAX, std::logic_error,
6692 "ABSMAX combine mode when the matrix has a dynamic graph is not yet "
6694 }
else if (combineMode ==
REPLACE) {
6695 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6696 !isStaticGraph() && combineMode ==
REPLACE, std::logic_error,
6697 "REPLACE combine mode when the matrix has a dynamic graph is not yet "
6700 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6701 true, std::logic_error,
6702 "Should never get here! Please report this "
6703 "bug to the Tpetra developers.");
6708template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6710 unpackAndCombine(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& importLIDs,
6711 Kokkos::DualView<char*, buffer_device_type> imports,
6712 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6713 const size_t constantNumPackets,
6719 const char tfecfFuncName[] =
"unpackAndCombine: ";
6720 ProfilingRegion regionUAC(
"Tpetra::CrsMatrix::unpackAndCombine");
6722 const bool debug = Behavior::debug(
"CrsMatrix");
6723 const bool verbose = Behavior::verbose(
"CrsMatrix");
6724 constexpr int numValidModes = 5;
6727 const char* validModeNames[numValidModes] =
6728 {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT",
"ZERO"};
6730 std::unique_ptr<std::string> prefix;
6732 prefix = this->createPrefix(
"CrsMatrix",
"unpackAndCombine");
6733 std::ostringstream os;
6734 os << *prefix <<
"Start:" << endl
6736 << dualViewStatusToString(importLIDs,
"importLIDs")
6739 << dualViewStatusToString(imports,
"imports")
6742 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6744 << *prefix <<
" constantNumPackets: " << constantNumPackets
6748 std::cerr << os.str();
6752 if (std::find(validModes, validModes + numValidModes, combineMode) ==
6753 validModes + numValidModes) {
6754 std::ostringstream os;
6755 os <<
"Invalid combine mode. Valid modes are {";
6756 for (
int k = 0; k < numValidModes; ++k) {
6757 os << validModeNames[k];
6758 if (k < numValidModes - 1) {
6763 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::invalid_argument, os.str());
6765 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(importLIDs.extent(0) != numPacketsPerLID.extent(0),
6766 std::invalid_argument,
"importLIDs.extent(0)=" << importLIDs.extent(0) <<
" != numPacketsPerLID.extent(0)=" << numPacketsPerLID.extent(0) <<
".");
6769 if (combineMode ==
ZERO) {
6774 using Teuchos::reduceAll;
6775 std::unique_ptr<std::ostringstream> msg(
new std::ostringstream());
6778 unpackAndCombineImpl(importLIDs, imports, numPacketsPerLID,
6779 constantNumPackets, combineMode,
6781 }
catch (std::exception& e) {
6786 const Teuchos::Comm<int>& comm = *(this->
getComm());
6787 reduceAll<int, int>(comm, Teuchos::REDUCE_MAX,
6788 lclBad, Teuchos::outArg(gblBad));
6794 std::ostringstream os;
6795 os <<
"Proc " << comm.getRank() <<
": " << msg->str() << endl;
6796 msg = std::unique_ptr<std::ostringstream>(
new std::ostringstream());
6798 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error, std::endl
6799 <<
"unpackAndCombineImpl "
6800 "threw an exception on one or more participating processes: "
6805 unpackAndCombineImpl(importLIDs, imports, numPacketsPerLID,
6806 constantNumPackets, combineMode,
6811 std::ostringstream os;
6812 os << *prefix <<
"Done!" << endl
6814 << dualViewStatusToString(importLIDs,
"importLIDs")
6817 << dualViewStatusToString(imports,
"imports")
6820 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6822 std::cerr << os.str();
6826template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6831 Kokkos::DualView<char*, buffer_device_type> imports,
6832 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6833 const size_t constantNumPackets,
6835 const bool verbose) {
6837 "Tpetra::CrsMatrix::unpackAndCombineImpl",
6840 const char tfecfFuncName[] =
"unpackAndCombineImpl";
6841 std::unique_ptr<std::string> prefix;
6843 prefix = this->createPrefix(
"CrsMatrix", tfecfFuncName);
6844 std::ostringstream os;
6845 os << *prefix <<
"isStaticGraph(): "
6846 << (isStaticGraph() ?
"true" :
"false")
6847 <<
", importLIDs.extent(0): "
6848 << importLIDs.extent(0)
6849 <<
", imports.extent(0): "
6850 << imports.extent(0)
6851 <<
", numPacketsPerLID.extent(0): "
6852 << numPacketsPerLID.extent(0)
6854 std::cerr << os.str();
6857 if (isStaticGraph()) {
6858 using Details::unpackCrsMatrixAndCombineNew;
6859 unpackCrsMatrixAndCombineNew(*
this, imports, numPacketsPerLID,
6860 importLIDs, constantNumPackets,
6864 using padding_type =
typename crs_graph_type::padding_type;
6865 std::unique_ptr<padding_type> padding;
6867 padding = myGraph_->computePaddingForCrsMatrixUnpack(
6868 importLIDs, imports, numPacketsPerLID, verbose);
6869 }
catch (std::exception& e) {
6870 const auto rowMap = getRowMap();
6871 const auto comm = rowMap.is_null() ? Teuchos::null : rowMap->
getComm();
6872 const int myRank = comm.is_null() ? -1 : comm->getRank();
6873 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"Proc " << myRank <<
": "
6874 "Tpetra::CrsGraph::computePaddingForCrsMatrixUnpack "
6875 "threw an exception: "
6879 std::ostringstream os;
6880 os << *prefix <<
"Call applyCrsPadding" << endl;
6881 std::cerr << os.str();
6883 applyCrsPadding(*padding, verbose);
6886 std::ostringstream os;
6887 os << *prefix <<
"Call unpackAndCombineImplNonStatic" << endl;
6888 std::cerr << os.str();
6890 unpackAndCombineImplNonStatic(importLIDs, imports,
6897 std::ostringstream os;
6898 os << *prefix <<
"Done" << endl;
6899 std::cerr << os.str();
6903template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6908 Kokkos::DualView<char*, buffer_device_type> imports,
6909 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6910 const size_t constantNumPackets,
6917 using Kokkos::MemoryUnmanaged;
6918 using Kokkos::subview;
6921 using LO = LocalOrdinal;
6922 using GO = GlobalOrdinal;
6924 using size_type =
typename Teuchos::ArrayView<LO>::size_type;
6926 typename View<int*, device_type>::host_mirror_type::execution_space;
6927 using pair_type = std::pair<typename View<int*, HES>::size_type,
6928 typename View<int*, HES>::size_type>;
6929 using gids_out_type = View<GO*, HES, MemoryUnmanaged>;
6930 using vals_out_type = View<ST*, HES, MemoryUnmanaged>;
6931 const char tfecfFuncName[] =
"unpackAndCombineImplNonStatic";
6933 const bool debug = Behavior::debug(
"CrsMatrix");
6934 const bool verbose = Behavior::verbose(
"CrsMatrix");
6935 std::unique_ptr<std::string> prefix;
6937 prefix = this->
createPrefix(
"CrsMatrix", tfecfFuncName);
6938 std::ostringstream os;
6939 os << *prefix << endl;
6940 std::cerr << os.str();
6942 const char*
const prefix_raw =
6943 verbose ? prefix.get()->c_str() :
nullptr;
6945 const size_type numImportLIDs = importLIDs.extent(0);
6946 if (combineMode ==
ZERO || numImportLIDs == 0) {
6951 "Tpetra::CrsMatrix::unpackAndCombineImplNonStatic",
6955 if (imports.need_sync_host()) {
6956 imports.sync_host();
6958 auto imports_h = imports.view_host();
6961 if (numPacketsPerLID.need_sync_host()) {
6962 numPacketsPerLID.sync_host();
6964 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
6966 TEUCHOS_ASSERT(!importLIDs.need_sync_host());
6967 auto importLIDs_h = importLIDs.view_host();
6969 size_t numBytesPerValue;
6980 numBytesPerValue = PackTraits<ST>::packValueCount(val);
6985 size_t maxRowNumEnt = 0;
6986 for (size_type i = 0; i < numImportLIDs; ++i) {
6987 const size_t numBytes = numPacketsPerLID_h[i];
6988 if (numBytes == 0) {
6993 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(offset + numBytes >
size_t(imports_h.extent(0)),
6994 std::logic_error,
": At local row index importLIDs_h[i=" << i <<
"]=" << importLIDs_h[i] <<
", offset (=" << offset <<
") + numBytes (=" << numBytes <<
") > "
6995 "imports_h.extent(0)="
6996 << imports_h.extent(0) <<
".");
7001 const size_t theNumBytes =
7002 PackTraits<LO>::packValueCount(numEntLO);
7003 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(theNumBytes > numBytes, std::logic_error,
": theNumBytes=" << theNumBytes <<
" > numBytes = " << numBytes <<
".");
7005 const char*
const inBuf = imports_h.data() + offset;
7006 const size_t actualNumBytes =
7007 PackTraits<LO>::unpackValue(numEntLO, inBuf);
7010 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(actualNumBytes > numBytes, std::logic_error,
": At i=" << i <<
", actualNumBytes=" << actualNumBytes <<
" > numBytes=" << numBytes <<
".");
7011 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numEntLO == 0, std::logic_error,
7012 ": At local row index "
7014 << i <<
"]=" << importLIDs_h[i] <<
", "
7015 "the number of entries read from the packed data is "
7017 << numEntLO <<
", but numBytes=" << numBytes
7021 maxRowNumEnt = std::max(
size_t(numEntLO), maxRowNumEnt);
7029 View<GO*, HES> gblColInds;
7030 View<LO*, HES> lclColInds;
7031 View<ST*, HES> vals;
7044 gblColInds = ScalarViewTraits<GO, HES>::allocateArray(
7045 gid, maxRowNumEnt,
"gids");
7046 lclColInds = ScalarViewTraits<LO, HES>::allocateArray(
7047 lid, maxRowNumEnt,
"lids");
7048 vals = ScalarViewTraits<ST, HES>::allocateArray(
7049 val, maxRowNumEnt,
"vals");
7053 for (size_type i = 0; i < numImportLIDs; ++i) {
7054 const size_t numBytes = numPacketsPerLID_h[i];
7055 if (numBytes == 0) {
7059 const char*
const inBuf = imports_h.data() + offset;
7060 (void)PackTraits<LO>::unpackValue(numEntLO, inBuf);
7062 const size_t numEnt =
static_cast<size_t>(numEntLO);
7064 const LO lclRow = importLIDs_h[i];
7066 gids_out_type gidsOut = subview(gblColInds, pair_type(0, numEnt));
7067 vals_out_type valsOut = subview(vals, pair_type(0, numEnt));
7069 const size_t numBytesOut =
7070 unpackRow(gidsOut.data(), valsOut.data(), imports_h.data(),
7071 offset, numBytes, numEnt, numBytesPerValue);
7072 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numBytes != numBytesOut, std::logic_error,
": At i=" << i <<
", numBytes=" << numBytes <<
" != numBytesOut=" << numBytesOut <<
".");
7074 const ST*
const valsRaw =
const_cast<const ST*
>(valsOut.data());
7075 const GO*
const gidsRaw =
const_cast<const GO*
>(gidsOut.data());
7076 combineGlobalValuesRaw(lclRow, numEnt, valsRaw, gidsRaw,
7077 combineMode, prefix_raw, debug, verbose);
7083 std::ostringstream os;
7084 os << *prefix <<
"Done" << endl;
7085 std::cerr << os.str();
7089template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7090Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7093 const bool force)
const {
7094 using Teuchos::null;
7098 TEUCHOS_TEST_FOR_EXCEPTION(
7100 "Tpetra::CrsMatrix::getColumn"
7101 "MapMultiVector: You may only call this method if the matrix has a "
7102 "column Map. If the matrix does not yet have a column Map, you should "
7103 "first call fillComplete (with domain and range Map if necessary).");
7107 TEUCHOS_TEST_FOR_EXCEPTION(
7110 "CrsMatrix::getColumnMapMultiVector: You may only call this method if "
7111 "this matrix's graph is fill complete.");
7114 RCP<const import_type> importer = this->
getGraph()->getImporter();
7115 RCP<const map_type> colMap = this->
getColMap();
7128 if (!importer.is_null() || force) {
7130 X_colMap = rcp(
new MV(colMap, numVecs));
7146template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7147Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7150 const bool force)
const {
7151 using Teuchos::null;
7157 TEUCHOS_TEST_FOR_EXCEPTION(
7160 "CrsMatrix::getRowMapMultiVector: You may only call this method if this "
7161 "matrix's graph is fill complete.");
7164 RCP<const export_type> exporter = this->
getGraph()->getExporter();
7168 RCP<const map_type> rowMap = this->
getRowMap();
7180 if (!exporter.is_null() || force) {
7182 Y_rowMap = rcp(
new MV(rowMap, numVecs));
7191template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7194 TEUCHOS_TEST_FOR_EXCEPTION(
7195 myGraph_.is_null(), std::logic_error,
7196 "Tpetra::CrsMatrix::"
7197 "removeEmptyProcessesInPlace: This method does not work when the matrix "
7198 "was created with a constant graph (that is, when it was created using "
7199 "the version of its constructor that takes an RCP<const CrsGraph>). "
7200 "This is because the matrix is not allowed to modify the graph in that "
7201 "case, but removing empty processes requires modifying the graph.");
7202 myGraph_->removeEmptyProcessesInPlace(newMap);
7210 staticGraph_ = Teuchos::rcp_const_cast<const Graph>(myGraph_);
7213template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7214Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7216 add(
const Scalar& alpha,
7219 const Teuchos::RCP<const map_type>& domainMap,
7220 const Teuchos::RCP<const map_type>& rangeMap,
7221 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
7223 using Teuchos::Array;
7224 using Teuchos::ArrayView;
7225 using Teuchos::ParameterList;
7228 using Teuchos::rcp_implicit_cast;
7229 using Teuchos::sublist;
7232 using crs_matrix_type =
7233 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
7234 const char errPfx[] =
"Tpetra::CrsMatrix::add: ";
7238 std::unique_ptr<std::string> prefix;
7240 prefix = this->createPrefix(
"CrsMatrix",
"add");
7241 std::ostringstream os;
7242 os << *prefix <<
"Start" << endl;
7243 std::cerr << os.str();
7246 const crs_matrix_type& B = *
this;
7247 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
7248 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
7256 RCP<const map_type> B_domainMap = B.getDomainMap();
7257 RCP<const map_type> B_rangeMap = B.getRangeMap();
7259 RCP<const map_type> theDomainMap = domainMap;
7260 RCP<const map_type> theRangeMap = rangeMap;
7262 if (domainMap.is_null()) {
7263 if (B_domainMap.is_null()) {
7264 TEUCHOS_TEST_FOR_EXCEPTION(
7265 A_domainMap.is_null(), std::invalid_argument,
7266 "Tpetra::CrsMatrix::add: If neither A nor B have a domain Map, "
7267 "then you must supply a nonnull domain Map to this method.");
7268 theDomainMap = A_domainMap;
7270 theDomainMap = B_domainMap;
7273 if (rangeMap.is_null()) {
7274 if (B_rangeMap.is_null()) {
7275 TEUCHOS_TEST_FOR_EXCEPTION(
7276 A_rangeMap.is_null(), std::invalid_argument,
7277 "Tpetra::CrsMatrix::add: If neither A nor B have a range Map, "
7278 "then you must supply a nonnull range Map to this method.");
7279 theRangeMap = A_rangeMap;
7281 theRangeMap = B_rangeMap;
7289 if (!A_domainMap.is_null() && !A_rangeMap.is_null()) {
7290 if (!B_domainMap.is_null() && !B_rangeMap.is_null()) {
7291 TEUCHOS_TEST_FOR_EXCEPTION(!B_domainMap->isSameAs(*A_domainMap),
7292 std::invalid_argument,
7293 errPfx <<
"The input RowMatrix A must have a domain Map "
7294 "which is the same as (isSameAs) this RowMatrix's "
7296 TEUCHOS_TEST_FOR_EXCEPTION(!B_rangeMap->isSameAs(*A_rangeMap), std::invalid_argument,
7297 errPfx <<
"The input RowMatrix A must have a range Map "
7298 "which is the same as (isSameAs) this RowMatrix's range "
7300 TEUCHOS_TEST_FOR_EXCEPTION(!domainMap.is_null() &&
7301 !domainMap->isSameAs(*B_domainMap),
7302 std::invalid_argument,
7303 errPfx <<
"The input domain Map must be the same as "
7304 "(isSameAs) this RowMatrix's domain Map.");
7305 TEUCHOS_TEST_FOR_EXCEPTION(!rangeMap.is_null() &&
7306 !rangeMap->isSameAs(*B_rangeMap),
7307 std::invalid_argument,
7308 errPfx <<
"The input range Map must be the same as "
7309 "(isSameAs) this RowMatrix's range Map.");
7311 }
else if (!B_domainMap.is_null() && !B_rangeMap.is_null()) {
7312 TEUCHOS_TEST_FOR_EXCEPTION(!domainMap.is_null() &&
7313 !domainMap->isSameAs(*B_domainMap),
7314 std::invalid_argument,
7315 errPfx <<
"The input domain Map must be the same as "
7316 "(isSameAs) this RowMatrix's domain Map.");
7317 TEUCHOS_TEST_FOR_EXCEPTION(!rangeMap.is_null() && !rangeMap->isSameAs(*B_rangeMap),
7318 std::invalid_argument,
7319 errPfx <<
"The input range Map must be the same as "
7320 "(isSameAs) this RowMatrix's range Map.");
7322 TEUCHOS_TEST_FOR_EXCEPTION(domainMap.is_null() || rangeMap.is_null(),
7323 std::invalid_argument, errPfx <<
"If neither A nor B "
7324 "have a domain and range Map, then you must supply a "
7325 "nonnull domain and range Map to this method.");
7332 bool callFillComplete =
true;
7333 RCP<ParameterList> constructorSublist;
7334 RCP<ParameterList> fillCompleteSublist;
7335 if (!params.is_null()) {
7337 params->get(
"Call fillComplete", callFillComplete);
7338 constructorSublist = sublist(params,
"Constructor parameters");
7339 fillCompleteSublist = sublist(params,
"fillComplete parameters");
7342 RCP<const map_type> A_rowMap = A.
getRowMap();
7343 RCP<const map_type> B_rowMap = B.getRowMap();
7344 RCP<const map_type> C_rowMap = B_rowMap;
7345 RCP<crs_matrix_type> C;
7351 if (A_rowMap->isSameAs(*B_rowMap)) {
7352 const LO localNumRows =
static_cast<LO
>(A_rowMap->getLocalNumElements());
7353 Array<size_t> C_maxNumEntriesPerRow(localNumRows, 0);
7356 if (alpha !=
ZERO) {
7357 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
7359 C_maxNumEntriesPerRow[localRow] += A_numEntries;
7364 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
7365 const size_t B_numEntries = B.getNumEntriesInLocalRow(localRow);
7366 C_maxNumEntriesPerRow[localRow] += B_numEntries;
7370 if (constructorSublist.is_null()) {
7371 C = rcp(
new crs_matrix_type(C_rowMap, C_maxNumEntriesPerRow()));
7373 C = rcp(
new crs_matrix_type(C_rowMap, C_maxNumEntriesPerRow(),
7374 constructorSublist));
7384 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, errPfx <<
"The row maps must "
7385 "be the same for statically allocated matrices, to ensure "
7386 "that there is sufficient space to do the addition.");
7389 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null(), std::logic_error,
7390 errPfx <<
"C should not be null at this point. "
7391 "Please report this bug to the Tpetra developers.");
7394 std::ostringstream os;
7395 os << *prefix <<
"Compute C = alpha*A + beta*B" << endl;
7396 std::cerr << os.str();
7398 using gids_type = nonconst_global_inds_host_view_type;
7399 using vals_type = nonconst_values_host_view_type;
7403 if (alpha !=
ZERO) {
7404 const LO A_localNumRows =
static_cast<LO
>(A_rowMap->getLocalNumElements());
7405 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
7407 const GO globalRow = A_rowMap->getGlobalElement(localRow);
7408 if (A_numEntries >
static_cast<size_t>(ind.size())) {
7409 Kokkos::resize(ind, A_numEntries);
7410 Kokkos::resize(val, A_numEntries);
7412 gids_type indView = Kokkos::subview(ind, std::make_pair((
size_t)0, A_numEntries));
7413 vals_type valView = Kokkos::subview(val, std::make_pair((
size_t)0, A_numEntries));
7417 for (
size_t k = 0; k < A_numEntries; ++k) {
7418 valView[k] *= alpha;
7421 C->insertGlobalValues(globalRow, A_numEntries,
7422 reinterpret_cast<Scalar*
>(valView.data()),
7428 const LO B_localNumRows =
static_cast<LO
>(B_rowMap->getLocalNumElements());
7429 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
7430 size_t B_numEntries = B.getNumEntriesInLocalRow(localRow);
7431 const GO globalRow = B_rowMap->getGlobalElement(localRow);
7432 if (B_numEntries >
static_cast<size_t>(ind.size())) {
7433 Kokkos::resize(ind, B_numEntries);
7434 Kokkos::resize(val, B_numEntries);
7436 gids_type indView = Kokkos::subview(ind, std::make_pair((
size_t)0, B_numEntries));
7437 vals_type valView = Kokkos::subview(val, std::make_pair((
size_t)0, B_numEntries));
7438 B.getGlobalRowCopy(globalRow, indView, valView, B_numEntries);
7441 for (
size_t k = 0; k < B_numEntries; ++k) {
7445 C->insertGlobalValues(globalRow, B_numEntries,
7446 reinterpret_cast<Scalar*
>(valView.data()),
7451 if (callFillComplete) {
7453 std::ostringstream os;
7454 os << *prefix <<
"Call fillComplete on C" << endl;
7455 std::cerr << os.str();
7457 if (fillCompleteSublist.is_null()) {
7458 C->fillComplete(theDomainMap, theRangeMap);
7460 C->fillComplete(theDomainMap, theRangeMap, fillCompleteSublist);
7462 }
else if (verbose) {
7463 std::ostringstream os;
7464 os << *prefix <<
"Do NOT call fillComplete on C" << endl;
7465 std::cerr << os.str();
7469 std::ostringstream os;
7470 os << *prefix <<
"Done" << endl;
7471 std::cerr << os.str();
7473 return rcp_implicit_cast<row_matrix_type>(C);
7476template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7479 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
7480 const Teuchos::RCP<const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>>& domainTransfer,
7481 const Teuchos::RCP<const map_type>& domainMap,
7482 const Teuchos::RCP<const map_type>& rangeMap,
7483 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
7490 using Teuchos::ArrayRCP;
7491 using Teuchos::ArrayView;
7492 using Teuchos::Comm;
7493 using Teuchos::ParameterList;
7495 typedef LocalOrdinal LO;
7496 typedef GlobalOrdinal GO;
7502 const bool debug = Behavior::debug(
"CrsMatrix");
7503 const bool verbose = Behavior::verbose(
"CrsMatrix");
7504 int MyPID = getComm()->getRank();
7506 std::unique_ptr<std::string> verbosePrefix;
7509 this->createPrefix(
"CrsMatrix",
"transferAndFillComplete");
7510 std::ostringstream os;
7511 os <<
"Start" << endl;
7512 std::cerr << os.str();
7519 bool reverseMode =
false;
7520 bool restrictComm =
false;
7522 int mm_optimization_core_count =
7523 Behavior::TAFC_OptimizationCoreCount();
7524 RCP<ParameterList> matrixparams;
7525 bool overrideAllreduce =
false;
7526 bool useKokkosPath =
false;
7527 if (!params.is_null()) {
7528 matrixparams = sublist(params,
"CrsMatrix");
7529 reverseMode = params->get(
"Reverse Mode", reverseMode);
7530 useKokkosPath = params->get(
"TAFC: use kokkos path", useKokkosPath);
7531 restrictComm = params->get(
"Restrict Communicator", restrictComm);
7532 auto& slist = params->sublist(
"matrixmatrix: kernel params",
false);
7533 isMM = slist.get(
"isMatrixMatrix_TransferAndFillComplete",
false);
7534 mm_optimization_core_count = slist.get(
"MM_TAFC_OptimizationCoreCount", mm_optimization_core_count);
7536 overrideAllreduce = slist.get(
"MM_TAFC_OverrideAllreduceCheck",
false);
7537 if (getComm()->getSize() < mm_optimization_core_count && isMM) isMM =
false;
7538 if (reverseMode) isMM =
false;
7542 std::shared_ptr<::Tpetra::Details::CommRequest> iallreduceRequest;
7544 int reduced_mismatch = 0;
7545 if (isMM && !overrideAllreduce) {
7547 const bool source_vals = !getGraph()->getImporter().is_null();
7548 const bool target_vals = !(rowTransfer.getExportLIDs().size() == 0 ||
7549 rowTransfer.getRemoteLIDs().size() == 0);
7550 mismatch = (source_vals != target_vals) ? 1 : 0;
7552 ::Tpetra::Details::iallreduce(mismatch, reduced_mismatch,
7553 Teuchos::REDUCE_MAX, *(getComm()));
7556#ifdef HAVE_TPETRA_MMM_TIMINGS
7557 using Teuchos::TimeMonitor;
7559 if (!params.is_null())
7560 label = params->get(
"Timer Label", label);
7561 std::string prefix = std::string(
"Tpetra ") + label + std::string(
": ");
7564 std::ostringstream os;
7572 Teuchos::TimeMonitor MMall(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC All") + tlstr));
7582 TEUCHOS_TEST_FOR_EXCEPTION(
7583 xferAsImport ==
nullptr && xferAsExport ==
nullptr, std::invalid_argument,
7584 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' input "
7585 "argument must be either an Import or an Export, and its template "
7586 "parameters must match the corresponding template parameters of the "
7594 Teuchos::RCP<const import_type> xferDomainAsImport = Teuchos::rcp_dynamic_cast<const import_type>(domainTransfer);
7595 Teuchos::RCP<const export_type> xferDomainAsExport = Teuchos::rcp_dynamic_cast<const export_type>(domainTransfer);
7597 if (!domainTransfer.is_null()) {
7598 TEUCHOS_TEST_FOR_EXCEPTION(
7599 (xferDomainAsImport.is_null() && xferDomainAsExport.is_null()), std::invalid_argument,
7600 "Tpetra::CrsMatrix::transferAndFillComplete: The 'domainTransfer' input "
7601 "argument must be either an Import or an Export, and its template "
7602 "parameters must match the corresponding template parameters of the "
7605 TEUCHOS_TEST_FOR_EXCEPTION(
7606 (xferAsImport !=
nullptr || !xferDomainAsImport.is_null()) &&
7607 ((xferAsImport !=
nullptr && xferDomainAsImport.is_null()) ||
7608 (xferAsImport ==
nullptr && !xferDomainAsImport.is_null())),
7609 std::invalid_argument,
7610 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
7611 "arguments must be of the same type (either Import or Export).");
7613 TEUCHOS_TEST_FOR_EXCEPTION(
7614 (xferAsExport !=
nullptr || !xferDomainAsExport.is_null()) &&
7615 ((xferAsExport !=
nullptr && xferDomainAsExport.is_null()) ||
7616 (xferAsExport ==
nullptr && !xferDomainAsExport.is_null())),
7617 std::invalid_argument,
7618 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
7619 "arguments must be of the same type (either Import or Export).");
7624 const bool communication_needed = rowTransfer.getSourceMap()->isDistributed();
7628 RCP<const map_type> MyRowMap = reverseMode ? rowTransfer.getSourceMap() : rowTransfer.getTargetMap();
7629 RCP<const map_type> MyColMap;
7630 RCP<const map_type> MyDomainMap = !domainMap.is_null() ? domainMap : getDomainMap();
7631 RCP<const map_type> MyRangeMap = !rangeMap.is_null() ? rangeMap : getRangeMap();
7632 RCP<const map_type> BaseRowMap = MyRowMap;
7633 RCP<const map_type> BaseDomainMap = MyDomainMap;
7641 if (!destMat.is_null()) {
7652 const bool NewFlag = !destMat->getGraph()->isLocallyIndexed() &&
7653 !destMat->getGraph()->isGloballyIndexed();
7654 TEUCHOS_TEST_FOR_EXCEPTION(
7655 !NewFlag, std::invalid_argument,
7656 "Tpetra::CrsMatrix::"
7657 "transferAndFillComplete: The input argument 'destMat' is only allowed "
7658 "to be nonnull, if its graph is empty (neither locally nor globally "
7667 TEUCHOS_TEST_FOR_EXCEPTION(
7668 !destMat->getRowMap()->isSameAs(*MyRowMap), std::invalid_argument,
7669 "Tpetra::CrsMatrix::transferAndFillComplete: The (row) Map of the "
7670 "input argument 'destMat' is not the same as the (row) Map specified "
7671 "by the input argument 'rowTransfer'.");
7672 TEUCHOS_TEST_FOR_EXCEPTION(
7673 !destMat->checkSizes(*
this), std::invalid_argument,
7674 "Tpetra::CrsMatrix::transferAndFillComplete: You provided a nonnull "
7675 "destination matrix, but checkSizes() indicates that it is not a legal "
7676 "legal target for redistribution from the source matrix (*this). This "
7677 "may mean that they do not have the same dimensions.");
7691 TEUCHOS_TEST_FOR_EXCEPTION(
7692 !(reverseMode || getRowMap()->isSameAs(*rowTransfer.getSourceMap())),
7693 std::invalid_argument,
7694 "Tpetra::CrsMatrix::transferAndFillComplete: "
7695 "rowTransfer->getSourceMap() must match this->getRowMap() in forward mode.");
7696 TEUCHOS_TEST_FOR_EXCEPTION(
7697 !(!reverseMode || getRowMap()->isSameAs(*rowTransfer.getTargetMap())),
7698 std::invalid_argument,
7699 "Tpetra::CrsMatrix::transferAndFillComplete: "
7700 "rowTransfer->getTargetMap() must match this->getRowMap() in reverse mode.");
7703 TEUCHOS_TEST_FOR_EXCEPTION(
7704 !xferDomainAsImport.is_null() && !xferDomainAsImport->getTargetMap()->isSameAs(*domainMap),
7705 std::invalid_argument,
7706 "Tpetra::CrsMatrix::transferAndFillComplete: The target map of the 'domainTransfer' input "
7707 "argument must be the same as the rebalanced domain map 'domainMap'");
7709 TEUCHOS_TEST_FOR_EXCEPTION(
7710 !xferDomainAsExport.is_null() && !xferDomainAsExport->getSourceMap()->isSameAs(*domainMap),
7711 std::invalid_argument,
7712 "Tpetra::CrsMatrix::transferAndFillComplete: The source map of the 'domainTransfer' input "
7713 "argument must be the same as the rebalanced domain map 'domainMap'");
7726 const size_t NumSameIDs = rowTransfer.getNumSameIDs();
7727 ArrayView<const LO> ExportLIDs = reverseMode ? rowTransfer.getRemoteLIDs() : rowTransfer.getExportLIDs();
7728 auto RemoteLIDs = reverseMode ? rowTransfer.getExportLIDs_dv() : rowTransfer.getRemoteLIDs_dv();
7729 auto PermuteToLIDs = reverseMode ? rowTransfer.getPermuteFromLIDs_dv() : rowTransfer.getPermuteToLIDs_dv();
7730 auto PermuteFromLIDs = reverseMode ? rowTransfer.getPermuteToLIDs_dv() : rowTransfer.getPermuteFromLIDs_dv();
7731 Distributor& Distor = rowTransfer.getDistributor();
7734 Teuchos::Array<int> SourcePids;
7737 RCP<const map_type> ReducedRowMap, ReducedColMap,
7738 ReducedDomainMap, ReducedRangeMap;
7739 RCP<const Comm<int>> ReducedComm;
7743 if (destMat.is_null()) {
7744 destMat = rcp(
new this_CRS_type(MyRowMap, 0, matrixparams));
7751#ifdef HAVE_TPETRA_MMM_TIMINGS
7752 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC restrictComm")));
7754 ReducedRowMap = MyRowMap->removeEmptyProcesses();
7755 ReducedComm = ReducedRowMap.is_null() ? Teuchos::null : ReducedRowMap->getComm();
7756 destMat->removeEmptyProcessesInPlace(ReducedRowMap);
7758 ReducedDomainMap = MyRowMap.getRawPtr() == MyDomainMap.getRawPtr() ? ReducedRowMap : MyDomainMap->replaceCommWithSubset(ReducedComm);
7759 ReducedRangeMap = MyRowMap.getRawPtr() == MyRangeMap.getRawPtr() ? ReducedRowMap : MyRangeMap->replaceCommWithSubset(ReducedComm);
7762 MyRowMap = ReducedRowMap;
7763 MyDomainMap = ReducedDomainMap;
7764 MyRangeMap = ReducedRangeMap;
7767 if (!ReducedComm.is_null()) {
7768 MyPID = ReducedComm->getRank();
7773 ReducedComm = MyRowMap->getComm();
7780 RCP<const import_type> MyImporter = getGraph()->getImporter();
7783 bool bSameDomainMap = BaseDomainMap->isSameAs(*getDomainMap());
7785 if (!restrictComm && !MyImporter.is_null() && bSameDomainMap) {
7786#ifdef HAVE_TPETRA_MMM_TIMINGS
7787 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs same map")));
7796 }
else if (restrictComm && !MyImporter.is_null() && bSameDomainMap) {
7799#ifdef HAVE_TPETRA_MMM_TIMINGS
7800 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs restricted comm")));
7802 IntVectorType SourceDomain_pids(getDomainMap(),
true);
7803 IntVectorType SourceCol_pids(getColMap());
7805 SourceDomain_pids.putScalar(MyPID);
7807 SourceCol_pids.doImport(SourceDomain_pids, *MyImporter,
INSERT);
7808 SourcePids.resize(getColMap()->getLocalNumElements());
7809 SourceCol_pids.get1dCopy(SourcePids());
7810 }
else if (MyImporter.is_null()) {
7812#ifdef HAVE_TPETRA_MMM_TIMINGS
7813 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs all local entries")));
7815 SourcePids.resize(getColMap()->getLocalNumElements());
7816 SourcePids.assign(getColMap()->getLocalNumElements(), MyPID);
7817 }
else if (!MyImporter.is_null() &&
7818 !domainTransfer.is_null()) {
7823#ifdef HAVE_TPETRA_MMM_TIMINGS
7824 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs rectangular case")));
7828 IntVectorType TargetDomain_pids(domainMap);
7829 TargetDomain_pids.putScalar(MyPID);
7832 IntVectorType SourceDomain_pids(getDomainMap());
7835 IntVectorType SourceCol_pids(getColMap());
7837 if (!reverseMode && !xferDomainAsImport.is_null()) {
7838 SourceDomain_pids.doExport(TargetDomain_pids, *xferDomainAsImport,
INSERT);
7839 }
else if (reverseMode && !xferDomainAsExport.is_null()) {
7840 SourceDomain_pids.doExport(TargetDomain_pids, *xferDomainAsExport,
INSERT);
7841 }
else if (!reverseMode && !xferDomainAsExport.is_null()) {
7842 SourceDomain_pids.doImport(TargetDomain_pids, *xferDomainAsExport,
INSERT);
7843 }
else if (reverseMode && !xferDomainAsImport.is_null()) {
7844 SourceDomain_pids.doImport(TargetDomain_pids, *xferDomainAsImport,
INSERT);
7846 TEUCHOS_TEST_FOR_EXCEPTION(
7847 true, std::logic_error,
7848 "Tpetra::CrsMatrix::"
7849 "transferAndFillComplete: Should never get here! "
7850 "Please report this bug to a Tpetra developer.");
7852 SourceCol_pids.doImport(SourceDomain_pids, *MyImporter,
INSERT);
7853 SourcePids.resize(getColMap()->getLocalNumElements());
7854 SourceCol_pids.get1dCopy(SourcePids());
7855 }
else if (!MyImporter.is_null() &&
7856 BaseDomainMap->isSameAs(*BaseRowMap) &&
7857 getDomainMap()->isSameAs(*getRowMap())) {
7859#ifdef HAVE_TPETRA_MMM_TIMINGS
7860 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs query import")));
7863 IntVectorType TargetRow_pids(domainMap);
7864 IntVectorType SourceRow_pids(getRowMap());
7865 IntVectorType SourceCol_pids(getColMap());
7867 TargetRow_pids.putScalar(MyPID);
7868 if (!reverseMode && xferAsImport !=
nullptr) {
7869 SourceRow_pids.doExport(TargetRow_pids, *xferAsImport,
INSERT);
7870 }
else if (reverseMode && xferAsExport !=
nullptr) {
7871 SourceRow_pids.doExport(TargetRow_pids, *xferAsExport,
INSERT);
7872 }
else if (!reverseMode && xferAsExport !=
nullptr) {
7873 SourceRow_pids.doImport(TargetRow_pids, *xferAsExport,
INSERT);
7874 }
else if (reverseMode && xferAsImport !=
nullptr) {
7875 SourceRow_pids.doImport(TargetRow_pids, *xferAsImport,
INSERT);
7877 TEUCHOS_TEST_FOR_EXCEPTION(
7878 true, std::logic_error,
7879 "Tpetra::CrsMatrix::"
7880 "transferAndFillComplete: Should never get here! "
7881 "Please report this bug to a Tpetra developer.");
7884 SourceCol_pids.doImport(SourceRow_pids, *MyImporter,
INSERT);
7885 SourcePids.resize(getColMap()->getLocalNumElements());
7886 SourceCol_pids.get1dCopy(SourcePids());
7888 TEUCHOS_TEST_FOR_EXCEPTION(
7889 true, std::invalid_argument,
7890 "Tpetra::CrsMatrix::"
7891 "transferAndFillComplete: This method only allows either domainMap == "
7892 "getDomainMap (), or (domainMap == rowTransfer.getTargetMap () and "
7893 "getDomainMap () == getRowMap ()).");
7897 size_t constantNumPackets = destMat->constantNumberOfPackets();
7899#ifdef HAVE_TPETRA_MMM_TIMINGS
7900 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC reallocate buffers")));
7902 if (constantNumPackets == 0) {
7903 destMat->reallocArraysForNumPacketsPerLid(ExportLIDs.size(),
7904 RemoteLIDs.view_host().size());
7910 const size_t rbufLen = RemoteLIDs.view_host().size() * constantNumPackets;
7911 destMat->reallocImportsIfNeeded(rbufLen,
false,
nullptr);
7917#ifdef HAVE_TPETRA_MMM_TIMINGS
7918 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC pack and prepare")));
7923 using Teuchos::outArg;
7924 using Teuchos::REDUCE_MAX;
7925 using Teuchos::reduceAll;
7926 RCP<const Teuchos::Comm<int>> comm = this->getComm();
7927 const int myRank = comm->getRank();
7929 std::ostringstream errStrm;
7933 Teuchos::ArrayView<size_t> numExportPacketsPerLID;
7936 destMat->numExportPacketsPerLID_.modify_host();
7937 numExportPacketsPerLID =
7939 }
catch (std::exception& e) {
7940 errStrm <<
"Proc " << myRank <<
": getArrayViewFromDualView threw: "
7941 << e.what() << std::endl;
7944 errStrm <<
"Proc " << myRank <<
": getArrayViewFromDualView threw "
7945 "an exception not a subclass of std::exception"
7950 if (!comm.is_null()) {
7951 reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
7955 TEUCHOS_TEST_FOR_EXCEPTION(
7956 true, std::runtime_error,
7957 "getArrayViewFromDualView threw an "
7958 "exception on at least one process.");
7962 std::ostringstream os;
7963 os << *verbosePrefix <<
"Calling packCrsMatrixWithOwningPIDs"
7965 std::cerr << os.str();
7970 numExportPacketsPerLID,
7973 constantNumPackets);
7974 }
catch (std::exception& e) {
7975 errStrm <<
"Proc " << myRank <<
": packCrsMatrixWithOwningPIDs threw: "
7976 << e.what() << std::endl;
7979 errStrm <<
"Proc " << myRank <<
": packCrsMatrixWithOwningPIDs threw "
7980 "an exception not a subclass of std::exception"
7986 std::ostringstream os;
7987 os << *verbosePrefix <<
"Done with packCrsMatrixWithOwningPIDs"
7989 std::cerr << os.str();
7992 if (!comm.is_null()) {
7993 reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
7997 TEUCHOS_TEST_FOR_EXCEPTION(
7998 true, std::runtime_error,
7999 "packCrsMatrixWithOwningPIDs threw an "
8000 "exception on at least one process.");
8004 destMat->numExportPacketsPerLID_.modify_host();
8005 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
8008 std::ostringstream os;
8009 os << *verbosePrefix <<
"Calling packCrsMatrixWithOwningPIDs"
8011 std::cerr << os.str();
8015 numExportPacketsPerLID,
8018 constantNumPackets);
8020 std::ostringstream os;
8021 os << *verbosePrefix <<
"Done with packCrsMatrixWithOwningPIDs"
8023 std::cerr << os.str();
8030#ifdef HAVE_TPETRA_MMM_TIMINGS
8031 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC getOwningPIDs exchange remote data")));
8033 if (!communication_needed) {
8035 std::ostringstream os;
8036 os << *verbosePrefix <<
"Communication not needed" << std::endl;
8037 std::cerr << os.str();
8041 if (constantNumPackets == 0) {
8043 std::ostringstream os;
8044 os << *verbosePrefix <<
"Reverse mode, variable # packets / LID"
8046 std::cerr << os.str();
8051 destMat->numExportPacketsPerLID_.sync_host();
8052 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8054 destMat->numImportPacketsPerLID_.sync_host();
8055 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8059 std::ostringstream os;
8060 os << *verbosePrefix <<
"Calling 3-arg doReversePostsAndWaits"
8062 std::cerr << os.str();
8064 Distor.doReversePostsAndWaits(destMat->numExportPacketsPerLID_.view_host(), 1,
8065 destMat->numImportPacketsPerLID_.view_host());
8067 std::ostringstream os;
8068 os << *verbosePrefix <<
"Finished 3-arg doReversePostsAndWaits"
8070 std::cerr << os.str();
8073 size_t totalImportPackets = 0;
8075 totalImportPackets += numImportPacketsPerLID[i];
8080 destMat->reallocImportsIfNeeded(totalImportPackets, verbose,
8081 verbosePrefix.get());
8082 destMat->imports_.modify_host();
8083 auto hostImports = destMat->imports_.view_host();
8086 destMat->exports_.sync_host();
8087 auto hostExports = destMat->exports_.view_host();
8089 std::ostringstream os;
8090 os << *verbosePrefix <<
"Calling 4-arg doReversePostsAndWaits"
8092 std::cerr << os.str();
8094 Distor.doReversePostsAndWaits(hostExports,
8095 numExportPacketsPerLID,
8097 numImportPacketsPerLID);
8099 std::ostringstream os;
8100 os << *verbosePrefix <<
"Finished 4-arg doReversePostsAndWaits"
8102 std::cerr << os.str();
8106 std::ostringstream os;
8107 os << *verbosePrefix <<
"Reverse mode, constant # packets / LID"
8109 std::cerr << os.str();
8111 destMat->imports_.modify_host();
8112 auto hostImports = destMat->imports_.view_host();
8115 destMat->exports_.sync_host();
8116 auto hostExports = destMat->exports_.view_host();
8118 std::ostringstream os;
8119 os << *verbosePrefix <<
"Calling 3-arg doReversePostsAndWaits"
8121 std::cerr << os.str();
8123 Distor.doReversePostsAndWaits(hostExports,
8127 std::ostringstream os;
8128 os << *verbosePrefix <<
"Finished 3-arg doReversePostsAndWaits"
8130 std::cerr << os.str();
8134 if (constantNumPackets == 0) {
8136 std::ostringstream os;
8137 os << *verbosePrefix <<
"Forward mode, variable # packets / LID"
8139 std::cerr << os.str();
8144 destMat->numExportPacketsPerLID_.sync_host();
8145 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8147 destMat->numImportPacketsPerLID_.sync_host();
8148 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8151 std::ostringstream os;
8152 os << *verbosePrefix <<
"Calling 3-arg doPostsAndWaits"
8154 std::cerr << os.str();
8156 Distor.doPostsAndWaits(destMat->numExportPacketsPerLID_.view_host(), 1,
8157 destMat->numImportPacketsPerLID_.view_host());
8159 std::ostringstream os;
8160 os << *verbosePrefix <<
"Finished 3-arg doPostsAndWaits"
8162 std::cerr << os.str();
8165 size_t totalImportPackets = 0;
8167 totalImportPackets += numImportPacketsPerLID[i];
8172 destMat->reallocImportsIfNeeded(totalImportPackets, verbose,
8173 verbosePrefix.get());
8174 destMat->imports_.modify_host();
8175 auto hostImports = destMat->imports_.view_host();
8178 destMat->exports_.sync_host();
8179 auto hostExports = destMat->exports_.view_host();
8181 std::ostringstream os;
8182 os << *verbosePrefix <<
"Calling 4-arg doPostsAndWaits"
8184 std::cerr << os.str();
8186 Distor.doPostsAndWaits(hostExports,
8187 numExportPacketsPerLID,
8189 numImportPacketsPerLID);
8191 std::ostringstream os;
8192 os << *verbosePrefix <<
"Finished 4-arg doPostsAndWaits"
8194 std::cerr << os.str();
8198 std::ostringstream os;
8199 os << *verbosePrefix <<
"Forward mode, constant # packets / LID"
8201 std::cerr << os.str();
8203 destMat->imports_.modify_host();
8204 auto hostImports = destMat->imports_.view_host();
8207 destMat->exports_.sync_host();
8208 auto hostExports = destMat->exports_.view_host();
8210 std::ostringstream os;
8211 os << *verbosePrefix <<
"Calling 3-arg doPostsAndWaits"
8213 std::cerr << os.str();
8215 Distor.doPostsAndWaits(hostExports,
8219 std::ostringstream os;
8220 os << *verbosePrefix <<
"Finished 3-arg doPostsAndWaits"
8222 std::cerr << os.str();
8233 bool runOnHost = std::is_same_v<typename device_type::memory_space, Kokkos::HostSpace> && !useKokkosPath;
8235 Teuchos::Array<int> RemotePids;
8237 Teuchos::Array<int> TargetPids;
8243 destMat->numImportPacketsPerLID_.modify_host();
8245#ifdef HAVE_TPETRA_MMM_TIMINGS
8246 RCP<TimeMonitor> tmCopySPRdata = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC unpack-count-resize + copy same-perm-remote data"))));
8248 ArrayRCP<size_t> CSR_rowptr;
8249 ArrayRCP<GO> CSR_colind_GID;
8250 ArrayRCP<LO> CSR_colind_LID;
8251 ArrayRCP<Scalar> CSR_vals;
8253 destMat->imports_.sync_device();
8254 destMat->numImportPacketsPerLID_.sync_device();
8256 size_t N = BaseRowMap->getLocalNumElements();
8258 auto RemoteLIDs_d = RemoteLIDs.view_device();
8259 auto PermuteToLIDs_d = PermuteToLIDs.view_device();
8260 auto PermuteFromLIDs_d = PermuteFromLIDs.view_device();
8265 destMat->imports_.view_device(),
8266 destMat->numImportPacketsPerLID_.view_device(),
8280 if (
typeid(LO) ==
typeid(GO)) {
8281 CSR_colind_LID = Teuchos::arcp_reinterpret_cast<LO>(CSR_colind_GID);
8283 CSR_colind_LID.resize(CSR_colind_GID.size());
8285 CSR_colind_LID.resize(CSR_colind_GID.size());
8290 for (
size_t i = 0; i < static_cast<size_t>(TargetPids.size()); i++) {
8291 if (TargetPids[i] == -1) TargetPids[i] = MyPID;
8293#ifdef HAVE_TPETRA_MMM_TIMINGS
8294 tmCopySPRdata = Teuchos::null;
8303 std::ostringstream os;
8304 os << *verbosePrefix <<
"Calling lowCommunicationMakeColMapAndReindex"
8306 std::cerr << os.str();
8309#ifdef HAVE_TPETRA_MMM_TIMINGS
8310 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC makeColMap")));
8312 Import_Util::lowCommunicationMakeColMapAndReindexSerial(CSR_rowptr(),
8322 std::ostringstream os;
8323 os << *verbosePrefix <<
"restrictComm="
8324 << (restrictComm ?
"true" :
"false") << std::endl;
8325 std::cerr << os.str();
8332#ifdef HAVE_TPETRA_MMM_TIMINGS
8333 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC restrict colmap")));
8336 ReducedColMap = (MyRowMap.getRawPtr() == MyColMap.getRawPtr()) ? ReducedRowMap : MyColMap->replaceCommWithSubset(ReducedComm);
8337 MyColMap = ReducedColMap;
8342 std::ostringstream os;
8343 os << *verbosePrefix <<
"Calling replaceColMap" << std::endl;
8344 std::cerr << os.str();
8346 destMat->replaceColMap(MyColMap);
8353 if (ReducedComm.is_null()) {
8355 std::ostringstream os;
8356 os << *verbosePrefix <<
"I am no longer in the communicator; "
8359 std::cerr << os.str();
8368 if ((!reverseMode && xferAsImport !=
nullptr) ||
8369 (reverseMode && xferAsExport !=
nullptr)) {
8371 std::ostringstream os;
8372 os << *verbosePrefix <<
"Calling sortCrsEntries" << endl;
8373 std::cerr << os.str();
8375#ifdef HAVE_TPETRA_MMM_TIMINGS
8376 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC sortCrsEntries")));
8381 }
else if ((!reverseMode && xferAsExport !=
nullptr) ||
8382 (reverseMode && xferAsImport !=
nullptr)) {
8384 std::ostringstream os;
8385 os << *verbosePrefix <<
"Calling sortAndMergeCrsEntries"
8387 std::cerr << os.str();
8389#ifdef HAVE_TPETRA_MMM_TIMINGS
8390 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC sortAndMergeCrsEntries")));
8395 if (CSR_rowptr[N] !=
static_cast<size_t>(CSR_vals.size())) {
8396 CSR_colind_LID.resize(CSR_rowptr[N]);
8397 CSR_vals.resize(CSR_rowptr[N]);
8400 TEUCHOS_TEST_FOR_EXCEPTION(
8401 true, std::logic_error,
8402 "Tpetra::CrsMatrix::"
8403 "transferAndFillComplete: Should never get here! "
8404 "Please report this bug to a Tpetra developer.");
8411 std::ostringstream os;
8412 os << *verbosePrefix <<
"Calling destMat->setAllValues" << endl;
8413 std::cerr << os.str();
8422#ifdef HAVE_TPETRA_MMM_TIMINGS
8423 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC setAllValues")));
8425 destMat->setAllValues(CSR_rowptr, CSR_colind_LID, CSR_vals);
8436 destMat->numImportPacketsPerLID_.modify_host();
8438#ifdef HAVE_TPETRA_MMM_TIMINGS
8439 RCP<TimeMonitor> tmCopySPRdata = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC unpack-count-resize + copy same-perm-remote data"))));
8441 ArrayRCP<size_t> CSR_rowptr;
8442 ArrayRCP<GO> CSR_colind_GID;
8443 ArrayRCP<LO> CSR_colind_LID;
8444 ArrayRCP<Scalar> CSR_vals;
8446 destMat->imports_.sync_device();
8447 destMat->numImportPacketsPerLID_.sync_device();
8449 size_t N = BaseRowMap->getLocalNumElements();
8451 auto RemoteLIDs_d = RemoteLIDs.view_device();
8452 auto PermuteToLIDs_d = PermuteToLIDs.view_device();
8453 auto PermuteFromLIDs_d = PermuteFromLIDs.view_device();
8455 Kokkos::View<size_t*, device_type> CSR_rowptr_d;
8456 Kokkos::View<GO*, device_type> CSR_colind_GID_d;
8457 Kokkos::View<LO*, device_type> CSR_colind_LID_d;
8458 Kokkos::View<impl_scalar_type*, device_type> CSR_vals_d;
8459 Kokkos::View<int*, device_type> TargetPids_d;
8464 destMat->imports_.view_device(),
8465 destMat->numImportPacketsPerLID_.view_device(),
8477 Kokkos::resize(CSR_colind_LID_d, CSR_colind_GID_d.size());
8479#ifdef HAVE_TPETRA_MMM_TIMINGS
8480 tmCopySPRdata = Teuchos::null;
8489 std::ostringstream os;
8490 os << *verbosePrefix <<
"Calling lowCommunicationMakeColMapAndReindex"
8492 std::cerr << os.str();
8495#ifdef HAVE_TPETRA_MMM_TIMINGS
8496 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC makeColMap")));
8508 std::ostringstream os;
8509 os << *verbosePrefix <<
"restrictComm="
8510 << (restrictComm ?
"true" :
"false") << std::endl;
8511 std::cerr << os.str();
8518#ifdef HAVE_TPETRA_MMM_TIMINGS
8519 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC restrict colmap")));
8522 ReducedColMap = (MyRowMap.getRawPtr() == MyColMap.getRawPtr()) ? ReducedRowMap : MyColMap->replaceCommWithSubset(ReducedComm);
8523 MyColMap = ReducedColMap;
8528 std::ostringstream os;
8529 os << *verbosePrefix <<
"Calling replaceColMap" << std::endl;
8530 std::cerr << os.str();
8532 destMat->replaceColMap(MyColMap);
8539 if (ReducedComm.is_null()) {
8541 std::ostringstream os;
8542 os << *verbosePrefix <<
"I am no longer in the communicator; "
8545 std::cerr << os.str();
8555 if ((!reverseMode && xferAsImport !=
nullptr) ||
8556 (reverseMode && xferAsExport !=
nullptr)) {
8558 std::ostringstream os;
8559 os << *verbosePrefix <<
"Calling sortCrsEntries" << endl;
8560 std::cerr << os.str();
8562#ifdef HAVE_TPETRA_MMM_TIMINGS
8563 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC sortCrsEntries")));
8568 }
else if ((!reverseMode && xferAsExport !=
nullptr) ||
8569 (reverseMode && xferAsImport !=
nullptr)) {
8571 std::ostringstream os;
8572 os << *verbosePrefix <<
"Calling sortAndMergeCrsEntries"
8574 std::cerr << os.str();
8576#ifdef HAVE_TPETRA_MMM_TIMINGS
8577 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC sortAndMergeCrsEntries")));
8583 TEUCHOS_TEST_FOR_EXCEPTION(
8584 true, std::logic_error,
8585 "Tpetra::CrsMatrix::"
8586 "transferAndFillComplete: Should never get here! "
8587 "Please report this bug to a Tpetra developer.");
8595 std::ostringstream os;
8596 os << *verbosePrefix <<
"Calling destMat->setAllValues" << endl;
8597 std::cerr << os.str();
8601#ifdef HAVE_TPETRA_MMM_TIMINGS
8602 Teuchos::TimeMonitor MMrc(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC setAllValues")));
8604 destMat->setAllValues(CSR_rowptr_d, CSR_colind_LID_d, CSR_vals_d);
8612#ifdef HAVE_TPETRA_MMM_TIMINGS
8613 RCP<TimeMonitor> tmIESFC = rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC build importer and esfc"))));
8616 Teuchos::ParameterList esfc_params;
8618 RCP<import_type> MyImport;
8621 if (iallreduceRequest.get() !=
nullptr) {
8623 std::ostringstream os;
8624 os << *verbosePrefix <<
"Calling iallreduceRequest->wait()"
8626 std::cerr << os.str();
8628 iallreduceRequest->wait();
8629 if (reduced_mismatch != 0) {
8635#ifdef HAVE_TPETRA_MMM_TIMINGS
8636 Teuchos::TimeMonitor MMisMM(*TimeMonitor::getNewTimer(prefix + std::string(
"isMM Block")));
8641 std::ostringstream os;
8642 os << *verbosePrefix <<
"Getting CRS pointers" << endl;
8643 std::cerr << os.str();
8646 Teuchos::ArrayRCP<LocalOrdinal> type3LIDs;
8647 Teuchos::ArrayRCP<int> type3PIDs;
8648 auto rowptr = getCrsGraph()->getLocalRowPtrsHost();
8649 auto colind = getCrsGraph()->getLocalIndicesHost();
8652 std::ostringstream os;
8653 os << *verbosePrefix <<
"Calling reverseNeighborDiscovery" << std::endl;
8654 std::cerr << os.str();
8658#ifdef HAVE_TPETRA_MMM_TIMINGS
8659 TimeMonitor tm_rnd(*TimeMonitor::getNewTimer(prefix + std::string(
"isMMrevNeighDis")));
8661 Import_Util::reverseNeighborDiscovery(*
this,
8673 std::ostringstream os;
8674 os << *verbosePrefix <<
"Done with reverseNeighborDiscovery" << std::endl;
8675 std::cerr << os.str();
8678 Teuchos::ArrayView<const int> EPID1 = MyImporter.is_null() ? Teuchos::ArrayView<const int>() : MyImporter->getExportPIDs();
8679 Teuchos::ArrayView<const LO> ELID1 = MyImporter.is_null() ? Teuchos::ArrayView<const LO>() : MyImporter->getExportLIDs();
8681 Teuchos::ArrayView<const int> TEPID2 = rowTransfer.getExportPIDs();
8682 Teuchos::ArrayView<const LO> TELID2 = rowTransfer.getExportLIDs();
8684 const int numCols = getGraph()->getColMap()->getLocalNumElements();
8686 std::vector<bool> IsOwned(numCols,
true);
8687 std::vector<int> SentTo(numCols, -1);
8688 if (!MyImporter.is_null()) {
8689 for (
auto&& rlid : MyImporter->getRemoteLIDs()) {
8690 IsOwned[rlid] =
false;
8694 std::vector<std::pair<int, GO>> usrtg;
8695 usrtg.reserve(TEPID2.size());
8698 const auto& colMap = *(this->getColMap());
8700 const LO row = TELID2[i];
8701 const int pid = TEPID2[i];
8702 for (
auto j = rowptr[row]; j < rowptr[row + 1]; ++j) {
8703 const int col = colind[j];
8704 if (IsOwned[col] && SentTo[col] != pid) {
8706 GO gid = colMap.getGlobalElement(col);
8707 usrtg.push_back(std::pair<int, GO>(pid, gid));
8714 std::sort(usrtg.begin(), usrtg.end());
8715 auto eopg = std ::unique(usrtg.begin(), usrtg.end());
8717 usrtg.erase(eopg, usrtg.end());
8720 Teuchos::ArrayRCP<int> EPID2 = Teuchos::arcp(
new int[type2_us_size], 0, type2_us_size,
true);
8721 Teuchos::ArrayRCP<LO> ELID2 = Teuchos::arcp(
new LO[type2_us_size], 0, type2_us_size,
true);
8724 for (
auto&& p : usrtg) {
8725 EPID2[pos] = p.first;
8726 ELID2[pos] = this->getDomainMap()->getLocalElement(p.second);
8730 Teuchos::ArrayView<int> EPID3 = type3PIDs();
8731 Teuchos::ArrayView<LO> ELID3 = type3LIDs();
8732 GO InfGID = std::numeric_limits<GO>::max();
8733 int InfPID = INT_MAX;
8737#define TPETRA_MIN3(x, y, z) ((x) < (y) ? (std::min(x, z)) : (std::min(y, z)))
8738 int i1 = 0, i2 = 0, i3 = 0;
8739 int Len1 = EPID1.size();
8740 int Len2 = EPID2.size();
8741 int Len3 = EPID3.size();
8743 int MyLen = Len1 + Len2 + Len3;
8744 Teuchos::ArrayRCP<LO> userExportLIDs = Teuchos::arcp(
new LO[MyLen], 0, MyLen,
true);
8745 Teuchos::ArrayRCP<int> userExportPIDs = Teuchos::arcp(
new int[MyLen], 0, MyLen,
true);
8748 while (i1 < Len1 || i2 < Len2 || i3 < Len3) {
8749 int PID1 = (i1 < Len1) ? (EPID1[i1]) : InfPID;
8750 int PID2 = (i2 < Len2) ? (EPID2[i2]) : InfPID;
8751 int PID3 = (i3 < Len3) ? (EPID3[i3]) : InfPID;
8753 GO GID1 = (i1 < Len1) ? getDomainMap()->getGlobalElement(ELID1[i1]) : InfGID;
8754 GO GID2 = (i2 < Len2) ? getDomainMap()->getGlobalElement(ELID2[i2]) : InfGID;
8755 GO GID3 = (i3 < Len3) ? getDomainMap()->getGlobalElement(ELID3[i3]) : InfGID;
8757 int MIN_PID = TPETRA_MIN3(PID1, PID2, PID3);
8758 GO MIN_GID = TPETRA_MIN3(((PID1 == MIN_PID) ? GID1 : InfGID), ((PID2 == MIN_PID) ? GID2 : InfGID), ((PID3 == MIN_PID) ? GID3 : InfGID));
8762 bool added_entry =
false;
8764 if (PID1 == MIN_PID && GID1 == MIN_GID) {
8765 userExportLIDs[iloc] = ELID1[i1];
8766 userExportPIDs[iloc] = EPID1[i1];
8771 if (PID2 == MIN_PID && GID2 == MIN_GID) {
8773 userExportLIDs[iloc] = ELID2[i2];
8774 userExportPIDs[iloc] = EPID2[i2];
8780 if (PID3 == MIN_PID && GID3 == MIN_GID) {
8782 userExportLIDs[iloc] = ELID3[i3];
8783 userExportPIDs[iloc] = EPID3[i3];
8791 std::ostringstream os;
8792 os << *verbosePrefix <<
"Create Import" << std::endl;
8793 std::cerr << os.str();
8796#ifdef HAVE_TPETRA_MMM_TIMINGS
8797 auto ismmIctor(*TimeMonitor::getNewTimer(prefix + std::string(
"isMMIportCtor")));
8799 Teuchos::RCP<Teuchos::ParameterList> plist = rcp(
new Teuchos::ParameterList());
8801 if ((MyDomainMap != MyColMap) && (!MyDomainMap->isSameAs(*MyColMap)))
8805 userExportLIDs.view(0, iloc).getConst(),
8806 userExportPIDs.view(0, iloc).getConst(),
8810 std::ostringstream os;
8811 os << *verbosePrefix <<
"Call expertStaticFillComplete" << std::endl;
8812 std::cerr << os.str();
8816#ifdef HAVE_TPETRA_MMM_TIMINGS
8817 TimeMonitor esfc(*TimeMonitor::getNewTimer(prefix + std::string(
"isMM::destMat->eSFC")));
8818 esfc_params.set(
"Timer Label", label + std::string(
"isMM eSFC"));
8820 if (!params.is_null())
8821 esfc_params.set(
"compute global constants", params->get(
"compute global constants",
true));
8822 destMat->expertStaticFillComplete(MyDomainMap, MyRangeMap, MyImport, Teuchos::null, rcp(
new Teuchos::ParameterList(esfc_params)));
8827#ifdef HAVE_TPETRA_MMM_TIMINGS
8828 TimeMonitor MMnotMMblock(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC notMMblock")));
8831 std::ostringstream os;
8832 os << *verbosePrefix <<
"Create Import" << std::endl;
8833 std::cerr << os.str();
8836#ifdef HAVE_TPETRA_MMM_TIMINGS
8837 TimeMonitor notMMIcTor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC notMMCreateImporter")));
8839 Teuchos::RCP<Teuchos::ParameterList> mypars = rcp(
new Teuchos::ParameterList);
8840 mypars->set(
"Timer Label",
"notMMFrom_tAFC");
8841 if ((MyDomainMap != MyColMap) && (!MyDomainMap->isSameAs(*MyColMap)))
8842 MyImport = rcp(
new import_type(MyDomainMap, MyColMap, RemotePids, mypars));
8845 std::ostringstream os;
8846 os << *verbosePrefix <<
"Call expertStaticFillComplete" << endl;
8847 std::cerr << os.str();
8850#ifdef HAVE_TPETRA_MMM_TIMINGS
8851 TimeMonitor esfcnotmm(*TimeMonitor::getNewTimer(prefix + std::string(
"notMMdestMat->expertStaticFillComplete")));
8852 esfc_params.set(
"Timer Label", prefix + std::string(
"notMM eSFC"));
8854 esfc_params.set(
"Timer Label", std::string(
"notMM eSFC"));
8857 if (!params.is_null()) {
8858 esfc_params.set(
"compute global constants",
8859 params->get(
"compute global constants",
true));
8861 destMat->expertStaticFillComplete(MyDomainMap, MyRangeMap,
8862 MyImport, Teuchos::null,
8863 rcp(
new Teuchos::ParameterList(esfc_params)));
8866#ifdef HAVE_TPETRA_MMM_TIMINGS
8867 tmIESFC = Teuchos::null;
8871 std::ostringstream os;
8872 os << *verbosePrefix <<
"Done" << endl;
8873 std::cerr << os.str();
8877template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8879 importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8881 const Teuchos::RCP<const map_type>& domainMap,
8882 const Teuchos::RCP<const map_type>& rangeMap,
8883 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8884 transferAndFillComplete(destMatrix, importer, Teuchos::null, domainMap, rangeMap, params);
8887template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8889 importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8892 const Teuchos::RCP<const map_type>& domainMap,
8893 const Teuchos::RCP<const map_type>& rangeMap,
8894 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8895 transferAndFillComplete(destMatrix, rowImporter, Teuchos::rcpFromRef(domainImporter), domainMap, rangeMap, params);
8898template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8900 exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8902 const Teuchos::RCP<const map_type>& domainMap,
8903 const Teuchos::RCP<const map_type>& rangeMap,
8904 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8905 transferAndFillComplete(destMatrix, exporter, Teuchos::null, domainMap, rangeMap, params);
8908template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8910 exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8913 const Teuchos::RCP<const map_type>& domainMap,
8914 const Teuchos::RCP<const map_type>& rangeMap,
8915 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8916 transferAndFillComplete(destMatrix, rowExporter, Teuchos::rcpFromRef(domainExporter), domainMap, rangeMap, params);
8919template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8920void copyAndPermuteStaticGraphNew(
8923 const size_t numSameIDs,
8924 const LocalOrdinal permuteToLIDs[],
8925 const LocalOrdinal permuteFromLIDs[],
8926 const size_t numPermutes) {
8929 using Teuchos::Array;
8930 using LO = LocalOrdinal;
8931 using GO = GlobalOrdinal;
8933#if KOKKOS_VERSION >= 40799
8934 using impl_scalar_type =
typename KokkosKernels::ArithTraits<Scalar>::val_type;
8941 typedef typename crs_matrix_type::local_inds_device_view_type::non_const_value_type local_inds_device_value_t;
8943 typedef typename Node::execution_space exec_space;
8944 typedef Kokkos::RangePolicy<exec_space, LO> range_type;
8946 const LocalOrdinal LINV = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
8948 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteStaticGraphNew");
8950 const crs_matrix_type* srcMatCrsPtr =
dynamic_cast<const crs_matrix_type*
>(&srcMat);
8951 TEUCHOS_TEST_FOR_EXCEPTION(srcMatCrsPtr ==
nullptr, std::runtime_error,
"bad srcMatCrsPtr");
8952 const crs_matrix_type& srcMatCrs = *srcMatCrsPtr;
8954 crs_matrix_type* tgtMatCrsPtr =
dynamic_cast<crs_matrix_type*
>(&tgtMat);
8955 TEUCHOS_TEST_FOR_EXCEPTION(tgtMatCrsPtr ==
nullptr, std::runtime_error,
"bad tgtMatCrsPtr");
8956 crs_matrix_type& tgtMatCrs = *tgtMatCrsPtr;
8964 const auto& srcRowMap = *(srcMat.
getRowMap());
8965 auto comm = srcRowMap.getComm();
8967 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
8969 auto my_replaceGlobalValuesImpl_scalar = KOKKOS_LAMBDA(
8970 const bool sorted,
const bool atomic,
size_t hint[],
8971 const size_t numInTgtRow,
const local_inds_device_value_t tgtColInds[],
impl_scalar_type tgtRowVals[],
8977 if (lclColInd != LINV) {
8978 const size_t offset = KokkosSparse::findRelOffset(tgtColInds, numInTgtRow, lclColInd, hint[0], sorted);
8979 if (offset != numInTgtRow) {
8980 Kokkos::atomic_store(&tgtRowVals[offset], newVals);
8981 hint[0] = offset + 1;
8986 if (lclColInd != LINV) {
8987 const size_t offset = KokkosSparse::findRelOffset(tgtColInds, numInTgtRow, lclColInd, hint[0], sorted);
8988 if (offset != numInTgtRow) {
8989 tgtRowVals[offset] = newVals;
8990 hint[0] = offset + 1;
8998 if (sourceIsLocallyIndexed) {
8999 typename crs_matrix_type::row_ptrs_device_view_type tgtLocalRowPtrsDevice = tgtMatCrs.getLocalRowPtrsDevice();
9000 typename crs_matrix_type::local_inds_device_view_type tgtLocalColIndsDevice = tgtMatCrs.getLocalIndicesDevice();
9001 typename crs_matrix_type::row_ptrs_host_view_type srcLocalRowPtrsHost = srcMatCrs.getLocalRowPtrsHost();
9002 typename crs_matrix_type::row_ptrs_device_view_type srcLocalRowPtrsDevice = srcMatCrs.getLocalRowPtrsDevice();
9003 typename crs_matrix_type::local_inds_device_view_type srcLocalColIndsDevice = srcMatCrs.getLocalIndicesDevice();
9005 bool tgtMatIsSorted = tgtMatCrs.getCrsGraph()->isSorted();
9007 using local_map_type =
typename crs_matrix_type::map_type::local_map_type;
9009 local_map_type local_map = srcMat.
getRowMap()->getLocalMap();
9010 local_map_type local_col_map = srcMat.
getColMap()->getLocalMap();
9011 local_map_type tgt_local_map = tgtMatCrs.getRowMap()->getLocalMap();
9012 local_map_type tgt_local_col_map = tgtMatCrs.getColMap()->getLocalMap();
9014 auto vals = srcMatCrs.getLocalValuesDevice(Access::ReadOnly);
9015 auto tvals = tgtMatCrs.getLocalValuesDevice(Access::ReadWrite);
9017 Kokkos::parallel_for(
9018 "Tpetra_CrsMatrix::copyAndPermuteStaticGraph",
9019 range_type(0, numSameIDs_as_LID),
9020 KOKKOS_LAMBDA(
const LO sourceLID) {
9021 local_inds_device_value_t
start = srcLocalRowPtrsDevice(sourceLID);
9022 local_inds_device_value_t end = srcLocalRowPtrsDevice(sourceLID + 1);
9023 local_inds_device_value_t rowLength = (end -
start);
9025 local_inds_device_value_t tstart = tgtLocalRowPtrsDevice(sourceLID);
9026 local_inds_device_value_t tend = tgtLocalRowPtrsDevice(sourceLID + 1);
9027 local_inds_device_value_t numInTgtRow = (tend - tstart);
9029 KOKKOS_ASSERT(
static_cast<size_t>(tstart) < tvals.extent(0));
9031 const local_inds_device_value_t* tgtColInds = &tgtLocalColIndsDevice(tstart);
9034 for (LO j = 0; j < rowLength; j++) {
9035 local_inds_device_value_t ci = srcLocalColIndsDevice(start + j);
9036 GO gi = local_col_map.getGlobalElement(ci);
9037 const local_inds_device_value_t lclColInd = tgt_local_col_map.getLocalElement(gi);
9038 my_replaceGlobalValuesImpl_scalar(
9039 tgtMatIsSorted,
false, &hint, numInTgtRow, tgtColInds, tgtRowVals, lclColInd, vals(start + j));
9043 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
9047 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
9048 const GO targetGID = sourceGID;
9050 Teuchos::ArrayView<const GO> rowIndsConstView;
9051 Teuchos::ArrayView<const Scalar> rowValsConstView;
9053 typename crs_matrix_type::global_inds_host_view_type rowIndsView;
9054 typename crs_matrix_type::values_host_view_type rowValsView;
9060 rowIndsConstView = Teuchos::ArrayView<const GO>(
9061 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9062 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9063 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9064 rowValsView.extent(0),
9065 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9072 tgtMatCrs.replaceGlobalValues(targetGID, rowIndsConstView, rowValsConstView);
9080 typename crs_matrix_type::nonconst_global_inds_host_view_type rowInds;
9081 typename crs_matrix_type::nonconst_values_host_view_type rowVals;
9083 const auto& tgtRowMap = *(tgtMat.
getRowMap());
9084 for (
size_t p = 0; p < numPermutes; ++p) {
9085 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
9086 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
9088 Teuchos::ArrayView<const GO> rowIndsConstView;
9089 Teuchos::ArrayView<const Scalar> rowValsConstView;
9091 if (sourceIsLocallyIndexed) {
9093 if (rowLength >
static_cast<size_t>(rowInds.size())) {
9094 Kokkos::resize(rowInds, rowLength);
9095 Kokkos::resize(rowVals, rowLength);
9099 typename crs_matrix_type::nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(
9100 rowInds, std::make_pair((
size_t)0, rowLength));
9101 typename crs_matrix_type::nonconst_values_host_view_type rowValsView = Kokkos::subview(
9102 rowVals, std::make_pair((
size_t)0, rowLength));
9107 size_t checkRowLength = 0;
9108 srcMat.
getGlobalRowCopy(sourceGID, rowIndsView, rowValsView, checkRowLength);
9114 rowIndsConstView = Teuchos::ArrayView<const GO>(
9115 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9116 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9117 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9118 rowValsView.extent(0),
9119 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9123 typename crs_matrix_type::global_inds_host_view_type rowIndsView;
9124 typename crs_matrix_type::values_host_view_type rowValsView;
9130 rowIndsConstView = Teuchos::ArrayView<const GO>(
9131 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9132 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9133 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9134 rowValsView.extent(0),
9135 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9140 tgtMatCrs.replaceGlobalValues(targetGID, rowIndsConstView, rowValsConstView);
9152#define TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9154 template class CrsMatrix<SCALAR, LO, GO, NODE>;
9156#define TPETRA_CRSMATRIX_CONVERT_INSTANT(SO, SI, LO, GO, NODE) \
9158 template Teuchos::RCP<CrsMatrix<SO, LO, GO, NODE>> \
9159 CrsMatrix<SI, LO, GO, NODE>::convert<SO>() const;
9161#define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9163 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9164 importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9165 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9166 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9167 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& importer, \
9168 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9169 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9170 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9171 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9172 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9173 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9174 const Teuchos::RCP<Teuchos::ParameterList>& params);
9176#define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9178 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9179 importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9180 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9181 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9182 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowImporter, \
9183 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9184 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9185 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainImporter, \
9186 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9187 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9188 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9189 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9190 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9191 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9192 const Teuchos::RCP<Teuchos::ParameterList>& params);
9194#define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9196 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9197 exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9198 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9199 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9200 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& exporter, \
9201 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9202 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9203 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9204 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9205 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9206 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9207 const Teuchos::RCP<Teuchos::ParameterList>& params);
9209#define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9211 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9212 exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9213 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9214 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9215 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowExporter, \
9216 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9217 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9218 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainExporter, \
9219 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9220 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9221 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9222 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9223 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9224 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9225 const Teuchos::RCP<Teuchos::ParameterList>& params);
9227#define TPETRA_CRSMATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9228 TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9229 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9230 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9231 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9232 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE)