diff --git a/packages/tpetra/core/src/Tpetra_CrsMatrix_def.hpp b/packages/tpetra/core/src/Tpetra_CrsMatrix_def.hpp index f0eef6b3b32e..132fd6c3c0f6 100644 --- a/packages/tpetra/core/src/Tpetra_CrsMatrix_def.hpp +++ b/packages/tpetra/core/src/Tpetra_CrsMatrix_def.hpp @@ -3271,22 +3271,22 @@ CrsMatrix:: "a view with global column indices by calling getGlobalRowCopy()."); const RowInfo rowInfo = staticGraph_->getRowInfo (localRow); - if (rowInfo.localRow != Teuchos::OrdinalTraits::invalid () && - rowInfo.numEntries > 0) { - indices = staticGraph_->lclIndsUnpacked_wdv.getHostSubview( - rowInfo.offset1D, - rowInfo.numEntries, - Access::ReadOnly); - values = valuesUnpacked_wdv.getHostSubview(rowInfo.offset1D, - rowInfo.numEntries, - Access::ReadOnly); - } - else { - // This does the right thing (reports an empty row) if the input - // row is invalid. - indices = local_inds_host_view_type(); - values = values_host_view_type(); - } + // if (rowInfo.localRow != Teuchos::OrdinalTraits::invalid () && + // rowInfo.numEntries > 0) { + // indices = staticGraph_->lclIndsUnpacked_wdv.getHostSubview( + // rowInfo.offset1D, + // rowInfo.numEntries, + // Access::ReadOnly); + // values = valuesUnpacked_wdv.getHostSubview(rowInfo.offset1D, + // rowInfo.numEntries, + // Access::ReadOnly); + // } + // else { + // // This does the right thing (reports an empty row) if the input + // // row is invalid. + // indices = local_inds_host_view_type(); + // values = values_host_view_type(); + // } #ifdef HAVE_TPETRA_DEBUG const char suffix[] = ". This should never happen. Please report this " @@ -3602,23 +3602,17 @@ CrsMatrix:: "diag.getMap ()->isCompatible (A.getRowMap ());"); #endif // HAVE_TPETRA_DEBUG - if (this->isFillComplete ()) { - const auto D_lcl = diag.getLocalViewDevice(Access::OverwriteAll); - // 1-D subview of the first (and only) column of D_lcl. - const auto D_lcl_1d = - Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0); + const auto D_lcl = diag.getLocalViewDevice(Access::OverwriteAll); + // 1-D subview of the first (and only) column of D_lcl. + const auto D_lcl_1d = + Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0); - const auto lclRowMap = rowMap.getLocalMap (); - const auto lclColMap = colMap.getLocalMap (); - using ::Tpetra::Details::getDiagCopyWithoutOffsets; - (void) getDiagCopyWithoutOffsets (D_lcl_1d, lclRowMap, - lclColMap, - getLocalMatrixDevice ()); - } - else { - using ::Tpetra::Details::getLocalDiagCopyWithoutOffsetsNotFillComplete; - (void) getLocalDiagCopyWithoutOffsetsNotFillComplete (diag, *this); - } + const auto lclRowMap = rowMap.getLocalMap (); + const auto lclColMap = colMap.getLocalMap (); + using ::Tpetra::Details::getDiagCopyWithoutOffsets; + (void) getDiagCopyWithoutOffsets (D_lcl_1d, lclRowMap, + lclColMap, + getLocalMatrixDevice ()); } template diff --git a/packages/tpetra/core/src/Tpetra_Details_getDiagCopyWithoutOffsets_def.hpp b/packages/tpetra/core/src/Tpetra_Details_getDiagCopyWithoutOffsets_def.hpp index c054dcc7a69d..b5f9281d8284 100644 --- a/packages/tpetra/core/src/Tpetra_Details_getDiagCopyWithoutOffsets_def.hpp +++ b/packages/tpetra/core/src/Tpetra_Details_getDiagCopyWithoutOffsets_def.hpp @@ -38,23 +38,23 @@ namespace Details { template class GetLocalDiagCopyWithoutOffsetsNotFillCompleteFunctor { public: - typedef ::Tpetra::RowMatrix row_matrix_type; - typedef ::Tpetra::Vector vec_type; + using row_matrix_type = ::Tpetra::RowMatrix; + using vec_type = ::Tpetra::Vector; - typedef typename vec_type::impl_scalar_type IST; + using IST = typename vec_type::impl_scalar_type; // The output Vector determines the execution space. + using host_execution_space = typename vec_type::dual_view_type::t_host::execution_space; private: - typedef typename vec_type::dual_view_type::t_host::execution_space host_execution_space; - typedef typename vec_type::map_type map_type; + using map_type = typename vec_type::map_type; static bool graphIsSorted (const row_matrix_type& A) { using Teuchos::RCP; using Teuchos::rcp_dynamic_cast; - typedef Tpetra::CrsGraph crs_graph_type; - typedef Tpetra::RowGraph row_graph_type; + using crs_graph_type = Tpetra::CrsGraph; + using row_graph_type = Tpetra::RowGraph; // We conservatively assume not sorted. RowGraph lacks an // "isSorted" predicate, so we can't know for sure unless the cast @@ -108,6 +108,7 @@ class GetLocalDiagCopyWithoutOffsetsNotFillCompleteFunctor { void operator () (const LO& lclRowInd, LO& errCount) const { using KokkosSparse::findRelOffset; + Kokkos::printf("Not hanging yet [0]\n"); D_lcl_1d_(lclRowInd) = Kokkos::ArithTraits::zero (); const GO gblInd = lclRowMap_.getGlobalElement (lclRowInd); @@ -119,19 +120,23 @@ class GetLocalDiagCopyWithoutOffsetsNotFillCompleteFunctor { else { // row index is also in the column Map on this process typename row_matrix_type::local_inds_host_view_type lclColInds; typename row_matrix_type::values_host_view_type curVals; + Kokkos::printf("Not hanging yet [1]\n"); A_.getLocalRowView(lclRowInd, lclColInds, curVals); + Kokkos::printf("Not hanging yet [2]\n"); LO numEnt = lclColInds.extent(0); // The search hint is always zero, since we only call this // once per row of the matrix. const LO hint = 0; const LO offset = findRelOffset (lclColInds, numEnt, lclColInd, hint, sorted_); + Kokkos::printf("Not hanging yet [3]\n"); if (offset == numEnt) { // didn't find the diagonal column index errCount++; } else { D_lcl_1d_(lclRowInd) = curVals[offset]; } + Kokkos::printf("Not hanging yet [4]\n"); } } @@ -155,8 +160,7 @@ getLocalDiagCopyWithoutOffsetsNotFillComplete ( ::Tpetra::Vector using Teuchos::outArg; using Teuchos::REDUCE_MIN; using Teuchos::reduceAll; - typedef GetLocalDiagCopyWithoutOffsetsNotFillCompleteFunctor functor_type; + using functor_type = GetLocalDiagCopyWithoutOffsetsNotFillCompleteFunctor; // The functor's constructor does error checking and executes the // thread-parallel kernel. @@ -212,6 +216,8 @@ getLocalDiagCopyWithoutOffsetsNotFillComplete ( ::Tpetra::Vector } } else { // ! debug + Kokkos::printf("\n\n\nRunning on Kokkos::Serial: %s\n", std::is_same_v ? "true" : "false"); + functor_type functor (lclNumErrs, diag, A); } diff --git a/packages/tpetra/core/test/CrsMatrix/CMakeLists.txt b/packages/tpetra/core/test/CrsMatrix/CMakeLists.txt index 914cf4910d57..e46a3319bac4 100644 --- a/packages/tpetra/core/test/CrsMatrix/CMakeLists.txt +++ b/packages/tpetra/core/test/CrsMatrix/CMakeLists.txt @@ -524,6 +524,15 @@ if ( ) endif() +TRIBITS_ADD_EXECUTABLE_AND_TEST( + CrsMatrix_getLocalDiagCopy + SOURCES + CrsMatrix_getLocalDiagCopy.cpp + ${TEUCHOS_STD_UNIT_TEST_MAIN} + COMM mpi + STANDARD_PASS_OUTPUT + ) + SET(TIMING_INSTALLS "") INSTALL(TARGETS ${TIMING_INSTALLS} diff --git a/packages/tpetra/core/test/CrsMatrix/CrsMatrix_getLocalDiagCopy.cpp b/packages/tpetra/core/test/CrsMatrix/CrsMatrix_getLocalDiagCopy.cpp new file mode 100644 index 000000000000..1738c506d059 --- /dev/null +++ b/packages/tpetra/core/test/CrsMatrix/CrsMatrix_getLocalDiagCopy.cpp @@ -0,0 +1,112 @@ +// @HEADER +// ***************************************************************************** +// Tpetra: Templated Linear Algebra Services Package +// +// Copyright 2008 NTESS and the Tpetra contributors. +// SPDX-License-Identifier: BSD-3-Clause +// ***************************************************************************** +// @HEADER + +#include +#include + +#include +#include +#include +#include + +namespace { // (anonymous) + + template + void getLocalDiagCopyTest(Teuchos::FancyOStream& out, bool& success) { + using Teuchos::RCP; + + using map_type = Tpetra::Map; + using vec_type = Tpetra::Vector; + using crs_matrix_type = Tpetra::CrsMatrix; + using row_matrix_type = Tpetra::RowMatrix; + + using STS = Teuchos::ScalarTraits; + using MT = typename STS::magnitudeType; + const Scalar SC_ONE = STS::one(); + + using LOT = Teuchos::OrdinalTraits; + const LO LO_INVALID = LOT::invalid(); + const LO LO_ONE = LOT::one(); + const GO GO_ONE = Teuchos::OrdinalTraits::one(); + + int lclSuccess = 0; + int gblSuccess = 0; + + RCP > comm = Tpetra::getDefaultComm (); + const size_t numProc = comm->getSize(); + const size_t myProc = comm->getRank(); + + // create a Map + RCP map = Tpetra::createContigMapWithNode (LO_INVALID, + LO_ONE + LO_ONE, + comm); + + // Create a matrix with at most 3 entries per row + RCP matrix = Teuchos::rcp (new crs_matrix_type (map, 3)); + const Scalar rankAsScalar = static_cast(static_cast(comm->getRank())); + + Teuchos::Array vals = {{SC_ONE, rankAsScalar + SC_ONE, SC_ONE}}; + for(size_t lclRowIdx = 0; lclRowIdx < 2; ++lclRowIdx) { + const GO gblRowIdx = Teuchos::as(2*myProc + lclRowIdx); + Teuchos::Array cols = {{gblRowIdx - GO_ONE, gblRowIdx, gblRowIdx + GO_ONE}}; + + if((myProc == 0) && (lclRowIdx == 0)) { // First row of the matrix + matrix->insertGlobalValues(gblRowIdx, cols(1, 2), vals(1, 2)); + } else if((myProc == numProc - 1) && (lclRowIdx == 1)) { // Last row of the matrix + matrix->insertGlobalValues(gblRowIdx, cols(0, 2), vals(0, 2)); + } else { + matrix->insertGlobalValues(gblRowIdx, cols(), vals()); + } + } + + matrix->fillComplete(); + + // Make sure that all processes got this far. + { + lclSuccess = success ? 1 : 0; + gblSuccess = 0; + Teuchos::reduceAll (*comm, Teuchos::REDUCE_MIN, lclSuccess, Teuchos::outArg (gblSuccess)); + success = success && (gblSuccess == 1); + TEST_EQUALITY_CONST( gblSuccess, 1 ); + } + + RCP diag = Teuchos::rcp(new vec_type(map)); + diag->putScalar(-1); + + if constexpr (Tag == 0) { + matrix->resumeFill(); + } + matrix->getLocalDiagCopy(*diag); + } + + // Unit test of getLocalDiagCopy + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL( CrsMatrix, getLocalDiagCopy, Scalar, LO, GO, Node ) + { + getLocalDiagCopyTest(out, success); + } + + // Unit test of getLocalDiagCopy + TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL( CrsMatrix, getLocalDiagCopyFillActive, Scalar, LO, GO, Node ) + { + getLocalDiagCopyTest(out, success); + } + + // + // INSTANTIATIONS + // + +#define UNIT_TEST_GROUP( SCALAR, LO, GO, NODE ) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT( CrsMatrix, getLocalDiagCopy, SCALAR, LO, GO, NODE ) \ + TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT( CrsMatrix, getLocalDiagCopyFillActive, SCALAR, LO, GO, NODE ) + + TPETRA_ETI_MANGLING_TYPEDEFS() + + TPETRA_INSTANTIATE_SLGN( UNIT_TEST_GROUP ) + +} // namespace (anonymous)