3992 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillActive(), std::runtime_error,
3993 "Fill must be active before "
3994 "you may call this method.");
3996 const size_t myNumNonlocalRows =
nonlocals_.size();
4003 const int iHaveNonlocalRows = (myNumNonlocalRows == 0) ? 0 : 1;
4004 int someoneHasNonlocalRows = 0;
4005 reduceAll<int, int>(*comm, REDUCE_MAX, iHaveNonlocalRows,
4006 outArg(someoneHasNonlocalRows));
4007 if (someoneHasNonlocalRows == 0) {
4021 RCP<const map_type> nonlocalRowMap;
4022 Teuchos::Array<size_t> numEntPerNonlocalRow(myNumNonlocalRows);
4024 Teuchos::Array<GO> myNonlocalGblRows(myNumNonlocalRows);
4025 size_type curPos = 0;
4026 for (
auto mapIter = nonlocals_.begin(); mapIter != nonlocals_.end();
4027 ++mapIter, ++curPos) {
4028 myNonlocalGblRows[curPos] = mapIter->first;
4031 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4032 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4039 sort2(gblCols.begin(), gblCols.end(), vals.begin());
4040 typename Teuchos::Array<GO>::iterator gblCols_newEnd;
4041 typename Teuchos::Array<Scalar>::iterator vals_newEnd;
4042 merge2(gblCols_newEnd, vals_newEnd,
4043 gblCols.begin(), gblCols.end(),
4044 vals.begin(), vals.end());
4045 gblCols.erase(gblCols_newEnd, gblCols.end());
4046 vals.erase(vals_newEnd, vals.end());
4047 numEntPerNonlocalRow[curPos] = gblCols.size();
4058 GO myMinNonlocalGblRow = std::numeric_limits<GO>::max();
4060 auto iter = std::min_element(myNonlocalGblRows.begin(),
4061 myNonlocalGblRows.end());
4062 if (iter != myNonlocalGblRows.end()) {
4063 myMinNonlocalGblRow = *iter;
4066 GO gblMinNonlocalGblRow = 0;
4067 reduceAll<int, GO>(*comm, REDUCE_MIN, myMinNonlocalGblRow,
4068 outArg(gblMinNonlocalGblRow));
4069 const GO indexBase = gblMinNonlocalGblRow;
4070 const global_size_t INV = Teuchos::OrdinalTraits<global_size_t>::invalid();
4071 nonlocalRowMap = rcp(
new map_type(INV, myNonlocalGblRows(), indexBase, comm));
4080 std::ostringstream os;
4081 os << *prefix <<
"Create nonlocal matrix" << endl;
4082 std::cerr << os.str();
4084 RCP<crs_matrix_type> nonlocalMatrix =
4085 rcp(
new crs_matrix_type(nonlocalRowMap, numEntPerNonlocalRow()));
4087 size_type curPos = 0;
4088 for (
auto mapIter = nonlocals_.begin(); mapIter != nonlocals_.end();
4089 ++mapIter, ++curPos) {
4090 const GO gblRow = mapIter->first;
4092 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4093 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4095 nonlocalMatrix->insertGlobalValues(gblRow, gblCols(), vals());
4107 auto origRowMap = this->getRowMap();
4108 const bool origRowMapIsOneToOne = origRowMap->isOneToOne();
4110 int isLocallyComplete = 1;
4112 if (origRowMapIsOneToOne) {
4114 std::ostringstream os;
4115 os << *prefix <<
"Original row Map is 1-to-1" << endl;
4116 std::cerr << os.str();
4118 export_type exportToOrig(nonlocalRowMap, origRowMap);
4119 if (!exportToOrig.isLocallyComplete()) {
4120 isLocallyComplete = 0;
4123 std::ostringstream os;
4124 os << *prefix <<
"doExport from nonlocalMatrix" << endl;
4125 std::cerr << os.str();
4127 this->doExport(*nonlocalMatrix, exportToOrig,
Tpetra::ADD);
4131 std::ostringstream os;
4132 os << *prefix <<
"Original row Map is NOT 1-to-1" << endl;
4133 std::cerr << os.str();
4140 export_type exportToOneToOne(nonlocalRowMap, oneToOneRowMap);
4141 if (!exportToOneToOne.isLocallyComplete()) {
4142 isLocallyComplete = 0;
4150 std::ostringstream os;
4151 os << *prefix <<
"Create & doExport into 1-to-1 matrix"
4153 std::cerr << os.str();
4155 crs_matrix_type oneToOneMatrix(oneToOneRowMap, 0);
4157 oneToOneMatrix.doExport(*nonlocalMatrix, exportToOneToOne,
4163 std::ostringstream os;
4164 os << *prefix <<
"Free nonlocalMatrix" << endl;
4165 std::cerr << os.str();
4167 nonlocalMatrix = Teuchos::null;
4171 std::ostringstream os;
4172 os << *prefix <<
"doImport from 1-to-1 matrix" << endl;
4173 std::cerr << os.str();
4175 import_type importToOrig(oneToOneRowMap, origRowMap);
4176 this->doImport(oneToOneMatrix, importToOrig,
Tpetra::ADD);
4184 std::ostringstream os;
4185 os << *prefix <<
"Free nonlocals_ (std::map)" << endl;
4186 std::cerr << os.str();
4188 decltype(nonlocals_) newNonlocals;
4189 std::swap(nonlocals_, newNonlocals);
4198 int isGloballyComplete = 0;
4199 reduceAll<int, int>(*comm, REDUCE_MIN, isLocallyComplete,
4200 outArg(isGloballyComplete));
4201 TEUCHOS_TEST_FOR_EXCEPTION(isGloballyComplete != 1, std::runtime_error,
4202 "On at least one process, "
4203 "you called insertGlobalValues with a global row index which is not in "
4204 "the matrix's row Map on any process in its communicator.");
4207template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4209 resumeFill(
const Teuchos::RCP<Teuchos::ParameterList>& params) {
4211 myGraph_->resumeFill(params);
4214 applyHelper.reset();
4218template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4221 return getCrsGraphRef().haveGlobalConstants();
4224template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4226 fillComplete(
const Teuchos::RCP<Teuchos::ParameterList>& params) {
4227 const char tfecfFuncName[] =
"fillComplete(params): ";
4229 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->
getCrsGraph().is_null(), std::logic_error,
4230 "getCrsGraph() returns null. This should not happen at this point. "
4231 "Please report this bug to the Tpetra developers.");
4240 Teuchos::RCP<const map_type> rangeMap = graph.
getRowMap();
4241 Teuchos::RCP<const map_type> domainMap = rangeMap;
4246template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4248 fillComplete(
const Teuchos::RCP<const map_type>& domainMap,
4249 const Teuchos::RCP<const map_type>& rangeMap,
4250 const Teuchos::RCP<Teuchos::ParameterList>& params) {
4254 using Teuchos::ArrayRCP;
4257 const char tfecfFuncName[] =
"fillComplete: ";
4258 ProfilingRegion regionFillComplete(
"Tpetra::CrsMatrix::fillComplete");
4259 const bool verbose = Behavior::verbose(
"CrsMatrix");
4260 std::unique_ptr<std::string> prefix;
4262 prefix = this->createPrefix(
"CrsMatrix",
"fillComplete(dom,ran,p)");
4263 std::ostringstream os;
4264 os << *prefix << endl;
4265 std::cerr << os.str();
4268 "Tpetra::CrsMatrix::fillComplete",
4272 "Matrix fill state must be active (isFillActive() "
4273 "must be true) before you may call fillComplete().");
4274 const int numProcs = this->
getComm()->getSize();
4284 bool assertNoNonlocalInserts =
false;
4287 bool sortGhosts =
true;
4289 if (!params.is_null()) {
4290 assertNoNonlocalInserts = params->get(
"No Nonlocal Changes",
4291 assertNoNonlocalInserts);
4292 if (params->isParameter(
"sort column map ghost gids")) {
4293 sortGhosts = params->get(
"sort column map ghost gids", sortGhosts);
4294 }
else if (params->isParameter(
"Sort column Map ghost GIDs")) {
4295 sortGhosts = params->get(
"Sort column Map ghost GIDs", sortGhosts);
4300 const bool needGlobalAssemble = !assertNoNonlocalInserts && numProcs > 1;
4302 if (!this->myGraph_.is_null()) {
4303 this->myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
4306 if (!this->getCrsGraphRef().indicesAreAllocated()) {
4315 if (needGlobalAssemble) {
4318 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numProcs == 1 &&
nonlocals_.size() > 0,
4320 "Cannot have nonlocal entries on a serial run. "
4321 "An invalid entry (i.e., with row index not in the row Map) must have "
4322 "been submitted to the CrsMatrix.");
4333#ifdef HAVE_TPETRA_DEBUG
4351 const bool domainMapsMatch =
4352 this->staticGraph_->getDomainMap()->isSameAs(*domainMap);
4353 const bool rangeMapsMatch =
4354 this->staticGraph_->getRangeMap()->isSameAs(*rangeMap);
4356 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!domainMapsMatch, std::runtime_error,
4357 "The CrsMatrix's domain Map does not match the graph's domain Map. "
4358 "The graph cannot be changed because it was given to the CrsMatrix "
4359 "constructor as const. You can fix this by passing in the graph's "
4360 "domain Map and range Map to the matrix's fillComplete call.");
4362 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!rangeMapsMatch, std::runtime_error,
4363 "The CrsMatrix's range Map does not match the graph's range Map. "
4364 "The graph cannot be changed because it was given to the CrsMatrix "
4365 "constructor as const. You can fix this by passing in the graph's "
4366 "domain Map and range Map to the matrix's fillComplete call.");
4378 this->myGraph_->setDomainRangeMaps(domainMap, rangeMap);
4381 Teuchos::Array<int> remotePIDs(0);
4382 const bool mustBuildColMap = !this->
hasColMap();
4383 if (mustBuildColMap) {
4384 this->myGraph_->makeColMap(remotePIDs);
4389 const std::pair<size_t, std::string> makeIndicesLocalResult =
4390 this->myGraph_->makeIndicesLocal(verbose);
4395 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(makeIndicesLocalResult.first != 0, std::runtime_error,
4396 makeIndicesLocalResult.second);
4399 const bool sorted = this->myGraph_->isSorted();
4400 const bool merged = this->myGraph_->isMerged();
4407 this->myGraph_->makeImportExport(remotePIDs, mustBuildColMap);
4413 const bool callGraphComputeGlobalConstants = params.get() ==
nullptr ||
4414 params->get(
"compute global constants",
true);
4415 if (callGraphComputeGlobalConstants) {
4416 this->myGraph_->computeGlobalConstants();
4418 this->myGraph_->computeLocalConstants();
4420 this->myGraph_->fillComplete_ =
true;
4421 this->myGraph_->checkInternalState();
4431template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4434 const Teuchos::RCP<const map_type>& rangeMap,
4435 const Teuchos::RCP<const import_type>& importer,
4436 const Teuchos::RCP<const export_type>& exporter,
4437 const Teuchos::RCP<Teuchos::ParameterList>& params) {
4440 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
4443 "Matrix fill state must be active (isFillActive() "
4444 "must be true) before calling fillComplete().");
4445 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4446 myGraph_.is_null(), std::logic_error,
"myGraph_ is null. This is not allowed.");
4451 myGraph_->expertStaticFillComplete(domainMap, rangeMap, importer, exporter, params);
4465#ifdef HAVE_TPETRA_DEBUG
4466 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillActive(), std::logic_error,
4467 ": We're at the end of fillComplete(), but isFillActive() is true. "
4468 "Please report this bug to the Tpetra developers.");
4469 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete(), std::logic_error,
4470 ": We're at the end of fillComplete(), but isFillActive() is true. "
4471 "Please report this bug to the Tpetra developers.");
4477template <
class execution_space,
class LO,
class rowptr_type,
class colinds_type,
class numRowEntries_type,
class values_type>
4478void prepareSortMergeUnpackedMatrix(rowptr_type rowptr, colinds_type colinds, numRowEntries_type numRowEntries, values_type values) {
4479 using ATS = KokkosKernels::ArithTraits<LO>;
4480 using scalarATS = KokkosKernels::ArithTraits<typename values_type::value_type>;
4481 const auto unused = ATS::max();
4482 const auto zero = scalarATS::zero();
4484 auto numRows = rowptr.extent(0) - 1;
4487 Kokkos::parallel_for(
4488 "flag_unused_entries", Kokkos::RangePolicy<execution_space, LO>(0, numRows), KOKKOS_LAMBDA(
const LO rlid) {
4489 for (
size_t jj = rowptr(rlid) + numRowEntries(rlid); jj < rowptr(rlid + 1); ++jj) {
4490 colinds(jj) = unused;
4496template <
class execution_space,
class LO,
class rowptr_type,
class colinds_type,
class numRowEntries_type,
class values_type>
4497void mergeUnpackedMatrix(rowptr_type rowptr, colinds_type colinds, numRowEntries_type numRowEntries, values_type values) {
4498 auto numRows = rowptr.extent(0) - 1;
4504 Kokkos::parallel_for(
4505 "merge_entries", Kokkos::RangePolicy<execution_space>(0, numRows), KOKKOS_LAMBDA(
const LO rlid) {
4506 auto rowNNZ = numRowEntries(rlid);
4510 auto rowBegin = rowptr(rlid);
4511 auto pos = rowBegin;
4512 for (
size_t offset = rowBegin + 1; offset < rowBegin + rowNNZ; ++offset) {
4513 if ((colinds(offset) != colinds(pos))) {
4515 colinds(pos) = colinds(offset);
4516 values(pos) = values(offset);
4518 values(pos) += values(offset);
4521 numRowEntries(rlid) = pos + 1 - rowBegin;
4525template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4530 const char tfecfFuncName[] =
"sortAndMergeIndicesAndValues: ";
4532 if (!sorted || !merged) {
4533 ProfilingRegion regionSAM(
"Tpetra::CrsMatrix::sortAndMergeIndicesAndValues");
4535 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->
isStaticGraph(), std::runtime_error,
4536 "Cannot sort or merge with "
4537 "\"static\" (const) graph, since the matrix does not own the graph.");
4538 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(this->myGraph_.is_null(), std::logic_error,
4539 "myGraph_ is null, but "
4540 "this matrix claims ! isStaticGraph(). "
4541 "Please report this bug to the Tpetra developers.");
4543 "It is invalid to call "
4544 "this method if the graph's storage has already been optimized. "
4545 "Please report this bug to the Tpetra developers.");
4548 auto rowptr = graph.rowPtrsUnpacked_dev_;
4550 auto values = valuesUnpacked_wdv.getDeviceView(Access::ReadWrite);
4555 prepareSortMergeUnpackedMatrix<execution_space, LocalOrdinal>(rowptr, colinds, k_numRowEntries_d, values);
4564 mergeUnpackedMatrix<execution_space, LocalOrdinal>(rowptr, colinds, k_numRowEntries_d, values);
4571template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4576 Scalar beta)
const {
4579 using Teuchos::rcp_const_cast;
4580 using Teuchos::rcpFromRef;
4582 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
4583 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
4589 if (alpha ==
ZERO) {
4592 }
else if (beta != ONE) {
4606 RCP<const import_type> importer = this->
getGraph()->getImporter();
4607 RCP<const export_type> exporter = this->
getGraph()->getExporter();
4613 const bool Y_is_overwritten = (beta ==
ZERO);
4616 const bool Y_is_replicated =
4625 if (Y_is_replicated && this->
getComm()->getRank() > 0) {
4632 RCP<const MV> X_colMap;
4633 if (importer.is_null()) {
4643 X_colMap = rcp_const_cast<const MV>(X_colMapNonConst);
4647 X_colMap = rcpFromRef(X_in);
4650 ProfilingRegion regionImport(
"Tpetra::CrsMatrix::apply: Import");
4659 X_colMapNonConst->doImport(X_in, *importer,
INSERT);
4660 X_colMap = rcp_const_cast<const MV>(X_colMapNonConst);
4674 if (!exporter.is_null()) {
4675 this->
localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha,
ZERO);
4677 ProfilingRegion regionExport(
"Tpetra::CrsMatrix::apply: Export");
4683 if (Y_is_overwritten) {
4714 this->
localApply(*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
4717 this->
localApply(*X_colMap, Y_in, Teuchos::NO_TRANS, alpha, beta);
4725 if (Y_is_replicated) {
4726 ProfilingRegion regionReduce(
"Tpetra::CrsMatrix::apply: Reduce Y");
4731template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4735 const Teuchos::ETransp mode,
4737 Scalar beta)
const {
4738 using Teuchos::null;
4741 using Teuchos::rcp_const_cast;
4742 using Teuchos::rcpFromRef;
4744 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
4747 if (alpha ==
ZERO) {
4769 RCP<const import_type> importer = this->
getGraph()->getImporter();
4770 RCP<const export_type> exporter = this->
getGraph()->getExporter();
4776 const bool Y_is_overwritten = (beta ==
ZERO);
4777 if (Y_is_replicated && this->
getComm()->getRank() > 0) {
4783 X = rcp(
new MV(X_in, Teuchos::Copy));
4785 X = rcpFromRef(X_in);
4789 if (importer != Teuchos::null) {
4797 if (exporter != Teuchos::null) {
4808 if (!exporter.is_null()) {
4809 ProfilingRegion regionImport(
"Tpetra::CrsMatrix::apply (transpose): Import");
4817 if (importer != Teuchos::null) {
4818 ProfilingRegion regionExport(
"Tpetra::CrsMatrix::apply (transpose): Export");
4829 if (Y_is_overwritten) {
4846 MV Y(Y_in, Teuchos::Copy);
4850 this->
localApply(*X, Y_in, mode, alpha, beta);
4857 if (Y_is_replicated) {
4858 ProfilingRegion regionReduce(
"Tpetra::CrsMatrix::apply (transpose): Reduce Y");
4863template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4867 const Teuchos::ETransp mode,
4868 const Scalar& alpha,
4869 const Scalar& beta)
const {
4870 using Teuchos::NO_TRANS;
4872 ProfilingRegion regionLocalApply(
"Tpetra::CrsMatrix::localApply");
4879 const char tfecfFuncName[] =
"localApply: ";
4882 "Y.getNumVectors() = "
4884 const bool transpose = (mode != Teuchos::NO_TRANS);
4885 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!transpose && X.
getLocalLength() !=
4888 "NO_TRANS case: X has the wrong number of local rows. "
4889 "X.getLocalLength() = "
4891 "getColMap()->getLocalNumElements() = "
4892 <<
getColMap()->getLocalNumElements() <<
".");
4893 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!transpose && Y.
getLocalLength() !=
4896 "NO_TRANS case: Y has the wrong number of local rows. "
4897 "Y.getLocalLength() = "
4899 "getRowMap()->getLocalNumElements() = "
4900 <<
getRowMap()->getLocalNumElements() <<
".");
4901 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(transpose && X.
getLocalLength() !=
4904 "TRANS or CONJ_TRANS case: X has the wrong number of local "
4905 "rows. X.getLocalLength() = "
4907 <<
" != getRowMap()->getLocalNumElements() = "
4908 <<
getRowMap()->getLocalNumElements() <<
".");
4909 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(transpose && Y.
getLocalLength() !=
4912 "TRANS or CONJ_TRANS case: X has the wrong number of local "
4913 "rows. Y.getLocalLength() = "
4915 <<
" != getColMap()->getLocalNumElements() = "
4916 <<
getColMap()->getLocalNumElements() <<
".");
4917 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete(), std::runtime_error,
4918 "The matrix is not "
4919 "fill complete. You must call fillComplete() (possibly with "
4920 "domain and range Map arguments) without an intervening "
4921 "resumeFill() call before you may call this method.");
4923 std::runtime_error,
"X and Y must be constant stride.");
4928 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(X_lcl.data() == Y_lcl.data() && X_lcl.data() !=
nullptr && X_lcl.extent(0) != 0,
4929 std::runtime_error,
"X and Y may not alias one another.");
4934 if (!applyHelper.get()) {
4937 bool useMergePath =
false;
4938#ifdef KOKKOSKERNELS_ENABLE_TPL_CUSPARSE
4944 if constexpr (std::is_same_v<execution_space, Kokkos::Cuda>) {
4946 LocalOrdinal maxRowImbalance = 0;
4951 useMergePath =
true;
4954 applyHelper = std::make_shared<ApplyHelper>(A_lcl.nnz(), A_lcl.graph.row_map,
4955 useMergePath ? KokkosSparse::SPMV_MERGE_PATH : KokkosSparse::SPMV_DEFAULT);
4959 const char* modeKK =
nullptr;
4961 case Teuchos::NO_TRANS:
4962 modeKK = KokkosSparse::NoTranspose;
4964 case Teuchos::TRANS:
4965 modeKK = KokkosSparse::Transpose;
4967 case Teuchos::CONJ_TRANS:
4968 modeKK = KokkosSparse::ConjugateTranspose;
4971 throw std::invalid_argument(
"Tpetra::CrsMatrix::localApply: invalid mode");
4974 if (applyHelper->shouldUseIntRowptrs()) {
4975 auto A_lcl_int_rowptrs = applyHelper->getIntRowptrMatrix(A_lcl);
4977 &applyHelper->handle_int, modeKK,
4981 &applyHelper->handle, modeKK,
4986template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4990 Teuchos::ETransp mode,
4992 Scalar beta)
const {
4994 const char fnName[] =
"Tpetra::CrsMatrix::apply";
4996 TEUCHOS_TEST_FOR_EXCEPTION(!
isFillComplete(), std::runtime_error,
4997 fnName <<
": Cannot call apply() until fillComplete() "
4998 "has been called.");
5000 if (mode == Teuchos::NO_TRANS) {
5001 ProfilingRegion regionNonTranspose(fnName);
5004 ProfilingRegion regionTranspose(
"Tpetra::CrsMatrix::apply (transpose)");
5009template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5011Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node>>
5015 typedef CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node> output_matrix_type;
5016 const char tfecfFuncName[] =
"convert: ";
5018 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!this->
isFillComplete(), std::runtime_error,
5019 "This matrix (the source "
5020 "of the conversion) is not fill complete. You must first call "
5021 "fillComplete() (possibly with the domain and range Map) without an "
5022 "intervening call to resumeFill(), before you may call this method.");
5024 RCP<output_matrix_type> newMatrix(
new output_matrix_type(this->
getCrsGraph()));
5028 copyConvert(newMatrix->getLocalMatrixDevice().values,
5029 this->getLocalMatrixDevice().values);
5038template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5044 const char tfecfFuncName[] =
"checkInternalState: ";
5046 "Internal state is not consistent. "
5047 "Please report this bug to the Tpetra developers.";
5048 ProfilingRegion(
"Tpetra::CrsMatrix::checkInternalState");
5052 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(staticGraph_.is_null(), std::logic_error, err);
5056 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!myGraph_.is_null() && myGraph_ != staticGraph_,
5057 std::logic_error, err);
5059 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillComplete() && !staticGraph_->isFillComplete(),
5060 std::logic_error, err <<
" Specifically, the matrix is fill complete, "
5061 "but its graph is NOT fill complete.");
5064 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(staticGraph_->indicesAreAllocated() &&
5065 staticGraph_->getLocalAllocationSize() > 0 &&
5066 staticGraph_->getLocalNumRows() > 0 &&
5067 valuesUnpacked_wdv.extent(0) == 0,
5068 std::logic_error, err);
5072template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5076 std::ostringstream os;
5078 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
5079 if (this->getObjectLabel() !=
"") {
5080 os <<
"Label: \"" << this->getObjectLabel() <<
"\", ";
5083 os <<
"isFillComplete: true"
5089 os <<
"isFillComplete: false"
5096template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5098 describe(Teuchos::FancyOStream& out,
5099 const Teuchos::EVerbosityLevel verbLevel)
const {
5102 using Teuchos::ArrayView;
5103 using Teuchos::Comm;
5105 using Teuchos::TypeNameTraits;
5106 using Teuchos::VERB_DEFAULT;
5107 using Teuchos::VERB_EXTREME;
5108 using Teuchos::VERB_HIGH;
5109 using Teuchos::VERB_LOW;
5110 using Teuchos::VERB_MEDIUM;
5111 using Teuchos::VERB_NONE;
5113 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
5115 if (vl == VERB_NONE) {
5120 Teuchos::OSTab tab0(out);
5122 RCP<const Comm<int>> comm = this->
getComm();
5123 const int myRank = comm->getRank();
5124 const int numProcs = comm->getSize();
5129 width = std::max<size_t>(width,
static_cast<size_t>(11)) + 2;
5139 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
5141 Teuchos::OSTab tab1(out);
5144 if (this->getObjectLabel() !=
"") {
5145 out <<
"Label: \"" << this->getObjectLabel() <<
"\", ";
5148 out <<
"Template parameters:" << endl;
5149 Teuchos::OSTab tab2(out);
5150 out <<
"Scalar: " << TypeNameTraits<Scalar>::name() << endl
5151 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name() << endl
5152 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name() << endl
5153 <<
"Node: " << TypeNameTraits<Node>::name() << endl;
5156 out <<
"isFillComplete: true" << endl
5161 <<
"Global max number of entries in a row: "
5164 out <<
"isFillComplete: false" << endl
5170 if (vl < VERB_MEDIUM) {
5177 <<
"Row Map:" << endl;
5181 out <<
"null" << endl;
5192 out <<
"Column Map: ";
5196 out <<
"null" << endl;
5200 out <<
"same as row Map" << endl;
5211 out <<
"Domain Map: ";
5215 out <<
"null" << endl;
5219 out <<
"same as row Map" << endl;
5223 out <<
"same as column Map" << endl;
5234 out <<
"Range Map: ";
5238 out <<
"null" << endl;
5242 out <<
"same as domain Map" << endl;
5246 out <<
"same as row Map" << endl;
5256 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5257 if (myRank == curRank) {
5258 out <<
"Process rank: " << curRank << endl;
5259 Teuchos::OSTab tab2(out);
5260 if (!staticGraph_->indicesAreAllocated()) {
5261 out <<
"Graph indices not allocated" << endl;
5263 out <<
"Number of allocated entries: "
5264 << staticGraph_->getLocalAllocationSize() << endl;
5276 if (vl < VERB_HIGH) {
5281 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5282 if (myRank == curRank) {
5283 out << std::setw(width) <<
"Proc Rank"
5284 << std::setw(width) <<
"Global Row"
5285 << std::setw(width) <<
"Num Entries";
5286 if (vl == VERB_EXTREME) {
5287 out << std::setw(width) <<
"(Index,Value)";
5292 GlobalOrdinal gid =
getRowMap()->getGlobalElement(r);
5293 out << std::setw(width) << myRank
5294 << std::setw(width) << gid
5295 << std::setw(width) << nE;
5296 if (vl == VERB_EXTREME) {
5298 global_inds_host_view_type rowinds;
5299 values_host_view_type rowvals;
5301 for (
size_t j = 0; j < nE; ++j) {
5302 out <<
" (" << rowinds[j]
5303 <<
", " << rowvals[j]
5307 local_inds_host_view_type rowinds;
5308 values_host_view_type rowvals;
5310 for (
size_t j = 0; j < nE; ++j) {
5311 out <<
" (" <<
getColMap()->getGlobalElement(rowinds[j])
5312 <<
", " << rowvals[j]
5328template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5340 return (srcRowMat !=
nullptr);
5343template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5346 const typename crs_graph_type::padding_type& padding,
5347 const bool verbose) {
5352 using row_ptrs_type =
5353 typename local_graph_device_type::row_map_type::non_const_type;
5354 using range_policy =
5355 Kokkos::RangePolicy<execution_space, Kokkos::IndexType<LO>>;
5356 const char tfecfFuncName[] =
"applyCrsPadding";
5357 const char suffix[] =
5358 ". Please report this bug to the Tpetra developers.";
5359 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::applyCrsPadding");
5361 std::unique_ptr<std::string> prefix;
5363 prefix = this->createPrefix(
"CrsMatrix", tfecfFuncName);
5364 std::ostringstream os;
5365 os << *prefix <<
"padding: ";
5368 std::cerr << os.str();
5370 const int myRank = !verbose ? -1 : [&]() {
5371 auto map = this->getMap();
5372 if (map.is_null()) {
5375 auto comm = map->getComm();
5376 if (comm.is_null()) {
5379 return comm->getRank();
5383 if (!myGraph_->indicesAreAllocated()) {
5385 std::ostringstream os;
5386 os << *prefix <<
"Call allocateIndices" << endl;
5387 std::cerr << os.str();
5389 allocateValues(GlobalIndices, GraphNotYetAllocated, verbose);
5401 std::ostringstream os;
5402 os << *prefix <<
"Allocate row_ptrs_beg: "
5403 << myGraph_->getRowPtrsUnpackedHost().extent(0) << endl;
5404 std::cerr << os.str();
5406 using Kokkos::view_alloc;
5407 using Kokkos::WithoutInitializing;
5408 row_ptrs_type row_ptr_beg(view_alloc(
"row_ptr_beg", WithoutInitializing),
5409 myGraph_->rowPtrsUnpacked_dev_.extent(0));
5411 Kokkos::deep_copy(
execution_space(), row_ptr_beg, myGraph_->rowPtrsUnpacked_dev_);
5413 const size_t N = row_ptr_beg.extent(0) == 0 ? size_t(0) : size_t(row_ptr_beg.extent(0) - 1);
5415 std::ostringstream os;
5416 os << *prefix <<
"Allocate row_ptrs_end: " << N << endl;
5417 std::cerr << os.str();
5419 row_ptrs_type row_ptr_end(
5420 view_alloc(
"row_ptr_end", WithoutInitializing), N);
5422 row_ptrs_type num_row_entries_d;
5424 const bool refill_num_row_entries =
5425 myGraph_->k_numRowEntries_.extent(0) != 0;
5427 if (refill_num_row_entries) {
5430 num_row_entries_d = create_mirror_view_and_copy(memory_space(),
5431 myGraph_->k_numRowEntries_);
5432 Kokkos::parallel_for(
5433 "Fill end row pointers", range_policy(0, N),
5434 KOKKOS_LAMBDA(
const size_t i) {
5435 row_ptr_end(i) = row_ptr_beg(i) + num_row_entries_d(i);
5441 Kokkos::parallel_for(
5442 "Fill end row pointers", range_policy(0, N),
5443 KOKKOS_LAMBDA(
const size_t i) {
5444 row_ptr_end(i) = row_ptr_beg(i + 1);
5448 if (myGraph_->isGloballyIndexed()) {
5450 myGraph_->gblInds_wdv,
5451 valuesUnpacked_wdv, padding, myRank, verbose);
5452 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5453 const auto newColIndsLen = myGraph_->gblInds_wdv.extent(0);
5454 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(newValuesLen != newColIndsLen, std::logic_error,
5455 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5456 <<
" != myGraph_->gblInds_wdv.extent(0)=" << newColIndsLen
5460 myGraph_->lclIndsUnpacked_wdv,
5461 valuesUnpacked_wdv, padding, myRank, verbose);
5462 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5463 const auto newColIndsLen = myGraph_->lclIndsUnpacked_wdv.extent(0);
5464 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(newValuesLen != newColIndsLen, std::logic_error,
5465 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5466 <<
" != myGraph_->lclIndsUnpacked_wdv.extent(0)=" << newColIndsLen
5470 if (refill_num_row_entries) {
5471 Kokkos::parallel_for(
5472 "Fill num entries", range_policy(0, N),
5473 KOKKOS_LAMBDA(
const size_t i) {
5474 num_row_entries_d(i) = row_ptr_end(i) - row_ptr_beg(i);
5476 Kokkos::deep_copy(myGraph_->k_numRowEntries_, num_row_entries_d);
5480 std::ostringstream os;
5481 os << *prefix <<
"Assign myGraph_->rowPtrsUnpacked_; "
5482 <<
"old size: " << myGraph_->rowPtrsUnpacked_host_.extent(0)
5483 <<
", new size: " << row_ptr_beg.extent(0) << endl;
5484 std::cerr << os.str();
5485 TEUCHOS_ASSERT(myGraph_->getRowPtrsUnpackedHost().extent(0) ==
5486 row_ptr_beg.extent(0));
5488 myGraph_->setRowPtrsUnpacked(row_ptr_beg);
5491template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5492void copyAndPermuteStaticGraphNew(
5495 const size_t numSameIDs,
5496 const LocalOrdinal permuteToLIDs[],
5497 const LocalOrdinal permuteFromLIDs[],
5498 const size_t numPermutes);
5500template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5504 const size_t numSameIDs,
5505 const LocalOrdinal permuteToLIDs[],
5506 const LocalOrdinal permuteFromLIDs[],
5507 const size_t numPermutes) {
5510 using Teuchos::Array;
5511 using Teuchos::ArrayView;
5512 using LO = LocalOrdinal;
5513 using GO = GlobalOrdinal;
5514 const char tfecfFuncName[] =
"copyAndPermuteStaticGraph";
5515 const char suffix[] =
5516 " Please report this bug to the Tpetra developers.";
5517 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteStaticGraph");
5521 std::unique_ptr<std::string> prefix;
5523 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5524 std::ostringstream os;
5525 os << *prefix <<
"Start" << endl;
5527 const char*
const prefix_raw =
5528 verbose ? prefix.get()->c_str() :
nullptr;
5530 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed();
5535 const auto& srcRowMap = *(srcMat.getRowMap());
5536 nonconst_global_inds_host_view_type rowInds;
5537 nonconst_values_host_view_type rowVals;
5538 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
5539 if (sourceIsLocallyIndexed) {
5540 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5544 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5545 const GO targetGID = sourceGID;
5547 ArrayView<const GO> rowIndsConstView;
5548 ArrayView<const Scalar> rowValsConstView;
5550 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5551 if (rowLength >
static_cast<size_t>(rowInds.size())) {
5552 Kokkos::resize(rowInds, rowLength);
5553 Kokkos::resize(rowVals, rowLength);
5557 nonconst_global_inds_host_view_type rowIndsView;
5558 nonconst_values_host_view_type rowValsView;
5563 size_t checkRowLength = 0;
5566 const crs_matrix_type* srcMatCrsPtr =
dynamic_cast<const crs_matrix_type*
>(&srcMat);
5567 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(srcMatCrsPtr ==
nullptr, std::runtime_error,
"bad srcMatCrsPtr");
5568 const crs_matrix_type& srcMatCrs = *srcMatCrsPtr;
5570 auto globalRow = sourceGID;
5571 auto StaticGraphRCP = srcMatCrs.getGraph();
5573 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(StaticGraphPtr ==
nullptr, std::runtime_error,
"bad StaticGraphPtr");
5575 const RowInfo rowinfo = StaticGraph.getRowInfoFromGlobalRowIndex(globalRow);
5576 const size_t theNumEntries = rowinfo.numEntries;
5577 checkRowLength = theNumEntries;
5578 auto numEntries = theNumEntries;
5580 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid()) {
5581 if (StaticGraph.isLocallyIndexed()) {
5582 const map_type& colMap = *(StaticGraph.getColMap());
5583 auto curLclInds = StaticGraph.getLocalIndsViewHost(rowinfo);
5584 auto rowValsViewLocal = srcMatCrs.getValuesViewHost(rowinfo);
5585 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5586 reinterpret_cast<const Scalar*
>(rowValsViewLocal.data()),
5587 rowValsViewLocal.extent(0),
5588 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5589 auto rowIndsViewLocal = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5590 rowIndsConstView = Teuchos::ArrayView<const GO>(
5591 rowIndsViewLocal.data(), rowIndsViewLocal.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
5592 bool err = colMap.getGlobalElements(curLclInds.data(), numEntries, rowIndsViewLocal.data());
5593 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(err, std::runtime_error,
"getGlobalElements error");
5594 }
else if (StaticGraph.isGloballyIndexed()) {
5595 auto rowIndsViewLocal = StaticGraph.getGlobalIndsViewHost(rowinfo);
5596 rowIndsConstView = Teuchos::ArrayView<const GO>(
5597 rowIndsViewLocal.data(), rowIndsViewLocal.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
5598 auto rowValsViewLocal = srcMatCrs.getValuesViewHost(rowinfo);
5599 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5600 reinterpret_cast<const Scalar*
>(rowValsViewLocal.data()),
5601 rowValsViewLocal.extent(0),
5602 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5607 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5608 rowLength != checkRowLength,
5610 "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);
5613 combineGlobalValues(
5614 targetGID, rowIndsConstView, rowValsConstView,
REPLACE, prefix_raw, debug, verbose);
5617 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5621 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5622 const GO targetGID = sourceGID;
5624 ArrayView<const GO> rowIndsConstView;
5625 ArrayView<const Scalar> rowValsConstView;
5627 global_inds_host_view_type rowIndsView;
5628 values_host_view_type rowValsView;
5629 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5634 rowIndsConstView = Teuchos::ArrayView<const GO>(
5635 rowIndsView.data(), rowIndsView.extent(0),
5636 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5637 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5638 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5639 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5644 combineGlobalValues(targetGID, rowIndsConstView,
5646 prefix_raw, debug, verbose);
5651 std::ostringstream os;
5652 os << *prefix <<
"Do permutes" << endl;
5659 const map_type& tgtRowMap = *(this->getRowMap());
5660 for (
size_t p = 0; p < numPermutes; ++p) {
5661 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
5662 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
5664 ArrayView<const GO> rowIndsConstView;
5665 ArrayView<const Scalar> rowValsConstView;
5667 if (sourceIsLocallyIndexed) {
5668 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5669 if (rowLength >
static_cast<size_t>(rowInds.size())) {
5670 Kokkos::resize(rowInds, rowLength);
5671 Kokkos::resize(rowVals, rowLength);
5675 nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5676 nonconst_values_host_view_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5681 size_t checkRowLength = 0;
5682 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
5683 rowValsView, checkRowLength);
5685 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5687 "source matrix global row index "
5688 << sourceGID <<
", "
5689 "getNumEntriesInGlobalRow returns a row length of "
5690 << rowLength <<
", but getGlobalRowCopy a row length of "
5691 << checkRowLength <<
"." << suffix);
5698 rowIndsConstView = Teuchos::ArrayView<const GO>(
5699 rowIndsView.data(), rowIndsView.extent(0),
5700 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5701 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5702 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5703 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5707 global_inds_host_view_type rowIndsView;
5708 values_host_view_type rowValsView;
5709 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5714 rowIndsConstView = Teuchos::ArrayView<const GO>(
5715 rowIndsView.data(), rowIndsView.extent(0),
5716 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5717 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5718 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5719 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5724 combineGlobalValues(targetGID, rowIndsConstView,
5726 prefix_raw, debug, verbose);
5730 std::ostringstream os;
5731 os << *prefix <<
"Done" << endl;
5735template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5739 const size_t numSameIDs,
5740 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs_dv,
5741 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs_dv,
5742 const size_t numPermutes) {
5745 using Teuchos::Array;
5746 using Teuchos::ArrayView;
5747 using LO = LocalOrdinal;
5748 using GO = GlobalOrdinal;
5749 const char tfecfFuncName[] =
"copyAndPermuteNonStaticGraph";
5750 const char suffix[] =
5751 " Please report this bug to the Tpetra developers.";
5752 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteNonStaticGraph");
5756 std::unique_ptr<std::string> prefix;
5758 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5759 std::ostringstream os;
5760 os << *prefix <<
"Start" << endl;
5762 const char*
const prefix_raw =
5763 verbose ? prefix.get()->c_str() :
nullptr;
5767 const row_graph_type& srcGraph = *(srcMat.getGraph());
5769 myGraph_->computeCrsPadding(srcGraph, numSameIDs,
5770 permuteToLIDs_dv, permuteFromLIDs_dv, verbose);
5771 applyCrsPadding(*padding, verbose);
5773 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed();
5778 const map_type& srcRowMap = *(srcMat.getRowMap());
5779 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
5780 using gids_type = nonconst_global_inds_host_view_type;
5781 using vals_type = nonconst_values_host_view_type;
5784 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5788 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
5789 const GO targetGID = sourceGID;
5791 ArrayView<const GO> rowIndsConstView;
5792 ArrayView<const Scalar> rowValsConstView;
5794 if (sourceIsLocallyIndexed) {
5795 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5796 if (rowLength >
static_cast<size_t>(rowInds.extent(0))) {
5797 Kokkos::resize(rowInds, rowLength);
5798 Kokkos::resize(rowVals, rowLength);
5802 gids_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5803 vals_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5808 size_t checkRowLength = 0;
5809 srcMat.getGlobalRowCopy(sourceGID, rowIndsView, rowValsView,
5812 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5815 << sourceGID <<
", the source "
5816 "matrix's getNumEntriesInGlobalRow returns a row length "
5818 << rowLength <<
", but getGlobalRowCopy reports "
5820 << checkRowLength <<
"." << suffix);
5822 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
5823 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar*
>(rowValsView.data()), rowLength);
5825 global_inds_host_view_type rowIndsView;
5826 values_host_view_type rowValsView;
5827 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5833 rowIndsConstView = Teuchos::ArrayView<const GO>(
5834 rowIndsView.data(), rowIndsView.extent(0),
5835 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5836 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5837 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5838 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5844 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
5845 rowValsConstView, prefix_raw, debug, verbose);
5849 std::ostringstream os;
5850 os << *prefix <<
"Do permutes" << endl;
5852 const LO*
const permuteFromLIDs = permuteFromLIDs_dv.view_host().data();
5853 const LO*
const permuteToLIDs = permuteToLIDs_dv.view_host().data();
5855 const map_type& tgtRowMap = *(this->getRowMap());
5856 for (
size_t p = 0; p < numPermutes; ++p) {
5857 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
5858 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
5860 ArrayView<const GO> rowIndsConstView;
5861 ArrayView<const Scalar> rowValsConstView;
5863 if (sourceIsLocallyIndexed) {
5864 const size_t rowLength = srcMat.getNumEntriesInGlobalRow(sourceGID);
5865 if (rowLength >
static_cast<size_t>(rowInds.extent(0))) {
5866 Kokkos::resize(rowInds, rowLength);
5867 Kokkos::resize(rowVals, rowLength);
5871 gids_type rowIndsView = Kokkos::subview(rowInds, std::make_pair((
size_t)0, rowLength));
5872 vals_type rowValsView = Kokkos::subview(rowVals, std::make_pair((
size_t)0, rowLength));
5877 size_t checkRowLength = 0;
5878 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
5879 rowValsView, checkRowLength);
5881 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength, std::logic_error,
5883 "source matrix global row index "
5884 << sourceGID <<
", "
5885 "getNumEntriesInGlobalRow returns a row length of "
5886 << rowLength <<
", but getGlobalRowCopy a row length of "
5887 << checkRowLength <<
"." << suffix);
5889 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
5890 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar*
>(rowValsView.data()), rowLength);
5892 global_inds_host_view_type rowIndsView;
5893 values_host_view_type rowValsView;
5894 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5900 rowIndsConstView = Teuchos::ArrayView<const GO>(
5901 rowIndsView.data(), rowIndsView.extent(0),
5902 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5903 rowValsConstView = Teuchos::ArrayView<const Scalar>(
5904 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5905 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5911 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
5912 rowValsConstView, prefix_raw, debug, verbose);
5916 std::ostringstream os;
5917 os << *prefix <<
"Done" << endl;
5921template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5925 const size_t numSameIDs,
5926 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs,
5927 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs,
5935 const char tfecfFuncName[] =
"copyAndPermute: ";
5936 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermute");
5938 const bool verbose = Behavior::verbose(
"CrsMatrix");
5939 std::unique_ptr<std::string> prefix;
5941 prefix = this->createPrefix(
"CrsMatrix",
"copyAndPermute");
5942 std::ostringstream os;
5943 os << *prefix << endl
5944 << *prefix <<
" numSameIDs: " << numSameIDs << endl
5945 << *prefix <<
" numPermute: " << permuteToLIDs.extent(0)
5948 << dualViewStatusToString(permuteToLIDs,
"permuteToLIDs")
5951 << dualViewStatusToString(permuteFromLIDs,
"permuteFromLIDs")
5954 <<
"isStaticGraph: " << (
isStaticGraph() ?
"true" :
"false")
5956 std::cerr << os.str();
5959 const auto numPermute = permuteToLIDs.extent(0);
5960 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numPermute != permuteFromLIDs.extent(0),
5961 std::invalid_argument,
"permuteToLIDs.extent(0) = " << numPermute <<
"!= permuteFromLIDs.extent(0) = " << permuteFromLIDs.extent(0) <<
".");
5966 const RMT& srcMat =
dynamic_cast<const RMT&
>(srcObj);
5969 TEUCHOS_ASSERT(!permuteToLIDs.need_sync_device());
5970 auto permuteToLIDs_d = permuteToLIDs.view_device();
5971 TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_device());
5972 auto permuteFromLIDs_d = permuteFromLIDs.view_device();
5973 copyAndPermuteStaticGraphNew(
5974 srcMat, *
this, numSameIDs, permuteToLIDs_d.data(), permuteFromLIDs_d.data(), numPermute);
5977 TEUCHOS_ASSERT(!permuteToLIDs.need_sync_host());
5978 auto permuteToLIDs_h = permuteToLIDs.view_host();
5979 TEUCHOS_ASSERT(!permuteFromLIDs.need_sync_host());
5980 auto permuteFromLIDs_h = permuteFromLIDs.view_host();
5982 copyAndPermuteStaticGraph(srcMat, numSameIDs,
5983 permuteToLIDs_h.data(),
5984 permuteFromLIDs_h.data(),
5988 copyAndPermuteNonStaticGraph(srcMat, numSameIDs, permuteToLIDs,
5989 permuteFromLIDs, numPermute);
5993 std::ostringstream os;
5994 os << *prefix <<
"Done" << endl;
5995 std::cerr << os.str();
5999template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6002 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6003 Kokkos::DualView<char*, buffer_device_type>& exports,
6004 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6005 size_t& constantNumPackets) {
6010 using Teuchos::outArg;
6011 using Teuchos::REDUCE_MAX;
6012 using Teuchos::reduceAll;
6013 typedef LocalOrdinal LO;
6014 typedef GlobalOrdinal GO;
6015 const char tfecfFuncName[] =
"packAndPrepare: ";
6016 ProfilingRegion regionPAP(
"Tpetra::CrsMatrix::packAndPrepare");
6018 const bool debug = Behavior::debug(
"CrsMatrix");
6019 const bool verbose = Behavior::verbose(
"CrsMatrix");
6022 Teuchos::RCP<const Teuchos::Comm<int>> pComm = this->getComm();
6023 if (pComm.is_null()) {
6026 const Teuchos::Comm<int>& comm = *pComm;
6027 const int myRank = comm.getSize();
6029 std::unique_ptr<std::string> prefix;
6031 prefix = this->createPrefix(
"CrsMatrix",
"packAndPrepare");
6032 std::ostringstream os;
6033 os << *prefix <<
"Start" << endl
6035 << dualViewStatusToString(exportLIDs,
"exportLIDs")
6038 << dualViewStatusToString(exports,
"exports")
6041 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6043 std::cerr << os.str();
6066 std::ostringstream msg;
6070 const crs_matrix_type* srcCrsMat =
6071 dynamic_cast<const crs_matrix_type*
>(&source);
6072 if (srcCrsMat !=
nullptr) {
6074 std::ostringstream os;
6075 os << *prefix <<
"Source matrix same (CrsMatrix) type as target; "
6078 std::cerr << os.str();
6081 srcCrsMat->packNew(exportLIDs, exports, numPacketsPerLID,
6082 constantNumPackets);
6083 }
catch (std::exception& e) {
6085 msg <<
"Proc " << myRank <<
": " << e.what() << std::endl;
6088 using Kokkos::HostSpace;
6089 using Kokkos::subview;
6090 using exports_type = Kokkos::DualView<char*, buffer_device_type>;
6091 using range_type = Kokkos::pair<size_t, size_t>;
6094 std::ostringstream os;
6095 os << *prefix <<
"Source matrix NOT same (CrsMatrix) type as target"
6097 std::cerr << os.str();
6100 const row_matrix_type* srcRowMat =
6101 dynamic_cast<const row_matrix_type*
>(&source);
6102 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(srcRowMat ==
nullptr, std::invalid_argument,
6103 "The source object of the Import or Export operation is neither a "
6104 "CrsMatrix (with the same template parameters as the target object), "
6105 "nor a RowMatrix (with the same first four template parameters as the "
6116 TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
6117 auto exportLIDs_h = exportLIDs.view_host();
6118 Teuchos::ArrayView<const LO> exportLIDs_av(exportLIDs_h.data(),
6119 exportLIDs_h.size());
6123 Teuchos::Array<char> exports_a;
6129 numPacketsPerLID.clear_sync_state();
6130 numPacketsPerLID.modify_host();
6131 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
6132 Teuchos::ArrayView<size_t> numPacketsPerLID_av(numPacketsPerLID_h.data(),
6133 numPacketsPerLID_h.size());
6138 srcRowMat->pack(exportLIDs_av, exports_a, numPacketsPerLID_av,
6139 constantNumPackets);
6140 }
catch (std::exception& e) {
6142 msg <<
"Proc " << myRank <<
": " << e.what() << std::endl;
6146 const size_t newAllocSize =
static_cast<size_t>(exports_a.size());
6147 if (
static_cast<size_t>(exports.extent(0)) < newAllocSize) {
6148 const std::string oldLabel = exports.view_device().label();
6149 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6150 exports = exports_type(newLabel, newAllocSize);
6155 exports.modify_host();
6157 auto exports_h = exports.view_host();
6158 auto exports_h_sub = subview(exports_h, range_type(0, newAllocSize));
6162 typedef typename exports_type::t_host::execution_space HES;
6163 typedef Kokkos::Device<HES, HostSpace> host_device_type;
6164 Kokkos::View<const char*, host_device_type>
6165 exports_a_kv(exports_a.getRawPtr(), newAllocSize);
6167 Kokkos::deep_copy(exports_h_sub, exports_a_kv);
6172 reduceAll<int, int>(comm, REDUCE_MAX, lclBad, outArg(gblBad));
6175 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error,
6176 "packNew() or pack() threw an exception on "
6177 "one or more participating processes.");
6180 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(lclBad != 0, std::logic_error,
6181 "packNew threw an exception on one "
6182 "or more participating processes. Here is this process' error "
6188 std::ostringstream os;
6189 os << *prefix <<
"packAndPrepare: Done!" << endl
6199 std::cerr << os.str();
6203template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6207 const size_t offset,
6208 const size_t numEnt,
6209 const GlobalOrdinal gidsIn[],
6211 const size_t numBytesPerValue)
const {
6212 using Kokkos::subview;
6214 using Tpetra::Details::PackTraits;
6215 typedef LocalOrdinal LO;
6216 typedef GlobalOrdinal GO;
6225 const LO numEntLO =
static_cast<size_t>(numEnt);
6227 const size_t numEntBeg = offset;
6228 const size_t numEntLen = PackTraits<LO>::packValueCount(numEntLO);
6229 const size_t gidsBeg = numEntBeg + numEntLen;
6230 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
6231 const size_t valsBeg = gidsBeg + gidsLen;
6232 const size_t valsLen = numEnt * numBytesPerValue;
6234 char*
const numEntOut = exports + numEntBeg;
6235 char*
const gidsOut = exports + gidsBeg;
6236 char*
const valsOut = exports + valsBeg;
6238 size_t numBytesOut = 0;
6240 numBytesOut += PackTraits<LO>::packValue(numEntOut, numEntLO);
6243 Kokkos::pair<int, size_t> p;
6244 p = PackTraits<GO>::packArray(gidsOut, gidsIn, numEnt);
6245 errorCode += p.first;
6246 numBytesOut += p.second;
6248 p = PackTraits<ST>::packArray(valsOut, valsIn, numEnt);
6249 errorCode += p.first;
6250 numBytesOut += p.second;
6253 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6254 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
6257 << numBytesOut <<
" != expectedNumBytes = "
6258 << expectedNumBytes <<
".");
6259 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
6261 "PackTraits::packArray returned a nonzero error code");
6266template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6271 const char imports[],
6272 const size_t offset,
6273 const size_t numBytes,
6274 const size_t numEnt,
6275 const size_t numBytesPerValue) {
6276 using Kokkos::subview;
6278 using Tpetra::Details::PackTraits;
6279 typedef LocalOrdinal LO;
6280 typedef GlobalOrdinal GO;
6284 "Tpetra::CrsMatrix::unpackRow",
6287 if (numBytes == 0) {
6290 const int myRank = this->getMap()->getComm()->getRank();
6291 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6292 "unpackRow: The number of bytes to unpack numBytes=0, but the "
6293 "number of entries to unpack (as reported by numPacketsPerLID) "
6294 "for this row numEnt="
6295 << numEnt <<
" != 0.");
6300 if (numEnt == 0 && numBytes != 0) {
6301 const int myRank = this->getMap()->getComm()->getRank();
6302 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6303 "unpackRow: The number of entries to unpack (as reported by "
6304 "numPacketsPerLID) numEnt=0, but the number of bytes to unpack "
6306 << numBytes <<
" != 0.");
6312 const size_t numEntBeg = offset;
6313 const size_t numEntLen = PackTraits<LO>::packValueCount(lid);
6314 const size_t gidsBeg = numEntBeg + numEntLen;
6315 const size_t gidsLen = numEnt * PackTraits<GO>::packValueCount(gid);
6316 const size_t valsBeg = gidsBeg + gidsLen;
6317 const size_t valsLen = numEnt * numBytesPerValue;
6319 const char*
const numEntIn = imports + numEntBeg;
6320 const char*
const gidsIn = imports + gidsBeg;
6321 const char*
const valsIn = imports + valsBeg;
6323 size_t numBytesOut = 0;
6326 numBytesOut += PackTraits<LO>::unpackValue(numEntOut, numEntIn);
6327 if (
static_cast<size_t>(numEntOut) != numEnt ||
6328 numEntOut ==
static_cast<LO
>(0)) {
6329 const int myRank = this->getMap()->getComm()->getRank();
6330 std::ostringstream os;
6331 os <<
"(Proc " << myRank <<
") CrsMatrix::unpackRow: ";
6332 bool firstErrorCondition =
false;
6333 if (
static_cast<size_t>(numEntOut) != numEnt) {
6334 os <<
"Number of entries from numPacketsPerLID numEnt=" << numEnt
6335 <<
" does not equal number of entries unpacked from imports "
6337 << numEntOut <<
".";
6338 firstErrorCondition =
true;
6340 if (numEntOut ==
static_cast<LO
>(0)) {
6341 if (firstErrorCondition) {
6344 os <<
"Number of entries unpacked from imports buffer numEntOut=0, "
6345 "but number of bytes to unpack for this row numBytes="
6347 <<
" != 0. This should never happen, since packRow should only "
6348 "ever pack rows with a nonzero number of entries. In this case, "
6349 "the number of entries from numPacketsPerLID is numEnt="
6353 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error, os.str());
6357 Kokkos::pair<int, size_t> p;
6358 p = PackTraits<GO>::unpackArray(gidsOut, gidsIn, numEnt);
6359 errorCode += p.first;
6360 numBytesOut += p.second;
6362 p = PackTraits<ST>::unpackArray(valsOut, valsIn, numEnt);
6363 errorCode += p.first;
6364 numBytesOut += p.second;
6367 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != numBytes, std::logic_error,
"unpackRow: numBytesOut = " << numBytesOut <<
" != numBytes = " << numBytes <<
".");
6369 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6370 TEUCHOS_TEST_FOR_EXCEPTION(numBytesOut != expectedNumBytes, std::logic_error,
6373 << numBytesOut <<
" != expectedNumBytes = "
6374 << expectedNumBytes <<
".");
6376 TEUCHOS_TEST_FOR_EXCEPTION(errorCode != 0, std::runtime_error,
6378 "PackTraits::unpackArray returned a nonzero error code");
6383template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6386 size_t& totalNumEntries,
6387 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs)
const {
6392 typedef LocalOrdinal LO;
6393 typedef GlobalOrdinal GO;
6399 const bool verbose = Behavior::verbose(
"CrsMatrix");
6400 std::unique_ptr<std::string> prefix;
6402 prefix = this->
createPrefix(
"CrsMatrix",
"allocatePackSpaceNew");
6403 std::ostringstream os;
6404 os << *prefix <<
"Before:"
6412 std::cerr << os.str();
6417 const LO numExportLIDs =
static_cast<LO
>(exportLIDs.extent(0));
6419 TEUCHOS_ASSERT(!exportLIDs.need_sync_host());
6420 auto exportLIDs_h = exportLIDs.view_host();
6423 totalNumEntries = 0;
6424 for (LO i = 0; i < numExportLIDs; ++i) {
6425 const LO lclRow = exportLIDs_h[i];
6426 size_t curNumEntries = this->getNumEntriesInLocalRow(lclRow);
6429 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid()) {
6432 totalNumEntries += curNumEntries;
6443 const size_t allocSize =
6444 static_cast<size_t>(numExportLIDs) *
sizeof(LO) +
6445 totalNumEntries * (
sizeof(IST) +
sizeof(GO));
6446 if (
static_cast<size_t>(exports.extent(0)) < allocSize) {
6447 using exports_type = Kokkos::DualView<char*, buffer_device_type>;
6449 const std::string oldLabel = exports.view_device().label();
6450 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6451 exports = exports_type(newLabel, allocSize);
6455 std::ostringstream os;
6456 os << *prefix <<
"After:"
6464 std::cerr << os.str();
6468template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6470 packNew(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6471 Kokkos::DualView<char*, buffer_device_type>& exports,
6472 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
6473 size_t& constantNumPackets)
const {
6478 packCrsMatrixNew(*
this, exports, numPacketsPerLID, exportLIDs,
6479 constantNumPackets);
6481 this->packNonStaticNew(exportLIDs, exports, numPacketsPerLID,
6482 constantNumPackets);
6486template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6488 packNonStaticNew(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6489 Kokkos::DualView<char*, buffer_device_type>& exports,
6490 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
6491 size_t& constantNumPackets)
const {
6498 using LO = LocalOrdinal;
6499 using GO = GlobalOrdinal;
6501 const char tfecfFuncName[] =
"packNonStaticNew: ";
6503 const bool verbose = Behavior::verbose(
"CrsMatrix");
6504 std::unique_ptr<std::string> prefix;
6506 prefix = this->createPrefix(
"CrsMatrix",
"packNonStaticNew");
6507 std::ostringstream os;
6508 os << *prefix <<
"Start" << endl;
6509 std::cerr << os.str();
6512 const size_t numExportLIDs =
static_cast<size_t>(exportLIDs.extent(0));
6513 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numExportLIDs !=
static_cast<size_t>(numPacketsPerLID.extent(0)),
6514 std::invalid_argument,
"exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size() = " << numPacketsPerLID.extent(0) <<
".");
6519 constantNumPackets = 0;
6524 size_t totalNumEntries = 0;
6525 this->allocatePackSpaceNew(exports, totalNumEntries, exportLIDs);
6526 const size_t bufSize =
static_cast<size_t>(exports.extent(0));
6529 exports.clear_sync_state();
6530 exports.modify_host();
6531 auto exports_h = exports.view_host();
6533 std::ostringstream os;
6534 os << *prefix <<
"After marking exports as modified on host, "
6535 << dualViewStatusToString(exports,
"exports") << endl;
6536 std::cerr << os.str();
6540 auto exportLIDs_h = exportLIDs.view_host();
6543 const_cast<Kokkos::DualView<size_t*, buffer_device_type>*
>(&numPacketsPerLID)->clear_sync_state();
6544 const_cast<Kokkos::DualView<size_t*, buffer_device_type>*
>(&numPacketsPerLID)->modify_host();
6545 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
6550 auto maxRowNumEnt = this->getLocalMaxNumRowEntries();
6553 typename global_inds_host_view_type::non_const_type gidsIn_k;
6554 if (this->isLocallyIndexed()) {
6556 typename global_inds_host_view_type::non_const_type(
"packGids",
6561 for (
size_t i = 0; i < numExportLIDs; ++i) {
6562 const LO lclRow = exportLIDs_h[i];
6564 size_t numBytes = 0;
6565 size_t numEnt = this->getNumEntriesInLocalRow(lclRow);
6572 numPacketsPerLID_h[i] = 0;
6576 if (this->isLocallyIndexed()) {
6577 typename global_inds_host_view_type::non_const_type gidsIn;
6578 values_host_view_type valsIn;
6582 local_inds_host_view_type lidsIn;
6583 this->getLocalRowView(lclRow, lidsIn, valsIn);
6584 const map_type& colMap = *(this->getColMap());
6585 for (
size_t k = 0; k < numEnt; ++k) {
6586 gidsIn_k[k] = colMap.getGlobalElement(lidsIn[k]);
6588 gidsIn = Kokkos::subview(gidsIn_k, Kokkos::make_pair(GO(0), GO(numEnt)));
6590 const size_t numBytesPerValue =
6591 PackTraits<ST>::packValueCount(valsIn[0]);
6592 numBytes = this->
packRow(exports_h.data(), offset, numEnt,
6593 gidsIn.data(), valsIn.data(),
6595 }
else if (this->isGloballyIndexed()) {
6596 global_inds_host_view_type gidsIn;
6597 values_host_view_type valsIn;
6603 const map_type& rowMap = *(this->getRowMap());
6605 this->getGlobalRowView(gblRow, gidsIn, valsIn);
6607 const size_t numBytesPerValue =
6608 PackTraits<ST>::packValueCount(valsIn[0]);
6609 numBytes = this->
packRow(exports_h.data(), offset, numEnt,
6610 gidsIn.data(), valsIn.data(),
6617 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(offset > bufSize || offset + numBytes > bufSize, std::logic_error,
6618 "First invalid offset into 'exports' pack buffer at index i = " << i
6619 <<
". exportLIDs_h[i]: " << exportLIDs_h[i] <<
", bufSize: " << bufSize <<
", offset: " << offset <<
", numBytes: " << numBytes <<
".");
6623 numPacketsPerLID_h[i] = numBytes;
6628 std::ostringstream os;
6629 os << *prefix <<
"Tpetra::CrsMatrix::packNonStaticNew: After:" << endl
6636 std::cerr << os.str();
6640template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6644 const LocalOrdinal numEnt,
6646 const GlobalOrdinal cols[],
6648 const char*
const prefix,
6650 const bool verbose) {
6651 using GO = GlobalOrdinal;
6655 const GO gblRow = myGraph_->rowMap_->getGlobalElement(lclRow);
6656 Teuchos::ArrayView<const GO> cols_av(numEnt == 0 ?
nullptr : cols, numEnt);
6657 Teuchos::ArrayView<const Scalar> vals_av(numEnt == 0 ?
nullptr :
reinterpret_cast<const Scalar*
>(vals), numEnt);
6662 combineGlobalValues(gblRow, cols_av, vals_av, combMode,
6663 prefix, debug, verbose);
6667template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6670 const GlobalOrdinal globalRowIndex,
6671 const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
6672 const Teuchos::ArrayView<const Scalar>& values,
6674 const char*
const prefix,
6676 const bool verbose) {
6677 const char tfecfFuncName[] =
"combineGlobalValues: ";
6679 if (isStaticGraph()) {
6683 if (combineMode ==
ADD) {
6684 sumIntoGlobalValues(globalRowIndex, columnIndices, values);
6685 }
else if (combineMode ==
REPLACE) {
6686 replaceGlobalValues(globalRowIndex, columnIndices, values);
6687 }
else if (combineMode ==
ABSMAX) {
6688 using ::Tpetra::Details::AbsMax;
6690 this->
template transformGlobalValues<AbsMax<Scalar>>(globalRowIndex,
6693 }
else if (combineMode ==
INSERT) {
6694 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isStaticGraph() && combineMode ==
INSERT,
6695 std::invalid_argument,
6696 "INSERT combine mode is forbidden "
6697 "if the matrix has a static (const) graph (i.e., was "
6698 "constructed with the CrsMatrix constructor that takes a "
6699 "const CrsGraph pointer).");
6701 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error,
6702 "Invalid combine mode; should "
6704 "Please report this bug to the Tpetra developers.");
6707 if (combineMode ==
ADD || combineMode ==
INSERT) {
6714 insertGlobalValuesFilteredChecked(globalRowIndex,
6715 columnIndices, values, prefix, debug, verbose);
6726 else if (combineMode ==
ABSMAX) {
6727 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6728 !isStaticGraph() && combineMode ==
ABSMAX, std::logic_error,
6729 "ABSMAX combine mode when the matrix has a dynamic graph is not yet "
6731 }
else if (combineMode ==
REPLACE) {
6732 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6733 !isStaticGraph() && combineMode ==
REPLACE, std::logic_error,
6734 "REPLACE combine mode when the matrix has a dynamic graph is not yet "
6737 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6738 true, std::logic_error,
6739 "Should never get here! Please report this "
6740 "bug to the Tpetra developers.");
6745template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6747 unpackAndCombine(
const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& importLIDs,
6748 Kokkos::DualView<char*, buffer_device_type> imports,
6749 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6750 const size_t constantNumPackets,
6756 const char tfecfFuncName[] =
"unpackAndCombine: ";
6757 ProfilingRegion regionUAC(
"Tpetra::CrsMatrix::unpackAndCombine");
6759 const bool debug = Behavior::debug(
"CrsMatrix");
6760 const bool verbose = Behavior::verbose(
"CrsMatrix");
6761 constexpr int numValidModes = 5;
6764 const char* validModeNames[numValidModes] =
6765 {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT",
"ZERO"};
6767 std::unique_ptr<std::string> prefix;
6769 prefix = this->createPrefix(
"CrsMatrix",
"unpackAndCombine");
6770 std::ostringstream os;
6771 os << *prefix <<
"Start:" << endl
6773 << dualViewStatusToString(importLIDs,
"importLIDs")
6776 << dualViewStatusToString(imports,
"imports")
6779 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6781 << *prefix <<
" constantNumPackets: " << constantNumPackets
6785 std::cerr << os.str();
6789 if (std::find(validModes, validModes + numValidModes, combineMode) ==
6790 validModes + numValidModes) {
6791 std::ostringstream os;
6792 os <<
"Invalid combine mode. Valid modes are {";
6793 for (
int k = 0; k < numValidModes; ++k) {
6794 os << validModeNames[k];
6795 if (k < numValidModes - 1) {
6800 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::invalid_argument, os.str());
6802 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(importLIDs.extent(0) != numPacketsPerLID.extent(0),
6803 std::invalid_argument,
"importLIDs.extent(0)=" << importLIDs.extent(0) <<
" != numPacketsPerLID.extent(0)=" << numPacketsPerLID.extent(0) <<
".");
6806 if (combineMode ==
ZERO) {
6811 using Teuchos::reduceAll;
6812 std::unique_ptr<std::ostringstream> msg(
new std::ostringstream());
6815 unpackAndCombineImpl(importLIDs, imports, numPacketsPerLID,
6816 constantNumPackets, combineMode,
6818 }
catch (std::exception& e) {
6823 const Teuchos::Comm<int>& comm = *(this->
getComm());
6824 reduceAll<int, int>(comm, Teuchos::REDUCE_MAX,
6825 lclBad, Teuchos::outArg(gblBad));
6831 std::ostringstream os;
6832 os <<
"Proc " << comm.getRank() <<
": " << msg->str() << endl;
6833 msg = std::unique_ptr<std::ostringstream>(
new std::ostringstream());
6835 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::logic_error, std::endl
6836 <<
"unpackAndCombineImpl "
6837 "threw an exception on one or more participating processes: "
6842 unpackAndCombineImpl(importLIDs, imports, numPacketsPerLID,
6843 constantNumPackets, combineMode,
6848 std::ostringstream os;
6849 os << *prefix <<
"Done!" << endl
6851 << dualViewStatusToString(importLIDs,
"importLIDs")
6854 << dualViewStatusToString(imports,
"imports")
6857 << dualViewStatusToString(numPacketsPerLID,
"numPacketsPerLID")
6859 std::cerr << os.str();
6863template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6868 Kokkos::DualView<char*, buffer_device_type> imports,
6869 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6870 const size_t constantNumPackets,
6872 const bool verbose) {
6874 "Tpetra::CrsMatrix::unpackAndCombineImpl",
6877 const char tfecfFuncName[] =
"unpackAndCombineImpl";
6878 std::unique_ptr<std::string> prefix;
6880 prefix = this->createPrefix(
"CrsMatrix", tfecfFuncName);
6881 std::ostringstream os;
6882 os << *prefix <<
"isStaticGraph(): "
6883 << (isStaticGraph() ?
"true" :
"false")
6884 <<
", importLIDs.extent(0): "
6885 << importLIDs.extent(0)
6886 <<
", imports.extent(0): "
6887 << imports.extent(0)
6888 <<
", numPacketsPerLID.extent(0): "
6889 << numPacketsPerLID.extent(0)
6891 std::cerr << os.str();
6894 if (isStaticGraph()) {
6895 using Details::unpackCrsMatrixAndCombineNew;
6896 unpackCrsMatrixAndCombineNew(*
this, imports, numPacketsPerLID,
6897 importLIDs, constantNumPackets,
6901 using padding_type =
typename crs_graph_type::padding_type;
6902 std::unique_ptr<padding_type> padding;
6904 padding = myGraph_->computePaddingForCrsMatrixUnpack(
6905 importLIDs, imports, numPacketsPerLID, verbose);
6906 }
catch (std::exception& e) {
6907 const auto rowMap = getRowMap();
6908 const auto comm = rowMap.is_null() ? Teuchos::null : rowMap->
getComm();
6909 const int myRank = comm.is_null() ? -1 : comm->getRank();
6910 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"Proc " << myRank <<
": "
6911 "Tpetra::CrsGraph::computePaddingForCrsMatrixUnpack "
6912 "threw an exception: "
6916 std::ostringstream os;
6917 os << *prefix <<
"Call applyCrsPadding" << endl;
6918 std::cerr << os.str();
6920 applyCrsPadding(*padding, verbose);
6923 std::ostringstream os;
6924 os << *prefix <<
"Call unpackAndCombineImplNonStatic" << endl;
6925 std::cerr << os.str();
6927 unpackAndCombineImplNonStatic(importLIDs, imports,
6934 std::ostringstream os;
6935 os << *prefix <<
"Done" << endl;
6936 std::cerr << os.str();
6940template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6945 Kokkos::DualView<char*, buffer_device_type> imports,
6946 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6947 const size_t constantNumPackets,
6954 using Kokkos::MemoryUnmanaged;
6955 using Kokkos::subview;
6958 using LO = LocalOrdinal;
6959 using GO = GlobalOrdinal;
6961 using size_type =
typename Teuchos::ArrayView<LO>::size_type;
6963 typename View<int*, device_type>::host_mirror_type::execution_space;
6964 using pair_type = std::pair<typename View<int*, HES>::size_type,
6965 typename View<int*, HES>::size_type>;
6966 using gids_out_type = View<GO*, HES, MemoryUnmanaged>;
6967 using vals_out_type = View<ST*, HES, MemoryUnmanaged>;
6968 const char tfecfFuncName[] =
"unpackAndCombineImplNonStatic";
6970 const bool debug = Behavior::debug(
"CrsMatrix");
6971 const bool verbose = Behavior::verbose(
"CrsMatrix");
6972 std::unique_ptr<std::string> prefix;
6974 prefix = this->
createPrefix(
"CrsMatrix", tfecfFuncName);
6975 std::ostringstream os;
6976 os << *prefix << endl;
6977 std::cerr << os.str();
6979 const char*
const prefix_raw =
6980 verbose ? prefix.get()->c_str() :
nullptr;
6982 const size_type numImportLIDs = importLIDs.extent(0);
6983 if (combineMode ==
ZERO || numImportLIDs == 0) {
6988 "Tpetra::CrsMatrix::unpackAndCombineImplNonStatic",
6992 if (imports.need_sync_host()) {
6993 imports.sync_host();
6995 auto imports_h = imports.view_host();
6998 if (numPacketsPerLID.need_sync_host()) {
6999 numPacketsPerLID.sync_host();
7001 auto numPacketsPerLID_h = numPacketsPerLID.view_host();
7003 TEUCHOS_ASSERT(!importLIDs.need_sync_host());
7004 auto importLIDs_h = importLIDs.view_host();
7006 size_t numBytesPerValue;
7017 numBytesPerValue = PackTraits<ST>::packValueCount(val);
7022 size_t maxRowNumEnt = 0;
7023 for (size_type i = 0; i < numImportLIDs; ++i) {
7024 const size_t numBytes = numPacketsPerLID_h[i];
7025 if (numBytes == 0) {
7030 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(offset + numBytes >
size_t(imports_h.extent(0)),
7031 std::logic_error,
": At local row index importLIDs_h[i=" << i <<
"]=" << importLIDs_h[i] <<
", offset (=" << offset <<
") + numBytes (=" << numBytes <<
") > "
7032 "imports_h.extent(0)="
7033 << imports_h.extent(0) <<
".");
7038 const size_t theNumBytes =
7039 PackTraits<LO>::packValueCount(numEntLO);
7040 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(theNumBytes > numBytes, std::logic_error,
": theNumBytes=" << theNumBytes <<
" > numBytes = " << numBytes <<
".");
7042 const char*
const inBuf = imports_h.data() + offset;
7043 const size_t actualNumBytes =
7044 PackTraits<LO>::unpackValue(numEntLO, inBuf);
7047 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(actualNumBytes > numBytes, std::logic_error,
": At i=" << i <<
", actualNumBytes=" << actualNumBytes <<
" > numBytes=" << numBytes <<
".");
7048 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numEntLO == 0, std::logic_error,
7049 ": At local row index "
7051 << i <<
"]=" << importLIDs_h[i] <<
", "
7052 "the number of entries read from the packed data is "
7054 << numEntLO <<
", but numBytes=" << numBytes
7058 maxRowNumEnt = std::max(
size_t(numEntLO), maxRowNumEnt);
7066 View<GO*, HES> gblColInds;
7067 View<LO*, HES> lclColInds;
7068 View<ST*, HES> vals;
7081 gblColInds = ScalarViewTraits<GO, HES>::allocateArray(
7082 gid, maxRowNumEnt,
"gids");
7083 lclColInds = ScalarViewTraits<LO, HES>::allocateArray(
7084 lid, maxRowNumEnt,
"lids");
7085 vals = ScalarViewTraits<ST, HES>::allocateArray(
7086 val, maxRowNumEnt,
"vals");
7090 for (size_type i = 0; i < numImportLIDs; ++i) {
7091 const size_t numBytes = numPacketsPerLID_h[i];
7092 if (numBytes == 0) {
7096 const char*
const inBuf = imports_h.data() + offset;
7097 (void)PackTraits<LO>::unpackValue(numEntLO, inBuf);
7099 const size_t numEnt =
static_cast<size_t>(numEntLO);
7101 const LO lclRow = importLIDs_h[i];
7103 gids_out_type gidsOut = subview(gblColInds, pair_type(0, numEnt));
7104 vals_out_type valsOut = subview(vals, pair_type(0, numEnt));
7106 const size_t numBytesOut =
7107 unpackRow(gidsOut.data(), valsOut.data(), imports_h.data(),
7108 offset, numBytes, numEnt, numBytesPerValue);
7109 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(numBytes != numBytesOut, std::logic_error,
": At i=" << i <<
", numBytes=" << numBytes <<
" != numBytesOut=" << numBytesOut <<
".");
7111 const ST*
const valsRaw =
const_cast<const ST*
>(valsOut.data());
7112 const GO*
const gidsRaw =
const_cast<const GO*
>(gidsOut.data());
7113 combineGlobalValuesRaw(lclRow, numEnt, valsRaw, gidsRaw,
7114 combineMode, prefix_raw, debug, verbose);
7120 std::ostringstream os;
7121 os << *prefix <<
"Done" << endl;
7122 std::cerr << os.str();
7126template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7127Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7130 const bool force)
const {
7131 using Teuchos::null;
7135 TEUCHOS_TEST_FOR_EXCEPTION(
7137 "Tpetra::CrsMatrix::getColumn"
7138 "MapMultiVector: You may only call this method if the matrix has a "
7139 "column Map. If the matrix does not yet have a column Map, you should "
7140 "first call fillComplete (with domain and range Map if necessary).");
7144 TEUCHOS_TEST_FOR_EXCEPTION(
7147 "CrsMatrix::getColumnMapMultiVector: You may only call this method if "
7148 "this matrix's graph is fill complete.");
7151 RCP<const import_type> importer = this->
getGraph()->getImporter();
7152 RCP<const map_type> colMap = this->
getColMap();
7165 if (!importer.is_null() || force) {
7167 X_colMap = rcp(
new MV(colMap, numVecs));
7183template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7184Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7187 const bool force)
const {
7188 using Teuchos::null;
7194 TEUCHOS_TEST_FOR_EXCEPTION(
7197 "CrsMatrix::getRowMapMultiVector: You may only call this method if this "
7198 "matrix's graph is fill complete.");
7201 RCP<const export_type> exporter = this->
getGraph()->getExporter();
7205 RCP<const map_type> rowMap = this->
getRowMap();
7217 if (!exporter.is_null() || force) {
7219 Y_rowMap = rcp(
new MV(rowMap, numVecs));
7228template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7231 TEUCHOS_TEST_FOR_EXCEPTION(
7232 myGraph_.is_null(), std::logic_error,
7233 "Tpetra::CrsMatrix::"
7234 "removeEmptyProcessesInPlace: This method does not work when the matrix "
7235 "was created with a constant graph (that is, when it was created using "
7236 "the version of its constructor that takes an RCP<const CrsGraph>). "
7237 "This is because the matrix is not allowed to modify the graph in that "
7238 "case, but removing empty processes requires modifying the graph.");
7239 myGraph_->removeEmptyProcessesInPlace(newMap);
7247 staticGraph_ = Teuchos::rcp_const_cast<const Graph>(myGraph_);
7250template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7251Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
7253 add(
const Scalar& alpha,
7256 const Teuchos::RCP<const map_type>& domainMap,
7257 const Teuchos::RCP<const map_type>& rangeMap,
7258 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
7260 using Teuchos::Array;
7261 using Teuchos::ArrayView;
7262 using Teuchos::ParameterList;
7265 using Teuchos::rcp_implicit_cast;
7266 using Teuchos::sublist;
7269 using crs_matrix_type =
7270 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
7271 const char errPfx[] =
"Tpetra::CrsMatrix::add: ";
7275 std::unique_ptr<std::string> prefix;
7277 prefix = this->createPrefix(
"CrsMatrix",
"add");
7278 std::ostringstream os;
7279 os << *prefix <<
"Start" << endl;
7280 std::cerr << os.str();
7283 const crs_matrix_type& B = *
this;
7284 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero();
7285 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one();
7293 RCP<const map_type> B_domainMap = B.getDomainMap();
7294 RCP<const map_type> B_rangeMap = B.getRangeMap();
7296 RCP<const map_type> theDomainMap = domainMap;
7297 RCP<const map_type> theRangeMap = rangeMap;
7299 if (domainMap.is_null()) {
7300 if (B_domainMap.is_null()) {
7301 TEUCHOS_TEST_FOR_EXCEPTION(
7302 A_domainMap.is_null(), std::invalid_argument,
7303 "Tpetra::CrsMatrix::add: If neither A nor B have a domain Map, "
7304 "then you must supply a nonnull domain Map to this method.");
7305 theDomainMap = A_domainMap;
7307 theDomainMap = B_domainMap;
7310 if (rangeMap.is_null()) {
7311 if (B_rangeMap.is_null()) {
7312 TEUCHOS_TEST_FOR_EXCEPTION(
7313 A_rangeMap.is_null(), std::invalid_argument,
7314 "Tpetra::CrsMatrix::add: If neither A nor B have a range Map, "
7315 "then you must supply a nonnull range Map to this method.");
7316 theRangeMap = A_rangeMap;
7318 theRangeMap = B_rangeMap;
7326 if (!A_domainMap.is_null() && !A_rangeMap.is_null()) {
7327 if (!B_domainMap.is_null() && !B_rangeMap.is_null()) {
7328 TEUCHOS_TEST_FOR_EXCEPTION(!B_domainMap->isSameAs(*A_domainMap),
7329 std::invalid_argument,
7330 errPfx <<
"The input RowMatrix A must have a domain Map "
7331 "which is the same as (isSameAs) this RowMatrix's "
7333 TEUCHOS_TEST_FOR_EXCEPTION(!B_rangeMap->isSameAs(*A_rangeMap), std::invalid_argument,
7334 errPfx <<
"The input RowMatrix A must have a range Map "
7335 "which is the same as (isSameAs) this RowMatrix's range "
7337 TEUCHOS_TEST_FOR_EXCEPTION(!domainMap.is_null() &&
7338 !domainMap->isSameAs(*B_domainMap),
7339 std::invalid_argument,
7340 errPfx <<
"The input domain Map must be the same as "
7341 "(isSameAs) this RowMatrix's domain Map.");
7342 TEUCHOS_TEST_FOR_EXCEPTION(!rangeMap.is_null() &&
7343 !rangeMap->isSameAs(*B_rangeMap),
7344 std::invalid_argument,
7345 errPfx <<
"The input range Map must be the same as "
7346 "(isSameAs) this RowMatrix's range Map.");
7348 }
else if (!B_domainMap.is_null() && !B_rangeMap.is_null()) {
7349 TEUCHOS_TEST_FOR_EXCEPTION(!domainMap.is_null() &&
7350 !domainMap->isSameAs(*B_domainMap),
7351 std::invalid_argument,
7352 errPfx <<
"The input domain Map must be the same as "
7353 "(isSameAs) this RowMatrix's domain Map.");
7354 TEUCHOS_TEST_FOR_EXCEPTION(!rangeMap.is_null() && !rangeMap->isSameAs(*B_rangeMap),
7355 std::invalid_argument,
7356 errPfx <<
"The input range Map must be the same as "
7357 "(isSameAs) this RowMatrix's range Map.");
7359 TEUCHOS_TEST_FOR_EXCEPTION(domainMap.is_null() || rangeMap.is_null(),
7360 std::invalid_argument, errPfx <<
"If neither A nor B "
7361 "have a domain and range Map, then you must supply a "
7362 "nonnull domain and range Map to this method.");
7369 bool callFillComplete =
true;
7370 RCP<ParameterList> constructorSublist;
7371 RCP<ParameterList> fillCompleteSublist;
7372 if (!params.is_null()) {
7374 params->get(
"Call fillComplete", callFillComplete);
7375 constructorSublist = sublist(params,
"Constructor parameters");
7376 fillCompleteSublist = sublist(params,
"fillComplete parameters");
7379 RCP<const map_type> A_rowMap = A.
getRowMap();
7380 RCP<const map_type> B_rowMap = B.getRowMap();
7381 RCP<const map_type> C_rowMap = B_rowMap;
7382 RCP<crs_matrix_type> C;
7388 if (A_rowMap->isSameAs(*B_rowMap)) {
7389 const LO localNumRows =
static_cast<LO
>(A_rowMap->getLocalNumElements());
7390 Array<size_t> C_maxNumEntriesPerRow(localNumRows, 0);
7393 if (alpha !=
ZERO) {
7394 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
7396 C_maxNumEntriesPerRow[localRow] += A_numEntries;
7401 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
7402 const size_t B_numEntries = B.getNumEntriesInLocalRow(localRow);
7403 C_maxNumEntriesPerRow[localRow] += B_numEntries;
7407 if (constructorSublist.is_null()) {
7408 C = rcp(
new crs_matrix_type(C_rowMap, C_maxNumEntriesPerRow()));
7410 C = rcp(
new crs_matrix_type(C_rowMap, C_maxNumEntriesPerRow(),
7411 constructorSublist));
7421 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument, errPfx <<
"The row maps must "
7422 "be the same for statically allocated matrices, to ensure "
7423 "that there is sufficient space to do the addition.");
7426 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null(), std::logic_error,
7427 errPfx <<
"C should not be null at this point. "
7428 "Please report this bug to the Tpetra developers.");
7431 std::ostringstream os;
7432 os << *prefix <<
"Compute C = alpha*A + beta*B" << endl;
7433 std::cerr << os.str();
7435 using gids_type = nonconst_global_inds_host_view_type;
7436 using vals_type = nonconst_values_host_view_type;
7440 if (alpha !=
ZERO) {
7441 const LO A_localNumRows =
static_cast<LO
>(A_rowMap->getLocalNumElements());
7442 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
7444 const GO globalRow = A_rowMap->getGlobalElement(localRow);
7445 if (A_numEntries >
static_cast<size_t>(ind.size())) {
7446 Kokkos::resize(ind, A_numEntries);
7447 Kokkos::resize(val, A_numEntries);
7449 gids_type indView = Kokkos::subview(ind, std::make_pair((
size_t)0, A_numEntries));
7450 vals_type valView = Kokkos::subview(val, std::make_pair((
size_t)0, A_numEntries));
7454 for (
size_t k = 0; k < A_numEntries; ++k) {
7455 valView[k] *= alpha;
7458 C->insertGlobalValues(globalRow, A_numEntries,
7459 reinterpret_cast<Scalar*
>(valView.data()),
7465 const LO B_localNumRows =
static_cast<LO
>(B_rowMap->getLocalNumElements());
7466 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
7467 size_t B_numEntries = B.getNumEntriesInLocalRow(localRow);
7468 const GO globalRow = B_rowMap->getGlobalElement(localRow);
7469 if (B_numEntries >
static_cast<size_t>(ind.size())) {
7470 Kokkos::resize(ind, B_numEntries);
7471 Kokkos::resize(val, B_numEntries);
7473 gids_type indView = Kokkos::subview(ind, std::make_pair((
size_t)0, B_numEntries));
7474 vals_type valView = Kokkos::subview(val, std::make_pair((
size_t)0, B_numEntries));
7475 B.getGlobalRowCopy(globalRow, indView, valView, B_numEntries);
7478 for (
size_t k = 0; k < B_numEntries; ++k) {
7482 C->insertGlobalValues(globalRow, B_numEntries,
7483 reinterpret_cast<Scalar*
>(valView.data()),
7488 if (callFillComplete) {
7490 std::ostringstream os;
7491 os << *prefix <<
"Call fillComplete on C" << endl;
7492 std::cerr << os.str();
7494 if (fillCompleteSublist.is_null()) {
7495 C->fillComplete(theDomainMap, theRangeMap);
7497 C->fillComplete(theDomainMap, theRangeMap, fillCompleteSublist);
7499 }
else if (verbose) {
7500 std::ostringstream os;
7501 os << *prefix <<
"Do NOT call fillComplete on C" << endl;
7502 std::cerr << os.str();
7506 std::ostringstream os;
7507 os << *prefix <<
"Done" << endl;
7508 std::cerr << os.str();
7510 return rcp_implicit_cast<row_matrix_type>(C);
7513template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7516 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
7517 const Teuchos::RCP<const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>>& domainTransfer,
7518 const Teuchos::RCP<const map_type>& domainMap,
7519 const Teuchos::RCP<const map_type>& rangeMap,
7520 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
7527 using Teuchos::ArrayRCP;
7528 using Teuchos::ArrayView;
7529 using Teuchos::Comm;
7530 using Teuchos::ParameterList;
7532 typedef LocalOrdinal LO;
7533 typedef GlobalOrdinal GO;
7539 const bool debug = Behavior::debug(
"CrsMatrix");
7540 const bool verbose = Behavior::verbose(
"CrsMatrix");
7541 int MyPID = getComm()->getRank();
7543 std::unique_ptr<std::string> verbosePrefix;
7546 this->createPrefix(
"CrsMatrix",
"transferAndFillComplete");
7547 std::ostringstream os;
7548 os <<
"Start" << endl;
7549 std::cerr << os.str();
7556 bool reverseMode =
false;
7557 bool restrictComm =
false;
7559 int mm_optimization_core_count =
7560 Behavior::TAFC_OptimizationCoreCount();
7561 RCP<ParameterList> matrixparams;
7562 bool overrideAllreduce =
false;
7563 bool useKokkosPath =
false;
7564 if (!params.is_null()) {
7565 matrixparams = sublist(params,
"CrsMatrix");
7566 reverseMode = params->get(
"Reverse Mode", reverseMode);
7567 useKokkosPath = params->get(
"TAFC: use kokkos path", useKokkosPath);
7568 restrictComm = params->get(
"Restrict Communicator", restrictComm);
7569 auto& slist = params->sublist(
"matrixmatrix: kernel params",
false);
7570 isMM = slist.get(
"isMatrixMatrix_TransferAndFillComplete",
false);
7571 mm_optimization_core_count = slist.get(
"MM_TAFC_OptimizationCoreCount", mm_optimization_core_count);
7573 overrideAllreduce = slist.get(
"MM_TAFC_OverrideAllreduceCheck",
false);
7574 if (getComm()->getSize() < mm_optimization_core_count && isMM) isMM =
false;
7575 if (reverseMode) isMM =
false;
7579 std::shared_ptr<::Tpetra::Details::CommRequest> iallreduceRequest;
7581 int reduced_mismatch = 0;
7582 if (isMM && !overrideAllreduce) {
7584 const bool source_vals = !getGraph()->getImporter().is_null();
7585 const bool target_vals = !(rowTransfer.getExportLIDs().size() == 0 ||
7586 rowTransfer.getRemoteLIDs().size() == 0);
7587 mismatch = (source_vals != target_vals) ? 1 : 0;
7589 ::Tpetra::Details::iallreduce(mismatch, reduced_mismatch,
7590 Teuchos::REDUCE_MAX, *(getComm()));
7593 std::string tlstr = isMM ?
"Tpetra: TAFC All:MMOpt" :
"Tpetra: TAFC All:MMLegacy";
7594 Tpetra::Details::ProfilingRegion MMall(tlstr.c_str());
7603 TEUCHOS_TEST_FOR_EXCEPTION(
7604 xferAsImport ==
nullptr && xferAsExport ==
nullptr, std::invalid_argument,
7605 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' input "
7606 "argument must be either an Import or an Export, and its template "
7607 "parameters must match the corresponding template parameters of the "
7615 Teuchos::RCP<const import_type> xferDomainAsImport = Teuchos::rcp_dynamic_cast<const import_type>(domainTransfer);
7616 Teuchos::RCP<const export_type> xferDomainAsExport = Teuchos::rcp_dynamic_cast<const export_type>(domainTransfer);
7618 if (!domainTransfer.is_null()) {
7619 TEUCHOS_TEST_FOR_EXCEPTION(
7620 (xferDomainAsImport.is_null() && xferDomainAsExport.is_null()), std::invalid_argument,
7621 "Tpetra::CrsMatrix::transferAndFillComplete: The 'domainTransfer' input "
7622 "argument must be either an Import or an Export, and its template "
7623 "parameters must match the corresponding template parameters of the "
7626 TEUCHOS_TEST_FOR_EXCEPTION(
7627 (xferAsImport !=
nullptr || !xferDomainAsImport.is_null()) &&
7628 ((xferAsImport !=
nullptr && xferDomainAsImport.is_null()) ||
7629 (xferAsImport ==
nullptr && !xferDomainAsImport.is_null())),
7630 std::invalid_argument,
7631 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
7632 "arguments must be of the same type (either Import or Export).");
7634 TEUCHOS_TEST_FOR_EXCEPTION(
7635 (xferAsExport !=
nullptr || !xferDomainAsExport.is_null()) &&
7636 ((xferAsExport !=
nullptr && xferDomainAsExport.is_null()) ||
7637 (xferAsExport ==
nullptr && !xferDomainAsExport.is_null())),
7638 std::invalid_argument,
7639 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
7640 "arguments must be of the same type (either Import or Export).");
7645 const bool communication_needed = rowTransfer.getSourceMap()->isDistributed();
7649 RCP<const map_type> MyRowMap = reverseMode ? rowTransfer.getSourceMap() : rowTransfer.getTargetMap();
7650 RCP<const map_type> MyColMap;
7651 RCP<const map_type> MyDomainMap = !domainMap.is_null() ? domainMap : getDomainMap();
7652 RCP<const map_type> MyRangeMap = !rangeMap.is_null() ? rangeMap : getRangeMap();
7653 RCP<const map_type> BaseRowMap = MyRowMap;
7654 RCP<const map_type> BaseDomainMap = MyDomainMap;
7662 if (!destMat.is_null()) {
7673 const bool NewFlag = !destMat->getGraph()->isLocallyIndexed() &&
7674 !destMat->getGraph()->isGloballyIndexed();
7675 TEUCHOS_TEST_FOR_EXCEPTION(
7676 !NewFlag, std::invalid_argument,
7677 "Tpetra::CrsMatrix::"
7678 "transferAndFillComplete: The input argument 'destMat' is only allowed "
7679 "to be nonnull, if its graph is empty (neither locally nor globally "
7688 TEUCHOS_TEST_FOR_EXCEPTION(
7689 !destMat->getRowMap()->isSameAs(*MyRowMap), std::invalid_argument,
7690 "Tpetra::CrsMatrix::transferAndFillComplete: The (row) Map of the "
7691 "input argument 'destMat' is not the same as the (row) Map specified "
7692 "by the input argument 'rowTransfer'.");
7693 TEUCHOS_TEST_FOR_EXCEPTION(
7694 !destMat->checkSizes(*
this), std::invalid_argument,
7695 "Tpetra::CrsMatrix::transferAndFillComplete: You provided a nonnull "
7696 "destination matrix, but checkSizes() indicates that it is not a legal "
7697 "legal target for redistribution from the source matrix (*this). This "
7698 "may mean that they do not have the same dimensions.");
7712 TEUCHOS_TEST_FOR_EXCEPTION(
7713 !(reverseMode || getRowMap()->isSameAs(*rowTransfer.getSourceMap())),
7714 std::invalid_argument,
7715 "Tpetra::CrsMatrix::transferAndFillComplete: "
7716 "rowTransfer->getSourceMap() must match this->getRowMap() in forward mode.");
7717 TEUCHOS_TEST_FOR_EXCEPTION(
7718 !(!reverseMode || getRowMap()->isSameAs(*rowTransfer.getTargetMap())),
7719 std::invalid_argument,
7720 "Tpetra::CrsMatrix::transferAndFillComplete: "
7721 "rowTransfer->getTargetMap() must match this->getRowMap() in reverse mode.");
7724 TEUCHOS_TEST_FOR_EXCEPTION(
7725 !xferDomainAsImport.is_null() && !xferDomainAsImport->getTargetMap()->isSameAs(*domainMap),
7726 std::invalid_argument,
7727 "Tpetra::CrsMatrix::transferAndFillComplete: The target map of the 'domainTransfer' input "
7728 "argument must be the same as the rebalanced domain map 'domainMap'");
7730 TEUCHOS_TEST_FOR_EXCEPTION(
7731 !xferDomainAsExport.is_null() && !xferDomainAsExport->getSourceMap()->isSameAs(*domainMap),
7732 std::invalid_argument,
7733 "Tpetra::CrsMatrix::transferAndFillComplete: The source map of the 'domainTransfer' input "
7734 "argument must be the same as the rebalanced domain map 'domainMap'");
7747 const size_t NumSameIDs = rowTransfer.getNumSameIDs();
7748 ArrayView<const LO> ExportLIDs = reverseMode ? rowTransfer.getRemoteLIDs() : rowTransfer.getExportLIDs();
7749 auto RemoteLIDs = reverseMode ? rowTransfer.getExportLIDs_dv() : rowTransfer.getRemoteLIDs_dv();
7750 auto PermuteToLIDs = reverseMode ? rowTransfer.getPermuteFromLIDs_dv() : rowTransfer.getPermuteToLIDs_dv();
7751 auto PermuteFromLIDs = reverseMode ? rowTransfer.getPermuteToLIDs_dv() : rowTransfer.getPermuteFromLIDs_dv();
7752 Distributor& Distor = rowTransfer.getDistributor();
7755 Teuchos::Array<int> SourcePids;
7758 RCP<const map_type> ReducedRowMap, ReducedColMap,
7759 ReducedDomainMap, ReducedRangeMap;
7760 RCP<const Comm<int>> ReducedComm;
7764 if (destMat.is_null()) {
7765 destMat = rcp(
new this_CRS_type(MyRowMap, 0, matrixparams));
7772 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC restrictComm");
7773 ReducedRowMap = MyRowMap->removeEmptyProcesses();
7774 ReducedComm = ReducedRowMap.is_null() ? Teuchos::null : ReducedRowMap->getComm();
7775 destMat->removeEmptyProcessesInPlace(ReducedRowMap);
7777 ReducedDomainMap = MyRowMap.getRawPtr() == MyDomainMap.getRawPtr() ? ReducedRowMap : MyDomainMap->replaceCommWithSubset(ReducedComm);
7778 ReducedRangeMap = MyRowMap.getRawPtr() == MyRangeMap.getRawPtr() ? ReducedRowMap : MyRangeMap->replaceCommWithSubset(ReducedComm);
7781 MyRowMap = ReducedRowMap;
7782 MyDomainMap = ReducedDomainMap;
7783 MyRangeMap = ReducedRangeMap;
7786 if (!ReducedComm.is_null()) {
7787 MyPID = ReducedComm->getRank();
7792 ReducedComm = MyRowMap->getComm();
7799 RCP<const import_type> MyImporter = getGraph()->getImporter();
7802 bool bSameDomainMap = BaseDomainMap->isSameAs(*getDomainMap());
7804 if (!restrictComm && !MyImporter.is_null() && bSameDomainMap) {
7805 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs same map");
7813 }
else if (restrictComm && !MyImporter.is_null() && bSameDomainMap) {
7816 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs restricted comm");
7817 IntVectorType SourceDomain_pids(getDomainMap(),
true);
7818 IntVectorType SourceCol_pids(getColMap());
7820 SourceDomain_pids.putScalar(MyPID);
7822 SourceCol_pids.doImport(SourceDomain_pids, *MyImporter,
INSERT);
7823 SourcePids.resize(getColMap()->getLocalNumElements());
7824 SourceCol_pids.get1dCopy(SourcePids());
7825 }
else if (MyImporter.is_null()) {
7827 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs all local entries");
7828 SourcePids.resize(getColMap()->getLocalNumElements());
7829 SourcePids.assign(getColMap()->getLocalNumElements(), MyPID);
7830 }
else if (!MyImporter.is_null() &&
7831 !domainTransfer.is_null()) {
7836 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs rectangular case");
7839 IntVectorType TargetDomain_pids(domainMap);
7840 TargetDomain_pids.putScalar(MyPID);
7843 IntVectorType SourceDomain_pids(getDomainMap());
7846 IntVectorType SourceCol_pids(getColMap());
7848 if (!reverseMode && !xferDomainAsImport.is_null()) {
7849 SourceDomain_pids.doExport(TargetDomain_pids, *xferDomainAsImport,
INSERT);
7850 }
else if (reverseMode && !xferDomainAsExport.is_null()) {
7851 SourceDomain_pids.doExport(TargetDomain_pids, *xferDomainAsExport,
INSERT);
7852 }
else if (!reverseMode && !xferDomainAsExport.is_null()) {
7853 SourceDomain_pids.doImport(TargetDomain_pids, *xferDomainAsExport,
INSERT);
7854 }
else if (reverseMode && !xferDomainAsImport.is_null()) {
7855 SourceDomain_pids.doImport(TargetDomain_pids, *xferDomainAsImport,
INSERT);
7857 TEUCHOS_TEST_FOR_EXCEPTION(
7858 true, std::logic_error,
7859 "Tpetra::CrsMatrix::"
7860 "transferAndFillComplete: Should never get here! "
7861 "Please report this bug to a Tpetra developer.");
7863 SourceCol_pids.doImport(SourceDomain_pids, *MyImporter,
INSERT);
7864 SourcePids.resize(getColMap()->getLocalNumElements());
7865 SourceCol_pids.get1dCopy(SourcePids());
7866 }
else if (!MyImporter.is_null() &&
7867 BaseDomainMap->isSameAs(*BaseRowMap) &&
7868 getDomainMap()->isSameAs(*getRowMap())) {
7870 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs query import");
7872 IntVectorType TargetRow_pids(domainMap);
7873 IntVectorType SourceRow_pids(getRowMap());
7874 IntVectorType SourceCol_pids(getColMap());
7876 TargetRow_pids.putScalar(MyPID);
7877 if (!reverseMode && xferAsImport !=
nullptr) {
7878 SourceRow_pids.doExport(TargetRow_pids, *xferAsImport,
INSERT);
7879 }
else if (reverseMode && xferAsExport !=
nullptr) {
7880 SourceRow_pids.doExport(TargetRow_pids, *xferAsExport,
INSERT);
7881 }
else if (!reverseMode && xferAsExport !=
nullptr) {
7882 SourceRow_pids.doImport(TargetRow_pids, *xferAsExport,
INSERT);
7883 }
else if (reverseMode && xferAsImport !=
nullptr) {
7884 SourceRow_pids.doImport(TargetRow_pids, *xferAsImport,
INSERT);
7886 TEUCHOS_TEST_FOR_EXCEPTION(
7887 true, std::logic_error,
7888 "Tpetra::CrsMatrix::"
7889 "transferAndFillComplete: Should never get here! "
7890 "Please report this bug to a Tpetra developer.");
7893 SourceCol_pids.doImport(SourceRow_pids, *MyImporter,
INSERT);
7894 SourcePids.resize(getColMap()->getLocalNumElements());
7895 SourceCol_pids.get1dCopy(SourcePids());
7897 TEUCHOS_TEST_FOR_EXCEPTION(
7898 true, std::invalid_argument,
7899 "Tpetra::CrsMatrix::"
7900 "transferAndFillComplete: This method only allows either domainMap == "
7901 "getDomainMap (), or (domainMap == rowTransfer.getTargetMap () and "
7902 "getDomainMap () == getRowMap ()).");
7906 size_t constantNumPackets = destMat->constantNumberOfPackets();
7908 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC reallocate buffers");
7909 if (constantNumPackets == 0) {
7910 destMat->reallocArraysForNumPacketsPerLid(ExportLIDs.size(),
7911 RemoteLIDs.view_host().size());
7917 const size_t rbufLen = RemoteLIDs.view_host().size() * constantNumPackets;
7918 destMat->reallocImportsIfNeeded(rbufLen,
false,
nullptr);
7924 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC pack and prepare");
7928 using Teuchos::outArg;
7929 using Teuchos::REDUCE_MAX;
7930 using Teuchos::reduceAll;
7931 RCP<const Teuchos::Comm<int>> comm = this->getComm();
7932 const int myRank = comm->getRank();
7934 std::ostringstream errStrm;
7938 Teuchos::ArrayView<size_t> numExportPacketsPerLID;
7941 destMat->numExportPacketsPerLID_.modify_host();
7942 numExportPacketsPerLID =
7944 }
catch (std::exception& e) {
7945 errStrm <<
"Proc " << myRank <<
": getArrayViewFromDualView threw: "
7946 << e.what() << std::endl;
7949 errStrm <<
"Proc " << myRank <<
": getArrayViewFromDualView threw "
7950 "an exception not a subclass of std::exception"
7955 if (!comm.is_null()) {
7956 reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
7960 TEUCHOS_TEST_FOR_EXCEPTION(
7961 true, std::runtime_error,
7962 "getArrayViewFromDualView threw an "
7963 "exception on at least one process.");
7967 std::ostringstream os;
7968 os << *verbosePrefix <<
"Calling packCrsMatrixWithOwningPIDs"
7970 std::cerr << os.str();
7975 numExportPacketsPerLID,
7978 constantNumPackets);
7979 }
catch (std::exception& e) {
7980 errStrm <<
"Proc " << myRank <<
": packCrsMatrixWithOwningPIDs threw: "
7981 << e.what() << std::endl;
7984 errStrm <<
"Proc " << myRank <<
": packCrsMatrixWithOwningPIDs threw "
7985 "an exception not a subclass of std::exception"
7991 std::ostringstream os;
7992 os << *verbosePrefix <<
"Done with packCrsMatrixWithOwningPIDs"
7994 std::cerr << os.str();
7997 if (!comm.is_null()) {
7998 reduceAll<int, int>(*comm, REDUCE_MAX, lclErr, outArg(gblErr));
8002 TEUCHOS_TEST_FOR_EXCEPTION(
8003 true, std::runtime_error,
8004 "packCrsMatrixWithOwningPIDs threw an "
8005 "exception on at least one process.");
8009 destMat->numExportPacketsPerLID_.modify_host();
8010 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
8013 std::ostringstream os;
8014 os << *verbosePrefix <<
"Calling packCrsMatrixWithOwningPIDs"
8016 std::cerr << os.str();
8020 numExportPacketsPerLID,
8023 constantNumPackets);
8025 std::ostringstream os;
8026 os << *verbosePrefix <<
"Done with packCrsMatrixWithOwningPIDs"
8028 std::cerr << os.str();
8035 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC getOwningPIDs exchange remote data");
8036 if (!communication_needed) {
8038 std::ostringstream os;
8039 os << *verbosePrefix <<
"Communication not needed" << std::endl;
8040 std::cerr << os.str();
8044 if (constantNumPackets == 0) {
8046 std::ostringstream os;
8047 os << *verbosePrefix <<
"Reverse mode, variable # packets / LID"
8049 std::cerr << os.str();
8054 destMat->numExportPacketsPerLID_.sync_host();
8055 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8057 destMat->numImportPacketsPerLID_.sync_host();
8058 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8062 std::ostringstream os;
8063 os << *verbosePrefix <<
"Calling 3-arg doReversePostsAndWaits"
8065 std::cerr << os.str();
8067 Distor.doReversePostsAndWaits(destMat->numExportPacketsPerLID_.view_host(), 1,
8068 destMat->numImportPacketsPerLID_.view_host());
8070 std::ostringstream os;
8071 os << *verbosePrefix <<
"Finished 3-arg doReversePostsAndWaits"
8073 std::cerr << os.str();
8076 size_t totalImportPackets = 0;
8078 totalImportPackets += numImportPacketsPerLID[i];
8083 destMat->reallocImportsIfNeeded(totalImportPackets, verbose,
8084 verbosePrefix.get());
8085 destMat->imports_.modify_host();
8086 auto hostImports = destMat->imports_.view_host();
8089 destMat->exports_.sync_host();
8090 auto hostExports = destMat->exports_.view_host();
8092 std::ostringstream os;
8093 os << *verbosePrefix <<
"Calling 4-arg doReversePostsAndWaits"
8095 std::cerr << os.str();
8097 Distor.doReversePostsAndWaits(hostExports,
8098 numExportPacketsPerLID,
8100 numImportPacketsPerLID);
8102 std::ostringstream os;
8103 os << *verbosePrefix <<
"Finished 4-arg doReversePostsAndWaits"
8105 std::cerr << os.str();
8109 std::ostringstream os;
8110 os << *verbosePrefix <<
"Reverse mode, constant # packets / LID"
8112 std::cerr << os.str();
8114 destMat->imports_.modify_host();
8115 auto hostImports = destMat->imports_.view_host();
8118 destMat->exports_.sync_host();
8119 auto hostExports = destMat->exports_.view_host();
8121 std::ostringstream os;
8122 os << *verbosePrefix <<
"Calling 3-arg doReversePostsAndWaits"
8124 std::cerr << os.str();
8126 Distor.doReversePostsAndWaits(hostExports,
8130 std::ostringstream os;
8131 os << *verbosePrefix <<
"Finished 3-arg doReversePostsAndWaits"
8133 std::cerr << os.str();
8137 if (constantNumPackets == 0) {
8139 std::ostringstream os;
8140 os << *verbosePrefix <<
"Forward mode, variable # packets / LID"
8142 std::cerr << os.str();
8147 destMat->numExportPacketsPerLID_.sync_host();
8148 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8150 destMat->numImportPacketsPerLID_.sync_host();
8151 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8154 std::ostringstream os;
8155 os << *verbosePrefix <<
"Calling 3-arg doPostsAndWaits"
8157 std::cerr << os.str();
8159 Distor.doPostsAndWaits(destMat->numExportPacketsPerLID_.view_host(), 1,
8160 destMat->numImportPacketsPerLID_.view_host());
8162 std::ostringstream os;
8163 os << *verbosePrefix <<
"Finished 3-arg doPostsAndWaits"
8165 std::cerr << os.str();
8168 size_t totalImportPackets = 0;
8170 totalImportPackets += numImportPacketsPerLID[i];
8175 destMat->reallocImportsIfNeeded(totalImportPackets, verbose,
8176 verbosePrefix.get());
8177 destMat->imports_.modify_host();
8178 auto hostImports = destMat->imports_.view_host();
8181 destMat->exports_.sync_host();
8182 auto hostExports = destMat->exports_.view_host();
8184 std::ostringstream os;
8185 os << *verbosePrefix <<
"Calling 4-arg doPostsAndWaits"
8187 std::cerr << os.str();
8189 Distor.doPostsAndWaits(hostExports,
8190 numExportPacketsPerLID,
8192 numImportPacketsPerLID);
8194 std::ostringstream os;
8195 os << *verbosePrefix <<
"Finished 4-arg doPostsAndWaits"
8197 std::cerr << os.str();
8201 std::ostringstream os;
8202 os << *verbosePrefix <<
"Forward mode, constant # packets / LID"
8204 std::cerr << os.str();
8206 destMat->imports_.modify_host();
8207 auto hostImports = destMat->imports_.view_host();
8210 destMat->exports_.sync_host();
8211 auto hostExports = destMat->exports_.view_host();
8213 std::ostringstream os;
8214 os << *verbosePrefix <<
"Calling 3-arg doPostsAndWaits"
8216 std::cerr << os.str();
8218 Distor.doPostsAndWaits(hostExports,
8222 std::ostringstream os;
8223 os << *verbosePrefix <<
"Finished 3-arg doPostsAndWaits"
8225 std::cerr << os.str();
8236 bool runOnHost = std::is_same_v<typename device_type::memory_space, Kokkos::HostSpace> && !useKokkosPath;
8238 Teuchos::Array<int> RemotePids;
8240 Teuchos::Array<int> TargetPids;
8246 destMat->numImportPacketsPerLID_.modify_host();
8248 auto tmCopySPRdata = Teuchos::rcp(
new Tpetra::Details::ProfilingRegion(
"Tpetra TAFC unpack-count-resize + copy same-perm-remote data"));
8249 ArrayRCP<size_t> CSR_rowptr;
8250 ArrayRCP<GO> CSR_colind_GID;
8251 ArrayRCP<LO> CSR_colind_LID;
8252 ArrayRCP<Scalar> CSR_vals;
8254 destMat->imports_.sync_device();
8255 destMat->numImportPacketsPerLID_.sync_device();
8257 size_t N = BaseRowMap->getLocalNumElements();
8259 auto RemoteLIDs_d = RemoteLIDs.view_device();
8260 auto PermuteToLIDs_d = PermuteToLIDs.view_device();
8261 auto PermuteFromLIDs_d = PermuteFromLIDs.view_device();
8266 destMat->imports_.view_device(),
8267 destMat->numImportPacketsPerLID_.view_device(),
8281 if (
typeid(LO) ==
typeid(GO)) {
8282 CSR_colind_LID = Teuchos::arcp_reinterpret_cast<LO>(CSR_colind_GID);
8284 CSR_colind_LID.resize(CSR_colind_GID.size());
8286 CSR_colind_LID.resize(CSR_colind_GID.size());
8291 for (
size_t i = 0; i < static_cast<size_t>(TargetPids.size()); i++) {
8292 if (TargetPids[i] == -1) TargetPids[i] = MyPID;
8294 tmCopySPRdata = Teuchos::null;
8302 std::ostringstream os;
8303 os << *verbosePrefix <<
"Calling lowCommunicationMakeColMapAndReindex"
8305 std::cerr << os.str();
8308 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC makeColMap");
8309 Import_Util::lowCommunicationMakeColMapAndReindexSerial(CSR_rowptr(),
8320 std::ostringstream os;
8321 os << *verbosePrefix <<
"restrictComm="
8322 << (restrictComm ?
"true" :
"false") << std::endl;
8323 std::cerr << os.str();
8330 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC restrict colmap");
8332 ReducedColMap = (MyRowMap.getRawPtr() == MyColMap.getRawPtr()) ? ReducedRowMap : MyColMap->replaceCommWithSubset(ReducedComm);
8333 MyColMap = ReducedColMap;
8338 std::ostringstream os;
8339 os << *verbosePrefix <<
"Calling replaceColMap" << std::endl;
8340 std::cerr << os.str();
8342 destMat->replaceColMap(MyColMap);
8349 if (ReducedComm.is_null()) {
8351 std::ostringstream os;
8352 os << *verbosePrefix <<
"I am no longer in the communicator; "
8355 std::cerr << os.str();
8364 if ((!reverseMode && xferAsImport !=
nullptr) ||
8365 (reverseMode && xferAsExport !=
nullptr)) {
8367 std::ostringstream os;
8368 os << *verbosePrefix <<
"Calling sortCrsEntries" << endl;
8369 std::cerr << os.str();
8371 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC sortCrsEntries");
8375 }
else if ((!reverseMode && xferAsExport !=
nullptr) ||
8376 (reverseMode && xferAsImport !=
nullptr)) {
8378 std::ostringstream os;
8379 os << *verbosePrefix <<
"Calling sortAndMergeCrsEntries"
8381 std::cerr << os.str();
8383 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC sortAndMergeCrsEntries");
8387 if (CSR_rowptr[N] !=
static_cast<size_t>(CSR_vals.size())) {
8388 CSR_colind_LID.resize(CSR_rowptr[N]);
8389 CSR_vals.resize(CSR_rowptr[N]);
8392 TEUCHOS_TEST_FOR_EXCEPTION(
8393 true, std::logic_error,
8394 "Tpetra::CrsMatrix::"
8395 "transferAndFillComplete: Should never get here! "
8396 "Please report this bug to a Tpetra developer.");
8403 std::ostringstream os;
8404 os << *verbosePrefix <<
"Calling destMat->setAllValues" << endl;
8405 std::cerr << os.str();
8414 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC setAllValues");
8415 destMat->setAllValues(CSR_rowptr, CSR_colind_LID, CSR_vals);
8426 destMat->numImportPacketsPerLID_.modify_host();
8428 auto tmCopySPRdata = Teuchos::rcp(
new Tpetra::Details::ProfilingRegion(
"Tpetra TAFC unpack-count-resize + copy same-perm-remote data"));
8429 ArrayRCP<size_t> CSR_rowptr;
8430 ArrayRCP<GO> CSR_colind_GID;
8431 ArrayRCP<LO> CSR_colind_LID;
8432 ArrayRCP<Scalar> CSR_vals;
8434 destMat->imports_.sync_device();
8435 destMat->numImportPacketsPerLID_.sync_device();
8437 size_t N = BaseRowMap->getLocalNumElements();
8439 auto RemoteLIDs_d = RemoteLIDs.view_device();
8440 auto PermuteToLIDs_d = PermuteToLIDs.view_device();
8441 auto PermuteFromLIDs_d = PermuteFromLIDs.view_device();
8443 Kokkos::View<size_t*, device_type> CSR_rowptr_d;
8444 Kokkos::View<GO*, device_type> CSR_colind_GID_d;
8445 Kokkos::View<LO*, device_type> CSR_colind_LID_d;
8446 Kokkos::View<impl_scalar_type*, device_type> CSR_vals_d;
8447 Kokkos::View<int*, device_type> TargetPids_d;
8452 destMat->imports_.view_device(),
8453 destMat->numImportPacketsPerLID_.view_device(),
8465 Kokkos::resize(CSR_colind_LID_d, CSR_colind_GID_d.size());
8467 tmCopySPRdata = Teuchos::null;
8475 std::ostringstream os;
8476 os << *verbosePrefix <<
"Calling lowCommunicationMakeColMapAndReindex"
8478 std::cerr << os.str();
8481 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC makeColMap");
8493 std::ostringstream os;
8494 os << *verbosePrefix <<
"restrictComm="
8495 << (restrictComm ?
"true" :
"false") << std::endl;
8496 std::cerr << os.str();
8503 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC restrict colmap");
8505 ReducedColMap = (MyRowMap.getRawPtr() == MyColMap.getRawPtr()) ? ReducedRowMap : MyColMap->replaceCommWithSubset(ReducedComm);
8506 MyColMap = ReducedColMap;
8511 std::ostringstream os;
8512 os << *verbosePrefix <<
"Calling replaceColMap" << std::endl;
8513 std::cerr << os.str();
8515 destMat->replaceColMap(MyColMap);
8522 if (ReducedComm.is_null()) {
8524 std::ostringstream os;
8525 os << *verbosePrefix <<
"I am no longer in the communicator; "
8528 std::cerr << os.str();
8538 if ((!reverseMode && xferAsImport !=
nullptr) ||
8539 (reverseMode && xferAsExport !=
nullptr)) {
8541 std::ostringstream os;
8542 os << *verbosePrefix <<
"Calling sortCrsEntries" << endl;
8543 std::cerr << os.str();
8545 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC sortCrsEntries");
8549 }
else if ((!reverseMode && xferAsExport !=
nullptr) ||
8550 (reverseMode && xferAsImport !=
nullptr)) {
8552 std::ostringstream os;
8553 os << *verbosePrefix <<
"Calling sortAndMergeCrsEntries"
8555 std::cerr << os.str();
8557 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC sortAndMergeCrsEntries");
8561 ::KokkosSparse::SortAlgorithm::DEFAULT);
8563 TEUCHOS_TEST_FOR_EXCEPTION(
8564 true, std::logic_error,
8565 "Tpetra::CrsMatrix::"
8566 "transferAndFillComplete: Should never get here! "
8567 "Please report this bug to a Tpetra developer.");
8575 std::ostringstream os;
8576 os << *verbosePrefix <<
"Calling destMat->setAllValues" << endl;
8577 std::cerr << os.str();
8581 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra TAFC setAllValues");
8582 destMat->setAllValues(CSR_rowptr_d, CSR_colind_LID_d, CSR_vals_d);
8590 auto tmIESFC = Teuchos::rcp(
new Tpetra::Details::ProfilingRegion(
"Tpetra TAFC build importer and esfc"));
8592 Teuchos::ParameterList esfc_params;
8594 RCP<import_type> MyImport;
8597 if (iallreduceRequest.get() !=
nullptr) {
8599 std::ostringstream os;
8600 os << *verbosePrefix <<
"Calling iallreduceRequest->wait()"
8602 std::cerr << os.str();
8604 iallreduceRequest->wait();
8605 if (reduced_mismatch != 0) {
8611 Tpetra::Details::ProfilingRegion MMisMM(
"Tpetra isMM Block");
8615 std::ostringstream os;
8616 os << *verbosePrefix <<
"Getting CRS pointers" << endl;
8617 std::cerr << os.str();
8620 Teuchos::ArrayRCP<LocalOrdinal> type3LIDs;
8621 Teuchos::ArrayRCP<int> type3PIDs;
8622 auto rowptr = getCrsGraph()->getLocalRowPtrsHost();
8623 auto colind = getCrsGraph()->getLocalIndicesHost();
8626 std::ostringstream os;
8627 os << *verbosePrefix <<
"Calling reverseNeighborDiscovery" << std::endl;
8628 std::cerr << os.str();
8632 Tpetra::Details::ProfilingRegion MMrc(
"Tpetra isMMrevNeighDis");
8633 Import_Util::reverseNeighborDiscovery(*
this,
8645 std::ostringstream os;
8646 os << *verbosePrefix <<
"Done with reverseNeighborDiscovery" << std::endl;
8647 std::cerr << os.str();
8650 Teuchos::ArrayView<const int> EPID1 = MyImporter.is_null() ? Teuchos::ArrayView<const int>() : MyImporter->getExportPIDs();
8651 Teuchos::ArrayView<const LO> ELID1 = MyImporter.is_null() ? Teuchos::ArrayView<const LO>() : MyImporter->getExportLIDs();
8653 Teuchos::ArrayView<const int> TEPID2 = rowTransfer.getExportPIDs();
8654 Teuchos::ArrayView<const LO> TELID2 = rowTransfer.getExportLIDs();
8656 const int numCols = getGraph()->getColMap()->getLocalNumElements();
8658 std::vector<bool> IsOwned(numCols,
true);
8659 std::vector<int> SentTo(numCols, -1);
8660 if (!MyImporter.is_null()) {
8661 for (
auto&& rlid : MyImporter->getRemoteLIDs()) {
8662 IsOwned[rlid] =
false;
8666 std::vector<std::pair<int, GO>> usrtg;
8667 usrtg.reserve(TEPID2.size());
8670 const auto& colMap = *(this->getColMap());
8672 const LO row = TELID2[i];
8673 const int pid = TEPID2[i];
8674 for (
auto j = rowptr[row]; j < rowptr[row + 1]; ++j) {
8675 const int col = colind[j];
8676 if (IsOwned[col] && SentTo[col] != pid) {
8678 GO gid = colMap.getGlobalElement(col);
8679 usrtg.push_back(std::pair<int, GO>(pid, gid));
8686 std::sort(usrtg.begin(), usrtg.end());
8687 auto eopg = std ::unique(usrtg.begin(), usrtg.end());
8689 usrtg.erase(eopg, usrtg.end());
8692 Teuchos::ArrayRCP<int> EPID2 = Teuchos::arcp(
new int[type2_us_size], 0, type2_us_size,
true);
8693 Teuchos::ArrayRCP<LO> ELID2 = Teuchos::arcp(
new LO[type2_us_size], 0, type2_us_size,
true);
8696 for (
auto&& p : usrtg) {
8697 EPID2[pos] = p.first;
8698 ELID2[pos] = this->getDomainMap()->getLocalElement(p.second);
8702 Teuchos::ArrayView<int> EPID3 = type3PIDs();
8703 Teuchos::ArrayView<LO> ELID3 = type3LIDs();
8704 GO InfGID = std::numeric_limits<GO>::max();
8705 int InfPID = INT_MAX;
8709#define TPETRA_MIN3(x, y, z) ((x) < (y) ? (std::min(x, z)) : (std::min(y, z)))
8710 int i1 = 0, i2 = 0, i3 = 0;
8711 int Len1 = EPID1.size();
8712 int Len2 = EPID2.size();
8713 int Len3 = EPID3.size();
8715 int MyLen = Len1 + Len2 + Len3;
8716 Teuchos::ArrayRCP<LO> userExportLIDs = Teuchos::arcp(
new LO[MyLen], 0, MyLen,
true);
8717 Teuchos::ArrayRCP<int> userExportPIDs = Teuchos::arcp(
new int[MyLen], 0, MyLen,
true);
8720 while (i1 < Len1 || i2 < Len2 || i3 < Len3) {
8721 int PID1 = (i1 < Len1) ? (EPID1[i1]) : InfPID;
8722 int PID2 = (i2 < Len2) ? (EPID2[i2]) : InfPID;
8723 int PID3 = (i3 < Len3) ? (EPID3[i3]) : InfPID;
8725 GO GID1 = (i1 < Len1) ? getDomainMap()->getGlobalElement(ELID1[i1]) : InfGID;
8726 GO GID2 = (i2 < Len2) ? getDomainMap()->getGlobalElement(ELID2[i2]) : InfGID;
8727 GO GID3 = (i3 < Len3) ? getDomainMap()->getGlobalElement(ELID3[i3]) : InfGID;
8729 int MIN_PID = TPETRA_MIN3(PID1, PID2, PID3);
8730 GO MIN_GID = TPETRA_MIN3(((PID1 == MIN_PID) ? GID1 : InfGID), ((PID2 == MIN_PID) ? GID2 : InfGID), ((PID3 == MIN_PID) ? GID3 : InfGID));
8734 bool added_entry =
false;
8736 if (PID1 == MIN_PID && GID1 == MIN_GID) {
8737 userExportLIDs[iloc] = ELID1[i1];
8738 userExportPIDs[iloc] = EPID1[i1];
8743 if (PID2 == MIN_PID && GID2 == MIN_GID) {
8745 userExportLIDs[iloc] = ELID2[i2];
8746 userExportPIDs[iloc] = EPID2[i2];
8752 if (PID3 == MIN_PID && GID3 == MIN_GID) {
8754 userExportLIDs[iloc] = ELID3[i3];
8755 userExportPIDs[iloc] = EPID3[i3];
8763 std::ostringstream os;
8764 os << *verbosePrefix <<
"Create Import" << std::endl;
8765 std::cerr << os.str();
8768 Tpetra::Details::ProfilingRegion ismmIctor(
"Tpetra isMMIportCtor");
8769 Teuchos::RCP<Teuchos::ParameterList> plist = rcp(
new Teuchos::ParameterList());
8771 if ((MyDomainMap != MyColMap) && (!MyDomainMap->isSameAs(*MyColMap)))
8775 userExportLIDs.view(0, iloc).getConst(),
8776 userExportPIDs.view(0, iloc).getConst(),
8780 std::ostringstream os;
8781 os << *verbosePrefix <<
"Call expertStaticFillComplete" << std::endl;
8782 std::cerr << os.str();
8786 Tpetra::Details::ProfilingRegion esfc(
"Tpetra isMM::destMat->eSFC");
8787 if (!params.is_null())
8788 esfc_params.set(
"compute global constants", params->get(
"compute global constants",
true));
8789 destMat->expertStaticFillComplete(MyDomainMap, MyRangeMap, MyImport, Teuchos::null, rcp(
new Teuchos::ParameterList(esfc_params)));
8794 Tpetra::Details::ProfilingRegion MMnotMMblock(
"Tpetra TAFC notMMblock");
8796 std::ostringstream os;
8797 os << *verbosePrefix <<
"Create Import" << std::endl;
8798 std::cerr << os.str();
8801 Tpetra::Details::ProfilingRegion notMMIcTor(
"Tpetra TAFC notMMCreateImporter");
8802 if ((MyDomainMap != MyColMap) && (!MyDomainMap->isSameAs(*MyColMap)))
8803 MyImport = rcp(
new import_type(MyDomainMap, MyColMap, RemotePids, Teuchos::null));
8806 std::ostringstream os;
8807 os << *verbosePrefix <<
"Call expertStaticFillComplete" << endl;
8808 std::cerr << os.str();
8811 Tpetra::Details::ProfilingRegion esfcnotmm(
"Tpetra notMMdestMat->expertStaticFillComplete");
8813 if (!params.is_null()) {
8814 esfc_params.set(
"compute global constants",
8815 params->get(
"compute global constants",
true));
8817 destMat->expertStaticFillComplete(MyDomainMap, MyRangeMap,
8818 MyImport, Teuchos::null,
8819 rcp(
new Teuchos::ParameterList(esfc_params)));
8822 tmIESFC = Teuchos::null;
8825 std::ostringstream os;
8826 os << *verbosePrefix <<
"Done" << endl;
8827 std::cerr << os.str();
8831template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8833 importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8835 const Teuchos::RCP<const map_type>& domainMap,
8836 const Teuchos::RCP<const map_type>& rangeMap,
8837 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8838 transferAndFillComplete(destMatrix, importer, Teuchos::null, domainMap, rangeMap, params);
8841template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8843 importAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8846 const Teuchos::RCP<const map_type>& domainMap,
8847 const Teuchos::RCP<const map_type>& rangeMap,
8848 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8849 transferAndFillComplete(destMatrix, rowImporter, Teuchos::rcpFromRef(domainImporter), domainMap, rangeMap, params);
8852template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8854 exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8856 const Teuchos::RCP<const map_type>& domainMap,
8857 const Teuchos::RCP<const map_type>& rangeMap,
8858 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8859 transferAndFillComplete(destMatrix, exporter, Teuchos::null, domainMap, rangeMap, params);
8862template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8864 exportAndFillComplete(Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>& destMatrix,
8867 const Teuchos::RCP<const map_type>& domainMap,
8868 const Teuchos::RCP<const map_type>& rangeMap,
8869 const Teuchos::RCP<Teuchos::ParameterList>& params)
const {
8870 transferAndFillComplete(destMatrix, rowExporter, Teuchos::rcpFromRef(domainExporter), domainMap, rangeMap, params);
8873template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8874void copyAndPermuteStaticGraphNew(
8877 const size_t numSameIDs,
8878 const LocalOrdinal permuteToLIDs[],
8879 const LocalOrdinal permuteFromLIDs[],
8880 const size_t numPermutes) {
8883 using Teuchos::Array;
8884 using LO = LocalOrdinal;
8885 using GO = GlobalOrdinal;
8887 using impl_scalar_type =
typename KokkosKernels::ArithTraits<Scalar>::val_type;
8891 typedef typename crs_matrix_type::local_inds_device_view_type::non_const_value_type local_inds_device_value_t;
8893 typedef typename Node::execution_space exec_space;
8894 typedef Kokkos::RangePolicy<exec_space, LO> range_type;
8896 const LocalOrdinal LINV = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
8898 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermuteStaticGraphNew");
8900 const crs_matrix_type* srcMatCrsPtr =
dynamic_cast<const crs_matrix_type*
>(&srcMat);
8901 TEUCHOS_TEST_FOR_EXCEPTION(srcMatCrsPtr ==
nullptr, std::runtime_error,
"bad srcMatCrsPtr");
8902 const crs_matrix_type& srcMatCrs = *srcMatCrsPtr;
8904 crs_matrix_type* tgtMatCrsPtr =
dynamic_cast<crs_matrix_type*
>(&tgtMat);
8905 TEUCHOS_TEST_FOR_EXCEPTION(tgtMatCrsPtr ==
nullptr, std::runtime_error,
"bad tgtMatCrsPtr");
8906 crs_matrix_type& tgtMatCrs = *tgtMatCrsPtr;
8914 const auto& srcRowMap = *(srcMat.
getRowMap());
8915 auto comm = srcRowMap.getComm();
8917 const LO numSameIDs_as_LID =
static_cast<LO
>(numSameIDs);
8919 auto my_replaceGlobalValuesImpl_scalar = KOKKOS_LAMBDA(
8920 const bool sorted,
const bool atomic,
size_t hint[],
8921 const size_t numInTgtRow,
const local_inds_device_value_t tgtColInds[],
impl_scalar_type tgtRowVals[],
8927 if (lclColInd != LINV) {
8928 const size_t offset = KokkosSparse::findRelOffset(tgtColInds, numInTgtRow, lclColInd, hint[0], sorted);
8929 if (offset != numInTgtRow) {
8930 Kokkos::atomic_store(&tgtRowVals[offset], newVals);
8931 hint[0] = offset + 1;
8936 if (lclColInd != LINV) {
8937 const size_t offset = KokkosSparse::findRelOffset(tgtColInds, numInTgtRow, lclColInd, hint[0], sorted);
8938 if (offset != numInTgtRow) {
8939 tgtRowVals[offset] = newVals;
8940 hint[0] = offset + 1;
8948 if (sourceIsLocallyIndexed) {
8949 typename crs_matrix_type::row_ptrs_device_view_type tgtLocalRowPtrsDevice = tgtMatCrs.getLocalRowPtrsDevice();
8950 typename crs_matrix_type::local_inds_device_view_type tgtLocalColIndsDevice = tgtMatCrs.getLocalIndicesDevice();
8951 typename crs_matrix_type::row_ptrs_host_view_type srcLocalRowPtrsHost = srcMatCrs.getLocalRowPtrsHost();
8952 typename crs_matrix_type::row_ptrs_device_view_type srcLocalRowPtrsDevice = srcMatCrs.getLocalRowPtrsDevice();
8953 typename crs_matrix_type::local_inds_device_view_type srcLocalColIndsDevice = srcMatCrs.getLocalIndicesDevice();
8955 bool tgtMatIsSorted = tgtMatCrs.getCrsGraph()->isSorted();
8957 using local_map_type =
typename crs_matrix_type::map_type::local_map_type;
8959 local_map_type local_map = srcMat.
getRowMap()->getLocalMap();
8960 local_map_type local_col_map = srcMat.
getColMap()->getLocalMap();
8961 local_map_type tgt_local_map = tgtMatCrs.getRowMap()->getLocalMap();
8962 local_map_type tgt_local_col_map = tgtMatCrs.getColMap()->getLocalMap();
8964 auto vals = srcMatCrs.getLocalValuesDevice(Access::ReadOnly);
8965 auto tvals = tgtMatCrs.getLocalValuesDevice(Access::ReadWrite);
8967 Kokkos::parallel_for(
8968 "Tpetra_CrsMatrix::copyAndPermuteStaticGraph",
8969 range_type(0, numSameIDs_as_LID),
8970 KOKKOS_LAMBDA(
const LO sourceLID) {
8971 local_inds_device_value_t
start = srcLocalRowPtrsDevice(sourceLID);
8972 local_inds_device_value_t end = srcLocalRowPtrsDevice(sourceLID + 1);
8973 local_inds_device_value_t rowLength = (end -
start);
8975 local_inds_device_value_t tstart = tgtLocalRowPtrsDevice(sourceLID);
8976 local_inds_device_value_t tend = tgtLocalRowPtrsDevice(sourceLID + 1);
8977 local_inds_device_value_t numInTgtRow = (tend - tstart);
8979 KOKKOS_ASSERT(
static_cast<size_t>(tstart) < tvals.extent(0));
8981 const local_inds_device_value_t* tgtColInds = &tgtLocalColIndsDevice(tstart);
8984 for (LO j = 0; j < rowLength; j++) {
8985 local_inds_device_value_t ci = srcLocalColIndsDevice(start + j);
8986 GO gi = local_col_map.getGlobalElement(ci);
8987 const local_inds_device_value_t lclColInd = tgt_local_col_map.getLocalElement(gi);
8988 my_replaceGlobalValuesImpl_scalar(
8989 tgtMatIsSorted,
false, &hint, numInTgtRow, tgtColInds, tgtRowVals, lclColInd, vals(start + j));
8993 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
8997 const GO sourceGID = srcRowMap.getGlobalElement(sourceLID);
8998 const GO targetGID = sourceGID;
9000 Teuchos::ArrayView<const GO> rowIndsConstView;
9001 Teuchos::ArrayView<const Scalar> rowValsConstView;
9003 typename crs_matrix_type::global_inds_host_view_type rowIndsView;
9004 typename crs_matrix_type::values_host_view_type rowValsView;
9010 rowIndsConstView = Teuchos::ArrayView<const GO>(
9011 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9012 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9013 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9014 rowValsView.extent(0),
9015 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9022 tgtMatCrs.replaceGlobalValues(targetGID, rowIndsConstView, rowValsConstView);
9030 typename crs_matrix_type::nonconst_global_inds_host_view_type rowInds;
9031 typename crs_matrix_type::nonconst_values_host_view_type rowVals;
9033 const auto& tgtRowMap = *(tgtMat.
getRowMap());
9034 for (
size_t p = 0; p < numPermutes; ++p) {
9035 const GO sourceGID = srcRowMap.getGlobalElement(permuteFromLIDs[p]);
9036 const GO targetGID = tgtRowMap.getGlobalElement(permuteToLIDs[p]);
9038 Teuchos::ArrayView<const GO> rowIndsConstView;
9039 Teuchos::ArrayView<const Scalar> rowValsConstView;
9041 if (sourceIsLocallyIndexed) {
9043 if (rowLength >
static_cast<size_t>(rowInds.size())) {
9044 Kokkos::resize(rowInds, rowLength);
9045 Kokkos::resize(rowVals, rowLength);
9049 typename crs_matrix_type::nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(
9050 rowInds, std::make_pair((
size_t)0, rowLength));
9051 typename crs_matrix_type::nonconst_values_host_view_type rowValsView = Kokkos::subview(
9052 rowVals, std::make_pair((
size_t)0, rowLength));
9057 size_t checkRowLength = 0;
9058 srcMat.
getGlobalRowCopy(sourceGID, rowIndsView, rowValsView, checkRowLength);
9064 rowIndsConstView = Teuchos::ArrayView<const GO>(
9065 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9066 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9067 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9068 rowValsView.extent(0),
9069 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9073 typename crs_matrix_type::global_inds_host_view_type rowIndsView;
9074 typename crs_matrix_type::values_host_view_type rowValsView;
9080 rowIndsConstView = Teuchos::ArrayView<const GO>(
9081 rowIndsView.data(), rowIndsView.extent(0), Teuchos::RCP_DISABLE_NODE_LOOKUP);
9082 rowValsConstView = Teuchos::ArrayView<const Scalar>(
9083 reinterpret_cast<const Scalar*
>(rowValsView.data()),
9084 rowValsView.extent(0),
9085 Teuchos::RCP_DISABLE_NODE_LOOKUP);
9090 tgtMatCrs.replaceGlobalValues(targetGID, rowIndsConstView, rowValsConstView);
9102#define TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9104 template class CrsMatrix<SCALAR, LO, GO, NODE>;
9106#define TPETRA_CRSMATRIX_CONVERT_INSTANT(SO, SI, LO, GO, NODE) \
9108 template Teuchos::RCP<CrsMatrix<SO, LO, GO, NODE>> \
9109 CrsMatrix<SI, LO, GO, NODE>::convert<SO>() const;
9111#define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9113 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9114 importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9115 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9116 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9117 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& importer, \
9118 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9119 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9120 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9121 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9122 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9123 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9124 const Teuchos::RCP<Teuchos::ParameterList>& params);
9126#define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9128 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9129 importAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9130 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9131 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9132 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowImporter, \
9133 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9134 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9135 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainImporter, \
9136 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9137 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9138 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9139 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9140 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9141 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9142 const Teuchos::RCP<Teuchos::ParameterList>& params);
9144#define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9146 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9147 exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9148 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9149 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9150 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& exporter, \
9151 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9152 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9153 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9154 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9155 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9156 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9157 const Teuchos::RCP<Teuchos::ParameterList>& params);
9159#define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9161 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE>> \
9162 exportAndFillCompleteCrsMatrix(const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE>>& sourceMatrix, \
9163 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9164 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9165 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowExporter, \
9166 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9167 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9168 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainExporter, \
9169 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9170 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9171 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& domainMap, \
9172 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
9173 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
9174 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>>& rangeMap, \
9175 const Teuchos::RCP<Teuchos::ParameterList>& params);
9177#define TPETRA_CRSMATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9178 TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \
9179 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9180 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
9181 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
9182 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE)