42 #ifndef TPETRA_CRSMATRIX_DEF_HPP 43 #define TPETRA_CRSMATRIX_DEF_HPP 53 #include "Tpetra_RowMatrix.hpp" 58 #include "Tpetra_Details_getDiagCopyWithoutOffsets.hpp" 61 #include "Teuchos_SerialDenseMatrix.hpp" 62 #include "Kokkos_Sparse_getDiagCopy.hpp" 80 template<
class Scalar>
84 typedef Teuchos::ScalarTraits<Scalar> STS;
85 return std::max (STS::magnitude (x), STS::magnitude (y));
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
95 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
96 CrsMatrix (
const Teuchos::RCP<const map_type>& rowMap,
97 size_t maxNumEntriesPerRow,
99 const Teuchos::RCP<Teuchos::ParameterList>& params) :
104 fillComplete_ (false),
105 frobNorm_ (-STM::one ())
107 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, size_t, " 108 "ProfileType[, RCP<ParameterList>]): ";
109 Teuchos::RCP<crs_graph_type> graph;
111 graph = Teuchos::rcp (
new crs_graph_type (rowMap, maxNumEntriesPerRow,
114 catch (std::exception& e) {
115 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
116 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 117 "size_t, ProfileType[, RCP<ParameterList>]) threw an exception: " 124 staticGraph_ = myGraph_;
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
132 const Teuchos::ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
134 const Teuchos::RCP<Teuchos::ParameterList>& params) :
139 fillComplete_ (false),
140 frobNorm_ (-STM::one ())
142 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, " 143 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]): ";
144 Teuchos::RCP<crs_graph_type> graph;
146 graph = Teuchos::rcp (
new crs_graph_type (rowMap, NumEntriesPerRowToAlloc,
149 catch (std::exception &e) {
150 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
151 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 152 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]) threw " 153 "an exception: " << e.what ());
159 staticGraph_ = graph;
164 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
167 const Teuchos::RCP<const map_type>& colMap,
168 size_t maxNumEntriesPerRow,
170 const Teuchos::RCP<Teuchos::ParameterList>& params) :
175 fillComplete_ (false),
176 frobNorm_ (-STM::one ())
178 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, RCP<const Map>, " 179 "size_t, ProfileType[, RCP<ParameterList>]): ";
181 #ifdef HAVE_TPETRA_DEBUG 183 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
184 (! staticGraph_.is_null (), std::logic_error,
185 "staticGraph_ is not null at the beginning of the constructor. " 186 "Please report this bug to the Tpetra developers.");
187 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
188 (! myGraph_.is_null (), std::logic_error,
189 "myGraph_ is not null at the beginning of the constructor. " 190 "Please report this bug to the Tpetra developers.");
191 #endif // HAVE_TPETRA_DEBUG 193 Teuchos::RCP<crs_graph_type> graph;
199 catch (std::exception &e) {
200 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
201 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 202 "RCP<const Map>, size_t, ProfileType[, RCP<ParameterList>]) threw an " 203 "exception: " << e.what ());
209 staticGraph_ = myGraph_;
214 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
217 const Teuchos::RCP<const map_type>& colMap,
218 const Teuchos::ArrayRCP<const size_t>& numEntPerRow,
220 const Teuchos::RCP<Teuchos::ParameterList>& params) :
225 fillComplete_ (false),
226 frobNorm_ (-STM::one ())
228 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, RCP<const Map>, " 229 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]): ";
230 Teuchos::RCP<crs_graph_type> graph;
232 graph = Teuchos::rcp (
new crs_graph_type (rowMap, colMap, numEntPerRow,
235 catch (std::exception &e) {
236 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
237 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 238 "RCP<const Map>, ArrayRCP<const size_t>, ProfileType[, " 239 "RCP<ParameterList>]) threw an exception: " << e.what ());
245 staticGraph_ = graph;
250 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
252 CrsMatrix (
const Teuchos::RCP<const crs_graph_type>& graph,
253 const Teuchos::RCP<Teuchos::ParameterList>& ) :
255 staticGraph_ (graph),
256 storageStatus_ (
Details::STORAGE_1D_PACKED),
257 fillComplete_ (false),
258 frobNorm_ (-STM::one ())
260 typedef typename local_matrix_type::values_type values_type;
261 const char tfecfFuncName[] =
"CrsMatrix(RCP<const CrsGraph>[, " 262 "RCP<ParameterList>]): ";
263 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
264 (graph.is_null (), std::runtime_error,
"Input graph is null.");
265 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
266 (! graph->isFillComplete (), std::runtime_error,
"Input graph is not " 267 "fill complete. You must call fillComplete on the graph before using " 268 "it to construct a CrsMatrix. Note that calling resumeFill on the " 269 "graph makes it not fill complete, even if you had previously called " 270 "fillComplete. In that case, you must call fillComplete on the graph " 279 const size_t numCols = graph->getColMap ()->getNodeNumElements ();
280 auto lclGraph = graph->getLocalGraph ();
281 const size_t numEnt = lclGraph.entries.dimension_0 ();
282 values_type val (
"Tpetra::CrsMatrix::val", numEnt);
285 numCols, val, lclGraph);
294 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
297 const Teuchos::RCP<const map_type>& colMap,
298 const typename local_matrix_type::row_map_type& rowPointers,
299 const typename local_graph_type::entries_type::non_const_type& columnIndices,
300 const typename local_matrix_type::values_type& values,
301 const Teuchos::RCP<Teuchos::ParameterList>& params) :
303 storageStatus_ (
Details::STORAGE_1D_PACKED),
304 fillComplete_ (false),
305 frobNorm_ (-STM::one ())
308 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, " 309 "RCP<const Map>, ptr, ind, val[, params]): ";
310 const char suffix[] =
". Please report this bug to the Tpetra developers.";
316 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
317 (values.dimension_0 () != columnIndices.dimension_0 (),
318 std::invalid_argument,
"Input arrays don't have matching dimensions. " 319 "values.dimension_0() = " << values.dimension_0 () <<
" != " 320 "columnIndices.dimension_0() = " << columnIndices.dimension_0 () <<
".");
321 #ifdef HAVE_TPETRA_DEBUG 322 if (rowPointers.dimension_0 () != 0) {
323 using Kokkos::subview;
327 auto ptr_last_d = subview (rowPointers, rowPointers.dimension_0 () - 1);
328 auto ptr_last_h = Kokkos::create_mirror_view (ptr_last_d);
330 const size_t numEnt =
static_cast<size_t> (ptr_last_h ());
331 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
332 (numEnt != static_cast<size_t> (columnIndices.dimension_0 ()) ||
333 numEnt != static_cast<size_t> (values.dimension_0 ()),
334 std::invalid_argument,
"Last entry of rowPointers says that the matrix" 335 " has " << numEnt <<
" entr" << (numEnt != 1 ?
"ies" :
"y") <<
", but " 336 "the dimensions of columnIndices and values don't match this. " 337 "columnIndices.dimension_0() = " << columnIndices.dimension_0 () <<
338 " and values.dimension_0() = " << values.dimension_0 () <<
".");
340 #endif // HAVE_TPETRA_DEBUG 342 RCP<crs_graph_type> graph;
344 graph = Teuchos::rcp (
new crs_graph_type (rowMap, colMap, rowPointers,
345 columnIndices, params));
347 catch (std::exception& e) {
348 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
349 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 350 "RCP<const Map>, ptr, ind[, params]) threw an exception: " 358 auto lclGraph = graph->getLocalGraph ();
359 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
360 (lclGraph.row_map.dimension_0 () != rowPointers.dimension_0 () ||
361 lclGraph.entries.dimension_0 () != columnIndices.dimension_0 (),
362 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, " 363 "ind[, params]) did not set the local graph correctly." << suffix);
364 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
365 (lclGraph.entries.dimension_0 () != values.dimension_0 (),
366 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, ind[, " 367 "params]) did not set the local graph correctly. " 368 "lclGraph.entries.dimension_0() = " << lclGraph.entries.dimension_0 ()
369 <<
" != values.dimension_0() = " << values.dimension_0 () << suffix);
375 staticGraph_ = graph;
384 const size_t numCols = graph->getColMap ()->getNodeNumElements ();
386 numCols, values, lclGraph);
387 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
388 (
lclMatrix_.values.dimension_0 () != values.dimension_0 (),
389 std::logic_error,
"Local matrix's constructor did not set the values " 390 "correctly. lclMatrix_.values.dimension_0() = " <<
391 lclMatrix_.values.dimension_0 () <<
" != values.dimension_0() = " <<
392 values.dimension_0 () << suffix);
402 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
405 const Teuchos::RCP<const map_type>& colMap,
406 const Teuchos::ArrayRCP<size_t>& ptr,
407 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
408 const Teuchos::ArrayRCP<Scalar>& val,
409 const Teuchos::RCP<Teuchos::ParameterList>& params) :
411 storageStatus_ (
Details::STORAGE_1D_PACKED),
412 fillComplete_ (false),
413 frobNorm_ (-STM::one ())
415 using Kokkos::Compat::getKokkosViewDeepCopy;
416 using Teuchos::av_reinterpret_cast;
418 typedef typename local_matrix_type::values_type values_type;
420 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, " 421 "RCP<const Map>, ptr, ind, val[, params]): ";
423 RCP<crs_graph_type> graph;
428 catch (std::exception& e) {
429 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
430 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 431 "RCP<const Map>, ArrayRCP<size_t>, ArrayRCP<LocalOrdinal>[, " 432 "RCP<ParameterList>]) threw an exception: " << e.what ());
438 staticGraph_ = graph;
451 auto lclGraph = staticGraph_->getLocalGraph ();
452 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
453 (static_cast<size_t> (lclGraph.row_map.dimension_0 ()) != static_cast<size_t> (ptr.size ()) ||
454 static_cast<size_t> (lclGraph.entries.dimension_0 ()) != static_cast<size_t> (ind.size ()),
455 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, " 456 "ind[, params]) did not set the local graph correctly. Please " 457 "report this bug to the Tpetra developers.");
459 const size_t numCols = staticGraph_->getColMap ()->getNodeNumElements ();
460 values_type valIn = getKokkosViewDeepCopy<device_type> (av_reinterpret_cast<IST> (val ()));
462 numCols, valIn, lclGraph);
471 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
474 const Teuchos::RCP<const map_type>& colMap,
476 const Teuchos::RCP<Teuchos::ParameterList>& params) :
478 lclMatrix_ (lclMatrix),
479 k_values1D_ (lclMatrix.values),
480 storageStatus_ (
Details::STORAGE_1D_PACKED),
481 fillComplete_ (true),
482 frobNorm_ (-STM::one ())
484 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, " 485 "RCP<const Map>, local_matrix_type[, RCP<ParameterList>]): ";
486 Teuchos::RCP<crs_graph_type> graph;
489 lclMatrix.graph, params));
491 catch (std::exception& e) {
492 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
493 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, " 494 "RCP<const Map>, local_graph_type[, RCP<ParameterList>]) threw an " 495 "exception: " << e.what ());
497 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
498 (!graph->isFillComplete (), std::logic_error,
"CrsGraph constructor (RCP" 499 "<const Map>, RCP<const Map>, local_graph_type[, RCP<ParameterList>]) " 500 "did not produce a fill-complete graph. Please report this bug to the " 501 "Tpetra developers.");
506 staticGraph_ = graph;
510 #ifdef HAVE_TPETRA_DEBUG 511 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillActive (), std::logic_error,
512 "We're at the end of fillComplete(), but isFillActive() is true. " 513 "Please report this bug to the Tpetra developers.");
514 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete (), std::logic_error,
515 "We're at the end of fillComplete(), but isFillComplete() is false. " 516 "Please report this bug to the Tpetra developers.");
517 #endif // HAVE_TPETRA_DEBUG 521 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
526 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
527 Teuchos::RCP<const Teuchos::Comm<int> >
530 return getCrsGraph ()->getComm ();
533 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
537 return getCrsGraph ()->getNode ();
540 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
544 return getCrsGraph ()->getProfileType ();
547 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
551 return fillComplete_;
554 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
558 return ! fillComplete_;
561 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
565 return getCrsGraph()->isStorageOptimized();
568 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
572 return getCrsGraph ()->isLocallyIndexed ();
575 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
579 return getCrsGraph ()->isGloballyIndexed ();
582 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
586 return getCrsGraph ()->hasColMap ();
589 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
593 return getCrsGraph ()->getGlobalNumEntries ();
596 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
600 return getCrsGraph ()->getNodeNumEntries ();
603 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
607 return getCrsGraph ()->getGlobalNumRows ();
610 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
614 return getCrsGraph ()->getGlobalNumCols ();
617 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
621 return getCrsGraph ()->getNodeNumRows ();
624 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
628 return getCrsGraph ()->getNodeNumCols ();
631 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
635 return getCrsGraph ()->getGlobalNumDiags ();
638 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
642 return getCrsGraph ()->getNodeNumDiags ();
645 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
649 return getCrsGraph ()->getNumEntriesInGlobalRow (globalRow);
652 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
656 return getCrsGraph ()->getNumEntriesInLocalRow (localRow);
659 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
663 return getCrsGraph ()->getGlobalMaxNumRowEntries ();
666 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
670 return getCrsGraph ()->getNodeMaxNumRowEntries ();
673 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
677 return getRowMap ()->getIndexBase ();
680 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
681 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
684 return getCrsGraph ()->getRowMap ();
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
688 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
691 return getCrsGraph ()->getColMap ();
694 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
695 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
698 return getCrsGraph ()->getDomainMap ();
701 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
702 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
705 return getCrsGraph()->getRangeMap();
708 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
709 Teuchos::RCP<const RowGraph<LocalOrdinal, GlobalOrdinal, Node> >
712 if (staticGraph_ != Teuchos::null) {
718 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
719 Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> >
722 if (staticGraph_ != Teuchos::null) {
728 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
732 return getCrsGraph ()->isLowerTriangular ();
735 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
739 return getCrsGraph ()->isUpperTriangular ();
742 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
746 return myGraph_.is_null ();
749 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
756 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
763 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
768 #ifdef HAVE_TPETRA_DEBUG 772 if ((gas == GraphAlreadyAllocated) != staticGraph_->indicesAreAllocated()) {
773 const std::string err1 (
"allocateValues: The caller has asserted that " 775 const std::string err2 (
"already allocated, but the static graph says " 776 "that its indices are ");
777 const std::string err3 (
"already allocated. Please report this bug to " 778 "the Tpetra developers.");
779 TEUCHOS_TEST_FOR_EXCEPTION(gas == GraphAlreadyAllocated && ! staticGraph_->indicesAreAllocated(),
780 std::logic_error, err1 << err2 <<
"not " << err3);
781 TEUCHOS_TEST_FOR_EXCEPTION(gas != GraphAlreadyAllocated && staticGraph_->indicesAreAllocated(),
782 std::logic_error, err1 <<
"not " << err2 << err3);
790 TEUCHOS_TEST_FOR_EXCEPTION(
791 ! staticGraph_->indicesAreAllocated() && myGraph_.is_null(),
793 "allocateValues: The static graph says that its indices are not " 794 "allocated, but the graph is not owned by the matrix. Please report " 795 "this bug to the Tpetra developers.");
796 #endif // HAVE_TPETRA_DEBUG 798 if (gas == GraphNotYetAllocated) {
799 myGraph_->allocateIndices (lg);
810 const size_t lclNumRows = staticGraph_->getNodeNumRows ();
811 typename Graph::local_graph_type::row_map_type k_ptrs =
812 staticGraph_->k_rowPtrs_;
813 TEUCHOS_TEST_FOR_EXCEPTION(
814 k_ptrs.dimension_0 () != lclNumRows+1, std::logic_error,
815 "Tpetra::CrsMatrix::allocateValues: With StaticProfile, row offsets " 816 "array has length " << k_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 817 << (lclNumRows+1) <<
".");
822 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, lclNumRows);
823 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
825 const size_t lclTotalNumEntries =
static_cast<size_t> (k_ptrs_ent_h ());
834 typedef typename local_matrix_type::values_type values_type;
835 k_values1D_ = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
843 values2D_ = staticGraph_->template allocateValues2D<impl_scalar_type> ();
847 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
850 getAllValues (Teuchos::ArrayRCP<const size_t>& rowPointers,
851 Teuchos::ArrayRCP<const LocalOrdinal>& columnIndices,
852 Teuchos::ArrayRCP<const Scalar>& values)
const 855 const char tfecfFuncName[] =
"getAllValues: ";
856 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
857 columnIndices.size () != values.size (), std::runtime_error,
858 "Requires that columnIndices and values are the same size.");
860 RCP<const crs_graph_type> relevantGraph = getCrsGraph ();
861 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
862 relevantGraph.is_null (), std::runtime_error,
863 "Requires that getCrsGraph() is not null.");
865 rowPointers = relevantGraph->getNodeRowPtrs ();
867 catch (std::exception &e) {
868 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
869 true, std::runtime_error,
870 "Caught exception while calling graph->getNodeRowPtrs(): " 874 columnIndices = relevantGraph->getNodePackedIndices ();
876 catch (std::exception &e) {
877 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
878 true, std::runtime_error,
879 "Caught exception while calling graph->getNodePackedIndices(): " 882 Teuchos::ArrayRCP<const impl_scalar_type> vals =
883 Kokkos::Compat::persistingView (k_values1D_);
884 values = Teuchos::arcp_reinterpret_cast<
const Scalar> (vals);
887 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
893 using Kokkos::create_mirror_view;
894 using Teuchos::arcp_const_cast;
895 using Teuchos::Array;
896 using Teuchos::ArrayRCP;
900 typedef typename local_matrix_type::row_map_type row_map_type;
901 typedef typename Graph::local_graph_type::entries_type::non_const_type lclinds_1d_type;
902 typedef typename local_matrix_type::values_type values_type;
903 #ifdef HAVE_TPETRA_DEBUG 904 const char tfecfFuncName[] =
"fillLocalGraphAndMatrix (called from " 905 "fillComplete or expertStaticFillComplete): ";
906 #endif // HAVE_TPETRA_DEBUG 908 #ifdef HAVE_TPETRA_DEBUG 911 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
912 (myGraph_.is_null (), std::logic_error,
"The nonconst graph (myGraph_) " 913 "is null. This means that the matrix has a const (a.k.a. \"static\") " 914 "graph. fillComplete or expertStaticFillComplete should never call " 915 "fillLocalGraphAndMatrix in that case. " 916 "Please report this bug to the Tpetra developers.");
917 #endif // HAVE_TPETRA_DEBUG 919 const size_t lclNumRows = this->getNodeNumRows ();
927 typename row_map_type::non_const_type k_ptrs;
928 row_map_type k_ptrs_const;
929 lclinds_1d_type k_inds;
935 lclinds_1d_type k_lclInds1D_ = myGraph_->k_lclInds1D_;
937 typedef decltype (myGraph_->k_numRowEntries_) row_entries_type;
958 typename row_entries_type::const_type numRowEnt_h =
959 myGraph_->k_numRowEntries_;
960 #ifdef HAVE_TPETRA_DEBUG 961 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
962 (static_cast<size_t> (numRowEnt_h.dimension_0 ()) != lclNumRows,
963 std::logic_error,
"(DynamicProfile branch) numRowEnt_h has the " 964 "wrong length. numRowEnt_h.dimension_0() = " 965 << numRowEnt_h.dimension_0 () <<
" != getNodeNumRows() = " 966 << lclNumRows <<
".");
967 #endif // HAVE_TPETRA_DEBUG 972 k_ptrs =
typename row_map_type::non_const_type (
"Tpetra::CrsGraph::ptr",
974 typename row_map_type::non_const_type::HostMirror h_ptrs =
975 create_mirror_view (k_ptrs);
983 const size_t lclTotalNumEntries =
985 #ifdef HAVE_TPETRA_DEBUG 986 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
987 (static_cast<size_t> (h_ptrs.dimension_0 ()) != lclNumRows + 1,
988 std::logic_error,
"(DynamicProfile branch) After packing h_ptrs, " 989 "h_ptrs.dimension_0() = " << h_ptrs.dimension_0 () <<
" != " 990 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
992 const size_t h_ptrs_lastEnt = h_ptrs(lclNumRows);
993 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
994 (h_ptrs_lastEnt != lclTotalNumEntries, std::logic_error,
995 "(DynamicProfile branch) After packing h_ptrs, h_ptrs(lclNumRows=" 996 << lclNumRows <<
") = " << h_ptrs_lastEnt <<
" != total number " 997 "of entries on the calling process = " << lclTotalNumEntries <<
".");
999 #endif // HAVE_TPETRA_DEBUG 1002 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
1003 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1006 typename lclinds_1d_type::HostMirror h_inds = create_mirror_view (k_inds);
1007 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
1010 ArrayRCP<Array<LocalOrdinal> > lclInds2D = myGraph_->lclInds2D_;
1011 for (
size_t row = 0; row < lclNumRows; ++row) {
1012 const size_t numEnt = numRowEnt_h(row);
1013 std::copy (lclInds2D[row].begin(),
1014 lclInds2D[row].begin() + numEnt,
1015 h_inds.ptr_on_device() + h_ptrs(row));
1016 std::copy (values2D_[row].begin(),
1017 values2D_[row].begin() + numEnt,
1018 h_vals.ptr_on_device() + h_ptrs(row));
1027 k_ptrs_const = k_ptrs;
1029 #ifdef HAVE_TPETRA_DEBUG 1031 if (k_ptrs.dimension_0 () != 0) {
1032 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
1033 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1034 (numOffsets != lclNumRows + 1, std::logic_error,
"(DynamicProfile " 1035 "branch) After copying into k_ptrs, k_ptrs.dimension_0() = " <<
1036 numOffsets <<
" != (lclNumRows+1) = " << (lclNumRows+1) <<
".");
1041 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, numOffsets-1);
1042 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
1044 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1045 (static_cast<size_t> (k_ptrs_ent_h ()) != k_vals.dimension_0 (),
1046 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs(" 1047 << (numOffsets-1) <<
") = " << k_ptrs_ent_h () <<
" != " 1048 "k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1049 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1050 (static_cast<size_t> (k_ptrs_ent_h ()) != k_inds.dimension_0 (),
1051 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs(" 1052 << (numOffsets-1) <<
") = " << k_ptrs_ent_h () <<
" != " 1053 "k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1055 #endif // HAVE_TPETRA_DEBUG 1064 typename Graph::local_graph_type::row_map_type curRowOffsets =
1065 myGraph_->k_rowPtrs_;
1067 #ifdef HAVE_TPETRA_DEBUG 1068 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1069 (curRowOffsets.dimension_0 () == 0, std::logic_error,
1070 "(StaticProfile branch) curRowOffsets.dimension_0() == 0.");
1071 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1072 (curRowOffsets.dimension_0 () != lclNumRows + 1, std::logic_error,
1073 "(StaticProfile branch) curRowOffsets.dimension_0() = " 1074 << curRowOffsets.dimension_0 () <<
" != lclNumRows + 1 = " 1075 << (lclNumRows + 1) <<
".")
1077 const size_t numOffsets = curRowOffsets.dimension_0 ();
1081 auto curRowOffsets_ent_d =
1082 Kokkos::subview (curRowOffsets, numOffsets - 1);
1083 auto curRowOffsets_ent_h = create_mirror_view (curRowOffsets_ent_d);
1085 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1087 myGraph_->k_lclInds1D_.dimension_0 () != curRowOffsets_ent_h (),
1088 std::logic_error,
"(StaticProfile branch) numOffsets = " <<
1089 numOffsets <<
" != 0 and myGraph_->k_lclInds1D_.dimension_0() = " 1090 << myGraph_->k_lclInds1D_.dimension_0 () <<
" != curRowOffsets(" 1091 << numOffsets <<
") = " << curRowOffsets_ent_h () <<
".");
1093 #endif // HAVE_TPETRA_DEBUG 1095 if (myGraph_->nodeNumEntries_ != myGraph_->nodeNumAllocated_) {
1102 #ifdef HAVE_TPETRA_DEBUG 1103 if (curRowOffsets.dimension_0 () != 0) {
1104 const size_t numOffsets =
1105 static_cast<size_t> (curRowOffsets.dimension_0 ());
1109 auto curRowOffsets_ent_d = Kokkos::subview (curRowOffsets, numOffsets-1);
1110 auto curRowOffsets_ent_h = create_mirror_view (curRowOffsets_ent_d);
1112 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1113 (static_cast<size_t> (curRowOffsets_ent_h ()) !=
1114 static_cast<size_t> (k_values1D_.dimension_0 ()),
1115 std::logic_error,
"(StaticProfile unpacked branch) Before " 1116 "allocating or packing, curRowOffsets(" << (numOffsets-1) <<
") = " 1117 << curRowOffsets_ent_h () <<
" != k_values1D_.dimension_0()" 1118 " = " << k_values1D_.dimension_0 () <<
".");
1119 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1120 (static_cast<size_t> (curRowOffsets_ent_h ()) !=
1121 static_cast<size_t> (myGraph_->k_lclInds1D_.dimension_0 ()),
1122 std::logic_error,
"(StaticProfile unpacked branch) Before " 1123 "allocating or packing, curRowOffsets(" << (numOffsets-1) <<
") = " 1124 << curRowOffsets_ent_h ()
1125 <<
" != myGraph_->k_lclInds1D_.dimension_0() = " 1126 << myGraph_->k_lclInds1D_.dimension_0 () <<
".");
1128 #endif // HAVE_TPETRA_DEBUG 1136 size_t lclTotalNumEntries = 0;
1138 typename row_map_type::non_const_type::HostMirror h_ptrs;
1143 typename row_map_type::non_const_type
1144 packedRowOffsets (
"Tpetra::CrsGraph::ptr", lclNumRows + 1);
1145 typename row_entries_type::const_type numRowEnt_h =
1146 myGraph_->k_numRowEntries_;
1149 lclTotalNumEntries =
1153 k_ptrs = packedRowOffsets;
1154 k_ptrs_const = k_ptrs;
1157 #ifdef HAVE_TPETRA_DEBUG 1158 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1159 (static_cast<size_t> (k_ptrs.dimension_0 ()) != lclNumRows + 1,
1161 "(StaticProfile unpacked branch) After packing k_ptrs, " 1162 "k_ptrs.dimension_0() = " << k_ptrs.dimension_0 () <<
" != " 1163 "lclNumRows+1 = " << (lclNumRows+1) <<
".");
1168 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, lclNumRows);
1169 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
1171 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1172 (k_ptrs_ent_h () != lclTotalNumEntries, std::logic_error,
1173 "(StaticProfile unpacked branch) After filling k_ptrs, " 1174 "k_ptrs(lclNumRows=" << lclNumRows <<
") = " << k_ptrs_ent_h ()
1175 <<
" != total number of entries on the calling process = " 1176 << lclTotalNumEntries <<
".");
1178 #endif // HAVE_TPETRA_DEBUG 1181 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
1182 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1194 typedef pack_functor<
typename Graph::local_graph_type::entries_type::non_const_type,
1195 typename Graph::local_graph_type::row_map_type>
1197 inds_packer_type indsPacker (k_inds, myGraph_->k_lclInds1D_,
1198 k_ptrs, curRowOffsets);
1199 Kokkos::parallel_for (lclNumRows, indsPacker);
1203 typedef pack_functor<values_type, row_map_type> vals_packer_type;
1204 vals_packer_type valsPacker (k_vals, this->k_values1D_,
1205 k_ptrs, curRowOffsets);
1206 Kokkos::parallel_for (lclNumRows, valsPacker);
1208 #ifdef HAVE_TPETRA_DEBUG 1209 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1210 (k_ptrs.dimension_0 () == 0, std::logic_error,
1211 "(StaticProfile \"Optimize Storage\" = " 1212 "true branch) After packing, k_ptrs.dimension_0() = 0. This " 1213 "probably means that k_rowPtrs_ was never allocated.");
1214 if (k_ptrs.dimension_0 () != 0) {
1215 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
1220 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, numOffsets - 1);
1221 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
1223 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1224 (static_cast<size_t> (k_ptrs_ent_h ()) != k_vals.dimension_0 (),
1226 "(StaticProfile \"Optimize Storage\"=true branch) After packing, " 1227 "k_ptrs(" << (numOffsets-1) <<
") = " << k_ptrs_ent_h () <<
1228 " != k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1229 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1230 (static_cast<size_t> (k_ptrs_ent_h ()) != k_inds.dimension_0 (),
1232 "(StaticProfile \"Optimize Storage\"=true branch) After packing, " 1233 "k_ptrs(" << (numOffsets-1) <<
") = " << k_ptrs_ent_h () <<
1234 " != k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1236 #endif // HAVE_TPETRA_DEBUG 1239 k_ptrs_const = myGraph_->k_rowPtrs_;
1240 k_inds = myGraph_->k_lclInds1D_;
1241 k_vals = this->k_values1D_;
1243 #ifdef HAVE_TPETRA_DEBUG 1244 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1245 (k_ptrs_const.dimension_0 () == 0, std::logic_error,
1246 "(StaticProfile \"Optimize Storage\"=false branch) " 1247 "k_ptrs_const.dimension_0() = 0. This probably means that " 1248 "k_rowPtrs_ was never allocated.");
1249 if (k_ptrs_const.dimension_0 () != 0) {
1250 const size_t numOffsets =
static_cast<size_t> (k_ptrs_const.dimension_0 ());
1255 auto k_ptrs_const_ent_d = Kokkos::subview (k_ptrs_const, numOffsets - 1);
1256 auto k_ptrs_const_ent_h = create_mirror_view (k_ptrs_const_ent_d);
1258 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1259 (static_cast<size_t> (k_ptrs_const_ent_h ()) != k_vals.dimension_0 (),
1261 "(StaticProfile \"Optimize Storage\"=false branch) " 1262 "k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const_ent_h ()
1263 <<
" != k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1264 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1265 (static_cast<size_t> (k_ptrs_const_ent_h ()) != k_inds.dimension_0 (),
1267 "(StaticProfile \"Optimize Storage\" = false branch) " 1268 "k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const_ent_h ()
1269 <<
" != k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1271 #endif // HAVE_TPETRA_DEBUG 1275 #ifdef HAVE_TPETRA_DEBUG 1277 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1278 (static_cast<size_t> (k_ptrs_const.dimension_0 ()) != lclNumRows + 1,
1279 std::logic_error,
"After packing, k_ptrs_const.dimension_0() = " <<
1280 k_ptrs_const.dimension_0 () <<
" != lclNumRows+1 = " << (lclNumRows+1)
1282 if (k_ptrs_const.dimension_0 () != 0) {
1283 const size_t numOffsets =
static_cast<size_t> (k_ptrs_const.dimension_0 ());
1288 auto k_ptrs_const_ent_d = Kokkos::subview (k_ptrs_const, numOffsets - 1);
1289 auto k_ptrs_const_ent_h = create_mirror_view (k_ptrs_const_ent_d);
1292 const size_t k_ptrs_const_numOffsetsMinus1 = k_ptrs_const_ent_h ();
1293 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1294 (k_ptrs_const_numOffsetsMinus1 != k_vals.dimension_0 (),
1295 std::logic_error,
"After packing, k_ptrs_const(" << (numOffsets-1) <<
1296 ") = " << k_ptrs_const_numOffsetsMinus1 <<
" != k_vals.dimension_0()" 1297 " = " << k_vals.dimension_0 () <<
".");
1298 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1299 (k_ptrs_const_numOffsetsMinus1 != k_inds.dimension_0 (),
1300 std::logic_error,
"After packing, k_ptrs_const(" << (numOffsets-1) <<
1301 ") = " << k_ptrs_const_numOffsetsMinus1 <<
" != k_inds.dimension_0()" 1302 " = " << k_inds.dimension_0 () <<
".");
1304 #endif // HAVE_TPETRA_DEBUG 1310 const bool defaultOptStorage =
1311 ! isStaticGraph () || staticGraph_->isStorageOptimized ();
1312 const bool requestOptimizedStorage =
1313 (! params.is_null () && params->get (
"Optimize Storage", defaultOptStorage)) ||
1314 (params.is_null () && defaultOptStorage);
1321 if (requestOptimizedStorage) {
1327 myGraph_->lclInds2D_ = null;
1328 myGraph_->k_numRowEntries_ = row_entries_type ();
1331 this->values2D_ = null;
1334 myGraph_->k_rowPtrs_ = k_ptrs_const;
1335 myGraph_->k_lclInds1D_ = k_inds;
1336 this->k_values1D_ = k_vals;
1340 myGraph_->nodeNumAllocated_ = myGraph_->nodeNumEntries_;
1344 myGraph_->storageStatus_ = Details::STORAGE_1D_PACKED;
1345 this->storageStatus_ = Details::STORAGE_1D_PACKED;
1356 myGraph_->lclGraph_ =
1361 getNodeNumCols (), k_vals,
1362 myGraph_->lclGraph_);
1365 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1370 using Kokkos::create_mirror_view;
1371 using Teuchos::ArrayRCP;
1372 using Teuchos::Array;
1373 using Teuchos::null;
1376 typedef LocalOrdinal LO;
1377 typedef typename Graph::local_graph_type::row_map_type row_map_type;
1378 typedef typename row_map_type::non_const_type non_const_row_map_type;
1379 typedef typename local_matrix_type::values_type values_type;
1380 #ifdef HAVE_TPETRA_DEBUG 1381 const char tfecfFuncName[] =
"fillLocalMatrix (called from fillComplete): ";
1382 #endif // HAVE_TPETRA_DEBUG 1384 const size_t lclNumRows = getNodeNumRows();
1385 const map_type& rowMap = * (getRowMap ());
1386 RCP<node_type> node = rowMap.
getNode ();
1397 ArrayRCP<Array<LO> > lclInds2D = staticGraph_->lclInds2D_;
1398 size_t nodeNumEntries = staticGraph_->nodeNumEntries_;
1399 size_t nodeNumAllocated = staticGraph_->nodeNumAllocated_;
1400 row_map_type k_rowPtrs_ = staticGraph_->lclGraph_.row_map;
1402 row_map_type k_ptrs;
1408 bool requestOptimizedStorage =
true;
1409 const bool default_OptimizeStorage =
1410 ! isStaticGraph () || staticGraph_->isStorageOptimized ();
1411 if (! params.is_null () && ! params->get (
"Optimize Storage", default_OptimizeStorage)) {
1412 requestOptimizedStorage =
false;
1419 if (! staticGraph_->isStorageOptimized () && requestOptimizedStorage) {
1421 "You requested optimized storage by setting the" 1422 "\"Optimize Storage\" flag to \"true\" in the parameter list, or by virtue" 1423 "of default behavior. However, the associated CrsGraph was filled separately" 1424 "and requested not to optimize storage. Therefore, the CrsMatrix cannot" 1425 "optimize storage.");
1426 requestOptimizedStorage =
false;
1429 typedef decltype (staticGraph_->k_numRowEntries_) row_entries_type;
1454 size_t lclTotalNumEntries = 0;
1456 typename non_const_row_map_type::HostMirror h_ptrs;
1458 typename row_entries_type::const_type numRowEnt_h =
1459 staticGraph_->k_numRowEntries_;
1461 non_const_row_map_type packedRowOffsets (
"Tpetra::CrsGraph::ptr",
1465 h_ptrs = create_mirror_view (packedRowOffsets);
1469 k_ptrs = packedRowOffsets;
1472 #ifdef HAVE_TPETRA_DEBUG 1473 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1474 (static_cast<size_t> (k_ptrs.dimension_0 ()) != lclNumRows + 1,
1475 std::logic_error,
"In DynamicProfile branch, after packing k_ptrs, " 1476 "k_ptrs.dimension_0() = " << k_ptrs.dimension_0 () <<
" != " 1477 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
1478 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1479 (static_cast<size_t> (h_ptrs.dimension_0 ()) != lclNumRows + 1,
1480 std::logic_error,
"In DynamicProfile branch, after packing h_ptrs, " 1481 "h_ptrs.dimension_0() = " << h_ptrs.dimension_0 () <<
" != " 1482 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
1487 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, lclNumRows);
1488 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
1490 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1491 (static_cast<size_t> (k_ptrs_ent_h ()) != lclTotalNumEntries,
1492 std::logic_error,
"(DynamicProfile branch) After packing k_ptrs, " 1493 "k_ptrs(lclNumRows = " << lclNumRows <<
") = " << k_ptrs_ent_h ()
1494 <<
" != total number of entries on the calling process = " 1495 << lclTotalNumEntries <<
".");
1497 #endif // HAVE_TPETRA_DEBUG 1500 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1502 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
1504 for (
size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) {
1505 const size_t numEnt = numRowEnt_h(lclRow);
1506 std::copy (values2D_[lclRow].begin(),
1507 values2D_[lclRow].begin() + numEnt,
1508 h_vals.ptr_on_device() + h_ptrs(lclRow));
1513 #ifdef HAVE_TPETRA_DEBUG 1515 if (k_ptrs.dimension_0 () != 0) {
1516 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
1521 auto k_ptrs_ent_d = Kokkos::subview (k_ptrs, numOffsets - 1);
1522 auto k_ptrs_ent_h = create_mirror_view (k_ptrs_ent_d);
1524 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1525 (static_cast<size_t> (k_ptrs_ent_h ()) != k_vals.dimension_0 (),
1526 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs(" 1527 << (numOffsets-1) <<
") = " << k_ptrs_ent_h () <<
" != " 1528 "k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1530 #endif // HAVE_TPETRA_DEBUG 1550 if (nodeNumEntries != nodeNumAllocated) {
1553 non_const_row_map_type tmpk_ptrs (
"Tpetra::CrsGraph::ptr",
1558 size_t lclTotalNumEntries = 0;
1561 typename row_entries_type::const_type numRowEnt_d =
1562 staticGraph_->k_numRowEntries_;
1570 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1573 typedef pack_functor<values_type, row_map_type> packer_type;
1574 packer_type valsPacker (k_vals, k_values1D_, tmpk_ptrs, k_rowPtrs_);
1575 Kokkos::parallel_for (lclNumRows, valsPacker);
1578 k_vals = k_values1D_;
1583 if (requestOptimizedStorage) {
1587 k_values1D_ = k_vals;
1588 this->storageStatus_ = Details::STORAGE_1D_PACKED;
1596 getColMap ()->getNodeNumElements (),
1598 staticGraph_->getLocalGraph ());
1601 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1605 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1606 const Teuchos::ArrayView<const Scalar>& values)
1608 using Teuchos::Array;
1609 using Teuchos::ArrayView;
1610 using Teuchos::av_reinterpret_cast;
1611 using Teuchos::toString;
1613 const char tfecfFuncName[] =
"insertLocalValues";
1615 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillActive (), std::runtime_error,
1616 ": Fill is not active. After calling fillComplete, you must call " 1617 "resumeFill before you may insert entries into the matrix again.");
1618 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isStaticGraph (), std::runtime_error,
1619 " cannot insert indices with static graph; use replaceLocalValues() instead.");
1620 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(myGraph_->isGloballyIndexed(),
1621 std::runtime_error,
": graph indices are global; use insertGlobalValues().");
1622 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! hasColMap (), std::runtime_error,
1623 " cannot insert local indices without a column map.");
1624 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
1625 std::runtime_error,
": values.size() must equal indices.size().");
1626 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1627 ! getRowMap()->isNodeLocalElement(localRow), std::runtime_error,
1628 ": Local row index " << localRow <<
" does not belong to this process.");
1630 if (! myGraph_->indicesAreAllocated ()) {
1632 allocateValues (LocalIndices, GraphNotYetAllocated);
1634 catch (std::exception& e) {
1635 TEUCHOS_TEST_FOR_EXCEPTION(
1636 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1637 "allocateValues(LocalIndices,GraphNotYetAllocated) threw an " 1638 "exception: " << e.what ());
1642 const size_t numEntriesToAdd =
static_cast<size_t> (indices.size ());
1643 #ifdef HAVE_TPETRA_DEBUG 1649 const map_type& colMap = * (getColMap ());
1650 Array<LocalOrdinal> badColInds;
1651 bool allInColMap =
true;
1652 for (
size_t k = 0; k < numEntriesToAdd; ++k) {
1654 allInColMap =
false;
1655 badColInds.push_back (indices[k]);
1658 if (! allInColMap) {
1659 std::ostringstream os;
1660 os <<
"Tpetra::CrsMatrix::insertLocalValues: You attempted to insert " 1661 "entries in owned row " << localRow <<
", at the following column " 1662 "indices: " << toString (indices) <<
"." << endl;
1663 os <<
"Of those, the following indices are not in the column Map on " 1664 "this process: " << toString (badColInds) <<
"." << endl <<
"Since " 1665 "the matrix has a column Map already, it is invalid to insert " 1666 "entries at those locations.";
1667 TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
1670 #endif // HAVE_TPETRA_DEBUG 1672 #ifdef HAVE_TPETRA_DEBUG 1675 rowInfo = myGraph_->getRowInfo (localRow);
1676 }
catch (std::exception& e) {
1677 TEUCHOS_TEST_FOR_EXCEPTION(
1678 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1679 "myGraph_->getRowInfo threw an exception: " << e.what ());
1682 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1683 #endif // HAVE_TPETRA_DEBUG 1685 const size_t curNumEntries = rowInfo.numEntries;
1686 const size_t newNumEntries = curNumEntries + numEntriesToAdd;
1687 if (newNumEntries > rowInfo.allocSize) {
1688 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1690 ": new indices exceed statically allocated graph structure.");
1694 rowInfo = myGraph_->template updateLocalAllocAndValues<impl_scalar_type> (rowInfo,
1696 values2D_[localRow]);
1697 }
catch (std::exception& e) {
1698 TEUCHOS_TEST_FOR_EXCEPTION(
1699 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1700 "myGraph_->updateGlobalAllocAndValues threw an exception: " 1704 typename Graph::SLocalGlobalViews indsView;
1705 indsView.linds = indices;
1707 #ifdef HAVE_TPETRA_DEBUG 1708 ArrayView<impl_scalar_type> valsView;
1710 valsView = this->getViewNonConst (rowInfo);
1711 }
catch (std::exception& e) {
1712 TEUCHOS_TEST_FOR_EXCEPTION(
1713 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1714 "getViewNonConst threw an exception: " << e.what ());
1717 ArrayView<impl_scalar_type> valsView = this->getViewNonConst (rowInfo);
1718 #endif // HAVE_TPETRA_DEBUG 1720 ArrayView<const impl_scalar_type> valsIn =
1723 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, indsView,
1727 }
catch (std::exception& e) {
1728 TEUCHOS_TEST_FOR_EXCEPTION(
1729 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1730 "myGraph_->insertIndicesAndValues threw an exception: " 1734 #ifdef HAVE_TPETRA_DEBUG 1735 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
1736 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1737 chkNewNumEntries != newNumEntries, std::logic_error,
1738 ": The row should have " << newNumEntries <<
" entries after insert, but " 1739 "instead has " << chkNewNumEntries <<
". Please report this bug to the " 1740 "Tpetra developers.");
1741 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isLocallyIndexed(), std::logic_error,
1742 ": At end of insertLocalValues(), this CrsMatrix is not locally indexed. " 1743 "Please report this bug to the Tpetra developers.");
1744 #endif // HAVE_TPETRA_DEBUG 1747 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1751 const LocalOrdinal numEnt,
1752 const Scalar vals[],
1753 const LocalOrdinal cols[])
1755 Teuchos::ArrayView<const LocalOrdinal> colsT (cols, numEnt);
1756 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
1757 this->insertLocalValues (localRow, colsT, valsT);
1760 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1764 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1765 const Teuchos::ArrayView<const Scalar>& values)
1767 using Teuchos::Array;
1768 using Teuchos::ArrayView;
1769 using Teuchos::av_reinterpret_cast;
1770 const char tfecfFuncName[] =
"insertLocalValues: ";
1772 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillActive (), std::runtime_error,
1773 "Requires that fill is active.");
1774 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isStaticGraph (), std::runtime_error,
1775 "Cannot insert indices with static graph; use replaceLocalValues() instead.");
1776 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(myGraph_->isGloballyIndexed(),
1777 std::runtime_error,
"Graph indices are global; use insertGlobalValues().");
1778 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1779 ! hasColMap (), std::runtime_error,
"The matrix has no column Map yet, " 1780 "so you cannot insert local indices. If you created the matrix without " 1781 "a column Map (or without a fill-complete graph), you must call " 1782 "fillComplete to create the column Map, before you may work with local " 1784 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1785 values.size () != indices.size (), std::runtime_error,
"values.size() = " 1786 << values.size () <<
" != indices.size() = " << indices.size ()<<
".");
1787 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1788 ! getRowMap()->isNodeLocalElement (localRow), std::runtime_error,
1789 "Local row index " << localRow <<
" does not belong to this process.");
1790 if (! myGraph_->indicesAreAllocated ()) {
1791 allocateValues (LocalIndices, GraphNotYetAllocated);
1795 Array<LocalOrdinal> f_inds (indices);
1796 ArrayView<const impl_scalar_type> valsIn =
1797 av_reinterpret_cast<
const impl_scalar_type> (values);
1798 Array<impl_scalar_type> f_vals (valsIn);
1799 const size_t numFilteredEntries =
1800 myGraph_->template filterLocalIndicesAndValues<impl_scalar_type> (f_inds (),
1802 if (numFilteredEntries > 0) {
1803 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1804 const size_t curNumEntries = rowInfo.numEntries;
1805 const size_t newNumEntries = curNumEntries + numFilteredEntries;
1806 if (newNumEntries > rowInfo.allocSize) {
1807 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1809 ": new indices exceed statically allocated graph structure. " 1810 "newNumEntries (" << newNumEntries <<
" > rowInfo.allocSize (" 1811 << rowInfo.allocSize <<
").");
1814 myGraph_->template updateLocalAllocAndValues<impl_scalar_type> (rowInfo,
1816 values2D_[localRow]);
1818 typename Graph::SLocalGlobalViews inds_view;
1819 inds_view.linds = f_inds (0, numFilteredEntries);
1820 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
1821 this->getViewNonConst (rowInfo),
1822 f_vals, LocalIndices,
1824 #ifdef HAVE_TPETRA_DEBUG 1825 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
1826 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
1827 std::logic_error,
": Internal logic error. Please contact Tpetra team.");
1828 #endif // HAVE_TPETRA_DEBUG 1830 #ifdef HAVE_TPETRA_DEBUG 1831 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isLocallyIndexed(), std::logic_error,
1832 ": At end of insertLocalValues(), this CrsMatrix is not locally indexed. " 1833 "Please report this bug to the Tpetra developers.");
1834 #endif // HAVE_TPETRA_DEBUG 1838 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1842 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
1843 const Teuchos::ArrayView<const Scalar>& values)
1845 using Teuchos::Array;
1846 using Teuchos::ArrayView;
1847 using Teuchos::av_reinterpret_cast;
1848 using Teuchos::toString;
1850 typedef LocalOrdinal LO;
1851 typedef GlobalOrdinal GO;
1852 typedef typename ArrayView<const GO>::size_type size_type;
1853 const char tfecfFuncName[] =
"insertGlobalValues: ";
1855 #ifdef HAVE_TPETRA_DEBUG 1856 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1857 values.size () != indices.size (), std::runtime_error,
1858 "values.size() = " << values.size() <<
" != indices.size() = " 1859 << indices.size() <<
".");
1860 #endif // HAVE_TPETRA_DEBUG 1862 const LO localRow = getRowMap ()->getLocalElement (globalRow);
1864 if (localRow == OTL::invalid ()) {
1865 insertNonownedGlobalValues (globalRow, indices, values);
1868 if (this->isStaticGraph ()) {
1870 std::ostringstream err;
1871 const int myRank = getRowMap ()->getComm ()->getRank ();
1872 const int numProcs = getRowMap ()->getComm ()->getSize ();
1874 err <<
"The matrix was constructed with a constant (\"static\") graph, " 1875 "yet the given global row index " << globalRow <<
" is in the row " 1876 "Map on the calling process (with rank " << myRank <<
", of " <<
1877 numProcs <<
" process(es)). In this case, you may not insert new " 1878 "entries into rows owned by the calling process.";
1880 if (! getRowMap ()->isNodeGlobalElement (globalRow)) {
1881 err <<
" Furthermore, GID->LID conversion with the row Map claims that " 1882 "the global row index is owned on the calling process, yet " 1883 "getRowMap()->isNodeGlobalElement(globalRow) returns false. That's" 1884 " weird! This might indicate a Map bug. Please report this to the" 1885 " Tpetra developers.";
1887 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1888 this->isStaticGraph (), std::runtime_error, err.str ());
1891 if (! myGraph_->indicesAreAllocated ()) {
1893 allocateValues (GlobalIndices, GraphNotYetAllocated);
1895 catch (std::exception& e) {
1896 TEUCHOS_TEST_FOR_EXCEPTION(
1897 true, std::runtime_error,
"Tpetra::CrsMatrix::insertGlobalValues: " 1898 "allocateValues(GlobalIndices,GraphNotYetAllocated) threw an " 1899 "exception: " << e.what ());
1903 const size_type numEntriesToInsert = indices.size ();
1911 const map_type& colMap = * (getColMap ());
1916 #ifdef HAVE_TPETRA_DEBUG 1917 Array<GO> badColInds;
1918 #endif // HAVE_TPETRA_DEBUG 1919 bool allInColMap =
true;
1920 for (size_type k = 0; k < numEntriesToInsert; ++k) {
1922 allInColMap =
false;
1923 #ifdef HAVE_TPETRA_DEBUG 1924 badColInds.push_back (indices[k]);
1927 #endif // HAVE_TPETRA_DEBUG 1930 if (! allInColMap) {
1931 std::ostringstream os;
1932 os <<
"You attempted to insert entries in owned row " << globalRow
1933 <<
", at the following column indices: " << toString (indices)
1935 #ifdef HAVE_TPETRA_DEBUG 1936 os <<
"Of those, the following indices are not in the column Map on " 1937 "this process: " << toString (badColInds) <<
"." << endl <<
"Since " 1938 "the matrix has a column Map already, it is invalid to insert " 1939 "entries at those locations.";
1941 os <<
"At least one of those indices is not in the column Map on this " 1942 "process." << endl <<
"It is invalid to insert into columns not in " 1943 "the column Map on the process that owns the row.";
1944 #endif // HAVE_TPETRA_DEBUG 1945 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1946 ! allInColMap, std::invalid_argument, os.str ());
1950 typename Graph::SLocalGlobalViews inds_view;
1951 ArrayView<const impl_scalar_type> vals_view;
1953 inds_view.ginds = indices;
1956 #ifdef HAVE_TPETRA_DEBUG 1959 rowInfo = myGraph_->getRowInfo (localRow);
1960 }
catch (std::exception& e) {
1961 TEUCHOS_TEST_FOR_EXCEPTION(
1962 true, std::runtime_error,
"myGraph_->getRowInfo(localRow=" << localRow
1963 <<
") threw an exception: " << e.what ());
1966 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1967 #endif // HAVE_TPETRA_DEBUG 1969 const size_t curNumEntries = rowInfo.numEntries;
1970 const size_t newNumEntries =
1971 curNumEntries +
static_cast<size_t> (numEntriesToInsert);
1972 if (newNumEntries > rowInfo.allocSize) {
1973 TEUCHOS_TEST_FOR_EXCEPTION(
1974 getProfileType () ==
StaticProfile && newNumEntries > rowInfo.allocSize,
1975 std::runtime_error,
"Tpetra::CrsMatrix::insertGlobalValues: new " 1976 "indices exceed statically allocated graph structure. curNumEntries" 1977 " (" << curNumEntries <<
") + numEntriesToInsert (" <<
1978 numEntriesToInsert <<
") > allocSize (" << rowInfo.allocSize <<
").");
1983 myGraph_->template updateGlobalAllocAndValues<impl_scalar_type> (rowInfo,
1985 values2D_[localRow]);
1986 }
catch (std::exception& e) {
1987 TEUCHOS_TEST_FOR_EXCEPTION(
1988 true, std::runtime_error,
"myGraph_->updateGlobalAllocAndValues" 1989 "(...) threw an exception: " << e.what ());
1993 if (isGloballyIndexed ()) {
1997 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
1998 this->getViewNonConst (rowInfo),
2000 GlobalIndices, GlobalIndices);
2006 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
2007 this->getViewNonConst (rowInfo),
2009 GlobalIndices, LocalIndices);
2012 catch (std::exception& e) {
2013 TEUCHOS_TEST_FOR_EXCEPTION(
2014 true, std::runtime_error,
"myGraph_->insertIndicesAndValues(...) " 2015 "threw an exception: " << e.what ());
2018 #ifdef HAVE_TPETRA_DEBUG 2019 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
2020 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
2021 std::logic_error,
": There should be a total of " << newNumEntries
2022 <<
" entries in the row, but the graph now reports " << chkNewNumEntries
2023 <<
" entries. Please report this bug to the Tpetra developers.");
2024 #endif // HAVE_TPETRA_DEBUG 2029 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2033 const LocalOrdinal numEnt,
2034 const Scalar vals[],
2035 const GlobalOrdinal inds[])
2037 Teuchos::ArrayView<const GlobalOrdinal> indsT (inds, numEnt);
2038 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
2039 this->insertGlobalValues (globalRow, indsT, valsT);
2043 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2047 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2048 const Teuchos::ArrayView<const Scalar>& values)
2050 using Teuchos::Array;
2051 using Teuchos::ArrayView;
2052 using Teuchos::av_reinterpret_cast;
2053 typedef LocalOrdinal LO;
2054 typedef GlobalOrdinal GO;
2055 typedef impl_scalar_type ST;
2056 const char tfecfFuncName[] =
"insertGlobalValuesFiltered: ";
2066 #ifdef HAVE_TPETRA_DEBUG 2067 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2068 values.size () != indices.size (), std::runtime_error,
2069 "values.size() = " << values.size() <<
" != indices.size() = " 2070 << indices.size() <<
".");
2071 #endif // HAVE_TPETRA_DEBUG 2073 ArrayView<const ST> valsIn = av_reinterpret_cast<
const ST> (values);
2074 const LO lrow = getRowMap ()->getLocalElement (globalRow);
2076 if (lrow != Teuchos::OrdinalTraits<LO>::invalid ()) {
2079 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2080 this->isStaticGraph (), std::runtime_error,
2081 "The matrix was constructed with a static graph. In that case, " 2082 "it is forbidden to insert new entries into rows owned by the " 2083 "calling process.");
2084 if (! myGraph_->indicesAreAllocated ()) {
2085 allocateValues (GlobalIndices, GraphNotYetAllocated);
2087 typename Graph::SLocalGlobalViews inds_view;
2088 ArrayView<const ST> vals_view;
2093 Array<GO> filtered_indices;
2094 Array<ST> filtered_values;
2098 filtered_indices.assign (indices.begin (), indices.end ());
2099 filtered_values.assign (valsIn.begin (), valsIn.end ());
2100 const size_t numFilteredEntries =
2101 myGraph_->template filterGlobalIndicesAndValues<ST> (filtered_indices (),
2102 filtered_values ());
2103 inds_view.ginds = filtered_indices (0, numFilteredEntries);
2104 vals_view = filtered_values (0, numFilteredEntries);
2107 inds_view.ginds = indices;
2110 const size_t numFilteredEntries = vals_view.size ();
2112 if (numFilteredEntries > 0) {
2113 RowInfo rowInfo = myGraph_->getRowInfo (lrow);
2114 const size_t curNumEntries = rowInfo.numEntries;
2115 const size_t newNumEntries = curNumEntries + numFilteredEntries;
2116 if (newNumEntries > rowInfo.allocSize) {
2117 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2119 "New indices exceed statically allocated graph structure.");
2122 rowInfo = myGraph_->template updateGlobalAllocAndValues<ST> (rowInfo,
2126 if (isGloballyIndexed ()) {
2130 myGraph_->template insertIndicesAndValues<ST> (rowInfo, inds_view,
2131 this->getViewNonConst (rowInfo),
2133 GlobalIndices, GlobalIndices);
2139 myGraph_->template insertIndicesAndValues<ST> (rowInfo, inds_view,
2140 this->getViewNonConst (rowInfo),
2142 GlobalIndices, LocalIndices);
2144 #ifdef HAVE_TPETRA_DEBUG 2146 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (lrow);
2147 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
2148 std::logic_error,
": There should be a total of " << newNumEntries
2149 <<
" entries in the row, but the graph now reports " << chkNewNumEntries
2150 <<
" entries. Please report this bug to the Tpetra developers.");
2152 #endif // HAVE_TPETRA_DEBUG 2156 insertNonownedGlobalValues (globalRow, indices, values);
2161 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2165 const Teuchos::ArrayView<const LocalOrdinal>& lclCols,
2166 const Teuchos::ArrayView<const Scalar>& vals)
const 2168 using Kokkos::MemoryUnmanaged;
2171 typedef LocalOrdinal LO;
2173 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2175 typedef View<const IST*, HD, MemoryUnmanaged> ISVT;
2176 typedef View<const LO*, HD, MemoryUnmanaged> LIVT;
2178 LIVT lclColsIn (lclCols.getRawPtr (), lclCols.size ());
2179 const IST* valsRaw =
reinterpret_cast<const IST*
> (vals.getRawPtr ());
2180 ISVT valsIn (valsRaw, vals.size ());
2181 return this->
template replaceLocalValues<LIVT, ISVT> (localRow,
2186 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2190 const LocalOrdinal numEnt,
2191 const Scalar inputVals[],
2192 const LocalOrdinal inputCols[])
const 2194 using Kokkos::MemoryUnmanaged;
2197 typedef LocalOrdinal LO;
2199 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2201 typedef View<const LO*, HD, MemoryUnmanaged> LIVT;
2202 typedef View<const IST*, HD, MemoryUnmanaged> ISVT;
2204 LIVT indsK (inputCols, numEnt);
2205 ISVT valsK (reinterpret_cast<const IST*> (inputVals), numEnt);
2206 return this->
template replaceLocalValues<LIVT, ISVT> (localRow,
2211 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2215 const Teuchos::ArrayView<const GlobalOrdinal>& inputInds,
2216 const Teuchos::ArrayView<const Scalar>& inputVals)
const 2218 using Kokkos::MemoryUnmanaged;
2221 typedef GlobalOrdinal GO;
2223 typedef typename View<GO*, DD>::HostMirror::device_type HD;
2225 typedef View<const GO*, HD, MemoryUnmanaged> GIVT;
2226 typedef View<const IST*, HD, MemoryUnmanaged> ISVT;
2228 const IST* inputValsRaw =
2229 reinterpret_cast<const IST*
> (inputVals.getRawPtr ());
2230 GIVT indsK (inputInds.getRawPtr (), inputInds.size ());
2231 ISVT valsK (inputValsRaw, inputVals.size ());
2232 return this->
template replaceGlobalValues<GIVT, ISVT> (globalRow,
2238 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2242 const LocalOrdinal numEnt,
2243 const Scalar inputVals[],
2244 const GlobalOrdinal inputCols[])
const 2246 using Kokkos::MemoryUnmanaged;
2249 typedef GlobalOrdinal GO;
2251 typedef typename View<GO*, DD>::HostMirror::device_type HD;
2253 typedef View<const GO*, HD, MemoryUnmanaged> GIVT;
2254 typedef View<const IST*, HD, MemoryUnmanaged> ISVT;
2256 GIVT indsK (inputCols, numEnt);
2257 ISVT valsK (reinterpret_cast<const IST*> (inputVals), numEnt);
2258 return this->
template replaceGlobalValues<GIVT, ISVT> (globalRow,
2264 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2268 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2269 const Teuchos::ArrayView<const Scalar>& values,
2272 using Kokkos::MemoryUnmanaged;
2275 typedef LocalOrdinal LO;
2276 typedef GlobalOrdinal GO;
2278 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2280 if (! isFillActive ()) {
2282 return Teuchos::OrdinalTraits<LO>::invalid ();
2289 const LO lrow = this->staticGraph_.is_null () ?
2290 myGraph_->rowMap_->getLocalElement (globalRow) :
2291 staticGraph_->rowMap_->getLocalElement (globalRow);
2294 if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) {
2299 this->insertNonownedGlobalValues (globalRow, indices, values);
2306 return static_cast<LO
> (indices.size ());
2309 if (staticGraph_.is_null ()) {
2310 return Teuchos::OrdinalTraits<LO>::invalid ();
2312 const RowInfo rowInfo = this->staticGraph_->getRowInfo (lrow);
2314 auto curVals = this->getRowViewNonConst (rowInfo);
2315 const ST* valsRaw =
reinterpret_cast<const ST*
> (values.getRawPtr ());
2316 View<const ST*, HD, MemoryUnmanaged> valsIn (valsRaw, values.size ());
2317 View<const GO*, HD, MemoryUnmanaged> indsIn (indices.getRawPtr (),
2319 return staticGraph_->template sumIntoGlobalValues<ST, HD, DD> (rowInfo,
2327 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2331 const LocalOrdinal numEnt,
2332 const Scalar vals[],
2333 const GlobalOrdinal cols[],
2336 Teuchos::ArrayView<const GlobalOrdinal> colsIn (cols, numEnt);
2337 Teuchos::ArrayView<const Scalar> valsIn (vals, numEnt);
2338 return this->sumIntoGlobalValues (globalRow, colsIn, valsIn, atomic);
2342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2346 const Teuchos::ArrayView<const LocalOrdinal>& indices,
2347 const Teuchos::ArrayView<const Scalar>& values,
2348 const bool atomic)
const 2350 using Kokkos::MemoryUnmanaged;
2353 typedef LocalOrdinal LO;
2355 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2356 typedef View<const IST*, HD, MemoryUnmanaged> IVT;
2357 typedef View<const LO*, HD, MemoryUnmanaged> IIT;
2359 const IST* valsRaw =
reinterpret_cast<const IST*
> (values.getRawPtr ());
2360 IVT valsIn (valsRaw, values.size ());
2361 IIT indsIn (indices.getRawPtr (), indices.size ());
2362 return this->
template sumIntoLocalValues<IIT, IVT> (localRow, indsIn,
2366 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2370 const LocalOrdinal numEnt,
2371 const Scalar vals[],
2372 const LocalOrdinal cols[],
2373 const bool atomic)
const 2375 using Kokkos::MemoryUnmanaged;
2378 typedef LocalOrdinal LO;
2380 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2381 typedef View<const IST*, HD, MemoryUnmanaged> IVT;
2382 typedef View<const LO*, HD, MemoryUnmanaged> IIT;
2384 const IST* valsRaw =
reinterpret_cast<const IST*
> (vals);
2385 IVT valsIn (valsRaw, numEnt);
2386 IIT indsIn (cols, numEnt);
2387 return this->
template sumIntoLocalValues<IIT, IVT> (localRow, indsIn,
2392 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2393 Teuchos::ArrayView<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type>
2397 using Kokkos::MemoryUnmanaged;
2399 using Teuchos::ArrayView;
2401 typedef std::pair<size_t, size_t> range_type;
2403 if (k_values1D_.dimension_0 () != 0 && rowinfo.allocSize > 0) {
2404 #ifdef HAVE_TPETRA_DEBUG 2405 TEUCHOS_TEST_FOR_EXCEPTION(
2406 rowinfo.offset1D + rowinfo.allocSize > k_values1D_.dimension_0 (),
2407 std::range_error,
"Tpetra::CrsMatrix::getView: Invalid access " 2408 "to 1-D storage of values." << std::endl <<
"rowinfo.offset1D (" <<
2409 rowinfo.offset1D <<
") + rowinfo.allocSize (" << rowinfo.allocSize <<
2410 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2411 #endif // HAVE_TPETRA_DEBUG 2412 range_type range (rowinfo.offset1D, rowinfo.offset1D + rowinfo.allocSize);
2413 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
2420 subview_type sv = Kokkos::subview (subview_type (k_values1D_), range);
2421 const ST*
const sv_raw = (rowinfo.allocSize == 0) ? NULL : sv.ptr_on_device ();
2422 return ArrayView<const ST> (sv_raw, rowinfo.allocSize);
2424 else if (values2D_ != Teuchos::null) {
2425 return values2D_[rowinfo.localRow] ();
2428 return ArrayView<impl_scalar_type> ();
2433 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2437 LocalOrdinal& numEnt,
2440 if (k_values1D_.dimension_0 () != 0 && rowinfo.allocSize > 0) {
2441 #ifdef HAVE_TPETRA_DEBUG 2442 if (rowinfo.offset1D + rowinfo.allocSize > k_values1D_.dimension_0 ()) {
2445 return Teuchos::OrdinalTraits<LocalOrdinal>::invalid ();
2447 #endif // HAVE_TPETRA_DEBUG 2448 vals = k_values1D_.ptr_on_device () + rowinfo.offset1D;
2449 numEnt = rowinfo.allocSize;
2451 else if (! values2D_.is_null ()) {
2452 #ifdef HAVE_TPETRA_DEBUG 2453 if (rowinfo.localRow >= static_cast<size_t> (values2D_.size ())) {
2456 return Teuchos::OrdinalTraits<LocalOrdinal>::invalid ();
2458 #endif // HAVE_TPETRA_DEBUG 2461 const auto& curRow = values2D_[rowinfo.localRow];
2462 vals = curRow.getRawPtr ();
2463 numEnt = curRow.size ();
2470 return static_cast<LocalOrdinal
> (0);
2473 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2477 LocalOrdinal& numEnt,
2481 const LocalOrdinal err = this->getViewRawConst (valsConst, numEnt, rowinfo);
2486 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2487 Kokkos::View<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type*,
2489 Kokkos::MemoryUnmanaged>
2493 using Kokkos::MemoryUnmanaged;
2495 typedef impl_scalar_type ST;
2496 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
2497 typedef std::pair<size_t, size_t> range_type;
2499 if (k_values1D_.dimension_0 () != 0 && rowInfo.allocSize > 0) {
2500 #ifdef HAVE_TPETRA_DEBUG 2501 TEUCHOS_TEST_FOR_EXCEPTION(
2502 rowInfo.offset1D + rowInfo.allocSize > k_values1D_.dimension_0 (),
2503 std::range_error,
"Tpetra::CrsMatrix::getRowView: Invalid access " 2504 "to 1-D storage of values." << std::endl <<
"rowInfo.offset1D (" <<
2505 rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize <<
2506 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2507 #endif // HAVE_TPETRA_DEBUG 2508 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
2515 return Kokkos::subview (subview_type (k_values1D_), range);
2517 else if (values2D_ != Teuchos::null) {
2518 Teuchos::ArrayView<const ST> rowView = values2D_[rowInfo.localRow] ();
2519 return subview_type (rowView.getRawPtr (), rowView.size ());
2522 return subview_type ();
2526 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2527 Kokkos::View<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type*,
2529 Kokkos::MemoryUnmanaged>
2530 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
2531 getRowViewNonConst (
const RowInfo& rowInfo)
const 2533 using Kokkos::MemoryUnmanaged;
2535 typedef impl_scalar_type ST;
2536 typedef View<ST*, execution_space, MemoryUnmanaged> subview_type;
2537 typedef std::pair<size_t, size_t> range_type;
2539 if (k_values1D_.dimension_0 () != 0 && rowInfo.allocSize > 0) {
2540 #ifdef HAVE_TPETRA_DEBUG 2541 TEUCHOS_TEST_FOR_EXCEPTION(
2542 rowInfo.offset1D + rowInfo.allocSize > k_values1D_.dimension_0 (),
2543 std::range_error,
"Tpetra::CrsMatrix::getRowViewNonConst: Invalid access " 2544 "to 1-D storage of values." << std::endl <<
"rowInfo.offset1D (" <<
2545 rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize <<
2546 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2547 #endif // HAVE_TPETRA_DEBUG 2548 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
2555 return Kokkos::subview (subview_type (k_values1D_), range);
2557 else if (values2D_ != Teuchos::null) {
2558 Teuchos::ArrayView<ST> rowView = values2D_[rowInfo.localRow] ();
2559 return subview_type (rowView.getRawPtr (), rowView.size ());
2562 return subview_type ();
2566 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2567 Teuchos::ArrayView<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type>
2574 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2578 const Teuchos::ArrayView<LocalOrdinal>& indices,
2579 const Teuchos::ArrayView<Scalar>& values,
2580 size_t& numEntries)
const 2582 using Teuchos::ArrayView;
2583 using Teuchos::av_reinterpret_cast;
2584 const char tfecfFuncName[] =
"getLocalRowCopy: ";
2586 TEUCHOS_TEST_FOR_EXCEPTION(
2587 isGloballyIndexed () && ! hasColMap (), std::runtime_error,
2588 "Tpetra::CrsMatrix::getLocalRowCopy: The matrix is globally indexed and " 2589 "does not have a column Map yet. That means we don't have local indices " 2590 "for columns yet, so it doesn't make sense to call this method. If the " 2591 "matrix doesn't have a column Map yet, you should call fillComplete on " 2593 #ifdef HAVE_TPETRA_DEBUG 2594 TEUCHOS_TEST_FOR_EXCEPTION(
2595 ! staticGraph_->hasRowInfo (), std::runtime_error,
2596 "Tpetra::CrsMatrix::getLocalRowCopy: The graph's row information was " 2597 "deleted at fillComplete().");
2598 #endif // HAVE_TPETRA_DEBUG 2600 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
2601 const size_t theNumEntries = rowinfo.numEntries;
2602 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2603 (static_cast<size_t> (indices.size ()) < theNumEntries ||
2604 static_cast<size_t> (values.size ()) < theNumEntries,
2605 std::runtime_error,
"Row with local index " << localRow <<
" has " <<
2606 theNumEntries <<
" entry/ies, but indices.size() = " <<
2607 indices.size () <<
" and values.size() = " << values.size () <<
".");
2608 numEntries = theNumEntries;
2610 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
2611 if (staticGraph_->isLocallyIndexed ()) {
2612 const LocalOrdinal* curLclInds;
2614 LocalOrdinal numSpots;
2619 #ifdef HAVE_TPETRA_DEBUG 2621 staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
2622 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2623 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2624 "staticGraph_->getLocalViewRawConst returned nonzero error code " 2626 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2627 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
2628 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
2630 const LocalOrdinal numSpotsBefore = numSpots;
2631 err = getViewRawConst (curVals, numSpots, rowinfo);
2632 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2633 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2634 "getViewRaw returned nonzero error code " << err <<
".");
2635 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2636 (numSpotsBefore != numSpots, std::logic_error,
2637 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = " 2638 << numSpots <<
".");
2640 (void) staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
2641 (void) getViewRawConst (curVals, numSpots, rowinfo);
2642 #endif // HAVE_TPETRA_DEBUG 2644 for (
size_t j = 0; j < theNumEntries; ++j) {
2645 values[j] = curVals[j];
2646 indices[j] = curLclInds[j];
2649 else if (staticGraph_->isGloballyIndexed ()) {
2650 const map_type& colMap = * (staticGraph_->colMap_);
2651 const GlobalOrdinal* curGblInds;
2653 LocalOrdinal numSpots;
2658 #ifdef HAVE_TPETRA_DEBUG 2660 staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
2661 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2662 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2663 "staticGraph_->getGlobalViewRawConst returned nonzero error code " 2665 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2666 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
2667 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
2669 const LocalOrdinal numSpotsBefore = numSpots;
2670 err = getViewRawConst (curVals, numSpots, rowinfo);
2671 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2672 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2673 "getViewRawConst returned nonzero error code " << err <<
".");
2674 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2675 (numSpotsBefore != numSpots, std::logic_error,
2676 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = " 2677 << numSpots <<
".");
2679 (void) staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
2680 (void) getViewRawConst (curVals, numSpots, rowinfo);
2681 #endif //HAVE_TPETRA_DEBUG 2683 for (
size_t j = 0; j < theNumEntries; ++j) {
2684 values[j] = curVals[j];
2691 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2695 const Teuchos::ArrayView<GlobalOrdinal>& indices,
2696 const Teuchos::ArrayView<Scalar>& values,
2697 size_t& numEntries)
const 2699 using Teuchos::ArrayView;
2700 using Teuchos::av_reinterpret_cast;
2701 const char tfecfFuncName[] =
"getGlobalRowCopy: ";
2704 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
2705 const size_t theNumEntries = rowinfo.numEntries;
2706 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2707 static_cast<size_t> (indices.size ()) < theNumEntries ||
2708 static_cast<size_t> (values.size ()) < theNumEntries,
2709 std::runtime_error,
"Row with global index " << globalRow <<
" has " 2710 << theNumEntries <<
" entry/ies, but indices.size() = " <<
2711 indices.size () <<
" and values.size() = " << values.size () <<
".");
2712 numEntries = theNumEntries;
2714 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
2715 if (staticGraph_->isLocallyIndexed ()) {
2716 const map_type& colMap = * (staticGraph_->colMap_);
2717 const LocalOrdinal* curLclInds;
2719 LocalOrdinal numSpots;
2724 #ifdef HAVE_TPETRA_DEBUG 2726 staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
2727 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2728 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2729 "staticGraph_->getLocalViewRawConst returned nonzero error code " 2731 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2732 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
2733 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
2735 const LocalOrdinal numSpotsBefore = numSpots;
2736 err = getViewRawConst (curVals, numSpots, rowinfo);
2737 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2738 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2739 "getViewRaw returned nonzero error code " << err <<
".");
2740 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2741 (numSpotsBefore != numSpots, std::logic_error,
2742 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = " 2743 << numSpots <<
".");
2745 (void) staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
2746 (void) getViewRawConst (curVals, numSpots, rowinfo);
2747 #endif //HAVE_TPETRA_DEBUG 2749 for (
size_t j = 0; j < theNumEntries; ++j) {
2750 values[j] = curVals[j];
2754 else if (staticGraph_->isGloballyIndexed ()) {
2755 const GlobalOrdinal* curGblInds;
2757 LocalOrdinal numSpots;
2762 #ifdef HAVE_TPETRA_DEBUG 2764 staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
2765 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2766 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2767 "staticGraph_->getGlobalViewRawConst returned nonzero error code " 2769 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2770 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
2771 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
2773 const LocalOrdinal numSpotsBefore = numSpots;
2774 err = getViewRawConst (curVals, numSpots, rowinfo);
2775 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2776 (err != static_cast<LocalOrdinal> (0), std::logic_error,
2777 "getViewRawConst returned nonzero error code " << err <<
".");
2778 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2779 (numSpotsBefore != numSpots, std::logic_error,
2780 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = " 2781 << numSpots <<
".");
2783 (void) staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
2784 (void) getViewRawConst (curVals, numSpots, rowinfo);
2785 #endif //HAVE_TPETRA_DEBUG 2787 for (
size_t j = 0; j < theNumEntries; ++j) {
2788 values[j] = curVals[j];
2789 indices[j] = curGblInds[j];
2795 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2799 Teuchos::ArrayView<const LocalOrdinal>& indices,
2800 Teuchos::ArrayView<const Scalar>& values)
const 2802 using Teuchos::ArrayView;
2803 using Teuchos::av_reinterpret_cast;
2804 typedef LocalOrdinal LO;
2805 const char tfecfFuncName[] =
"getLocalRowView: ";
2807 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2808 isGloballyIndexed (), std::runtime_error,
"The matrix currently stores " 2809 "its indices as global indices, so you cannot get a view with local " 2810 "column indices. If the matrix has a column Map, you may call " 2811 "getLocalRowCopy() to get local column indices; otherwise, you may get " 2812 "a view with global column indices by calling getGlobalRowCopy().");
2813 indices = Teuchos::null;
2814 values = Teuchos::null;
2815 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
2816 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
2817 rowinfo.numEntries > 0) {
2818 ArrayView<const LO> indTmp = staticGraph_->getLocalView (rowinfo);
2819 ArrayView<const Scalar> valTmp =
2820 av_reinterpret_cast<
const Scalar> (this->getView (rowinfo));
2821 indices = indTmp (0, rowinfo.numEntries);
2822 values = valTmp (0, rowinfo.numEntries);
2825 #ifdef HAVE_TPETRA_DEBUG 2826 const char suffix[] =
". This should never happen. Please report this " 2827 "bug to the Tpetra developers.";
2828 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2829 (static_cast<size_t> (indices.size ()) !=
2830 static_cast<size_t> (values.size ()), std::logic_error,
2831 "At the end of this method, for local row " << localRow <<
", " 2832 "indices.size() = " << indices.size () <<
" != values.size () = " 2833 << values.size () << suffix);
2834 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2835 (static_cast<size_t> (indices.size ()) !=
2836 static_cast<size_t> (rowinfo.numEntries), std::logic_error,
2837 "At the end of this method, for local row " << localRow <<
", " 2838 "indices.size() = " << indices.size () <<
" != rowinfo.numEntries = " 2839 << rowinfo.numEntries << suffix);
2840 const size_t expectedNumEntries = getNumEntriesInLocalRow (localRow);
2841 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2842 (rowinfo.numEntries != expectedNumEntries, std::logic_error,
"At the end " 2843 "of this method, for local row " << localRow <<
", rowinfo.numEntries = " 2844 << rowinfo.numEntries <<
" != getNumEntriesInLocalRow(localRow) = " <<
2845 expectedNumEntries << suffix);
2846 #endif // HAVE_TPETRA_DEBUG 2849 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2853 LocalOrdinal& numEnt,
2855 const LocalOrdinal*& ind)
const 2857 typedef LocalOrdinal LO;
2864 if (staticGraph_.is_null () || staticGraph_->isGloballyIndexed ()) {
2868 const RowInfo rowInfo = staticGraph_->getRowInfo (lclRow);
2874 return static_cast<LO
> (1);
2877 numEnt =
static_cast<LO
> (rowInfo.numEntries);
2878 auto lclColInds = staticGraph_->getLocalKokkosRowView (rowInfo);
2879 ind = lclColInds.ptr_on_device ();
2880 const LO err = this->getViewRawConst (val, numEnt, rowInfo);
2886 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2890 LocalOrdinal& numEnt,
2891 const LocalOrdinal*& lclColInds,
2892 const Scalar*& vals)
const 2895 const LocalOrdinal errCode =
2896 this->getLocalRowView (lclRow, numEnt, vals_ist, lclColInds);
2897 vals =
reinterpret_cast<const Scalar*
> (vals_ist);
2901 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2905 Teuchos::ArrayView<const GlobalOrdinal>& indices,
2906 Teuchos::ArrayView<const Scalar>& values)
const 2908 using Teuchos::ArrayView;
2909 using Teuchos::av_reinterpret_cast;
2910 typedef GlobalOrdinal GO;
2911 const char tfecfFuncName[] =
"getGlobalRowView: ";
2913 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2914 isLocallyIndexed (), std::runtime_error,
2915 "The matrix is locally indexed, so we cannot return a view of the row " 2916 "with global column indices. Use getGlobalRowCopy() instead.");
2917 indices = Teuchos::null;
2918 values = Teuchos::null;
2920 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
2921 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
2922 rowinfo.numEntries > 0) {
2923 ArrayView<const GO> indTmp = staticGraph_->getGlobalView (rowinfo);
2924 ArrayView<const Scalar> valTmp =
2925 av_reinterpret_cast<
const Scalar> (this->getView (rowinfo));
2926 indices = indTmp (0, rowinfo.numEntries);
2927 values = valTmp (0, rowinfo.numEntries);
2930 #ifdef HAVE_TPETRA_DEBUG 2931 const char suffix[] =
". This should never happen. Please report this " 2932 "bug to the Tpetra developers.";
2933 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2934 (static_cast<size_t> (indices.size ()) !=
2935 static_cast<size_t> (values.size ()), std::logic_error,
2936 "At the end of this method, for global row " << globalRow <<
", " 2937 "indices.size() = " << indices.size () <<
" != values.size () = " 2938 << values.size () << suffix);
2939 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2940 (static_cast<size_t> (indices.size ()) !=
2941 static_cast<size_t> (rowinfo.numEntries), std::logic_error,
2942 "At the end of this method, for global row " << globalRow <<
", " 2943 "indices.size() = " << indices.size () <<
" != rowinfo.numEntries = " 2944 << rowinfo.numEntries << suffix);
2945 const size_t expectedNumEntries = getNumEntriesInGlobalRow (globalRow);
2946 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2947 (rowinfo.numEntries != expectedNumEntries, std::logic_error,
"At the end " 2948 "of this method, for global row " << globalRow <<
", rowinfo.numEntries " 2949 "= " << rowinfo.numEntries <<
" != getNumEntriesInGlobalRow(globalRow) =" 2950 " " << expectedNumEntries << suffix);
2951 #endif // HAVE_TPETRA_DEBUG 2954 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2959 typedef LocalOrdinal LO;
2960 typedef typename Teuchos::Array<Scalar>::size_type size_type;
2961 const char tfecfFuncName[] =
"scale: ";
2964 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2965 ! isFillActive (), std::runtime_error,
2966 "Fill must be active before you may call this method. " 2967 "Please call resumeFill() to make fill active.");
2969 const size_t nlrs = staticGraph_->getNodeNumRows ();
2970 const size_t numAlloc = staticGraph_->getNodeAllocationSize ();
2971 const size_t numEntries = staticGraph_->getNodeNumEntries ();
2972 if (! staticGraph_->indicesAreAllocated () || nlrs == 0 ||
2973 numAlloc == 0 || numEntries == 0) {
2978 const LO lclNumRows = lclMatrix_.numRows ();
2979 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
2980 auto row_i = lclMatrix_.row (lclRow);
2981 for (LO k = 0; k < row_i.length; ++k) {
2983 row_i.value (k) *= theAlpha;
2988 for (
size_t row = 0; row < nlrs; ++row) {
2989 const size_type numEnt = getNumEntriesInLocalRow (row);
2990 Teuchos::ArrayView<impl_scalar_type> rowVals = values2D_[row] ();
2991 for (size_type k = 0; k < numEnt; ++k) {
2992 rowVals[k] *= theAlpha;
2999 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3004 const char tfecfFuncName[] =
"setAllToScalar: ";
3006 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3007 ! isFillActive (), std::runtime_error,
3008 "Fill must be active before you may call this method. " 3009 "Please call resumeFill() to make fill active.");
3015 const size_t nlrs = staticGraph_->getNodeNumRows(),
3016 numAlloc = staticGraph_->getNodeAllocationSize(),
3017 numEntries = staticGraph_->getNodeNumEntries();
3018 if (! staticGraph_->indicesAreAllocated () || numAlloc == 0 || numEntries == 0) {
3022 const ProfileType profType = staticGraph_->getProfileType ();
3030 for (
size_t row = 0; row < nlrs; ++row) {
3031 std::fill (values2D_[row].begin (), values2D_[row].end (), theAlpha);
3037 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3041 const typename local_graph_type::entries_type::non_const_type& columnIndices,
3042 const typename local_matrix_type::values_type& values)
3044 const char tfecfFuncName[] =
"setAllValues: ";
3045 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3046 (columnIndices.size () != values.size (), std::invalid_argument,
3047 "columnIndices.size() = " << columnIndices.size () <<
" != values.size()" 3048 " = " << values.size () <<
".");
3049 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3050 (myGraph_.is_null (), std::runtime_error,
"myGraph_ must not be null.");
3053 myGraph_->setAllIndices (rowPointers, columnIndices);
3055 catch (std::exception &e) {
3056 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3057 (
true, std::runtime_error,
"myGraph_->setAllIndices() threw an " 3058 "exception: " << e.what ());
3064 auto lclGraph = myGraph_->getLocalGraph ();
3065 const size_t numEnt = lclGraph.entries.dimension_0 ();
3066 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3067 (lclGraph.row_map.dimension_0 () != rowPointers.dimension_0 () ||
3068 numEnt !=
static_cast<size_t> (columnIndices.dimension_0 ()),
3069 std::logic_error,
"myGraph_->setAllIndices() did not correctly create " 3070 "local graph. Please report this bug to the Tpetra developers.");
3072 const size_t numCols = myGraph_->getColMap ()->getNodeNumElements ();
3074 numCols, values, lclGraph);
3078 this->k_values1D_ = this->lclMatrix_.values;
3080 checkInternalState ();
3083 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3087 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
3088 const Teuchos::ArrayRCP<Scalar>& val)
3090 using Kokkos::Compat::getKokkosViewDeepCopy;
3091 using Teuchos::ArrayRCP;
3092 using Teuchos::av_reinterpret_cast;
3095 typedef typename local_matrix_type::row_map_type row_map_type;
3097 const char tfecfFuncName[] =
"setAllValues(ArrayRCP<size_t>, ArrayRCP<LO>, ArrayRCP<Scalar>): ";
3103 typename row_map_type::non_const_type ptrNative (
"ptr", ptr.size ());
3104 Kokkos::View<
const size_t*,
3105 typename row_map_type::array_layout,
3107 Kokkos::MemoryUnmanaged> ptrSizeT (ptr.getRawPtr (), ptr.size ());
3110 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3111 (ptrNative.dimension_0 () != ptrSizeT.dimension_0 (),
3112 std::logic_error,
"ptrNative.dimension_0() = " <<
3113 ptrNative.dimension_0 () <<
" != ptrSizeT.dimension_0() = " 3114 << ptrSizeT.dimension_0 () <<
". Please report this bug to the " 3115 "Tpetra developers.");
3117 auto indIn = getKokkosViewDeepCopy<DT> (ind ());
3118 auto valIn = getKokkosViewDeepCopy<DT> (av_reinterpret_cast<IST> (val ()));
3119 this->setAllValues (ptrNative, indIn, valIn);
3122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3127 const char tfecfFuncName[] =
"getLocalDiagOffsets: ";
3128 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3129 (staticGraph_.is_null (), std::runtime_error,
"The matrix has no graph.");
3136 const size_t lclNumRows = staticGraph_->getNodeNumRows ();
3137 if (static_cast<size_t> (offsets.size ()) < lclNumRows) {
3138 offsets.resize (lclNumRows);
3144 typedef typename device_type::memory_space memory_space;
3145 if (std::is_same<memory_space, Kokkos::HostSpace>::value) {
3150 Kokkos::MemoryUnmanaged> output_type;
3151 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3152 staticGraph_->getLocalDiagOffsets (offsetsOut);
3155 Kokkos::View<size_t*, device_type> offsetsTmp (
"diagOffsets", lclNumRows);
3156 staticGraph_->getLocalDiagOffsets (offsetsTmp);
3157 typedef Kokkos::View<
size_t*, Kokkos::HostSpace,
3158 Kokkos::MemoryUnmanaged> output_type;
3159 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3164 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3169 using Teuchos::ArrayRCP;
3170 using Teuchos::ArrayView;
3171 using Teuchos::av_reinterpret_cast;
3172 const char tfecfFuncName[] =
"getLocalDiagCopy (1-arg): ";
3176 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3177 staticGraph_.is_null (), std::runtime_error,
3178 "This method requires that the matrix have a graph.");
3179 auto rowMapPtr = this->getRowMap ();
3180 if (rowMapPtr.is_null () || rowMapPtr->getComm ().is_null ()) {
3186 auto colMapPtr = this->getColMap ();
3187 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3188 (! this->hasColMap () || colMapPtr.is_null (), std::runtime_error,
3189 "This method requires that the matrix have a column Map.");
3190 const map_type& rowMap = * rowMapPtr;
3191 const map_type& colMap = * colMapPtr;
3192 const LO myNumRows =
static_cast<LO
> (this->getNodeNumRows ());
3194 #ifdef HAVE_TPETRA_DEBUG 3197 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3198 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3199 "The input Vector's Map must be compatible with the CrsMatrix's row " 3200 "Map. You may check this by using Map's isCompatible method: " 3201 "diag.getMap ()->isCompatible (A.getRowMap ());");
3202 #endif // HAVE_TPETRA_DEBUG 3204 #ifdef HAVE_TPETRA_DEBUG 3207 #endif // HAVE_TPETRA_DEBUG 3208 if (this->isFillComplete ()) {
3209 diag.template modify<device_type> ();
3210 const auto D_lcl = diag.template getLocalView<device_type> ();
3212 const auto D_lcl_1d =
3213 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3215 const auto lclRowMap = rowMap.getLocalMap ();
3217 const auto lclMatrix = this->lclMatrix_;
3219 #ifdef HAVE_TPETRA_DEBUG 3221 lclColMap, lclMatrix);
3224 lclColMap, lclMatrix);
3225 #endif // HAVE_TPETRA_DEBUG 3229 #ifdef HAVE_TPETRA_DEBUG 3233 #endif // HAVE_TPETRA_DEBUG 3236 #ifdef HAVE_TPETRA_DEBUG 3237 if (! this->getComm ().is_null ()) {
3238 using Teuchos::outArg;
3239 using Teuchos::REDUCE_SUM;
3240 using Teuchos::reduceAll;
3244 Teuchos::RCP<const Teuchos::Comm<int> > comm = this->getComm ();
3245 if (! comm.is_null ()) {
3246 reduceAll<int, GO> (*comm, REDUCE_SUM,
static_cast<GO
> (lclNumErrs),
3247 outArg (gblNumErrs));
3249 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3250 (gblNumErrs != 0, std::logic_error,
"Something went wrong on " 3251 << gblNumErrs <<
" out of " << this->getComm ()->getSize ()
3252 <<
" process(es).");
3254 #endif // HAVE_TPETRA_DEBUG 3257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3262 Kokkos::MemoryUnmanaged>& offsets)
const 3264 typedef LocalOrdinal LO;
3266 #ifdef HAVE_TPETRA_DEBUG 3267 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3268 const map_type& rowMap = * (this->getRowMap ());
3271 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3272 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3273 "The input Vector's Map must be compatible with (in the sense of Map::" 3274 "isCompatible) the CrsMatrix's row Map.");
3275 #endif // HAVE_TPETRA_DEBUG 3284 diag.template modify<device_type> ();
3285 auto D_lcl = diag.template getLocalView<device_type> ();
3286 const LO myNumRows =
static_cast<LO
> (this->getNodeNumRows ());
3289 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3291 KokkosSparse::getDiagCopy (D_lcl_1d, offsets, this->lclMatrix_);
3294 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3298 const Teuchos::ArrayView<const size_t>& offsets)
const 3300 typedef LocalOrdinal LO;
3303 typedef typename vec_type::dual_view_type dual_view_type;
3304 typedef typename dual_view_type::host_mirror_space::execution_space host_execution_space;
3306 #ifdef HAVE_TPETRA_DEBUG 3307 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3308 const map_type& rowMap = * (this->getRowMap ());
3311 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3312 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3313 "The input Vector's Map must be compatible with (in the sense of Map::" 3314 "isCompatible) the CrsMatrix's row Map.");
3315 #endif // HAVE_TPETRA_DEBUG 3320 diag.template modify<host_execution_space> ();
3321 auto lclVecHost = diag.template getLocalView<host_execution_space> ();
3323 auto lclVecHost1d = Kokkos::subview (lclVecHost, Kokkos::ALL (), 0);
3325 Kokkos::View<
const size_t*, Kokkos::HostSpace,
3326 Kokkos::MemoryTraits<Kokkos::Unmanaged> >
3327 h_offsets (offsets.getRawPtr (), offsets.size ());
3329 const LO myNumRows =
static_cast<LO
> (this->getNodeNumRows ());
3330 typedef Kokkos::RangePolicy<host_execution_space, LO> policy_type;
3333 Kokkos::parallel_for (policy_type (0, myNumRows), [&] (
const LO& lclRow) {
3334 lclVecHost1d(lclRow) = STS::zero ();
3335 if (h_offsets[lclRow] != INV) {
3336 auto curRow = lclMatrix_.rowConst (lclRow);
3337 lclVecHost1d(lclRow) =
static_cast<IST
> (curRow.value(h_offsets[lclRow]));
3340 diag.template sync<execution_space> ();
3344 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3349 using Teuchos::ArrayRCP;
3350 using Teuchos::ArrayView;
3351 using Teuchos::null;
3354 using Teuchos::rcpFromRef;
3356 const char tfecfFuncName[] =
"leftScale";
3360 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3361 ! isFillComplete (), std::runtime_error,
3362 ": matrix must be fill complete.");
3363 RCP<const vec_type> xp;
3365 if (getRangeMap ()->isSameAs (* (x.
getMap ()))){
3370 if (getCrsGraph ()->getExporter () != Teuchos::null) {
3371 RCP<vec_type> tempVec = rcp (
new vec_type (getRowMap ()));
3372 tempVec->doImport (x, * (getCrsGraph ()->getExporter ()),
INSERT);
3376 xp = rcpFromRef (x);
3379 else if (getRowMap ()->isSameAs (* (x.
getMap ()))) {
3380 xp = rcpFromRef (x);
3383 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::invalid_argument,
": The " 3384 "input scaling vector x's Map must be the same as either the row Map or " 3385 "the range Map of the CrsMatrix.");
3387 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
3388 ArrayView<impl_scalar_type> rowValues = null;
3390 const LocalOrdinal lclNumRows =
3391 static_cast<LocalOrdinal
> (this->getNodeNumRows ());
3392 for (LocalOrdinal i = 0; i < lclNumRows; ++i) {
3393 const RowInfo rowinfo = staticGraph_->getRowInfo (i);
3394 rowValues = this->getViewNonConst (rowinfo);
3396 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
3397 rowValues[j] *= scaleValue;
3402 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3407 using Teuchos::ArrayRCP;
3408 using Teuchos::ArrayView;
3409 using Teuchos::null;
3412 using Teuchos::rcpFromRef;
3414 const char tfecfFuncName[] =
"rightScale: ";
3418 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3419 ! isFillComplete (), std::runtime_error,
"Matrix must be fill complete.");
3420 RCP<const vec_type> xp;
3421 if (getDomainMap ()->isSameAs (* (x.
getMap ()))) {
3425 if (getCrsGraph ()->getImporter () != Teuchos::null) {
3426 RCP<vec_type> tempVec = rcp (
new vec_type (getColMap ()));
3427 tempVec->doImport (x, * (getCrsGraph ()->getImporter ()),
INSERT);
3431 xp = rcpFromRef (x);
3434 else if (getRowMap ()->isSameAs (* (x.
getMap ()))) {
3435 xp = rcpFromRef (x);
3437 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3438 true, std::runtime_error,
"The vector x must have the same Map as " 3439 "either the row Map or the range Map.");
3442 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
3443 ArrayView<impl_scalar_type> rowValues = null;
3445 const LocalOrdinal lclNumRows =
3446 static_cast<LocalOrdinal
> (this->getNodeNumRows ());
3447 for (LocalOrdinal i = 0; i < lclNumRows; ++i) {
3448 const RowInfo rowinfo = staticGraph_->getRowInfo (i);
3449 rowValues = this->getViewNonConst (rowinfo);
3450 ArrayView<const LocalOrdinal> colInds;
3451 getCrsGraph ()->getLocalRowView (i, colInds);
3452 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
3458 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3463 using Teuchos::ArrayView;
3464 using Teuchos::outArg;
3465 using Teuchos::REDUCE_SUM;
3466 using Teuchos::reduceAll;
3467 typedef typename Teuchos::ArrayRCP<const impl_scalar_type>::size_type size_type;
3475 if (frobNorm == -STM::one ()) {
3477 if (getNodeNumEntries() > 0) {
3478 if (isStorageOptimized ()) {
3481 const size_type numEntries =
3482 static_cast<size_type
> (getNodeNumEntries ());
3483 for (size_type k = 0; k < numEntries; ++k) {
3489 const mag_type val_abs = STS::abs (val);
3490 mySum += val_abs * val_abs;
3494 const LocalOrdinal numRows =
3495 static_cast<LocalOrdinal
> (this->getNodeNumRows ());
3496 for (LocalOrdinal r = 0; r < numRows; ++r) {
3497 const RowInfo rowInfo = myGraph_->getRowInfo (r);
3498 const size_type numEntries =
3499 static_cast<size_type
> (rowInfo.numEntries);
3500 ArrayView<const impl_scalar_type> A_r =
3501 this->getView (rowInfo).view (0, numEntries);
3502 for (size_type k = 0; k < numEntries; ++k) {
3504 const mag_type val_abs = STS::abs (val);
3505 mySum += val_abs * val_abs;
3511 reduceAll<int, mag_type> (* (getComm ()), REDUCE_SUM,
3512 mySum, outArg (totalSum));
3513 frobNorm = STM::sqrt (totalSum);
3515 if (isFillComplete ()) {
3519 frobNorm_ = frobNorm;
3524 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3529 const char tfecfFuncName[] =
"replaceColMap: ";
3533 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3534 myGraph_.is_null (), std::runtime_error,
3535 "This method does not work if the matrix has a const graph. The whole " 3536 "idea of a const graph is that you are not allowed to change it, but " 3537 "this method necessarily must modify the graph, since the graph owns " 3538 "the matrix's column Map.");
3539 myGraph_->replaceColMap (newColMap);
3542 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3546 const Teuchos::RCP<const map_type>& newColMap,
3547 const Teuchos::RCP<const import_type>& newImport,
3548 const bool sortEachRow)
3550 const char tfecfFuncName[] =
"reindexColumns: ";
3551 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3552 graph == NULL && myGraph_.is_null (), std::invalid_argument,
3553 "The input graph is NULL, but the matrix does not own its graph.");
3556 const bool sortGraph =
false;
3558 if (sortEachRow && theGraph.isLocallyIndexed () && ! theGraph.isSorted ()) {
3562 const LocalOrdinal lclNumRows =
3563 static_cast<LocalOrdinal
> (theGraph.getNodeNumRows ());
3564 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
3565 const RowInfo rowInfo = theGraph.getRowInfo (row);
3566 Teuchos::ArrayView<impl_scalar_type> rv = this->getViewNonConst (rowInfo);
3567 theGraph.template sortRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
3569 theGraph.indicesAreSorted_ =
true;
3573 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3577 Teuchos::RCP<const import_type>& newImporter)
3579 const char tfecfFuncName[] =
"replaceDomainMapAndImporter: ";
3580 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3581 myGraph_.is_null (), std::runtime_error,
3582 "This method does not work if the matrix has a const graph. The whole " 3583 "idea of a const graph is that you are not allowed to change it, but this" 3584 " method necessarily must modify the graph, since the graph owns the " 3585 "matrix's domain Map and Import objects.");
3586 myGraph_->replaceDomainMapAndImporter (newDomainMap, newImporter);
3589 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3593 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
3594 const Teuchos::ArrayView<const Scalar>& values)
3596 using Teuchos::Array;
3597 typedef GlobalOrdinal GO;
3598 typedef typename Array<GO>::size_type size_type;
3600 const size_type numToInsert = indices.size ();
3603 std::pair<Array<GO>, Array<Scalar> >& curRow = nonlocals_[globalRow];
3604 Array<GO>& curRowInds = curRow.first;
3605 Array<Scalar>& curRowVals = curRow.second;
3606 const size_type newCapacity = curRowInds.size () + numToInsert;
3607 curRowInds.reserve (newCapacity);
3608 curRowVals.reserve (newCapacity);
3609 for (size_type k = 0; k < numToInsert; ++k) {
3610 curRowInds.push_back (indices[k]);
3611 curRowVals.push_back (values[k]);
3615 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3620 using Teuchos::Comm;
3621 using Teuchos::outArg;
3624 using Teuchos::REDUCE_MAX;
3625 using Teuchos::REDUCE_MIN;
3626 using Teuchos::reduceAll;
3629 typedef GlobalOrdinal GO;
3630 typedef typename Teuchos::Array<GO>::size_type size_type;
3631 const char tfecfFuncName[] =
"globalAssemble: ";
3633 RCP<const Comm<int> > comm = getComm ();
3635 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3636 (! isFillActive (), std::runtime_error,
"Fill must be active before " 3637 "you may call this method.");
3639 const size_t myNumNonlocalRows = nonlocals_.size ();
3646 const int iHaveNonlocalRows = (myNumNonlocalRows == 0) ? 0 : 1;
3647 int someoneHasNonlocalRows = 0;
3648 reduceAll<int, int> (*comm, REDUCE_MAX, iHaveNonlocalRows,
3649 outArg (someoneHasNonlocalRows));
3650 if (someoneHasNonlocalRows == 0) {
3664 RCP<const map_type> nonlocalRowMap;
3666 Teuchos::ArrayRCP<size_t> numEntPerNonlocalRow (myNumNonlocalRows);
3668 Teuchos::Array<GO> myNonlocalGblRows (myNumNonlocalRows);
3669 size_type curPos = 0;
3670 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
3671 ++mapIter, ++curPos) {
3672 myNonlocalGblRows[curPos] = mapIter->first;
3675 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
3676 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
3683 sort2 (gblCols.begin (), gblCols.end (), vals.begin ());
3684 typename Teuchos::Array<GO>::iterator gblCols_newEnd;
3685 typename Teuchos::Array<Scalar>::iterator vals_newEnd;
3686 merge2 (gblCols_newEnd, vals_newEnd,
3687 gblCols.begin (), gblCols.end (),
3688 vals.begin (), vals.end ());
3689 gblCols.erase (gblCols_newEnd, gblCols.end ());
3690 vals.erase (vals_newEnd, vals.end ());
3691 numEntPerNonlocalRow[curPos] = gblCols.size ();
3702 GO myMinNonlocalGblRow = std::numeric_limits<GO>::max ();
3704 auto iter = std::min_element (myNonlocalGblRows.begin (),
3705 myNonlocalGblRows.end ());
3706 if (iter != myNonlocalGblRows.end ()) {
3707 myMinNonlocalGblRow = *iter;
3710 GO gblMinNonlocalGblRow = 0;
3711 reduceAll<int, GO> (*comm, REDUCE_MIN, myMinNonlocalGblRow,
3712 outArg (gblMinNonlocalGblRow));
3713 const GO indexBase = gblMinNonlocalGblRow;
3714 const global_size_t INV = Teuchos::OrdinalTraits<global_size_t>::invalid ();
3715 nonlocalRowMap = rcp (
new map_type (INV, myNonlocalGblRows (), indexBase, comm));
3723 RCP<crs_matrix_type> nonlocalMatrix =
3724 rcp (
new crs_matrix_type (nonlocalRowMap, numEntPerNonlocalRow,
3727 size_type curPos = 0;
3728 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
3729 ++mapIter, ++curPos) {
3730 const GO gblRow = mapIter->first;
3732 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
3733 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
3735 nonlocalMatrix->insertGlobalValues (gblRow, gblCols (), vals ());
3747 auto origRowMap = this->getRowMap ();
3748 const bool origRowMapIsOneToOne = origRowMap->isOneToOne ();
3750 int isLocallyComplete = 1;
3752 if (origRowMapIsOneToOne) {
3753 export_type exportToOrig (nonlocalRowMap, origRowMap);
3755 isLocallyComplete = 0;
3757 this->doExport (*nonlocalMatrix, exportToOrig,
Tpetra::ADD);
3766 export_type exportToOneToOne (nonlocalRowMap, oneToOneRowMap);
3768 isLocallyComplete = 0;
3775 crs_matrix_type oneToOneMatrix (oneToOneRowMap, 0);
3777 oneToOneMatrix.doExport (*nonlocalMatrix, exportToOneToOne,
Tpetra::ADD);
3781 nonlocalMatrix = Teuchos::null;
3784 import_type importToOrig (oneToOneRowMap, origRowMap);
3785 this->doImport (oneToOneMatrix, importToOrig,
Tpetra::ADD);
3792 decltype (nonlocals_) newNonlocals;
3793 std::swap (nonlocals_, newNonlocals);
3802 int isGloballyComplete = 0;
3803 reduceAll<int, int> (*comm, REDUCE_MIN, isLocallyComplete,
3804 outArg (isGloballyComplete));
3805 TEUCHOS_TEST_FOR_EXCEPTION
3806 (isGloballyComplete != 1, std::runtime_error,
"On at least one process, " 3807 "you called insertGlobalValues with a global row index which is not in " 3808 "the matrix's row Map on any process in its communicator.");
3811 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3816 if (! isStaticGraph ()) {
3817 myGraph_->resumeFill (params);
3819 clearGlobalConstants ();
3820 fillComplete_ =
false;
3823 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3837 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3849 frobNorm_ = -STM::one ();
3852 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3857 TEUCHOS_TEST_FOR_EXCEPTION(
3858 getCrsGraph ().is_null (), std::logic_error,
"Tpetra::CrsMatrix::" 3859 "fillComplete(params): getCrsGraph() returns null. " 3860 "This should not happen at this point. " 3861 "Please report this bug to the Tpetra developers.");
3863 if (isStaticGraph () && getCrsGraph ()->isFillComplete ()) {
3864 fillComplete (getCrsGraph ()->getDomainMap (),
3865 getCrsGraph ()->getRangeMap (), params);
3867 fillComplete (getRowMap (), getRowMap (), params);
3871 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3875 const Teuchos::RCP<const map_type>& rangeMap,
3876 const Teuchos::RCP<Teuchos::ParameterList>& params)
3878 using Teuchos::ArrayRCP;
3881 const char tfecfFuncName[] =
"fillComplete";
3883 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3884 ! isFillActive () || isFillComplete (),
3885 std::runtime_error,
": Matrix fill state must be active (isFillActive() " 3886 "must be true) before you may call fillComplete().");
3887 const int numProcs = getComm ()->getSize ();
3895 bool assertNoNonlocalInserts =
false;
3898 bool sortGhosts =
true;
3900 if (! params.is_null ()) {
3901 assertNoNonlocalInserts = params->get (
"No Nonlocal Changes",
3902 assertNoNonlocalInserts);
3903 if (params->isParameter (
"sort column map ghost gids")) {
3904 sortGhosts = params->get (
"sort column map ghost gids", sortGhosts);
3906 else if (params->isParameter (
"Sort column Map ghost GIDs")) {
3907 sortGhosts = params->get (
"Sort column Map ghost GIDs", sortGhosts);
3912 const bool needGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
3914 if (! myGraph_.is_null ()) {
3915 myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
3918 if (! getCrsGraph()->indicesAreAllocated()) {
3921 allocateValues (LocalIndices, GraphNotYetAllocated);
3924 allocateValues (GlobalIndices, GraphNotYetAllocated);
3929 if (needGlobalAssemble) {
3933 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3934 numProcs == 1 && nonlocals_.size() > 0,
3935 std::runtime_error,
": cannot have nonlocal entries on a serial run. " 3936 "An invalid entry (i.e., with row index not in the row Map) must have " 3937 "been submitted to the CrsMatrix.");
3940 if (isStaticGraph ()) {
3958 const bool domainMapsMatch = staticGraph_->getDomainMap ()->isSameAs (*domainMap);
3959 const bool rangeMapsMatch = staticGraph_->getRangeMap ()->isSameAs (*rangeMap);
3961 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3962 ! domainMapsMatch, std::runtime_error,
3963 ": The CrsMatrix's domain Map does not match the graph's domain Map. " 3964 "The graph cannot be changed because it was given to the CrsMatrix " 3965 "constructor as const. You can fix this by passing in the graph's " 3966 "domain Map and range Map to the matrix's fillComplete call.");
3968 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3969 ! rangeMapsMatch, std::runtime_error,
3970 ": The CrsMatrix's range Map does not match the graph's range Map. " 3971 "The graph cannot be changed because it was given to the CrsMatrix " 3972 "constructor as const. You can fix this by passing in the graph's " 3973 "domain Map and range Map to the matrix's fillComplete call.");
3980 myGraph_->setDomainRangeMaps (domainMap, rangeMap);
3983 if (! myGraph_->hasColMap ()) {
3984 myGraph_->makeColMap ();
3989 myGraph_->makeIndicesLocal ();
3991 if (! myGraph_->isSorted ()) {
3994 if (! myGraph_->isMerged ()) {
3995 mergeRedundantEntries ();
3998 myGraph_->makeImportExport ();
3999 myGraph_->computeGlobalConstants ();
4000 myGraph_->fillComplete_ =
true;
4001 myGraph_->checkInternalState ();
4003 computeGlobalConstants ();
4005 if (myGraph_.is_null ()) {
4008 fillLocalMatrix (params);
4012 fillLocalGraphAndMatrix (params);
4022 fillComplete_ =
true;
4023 checkInternalState ();
4026 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4030 const Teuchos::RCP<const map_type> & rangeMap,
4031 const Teuchos::RCP<const import_type>& importer,
4032 const Teuchos::RCP<const export_type>& exporter,
4033 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
4035 #ifdef HAVE_TPETRA_MMM_TIMINGS 4037 if(!params.is_null())
4038 label = params->get(
"Timer Label",label);
4039 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
4040 using Teuchos::TimeMonitor;
4041 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-Graph"))));
4044 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
4045 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
4046 std::runtime_error,
"Matrix fill state must be active (isFillActive() " 4047 "must be true) before calling fillComplete().");
4048 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4049 myGraph_.is_null (), std::logic_error,
"myGraph_ is null. This is not allowed.");
4053 myGraph_->expertStaticFillComplete (domainMap, rangeMap, importer, exporter,params);
4055 #ifdef HAVE_TPETRA_MMM_TIMINGS 4056 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cGC"))));
4059 computeGlobalConstants ();
4061 #ifdef HAVE_TPETRA_MMM_TIMINGS 4062 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-fLGAM"))));
4066 fillLocalGraphAndMatrix (params);
4071 fillComplete_ =
true;
4074 #ifdef HAVE_TPETRA_DEBUG 4075 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isFillActive(), std::logic_error,
4076 ": We're at the end of fillComplete(), but isFillActive() is true. " 4077 "Please report this bug to the Tpetra developers.");
4078 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillComplete(), std::logic_error,
4079 ": We're at the end of fillComplete(), but isFillActive() is true. " 4080 "Please report this bug to the Tpetra developers.");
4081 #endif // HAVE_TPETRA_DEBUG 4083 #ifdef HAVE_TPETRA_MMM_TIMINGS 4084 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cIS"))));
4087 checkInternalState();
4090 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4095 TEUCHOS_TEST_FOR_EXCEPTION(
4096 isStaticGraph (), std::runtime_error,
"Tpetra::CrsMatrix::sortEntries: " 4097 "Cannot sort with static graph.");
4098 if (! myGraph_->isSorted ()) {
4099 const LocalOrdinal lclNumRows =
4100 static_cast<LocalOrdinal
> (this->getNodeNumRows ());
4101 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
4102 const RowInfo rowInfo = myGraph_->getRowInfo (row);
4103 Teuchos::ArrayView<impl_scalar_type> rv = this->getViewNonConst (rowInfo);
4104 myGraph_->template sortRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
4107 myGraph_->indicesAreSorted_ =
true;
4111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4116 TEUCHOS_TEST_FOR_EXCEPTION(
4117 isStaticGraph (), std::runtime_error,
"Tpetra::CrsMatrix::" 4118 "mergeRedundantEntries: Cannot merge with static graph.");
4119 if (! myGraph_->isMerged ()) {
4120 const LocalOrdinal lclNumRows =
4121 static_cast<LocalOrdinal
> (this->getNodeNumRows ());
4122 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
4123 const RowInfo rowInfo = myGraph_->getRowInfo (row);
4124 Teuchos::ArrayView<impl_scalar_type> rv = this->getViewNonConst (rowInfo);
4125 myGraph_->template mergeRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
4127 myGraph_->noRedundancies_ =
true;
4131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4139 using Teuchos::null;
4142 using Teuchos::rcp_const_cast;
4143 using Teuchos::rcpFromRef;
4144 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4145 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
4151 if (alpha ==
ZERO) {
4154 }
else if (beta != ONE) {
4168 RCP<const import_type> importer = this->getGraph ()->getImporter ();
4169 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
4175 const bool Y_is_overwritten = (beta ==
ZERO);
4186 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
4193 RCP<const MV> X_colMap;
4194 if (importer.is_null ()) {
4202 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in,
true);
4204 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
4209 X_colMap = rcpFromRef (X_in);
4217 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in);
4220 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
4221 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
4226 RCP<MV> Y_rowMap = getRowMapMultiVector (Y_in);
4233 if (! exporter.is_null ()) {
4234 this->
template localMultiply<Scalar, Scalar> (*X_colMap, *Y_rowMap,
4241 if (Y_is_overwritten) {
4266 Y_rowMap = getRowMapMultiVector (Y_in,
true);
4273 this->
template localMultiply<Scalar, Scalar> (*X_colMap,
4280 this->
template localMultiply<Scalar, Scalar> (*X_colMap, Y_in,
4290 if (Y_is_replicated) {
4295 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4300 const Teuchos::ETransp mode,
4304 using Teuchos::null;
4307 using Teuchos::rcp_const_cast;
4308 using Teuchos::rcpFromRef;
4309 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4312 if (alpha ==
ZERO) {
4335 RCP<const import_type> importer = this->getGraph ()->getImporter ();
4336 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
4342 const bool Y_is_overwritten = (beta ==
ZERO);
4343 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
4349 X = rcp (
new MV (X_in, Teuchos::Copy));
4351 X = rcpFromRef (X_in);
4355 if (importer != Teuchos::null) {
4356 if (importMV_ != Teuchos::null && importMV_->getNumVectors() != numVectors) {
4359 if (importMV_ == null) {
4360 importMV_ = rcp (
new MV (this->getColMap (), numVectors));
4363 if (exporter != Teuchos::null) {
4364 if (exportMV_ != Teuchos::null && exportMV_->getNumVectors() != numVectors) {
4367 if (exportMV_ == null) {
4368 exportMV_ = rcp (
new MV (this->getRowMap (), numVectors));
4374 if (! exporter.is_null ()) {
4375 exportMV_->doImport (X_in, *exporter,
INSERT);
4382 if (importer != Teuchos::null) {
4388 importMV_->putScalar (
ZERO);
4390 this->
template localMultiply<Scalar, Scalar> (*X, *importMV_, mode,
4392 if (Y_is_overwritten) {
4409 MV Y (Y_in, Teuchos::Copy);
4410 this->
template localMultiply<Scalar, Scalar> (*X, Y, mode, alpha, beta);
4413 this->
template localMultiply<Scalar, Scalar> (*X, Y_in, mode, alpha, beta);
4420 if (Y_is_replicated) {
4425 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4430 Teuchos::ETransp mode,
4434 TEUCHOS_TEST_FOR_EXCEPTION(
4435 ! isFillComplete (), std::runtime_error,
4436 "Tpetra::CrsMatrix::apply(): Cannot call apply() until fillComplete() " 4437 "has been called.");
4439 if (mode == Teuchos::NO_TRANS) {
4440 applyNonTranspose (X, Y, alpha, beta);
4447 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4450 applyTranspose (X, Y, mode, alpha, beta);
4454 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4460 const Scalar& dampingFactor,
4462 const int numSweeps)
const 4464 reorderedGaussSeidel (B, X, D, Teuchos::null, dampingFactor, direction, numSweeps);
4467 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4473 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
4474 const Scalar& dampingFactor,
4476 const int numSweeps)
const 4478 using Teuchos::null;
4481 using Teuchos::rcp_const_cast;
4482 using Teuchos::rcpFromRef;
4485 TEUCHOS_TEST_FOR_EXCEPTION(
4486 isFillComplete() ==
false, std::runtime_error,
4487 "Tpetra::CrsMatrix::gaussSeidel: cannot call this method until " 4488 "fillComplete() has been called.");
4489 TEUCHOS_TEST_FOR_EXCEPTION(
4491 std::invalid_argument,
4492 "Tpetra::CrsMatrix::gaussSeidel: The number of sweeps must be , " 4493 "nonnegative but you provided numSweeps = " << numSweeps <<
" < 0.");
4497 KokkosClassic::ESweepDirection localDirection;
4498 if (direction == Forward) {
4499 localDirection = KokkosClassic::Forward;
4501 else if (direction == Backward) {
4502 localDirection = KokkosClassic::Backward;
4504 else if (direction == Symmetric) {
4506 localDirection = KokkosClassic::Forward;
4509 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
4510 "Tpetra::CrsMatrix::gaussSeidel: The 'direction' enum does not have " 4511 "any of its valid values: Forward, Backward, or Symmetric.");
4514 if (numSweeps == 0) {
4521 RCP<const import_type> importer = this->getGraph()->getImporter();
4522 RCP<const export_type> exporter = this->getGraph()->getExporter();
4523 TEUCHOS_TEST_FOR_EXCEPTION(
4524 ! exporter.is_null (), std::runtime_error,
4525 "Tpetra's gaussSeidel implementation requires that the row, domain, " 4526 "and range Maps be the same. This cannot be the case, because the " 4527 "matrix has a nontrivial Export object.");
4529 RCP<const map_type> domainMap = this->getDomainMap ();
4530 RCP<const map_type> rangeMap = this->getRangeMap ();
4531 RCP<const map_type> rowMap = this->getGraph ()->getRowMap ();
4532 RCP<const map_type> colMap = this->getGraph ()->getColMap ();
4534 #ifdef HAVE_TEUCHOS_DEBUG 4539 TEUCHOS_TEST_FOR_EXCEPTION(
4540 ! X.
getMap ()->isSameAs (*domainMap),
4542 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4543 "multivector X be in the domain Map of the matrix.");
4544 TEUCHOS_TEST_FOR_EXCEPTION(
4545 ! B.
getMap ()->isSameAs (*rangeMap),
4547 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4548 "B be in the range Map of the matrix.");
4549 TEUCHOS_TEST_FOR_EXCEPTION(
4550 ! D.
getMap ()->isSameAs (*rowMap),
4552 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4553 "D be in the row Map of the matrix.");
4554 TEUCHOS_TEST_FOR_EXCEPTION(
4555 ! rowMap->isSameAs (*rangeMap),
4557 "Tpetra::CrsMatrix::gaussSeidel requires that the row Map and the " 4558 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
4559 TEUCHOS_TEST_FOR_EXCEPTION(
4560 ! domainMap->isSameAs (*rangeMap),
4562 "Tpetra::CrsMatrix::gaussSeidel requires that the domain Map and " 4563 "the range Map of the matrix be the same.");
4569 #endif // HAVE_TEUCHOS_DEBUG 4579 B_in = rcpFromRef (B);
4586 RCP<MV> B_in_nonconst = getRowMapMultiVector (B,
true);
4588 B_in = rcp_const_cast<
const MV> (B_in_nonconst);
4593 "gaussSeidel: The current implementation of the Gauss-Seidel kernel " 4594 "requires that X and B both have constant stride. Since B does not " 4595 "have constant stride, we had to make a copy. This is a limitation of " 4596 "the current implementation and not your fault, but we still report it " 4597 "as an efficiency warning for your information.");
4606 RCP<MV> X_domainMap;
4608 bool copiedInput =
false;
4610 if (importer.is_null ()) {
4612 X_domainMap = rcpFromRef (X);
4613 X_colMap = X_domainMap;
4614 copiedInput =
false;
4620 X_colMap = getColumnMapMultiVector (X,
true);
4621 X_domainMap = X_colMap;
4626 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the " 4627 "Gauss-Seidel kernel requires that X and B both have constant " 4628 "stride. Since X does not have constant stride, we had to make a " 4629 "copy. This is a limitation of the current implementation and not " 4630 "your fault, but we still report it as an efficiency warning for " 4631 "your information.");
4636 X_domainMap = rcpFromRef (X);
4640 X_colMap = X_domainMap->offsetViewNonConst (colMap, 0);
4650 X_colMap->doImport (X, *importer,
INSERT);
4651 copiedInput =
false;
4658 X_colMap = getColumnMapMultiVector (X,
true);
4659 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
4660 X_colMap->doImport (X, *importer,
INSERT);
4664 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the " 4665 "Gauss-Seidel kernel requires that X and B both have constant stride. " 4666 "Since X does not have constant stride, we had to make a copy. " 4667 "This is a limitation of the current implementation and not your fault, " 4668 "but we still report it as an efficiency warning for your information.");
4672 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
4673 if (! importer.is_null () && sweep > 0) {
4675 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4679 if (direction != Symmetric) {
4680 if (rowIndices.is_null ()) {
4681 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4686 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4693 const bool doImportBetweenDirections =
false;
4694 if (rowIndices.is_null ()) {
4695 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4697 KokkosClassic::Forward);
4702 if (doImportBetweenDirections) {
4704 if (! importer.is_null ()) {
4705 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4708 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4710 KokkosClassic::Backward);
4713 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4716 KokkosClassic::Forward);
4717 if (doImportBetweenDirections) {
4719 if (! importer.is_null ()) {
4720 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4723 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4726 KokkosClassic::Backward);
4736 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4742 const Scalar& dampingFactor,
4744 const int numSweeps,
4745 const bool zeroInitialGuess)
const 4747 reorderedGaussSeidelCopy (X, B, D, Teuchos::null, dampingFactor, direction,
4748 numSweeps, zeroInitialGuess);
4751 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4757 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
4758 const Scalar& dampingFactor,
4760 const int numSweeps,
4761 const bool zeroInitialGuess)
const 4763 using Teuchos::null;
4766 using Teuchos::rcpFromRef;
4767 using Teuchos::rcp_const_cast;
4769 const char prefix[] =
"Tpetra::CrsMatrix::(reordered)gaussSeidelCopy: ";
4770 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4772 TEUCHOS_TEST_FOR_EXCEPTION(
4773 ! isFillComplete (), std::runtime_error,
4774 prefix <<
"The matrix is not fill complete.");
4775 TEUCHOS_TEST_FOR_EXCEPTION(
4776 numSweeps < 0, std::invalid_argument,
4777 prefix <<
"The number of sweeps must be nonnegative, " 4778 "but you provided numSweeps = " << numSweeps <<
" < 0.");
4782 KokkosClassic::ESweepDirection localDirection;
4783 if (direction == Forward) {
4784 localDirection = KokkosClassic::Forward;
4786 else if (direction == Backward) {
4787 localDirection = KokkosClassic::Backward;
4789 else if (direction == Symmetric) {
4791 localDirection = KokkosClassic::Forward;
4794 TEUCHOS_TEST_FOR_EXCEPTION(
4795 true, std::invalid_argument,
4796 prefix <<
"The 'direction' enum does not have any of its valid " 4797 "values: Forward, Backward, or Symmetric.");
4800 if (numSweeps == 0) {
4804 RCP<const import_type> importer = this->getGraph ()->getImporter ();
4805 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
4806 TEUCHOS_TEST_FOR_EXCEPTION(
4807 ! exporter.is_null (), std::runtime_error,
4808 "This method's implementation currently requires that the matrix's row, " 4809 "domain, and range Maps be the same. This cannot be the case, because " 4810 "the matrix has a nontrivial Export object.");
4812 RCP<const map_type> domainMap = this->getDomainMap ();
4813 RCP<const map_type> rangeMap = this->getRangeMap ();
4814 RCP<const map_type> rowMap = this->getGraph ()->getRowMap ();
4815 RCP<const map_type> colMap = this->getGraph ()->getColMap ();
4817 #ifdef HAVE_TEUCHOS_DEBUG 4822 TEUCHOS_TEST_FOR_EXCEPTION(
4823 ! X.
getMap ()->isSameAs (*domainMap), std::runtime_error,
4824 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4825 "multivector X be in the domain Map of the matrix.");
4826 TEUCHOS_TEST_FOR_EXCEPTION(
4827 ! B.
getMap ()->isSameAs (*rangeMap), std::runtime_error,
4828 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4829 "B be in the range Map of the matrix.");
4830 TEUCHOS_TEST_FOR_EXCEPTION(
4831 ! D.
getMap ()->isSameAs (*rowMap), std::runtime_error,
4832 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4833 "D be in the row Map of the matrix.");
4834 TEUCHOS_TEST_FOR_EXCEPTION(
4835 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
4836 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the " 4837 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
4838 TEUCHOS_TEST_FOR_EXCEPTION(
4839 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
4840 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and " 4841 "the range Map of the matrix be the same.");
4847 #endif // HAVE_TEUCHOS_DEBUG 4858 RCP<MV> X_domainMap;
4859 bool copyBackOutput =
false;
4860 if (importer.is_null ()) {
4862 X_colMap = rcpFromRef (X);
4863 X_domainMap = rcpFromRef (X);
4869 if (zeroInitialGuess) {
4870 X_colMap->putScalar (
ZERO);
4878 X_colMap = getColumnMapMultiVector (X,
true);
4882 X_domainMap = X_colMap;
4883 if (! zeroInitialGuess) {
4886 }
catch (std::exception& e) {
4887 std::ostringstream os;
4888 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: " 4889 "deep_copy(*X_domainMap, X) threw an exception: " 4890 << e.what () <<
".";
4891 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
4894 copyBackOutput =
true;
4898 "gaussSeidelCopy: The current implementation of the Gauss-Seidel " 4899 "kernel requires that X and B both have constant stride. Since X " 4900 "does not have constant stride, we had to make a copy. This is a " 4901 "limitation of the current implementation and not your fault, but we " 4902 "still report it as an efficiency warning for your information.");
4906 X_colMap = getColumnMapMultiVector (X);
4907 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
4909 #ifdef HAVE_TPETRA_DEBUG 4910 auto X_colMap_host_view =
4911 X_colMap->template getLocalView<Kokkos::HostSpace> ();
4912 auto X_domainMap_host_view =
4913 X_domainMap->template getLocalView<Kokkos::HostSpace> ();
4915 if (X_colMap->getLocalLength () != 0 && X_domainMap->getLocalLength ()) {
4916 TEUCHOS_TEST_FOR_EXCEPTION
4917 (X_colMap_host_view.ptr_on_device () != X_domainMap_host_view.ptr_on_device (),
4918 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: Pointer to " 4919 "start of column Map view of X is not equal to pointer to start of " 4920 "(domain Map view of) X. This may mean that Tpetra::MultiVector::" 4921 "offsetViewNonConst is broken. " 4922 "Please report this bug to the Tpetra developers.");
4925 TEUCHOS_TEST_FOR_EXCEPTION(
4926 X_colMap_host_view.dimension_0 () < X_domainMap_host_view.dimension_0 () ||
4927 X_colMap->getLocalLength () < X_domainMap->getLocalLength (),
4928 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: " 4929 "X_colMap has fewer local rows than X_domainMap. " 4930 "X_colMap_host_view.dimension_0() = " << X_colMap_host_view.dimension_0 ()
4931 <<
", X_domainMap_host_view.dimension_0() = " 4932 << X_domainMap_host_view.dimension_0 ()
4933 <<
", X_colMap->getLocalLength() = " << X_colMap->getLocalLength ()
4934 <<
", and X_domainMap->getLocalLength() = " 4935 << X_domainMap->getLocalLength ()
4936 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 4937 "is broken. Please report this bug to the Tpetra developers.");
4939 TEUCHOS_TEST_FOR_EXCEPTION(
4940 X_colMap->getNumVectors () != X_domainMap->getNumVectors (),
4941 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: " 4942 "X_colMap has a different number of columns than X_domainMap. " 4943 "X_colMap->getNumVectors() = " << X_colMap->getNumVectors ()
4944 <<
" != X_domainMap->getNumVectors() = " 4945 << X_domainMap->getNumVectors ()
4946 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 4947 "is broken. Please report this bug to the Tpetra developers.");
4948 #endif // HAVE_TPETRA_DEBUG 4950 if (zeroInitialGuess) {
4952 X_colMap->putScalar (
ZERO);
4962 X_colMap->doImport (X, *importer,
INSERT);
4964 copyBackOutput =
true;
4972 B_in = rcpFromRef (B);
4978 RCP<MV> B_in_nonconst = getRowMapMultiVector (B,
true);
4981 }
catch (std::exception& e) {
4982 std::ostringstream os;
4983 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: " 4984 "deep_copy(*B_in_nonconst, B) threw an exception: " 4985 << e.what () <<
".";
4986 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
4988 B_in = rcp_const_cast<
const MV> (B_in_nonconst);
4993 "gaussSeidelCopy: The current implementation requires that B have " 4994 "constant stride. Since B does not have constant stride, we had to " 4995 "copy it into a separate constant-stride multivector. This is a " 4996 "limitation of the current implementation and not your fault, but we " 4997 "still report it as an efficiency warning for your information.");
5000 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
5001 if (! importer.is_null () && sweep > 0) {
5004 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
5008 if (direction != Symmetric) {
5009 if (rowIndices.is_null ()) {
5010 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5015 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5022 if (rowIndices.is_null ()) {
5023 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5025 KokkosClassic::Forward);
5031 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5033 KokkosClassic::Backward);
5036 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5039 KokkosClassic::Forward);
5040 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5043 KokkosClassic::Backward);
5049 if (copyBackOutput) {
5052 }
catch (std::exception& e) {
5053 TEUCHOS_TEST_FOR_EXCEPTION(
5054 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) " 5055 "threw an exception: " << e.what ());
5060 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5062 Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node, classic> >
5066 using Teuchos::ArrayRCP;
5070 typedef typename out_mat_type::local_matrix_type out_lcl_mat_type;
5071 typedef typename out_lcl_mat_type::values_type out_vals_type;
5072 typedef ArrayRCP<size_t>::size_type size_type;
5073 const char tfecfFuncName[] =
"convert";
5075 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5076 ! isFillComplete (), std::runtime_error,
"This matrix (the source of " 5077 "the conversion) is not fill complete. You must first call " 5078 "fillComplete() (possibly with the domain and range Map) without an " 5079 "intervening call to resumeFill(), before you may call this method.");
5086 RCP<out_mat_type> newmat;
5088 if (this->isStaticGraph ()) {
5090 newmat = rcp (
new out_mat_type (this->getCrsGraph ()));
5094 const size_type numVals =
5095 static_cast<size_type
> (this->lclMatrix_.values.dimension_0 ());
5101 out_vals_type newVals1D (
"Tpetra::CrsMatrix::val", numVals);
5102 for (size_type k = 0; k < numVals; ++k) {
5103 newVals1D(k) =
static_cast<T
> (this->k_values1D_(k));
5105 newmat->lclMatrix_ =
5106 out_lcl_mat_type (
"Tpetra::CrsMatrix::lclMatrix_",
5107 this->lclMatrix_.numCols (), newVals1D,
5108 this->lclMatrix_.graph);
5109 newmat->k_values1D_ = newVals1D;
5113 newmat->fillComplete (this->getDomainMap (), this->getRangeMap ());
5129 ArrayRCP<const size_t> ptr;
5130 ArrayRCP<const LocalOrdinal> ind;
5131 ArrayRCP<const Scalar> oldVal;
5132 this->getAllValues (ptr, ind, oldVal);
5134 RCP<const map_type> rowMap = this->getRowMap ();
5135 RCP<const map_type> colMap = this->getColMap ();
5139 const size_type numLocalRows =
5140 static_cast<size_type
> (rowMap->getNodeNumElements ());
5141 ArrayRCP<size_t> numEntriesPerRow (numLocalRows);
5142 for (size_type localRow = 0; localRow < numLocalRows; ++localRow) {
5143 numEntriesPerRow[localRow] =
5144 static_cast<size_type
> (getNumEntriesInLocalRow (localRow));
5147 newmat = rcp (
new out_mat_type (rowMap, colMap, numEntriesPerRow,
5151 const size_type numVals = this->lclMatrix_.values.dimension_0 ();
5152 ArrayRCP<T> newVals1D (numVals);
5154 for (size_type k = 0; k < numVals; ++k) {
5155 newVals1D[k] =
static_cast<T
> (this->k_values1D_(k));
5162 ArrayRCP<size_t> newPtr (ptr.size ());
5163 std::copy (ptr.begin (), ptr.end (), newPtr.begin ());
5164 ArrayRCP<LocalOrdinal> newInd (ind.size ());
5165 std::copy (ind.begin (), ind.end (), newInd.begin ());
5166 newmat->setAllValues (newPtr, newInd, newVals1D);
5171 RCP<const map_type> domainMap = this->getDomainMap ();
5172 RCP<const map_type> rangeMap = this->getRangeMap ();
5173 RCP<const import_type> importer = this->getCrsGraph ()->getImporter ();
5174 RCP<const export_type> exporter = this->getCrsGraph ()->getExporter ();
5175 newmat->expertStaticFillComplete (domainMap, rangeMap, importer, exporter);
5182 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5187 #ifdef HAVE_TPETRA_DEBUG 5188 const char tfecfFuncName[] =
"checkInternalState: ";
5189 const char err[] =
"Internal state is not consistent. " 5190 "Please report this bug to the Tpetra developers.";
5194 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5195 staticGraph_.is_null (),
5196 std::logic_error, err);
5200 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5201 ! myGraph_.is_null () && myGraph_ != staticGraph_,
5202 std::logic_error, err);
5204 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5205 isFillComplete () && ! staticGraph_->isFillComplete (),
5206 std::logic_error, err <<
" Specifically, the matrix is fill complete, " 5207 "but its graph is NOT fill complete.");
5209 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5210 isStorageOptimized () && ! values2D_.is_null (),
5211 std::logic_error, err);
5213 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5214 getProfileType() ==
StaticProfile && values2D_ != Teuchos::null,
5215 std::logic_error, err);
5217 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5218 getProfileType() ==
DynamicProfile && k_values1D_.dimension_0 () > 0,
5219 std::logic_error, err);
5222 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5223 staticGraph_->indicesAreAllocated () &&
5224 staticGraph_->getNodeAllocationSize() > 0 &&
5225 staticGraph_->getNodeNumRows() > 0
5226 && values2D_.is_null () &&
5227 k_values1D_.dimension_0 () == 0,
5228 std::logic_error, err);
5230 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5231 k_values1D_.dimension_0 () > 0 && values2D_ != Teuchos::null,
5232 std::logic_error, err <<
" Specifically, k_values1D_ is allocated (has " 5233 "size " << k_values1D_.dimension_0 () <<
" > 0) and values2D_ is also " 5234 "allocated. CrsMatrix is not suppose to have both a 1-D and a 2-D " 5235 "allocation at the same time.");
5239 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5244 std::ostringstream os;
5246 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
5247 if (this->getObjectLabel () !=
"") {
5248 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
5250 if (isFillComplete ()) {
5251 os <<
"isFillComplete: true" 5252 <<
", global dimensions: [" << getGlobalNumRows () <<
", " 5253 << getGlobalNumCols () <<
"]" 5254 <<
", global number of entries: " << getGlobalNumEntries ()
5258 os <<
"isFillComplete: false" 5259 <<
", global dimensions: [" << getGlobalNumRows () <<
", " 5260 << getGlobalNumCols () <<
"]}";
5265 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5269 const Teuchos::EVerbosityLevel verbLevel)
const 5273 using Teuchos::ArrayView;
5274 using Teuchos::Comm;
5276 using Teuchos::TypeNameTraits;
5277 using Teuchos::VERB_DEFAULT;
5278 using Teuchos::VERB_NONE;
5279 using Teuchos::VERB_LOW;
5280 using Teuchos::VERB_MEDIUM;
5281 using Teuchos::VERB_HIGH;
5282 using Teuchos::VERB_EXTREME;
5284 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
5286 if (vl == VERB_NONE) {
5290 Teuchos::OSTab tab0 (out);
5292 RCP<const Comm<int> > comm = this->getComm();
5293 const int myRank = comm->getRank();
5294 const int numProcs = comm->getSize();
5296 for (
size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
5299 width = std::max<size_t> (width,
static_cast<size_t> (11)) + 2;
5309 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
5311 Teuchos::OSTab tab1 (out);
5314 if (this->getObjectLabel () !=
"") {
5315 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
5318 out <<
"Template parameters:" << endl;
5319 Teuchos::OSTab tab2 (out);
5320 out <<
"Scalar: " << TypeNameTraits<Scalar>::name () << endl
5321 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name () << endl
5322 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name () << endl
5323 <<
"Node: " << TypeNameTraits<Node>::name () << endl;
5325 if (isFillComplete()) {
5326 out <<
"isFillComplete: true" << endl
5327 <<
"Global dimensions: [" << getGlobalNumRows () <<
", " 5328 << getGlobalNumCols () <<
"]" << endl
5329 <<
"Global number of entries: " << getGlobalNumEntries () << endl
5330 <<
"Global number of diagonal entries: " << getGlobalNumDiags ()
5331 << endl <<
"Global max number of entries in a row: " 5332 << getGlobalMaxNumRowEntries () << endl;
5335 out <<
"isFillComplete: false" << endl
5336 <<
"Global dimensions: [" << getGlobalNumRows () <<
", " 5337 << getGlobalNumCols () <<
"]" << endl;
5341 if (vl < VERB_MEDIUM) {
5347 out << endl <<
"Row Map:" << endl;
5349 if (getRowMap ().is_null ()) {
5351 out <<
"null" << endl;
5358 getRowMap ()->describe (out, vl);
5363 out <<
"Column Map: ";
5365 if (getColMap ().is_null ()) {
5367 out <<
"null" << endl;
5369 }
else if (getColMap () == getRowMap ()) {
5371 out <<
"same as row Map" << endl;
5377 getColMap ()->describe (out, vl);
5382 out <<
"Domain Map: ";
5384 if (getDomainMap ().is_null ()) {
5386 out <<
"null" << endl;
5388 }
else if (getDomainMap () == getRowMap ()) {
5390 out <<
"same as row Map" << endl;
5392 }
else if (getDomainMap () == getColMap ()) {
5394 out <<
"same as column Map" << endl;
5400 getDomainMap ()->describe (out, vl);
5405 out <<
"Range Map: ";
5407 if (getRangeMap ().is_null ()) {
5409 out <<
"null" << endl;
5411 }
else if (getRangeMap () == getDomainMap ()) {
5413 out <<
"same as domain Map" << endl;
5415 }
else if (getRangeMap () == getRowMap ()) {
5417 out <<
"same as row Map" << endl;
5423 getRangeMap ()->describe (out, vl);
5427 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5428 if (myRank == curRank) {
5429 out <<
"Process rank: " << curRank << endl;
5430 Teuchos::OSTab tab2 (out);
5431 if (! staticGraph_->indicesAreAllocated ()) {
5432 out <<
"Graph indices not allocated" << endl;
5435 out <<
"Number of allocated entries: " 5436 << staticGraph_->getNodeAllocationSize () << endl;
5438 out <<
"Number of entries: " << getNodeNumEntries () << endl;
5439 if (isFillComplete ()) {
5440 out <<
"Number of diagonal entries: " << getNodeNumDiags () << endl;
5442 out <<
"Max number of entries per row: " << getNodeMaxNumRowEntries ()
5451 if (vl < VERB_HIGH) {
5456 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5457 if (myRank == curRank) {
5458 out << std::setw(width) <<
"Proc Rank" 5459 << std::setw(width) <<
"Global Row" 5460 << std::setw(width) <<
"Num Entries";
5461 if (vl == VERB_EXTREME) {
5462 out << std::setw(width) <<
"(Index,Value)";
5465 for (
size_t r = 0; r < getNodeNumRows (); ++r) {
5466 const size_t nE = getNumEntriesInLocalRow(r);
5467 GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
5468 out << std::setw(width) << myRank
5469 << std::setw(width) << gid
5470 << std::setw(width) << nE;
5471 if (vl == VERB_EXTREME) {
5472 if (isGloballyIndexed()) {
5473 ArrayView<const GlobalOrdinal> rowinds;
5474 ArrayView<const Scalar> rowvals;
5475 getGlobalRowView (gid, rowinds, rowvals);
5476 for (
size_t j = 0; j < nE; ++j) {
5477 out <<
" (" << rowinds[j]
5478 <<
", " << rowvals[j]
5482 else if (isLocallyIndexed()) {
5483 ArrayView<const LocalOrdinal> rowinds;
5484 ArrayView<const Scalar> rowvals;
5485 getLocalRowView (r, rowinds, rowvals);
5486 for (
size_t j=0; j < nE; ++j) {
5487 out <<
" (" << getColMap()->getGlobalElement(rowinds[j])
5488 <<
", " << rowvals[j]
5504 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5517 const row_matrix_type* srcRowMat =
5518 dynamic_cast<const row_matrix_type*
> (&source);
5519 return (srcRowMat != NULL);
5522 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5527 const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs,
5528 const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs)
5530 using Teuchos::Array;
5531 using Teuchos::ArrayView;
5532 typedef LocalOrdinal LO;
5533 typedef GlobalOrdinal GO;
5536 const char tfecfFuncName[] =
"copyAndPermute: ";
5538 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5539 permuteToLIDs.size() != permuteFromLIDs.size(),
5540 std::invalid_argument,
"permuteToLIDs.size() = " << permuteToLIDs.size()
5541 <<
"!= permuteFromLIDs.size() = " << permuteFromLIDs.size() <<
".");
5546 const row_matrix_type& srcMat =
dynamic_cast<const row_matrix_type&
> (source);
5553 const map_type& srcRowMap = * (srcMat.getRowMap ());
5555 Array<Scalar> rowVals;
5556 const LO numSameIDs_as_LID =
static_cast<LO
> (numSameIDs);
5557 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5562 const GO targetGID = sourceGID;
5565 ArrayView<const GO> rowIndsConstView;
5566 ArrayView<const Scalar> rowValsConstView;
5568 if (sourceIsLocallyIndexed) {
5569 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5570 if (rowLength > static_cast<size_t> (rowInds.size())) {
5571 rowInds.resize (rowLength);
5572 rowVals.resize (rowLength);
5576 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
5577 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
5582 size_t checkRowLength = 0;
5583 srcMat.getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
5585 #ifdef HAVE_TPETRA_DEBUG 5586 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
5587 std::logic_error,
"For global row index " << sourceGID <<
", the source" 5588 " matrix's getNumEntriesInGlobalRow() method returns a row length of " 5589 << rowLength <<
", but the getGlobalRowCopy() method reports that " 5590 "the row length is " << checkRowLength <<
". Please report this bug " 5591 "to the Tpetra developers.");
5592 #endif // HAVE_TPETRA_DEBUG 5594 rowIndsConstView = rowIndsView.view (0, rowLength);
5595 rowValsConstView = rowValsView.view (0, rowLength);
5598 srcMat.getGlobalRowView (sourceGID, rowIndsConstView, rowValsConstView);
5602 if (isStaticGraph()) {
5605 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
REPLACE);
5611 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
INSERT);
5618 const map_type& tgtRowMap = * (this->getRowMap ());
5619 const size_t numPermuteToLIDs =
static_cast<size_t> (permuteToLIDs.size ());
5620 for (
size_t p = 0; p < numPermuteToLIDs; ++p) {
5625 ArrayView<const GO> rowIndsConstView;
5626 ArrayView<const Scalar> rowValsConstView;
5628 if (sourceIsLocallyIndexed) {
5629 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5630 if (rowLength > static_cast<size_t> (rowInds.size ())) {
5631 rowInds.resize (rowLength);
5632 rowVals.resize (rowLength);
5636 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
5637 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
5642 size_t checkRowLength = 0;
5643 srcMat.getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
5645 #ifdef HAVE_TPETRA_DEBUG 5646 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
5647 std::logic_error,
"For the source matrix's global row index " 5648 << sourceGID <<
", the source matrix's getNumEntriesInGlobalRow() " 5649 "method returns a row length of " << rowLength <<
", but the " 5650 "getGlobalRowCopy() method reports that the row length is " 5651 << checkRowLength <<
". Please report this bug to the Tpetra " 5653 #endif // HAVE_TPETRA_DEBUG 5655 rowIndsConstView = rowIndsView.view (0, rowLength);
5656 rowValsConstView = rowValsView.view (0, rowLength);
5659 srcMat.getGlobalRowView (sourceGID, rowIndsConstView, rowValsConstView);
5663 if (isStaticGraph()) {
5664 this->combineGlobalValues (targetGID, rowIndsConstView,
5668 this->combineGlobalValues (targetGID, rowIndsConstView,
5669 rowValsConstView,
INSERT);
5674 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5678 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
5679 Teuchos::Array<char>& exports,
5680 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5681 size_t& constantNumPackets,
5684 using Teuchos::Array;
5685 using Teuchos::ArrayView;
5686 using Teuchos::av_reinterpret_cast;
5687 typedef LocalOrdinal LO;
5688 typedef GlobalOrdinal GO;
5689 const char tfecfFuncName[] =
"packAndPrepare: ";
5712 const row_matrix_type* srcRowMat =
5713 dynamic_cast<const row_matrix_type*
> (&source);
5714 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5715 srcRowMat == NULL, std::invalid_argument,
5716 "The source object of the Import or Export operation is neither a " 5717 "CrsMatrix (with the same template parameters as the target object), " 5718 "nor a RowMatrix (with the same first four template parameters as the " 5720 #ifdef HAVE_TPETRA_DEBUG 5722 using Teuchos::reduceAll;
5723 std::ostringstream msg;
5726 srcRowMat->pack (exportLIDs, exports, numPacketsPerLID,
5727 constantNumPackets, distor);
5728 }
catch (std::exception& e) {
5733 const Teuchos::Comm<int>& comm = * (this->getComm ());
5734 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
5735 lclBad, Teuchos::outArg (gblBad));
5737 const int myRank = comm.getRank ();
5738 const int numProcs = comm.getSize ();
5739 for (
int r = 0; r < numProcs; ++r) {
5740 if (r == myRank && lclBad != 0) {
5741 std::ostringstream os;
5742 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
5743 std::cerr << os.str ();
5749 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5750 true, std::logic_error,
"pack() threw an exception on one or " 5751 "more participating processes.");
5755 srcRowMat->pack (exportLIDs, exports, numPacketsPerLID,
5756 constantNumPackets, distor);
5757 #endif // HAVE_TPETRA_DEBUG 5760 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5762 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
5763 packRow (
char*
const numEntOut,
5766 const size_t numEnt,
5767 const LocalOrdinal lclRow)
const 5769 using Teuchos::ArrayView;
5770 typedef LocalOrdinal LO;
5771 typedef GlobalOrdinal GO;
5773 const LO numEntLO =
static_cast<LO
> (numEnt);
5774 memcpy (numEntOut, &numEntLO,
sizeof (LO));
5775 if (this->isLocallyIndexed ()) {
5779 ArrayView<const LO> indIn;
5780 ArrayView<const Scalar> valIn;
5781 this->getLocalRowView (lclRow, indIn, valIn);
5782 const map_type& colMap = * (this->getColMap ());
5785 for (
size_t k = 0; k < numEnt; ++k) {
5786 const GO gblIndIn = colMap.getGlobalElement (indIn[k]);
5787 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
5789 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
5791 else if (this->isGloballyIndexed ()) {
5797 ArrayView<const GO> indIn;
5798 ArrayView<const Scalar> valIn;
5799 const map_type& rowMap = * (this->getRowMap ());
5800 const GO gblRow = rowMap.getGlobalElement (lclRow);
5801 this->getGlobalRowView (gblRow, indIn, valIn);
5802 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
5803 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
5814 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5816 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
5817 unpackRow (Scalar*
const valInTmp,
5818 GlobalOrdinal*
const indInTmp,
5819 const size_t tmpSize,
5820 const char*
const valIn,
5821 const char*
const indIn,
5822 const size_t numEnt,
5823 const LocalOrdinal lclRow,
5826 if (tmpSize < numEnt || (numEnt != 0 && (valInTmp == NULL || indInTmp == NULL))) {
5829 memcpy (valInTmp, valIn, numEnt *
sizeof (Scalar));
5830 memcpy (indInTmp, indIn, numEnt *
sizeof (GlobalOrdinal));
5831 const GlobalOrdinal gblRow = this->getRowMap ()->getGlobalElement (lclRow);
5832 Teuchos::ArrayView<Scalar> val ((numEnt == 0) ? NULL : valInTmp, numEnt);
5833 Teuchos::ArrayView<GlobalOrdinal> ind ((numEnt == 0) ? NULL : indInTmp, numEnt);
5834 this->combineGlobalValues (gblRow, ind, val, combineMode);
5839 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5841 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
5842 allocatePackSpace (Teuchos::Array<char>& exports,
5843 size_t& totalNumEntries,
5844 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs)
const 5846 typedef LocalOrdinal LO;
5847 typedef GlobalOrdinal GO;
5848 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
5850 const size_type numExportLIDs = exportLIDs.size ();
5853 totalNumEntries = 0;
5854 for (size_type i = 0; i < numExportLIDs; ++i) {
5855 const LO lclRow = exportLIDs[i];
5856 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
5859 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
5862 totalNumEntries += curNumEntries;
5873 const size_t allocSize =
5874 static_cast<size_t> (numExportLIDs) *
sizeof (LO) +
5875 totalNumEntries * (
sizeof (Scalar) +
sizeof (GO));
5876 if (static_cast<size_t> (exports.size ()) < allocSize) {
5877 exports.resize (allocSize);
5881 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5884 pack (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
5885 Teuchos::Array<char>& exports,
5886 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5887 size_t& constantNumPackets,
5890 using Teuchos::Array;
5891 using Teuchos::ArrayView;
5892 using Teuchos::av_reinterpret_cast;
5894 typedef LocalOrdinal LO;
5895 typedef GlobalOrdinal GO;
5896 typedef typename ArrayView<const LO>::size_type size_type;
5897 const char tfecfFuncName[] =
"pack: ";
5899 const size_type numExportLIDs = exportLIDs.size ();
5900 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5901 numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
5902 "exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size()" 5903 " = " << numPacketsPerLID.size () <<
".");
5908 constantNumPackets = 0;
5913 size_t totalNumEntries = 0;
5914 allocatePackSpace (exports, totalNumEntries, exportLIDs);
5915 const size_t bufSize =
static_cast<size_t> (exports.size ());
5927 size_type firstBadIndex = 0;
5928 size_t firstBadOffset = 0;
5929 size_t firstBadNumBytes = 0;
5930 bool outOfBounds =
false;
5931 bool packErr =
false;
5933 char*
const exportsRawPtr = exports.getRawPtr ();
5935 for (size_type i = 0; i < numExportLIDs; ++i) {
5936 const LO lclRow = exportLIDs[i];
5937 const size_t numEnt = this->getNumEntriesInLocalRow (lclRow);
5941 numPacketsPerLID[i] = 0;
5944 char*
const numEntBeg = exportsRawPtr + offset;
5945 char*
const numEntEnd = numEntBeg +
sizeof (LO);
5946 char*
const valBeg = numEntEnd;
5947 char*
const valEnd = valBeg + numEnt *
sizeof (Scalar);
5948 char*
const indBeg = valEnd;
5949 const size_t numBytes =
sizeof (LO) +
5950 numEnt * (
sizeof (Scalar) +
sizeof (GO));
5951 if (offset > bufSize || offset + numBytes > bufSize) {
5953 firstBadOffset = offset;
5954 firstBadNumBytes = numBytes;
5958 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
5961 firstBadOffset = offset;
5962 firstBadNumBytes = numBytes;
5968 numPacketsPerLID[i] = numBytes;
5973 TEUCHOS_TEST_FOR_EXCEPTION(
5974 outOfBounds, std::logic_error,
"First invalid offset into 'exports' " 5975 "pack buffer at index i = " << firstBadIndex <<
". exportLIDs[i]: " 5976 << exportLIDs[firstBadIndex] <<
", bufSize: " << bufSize <<
", offset: " 5977 << firstBadOffset <<
", numBytes: " << firstBadNumBytes <<
".");
5978 TEUCHOS_TEST_FOR_EXCEPTION(
5979 packErr, std::logic_error,
"First error in packRow() at index i = " 5980 << firstBadIndex <<
". exportLIDs[i]: " << exportLIDs[firstBadIndex]
5981 <<
", bufSize: " << bufSize <<
", offset: " << firstBadOffset
5982 <<
", numBytes: " << firstBadNumBytes <<
".");
5985 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5989 const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
5990 const Teuchos::ArrayView<const Scalar>& values,
5993 const char tfecfFuncName[] =
"combineGlobalValues: ";
5995 if (isStaticGraph ()) {
5999 if (combineMode ==
ADD) {
6000 sumIntoGlobalValues (globalRowIndex, columnIndices, values);
6002 else if (combineMode ==
REPLACE) {
6003 replaceGlobalValues (globalRowIndex, columnIndices, values);
6005 else if (combineMode ==
ABSMAX) {
6006 using Details::AbsMax;
6008 this->
template transformGlobalValues<AbsMax<Scalar> > (globalRowIndex,
6012 else if (combineMode ==
INSERT) {
6013 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6014 isStaticGraph () && combineMode ==
INSERT, std::invalid_argument,
6015 "INSERT combine mode is not allowed if the matrix has a static graph " 6016 "(i.e., was constructed with the CrsMatrix constructor that takes a " 6017 "const CrsGraph pointer).");
6020 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6021 true, std::logic_error,
"Invalid combine mode; should never get " 6022 "here! Please report this bug to the Tpetra developers.");
6026 if (combineMode ==
ADD || combineMode ==
INSERT) {
6033 insertGlobalValuesFiltered (globalRowIndex, columnIndices, values);
6044 else if (combineMode ==
ABSMAX) {
6045 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6046 ! isStaticGraph () && combineMode ==
ABSMAX, std::logic_error,
6047 "ABSMAX combine mode when the matrix has a dynamic graph is not yet " 6050 else if (combineMode ==
REPLACE) {
6051 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6052 ! isStaticGraph () && combineMode ==
REPLACE, std::logic_error,
6053 "REPLACE combine mode when the matrix has a dynamic graph is not yet " 6057 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6058 true, std::logic_error,
"Should never get here! Please report this " 6059 "bug to the Tpetra developers.");
6065 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6069 const Teuchos::ArrayView<const char>& imports,
6070 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
6071 size_t constantNumPackets,
6075 #ifdef HAVE_TPETRA_DEBUG 6076 const char tfecfFuncName[] =
"unpackAndCombine: ";
6078 const char* validModeNames[4] = {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT"};
6079 const int numValidModes = 4;
6081 if (std::find (validModes, validModes+numValidModes, combineMode) ==
6082 validModes+numValidModes) {
6083 std::ostringstream os;
6084 os <<
"Invalid combine mode. Valid modes are {";
6085 for (
int k = 0; k < numValidModes; ++k) {
6086 os << validModeNames[k];
6087 if (k < numValidModes - 1) {
6092 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6093 true, std::invalid_argument, os.str ());
6097 using Teuchos::reduceAll;
6098 std::ostringstream msg;
6101 this->unpackAndCombineImpl (importLIDs, imports, numPacketsPerLID,
6102 constantNumPackets, distor, combineMode);
6103 }
catch (std::exception& e) {
6108 const Teuchos::Comm<int>& comm = * (this->getComm ());
6109 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
6110 lclBad, Teuchos::outArg (gblBad));
6112 const int myRank = comm.getRank ();
6113 const int numProcs = comm.getSize ();
6114 for (
int r = 0; r < numProcs; ++r) {
6115 if (r == myRank && lclBad != 0) {
6116 std::ostringstream os;
6117 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
6118 std::cerr << os.str ();
6124 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6125 true, std::logic_error,
"unpackAndCombineImpl() threw an " 6126 "exception on one or more participating processes.");
6130 this->unpackAndCombineImpl (importLIDs, imports, numPacketsPerLID,
6131 constantNumPackets, distor, combineMode);
6132 #endif // HAVE_TPETRA_DEBUG 6135 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6139 const Teuchos::ArrayView<const char>& imports,
6140 const Teuchos::ArrayView<const size_t>& numPacketsPerLID,
6141 size_t constantNumPackets,
6145 using Teuchos::Array;
6146 typedef LocalOrdinal LO;
6147 typedef GlobalOrdinal GO;
6148 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
6149 const char tfecfFuncName[] =
"unpackAndCombine: ";
6151 #ifdef HAVE_TPETRA_DEBUG 6153 const char* validModeNames[4] = {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT"};
6154 const int numValidModes = 4;
6156 if (std::find (validModes, validModes+numValidModes, combineMode) ==
6157 validModes+numValidModes) {
6158 std::ostringstream os;
6159 os <<
"Invalid combine mode. Valid modes are {";
6160 for (
int k = 0; k < numValidModes; ++k) {
6161 os << validModeNames[k];
6162 if (k < numValidModes - 1) {
6167 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6168 true, std::invalid_argument, os.str ());
6170 #endif // HAVE_TPETRA_DEBUG 6172 const size_type numImportLIDs = importLIDs.size ();
6173 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6174 numImportLIDs != numPacketsPerLID.size (), std::invalid_argument,
6175 "importLIDs.size() = " << numImportLIDs <<
" != numPacketsPerLID.size()" 6176 <<
" = " << numPacketsPerLID.size () <<
".");
6182 size_type firstBadIndex = 0;
6183 size_t firstBadOffset = 0;
6184 size_t firstBadExpectedNumBytes = 0;
6185 size_t firstBadNumBytes = 0;
6186 LO firstBadNumEnt = 0;
6196 bool wrongNumBytes =
false;
6197 bool unpackErr =
false;
6199 const size_t bufSize =
static_cast<size_t> (imports.size ());
6200 const char*
const importsRawPtr = imports.getRawPtr ();
6209 Array<Scalar> valInTmp;
6211 for (size_type i = 0; i < numImportLIDs; ++i) {
6212 const LO lclRow = importLIDs[i];
6213 const size_t numBytes = numPacketsPerLID[i];
6216 const char*
const numEntBeg = importsRawPtr + offset;
6217 const char*
const numEntEnd = numEntBeg +
sizeof (LO);
6222 memcpy (&numEnt, numEntBeg,
sizeof (LO));
6224 const char*
const valBeg = numEntEnd;
6225 const char*
const valEnd =
6226 valBeg +
static_cast<size_t> (numEnt) *
sizeof (Scalar);
6227 const char*
const indBeg = valEnd;
6228 const size_t expectedNumBytes =
sizeof (LO) +
6229 static_cast<size_t> (numEnt) * (
sizeof (Scalar) +
sizeof (GO));
6231 if (expectedNumBytes > numBytes) {
6233 firstBadOffset = offset;
6234 firstBadExpectedNumBytes = expectedNumBytes;
6235 firstBadNumBytes = numBytes;
6236 firstBadNumEnt = numEnt;
6237 wrongNumBytes =
true;
6240 if (offset > bufSize || offset + numBytes > bufSize) {
6242 firstBadOffset = offset;
6243 firstBadExpectedNumBytes = expectedNumBytes;
6244 firstBadNumBytes = numBytes;
6245 firstBadNumEnt = numEnt;
6249 size_t tmpNumEnt =
static_cast<size_t> (valInTmp.size ());
6250 if (tmpNumEnt < static_cast<size_t> (numEnt) ||
6251 static_cast<size_t> (indInTmp.size ()) < static_cast<size_t> (numEnt)) {
6253 tmpNumEnt = std::max (static_cast<size_t> (numEnt), tmpNumEnt * 2);
6254 valInTmp.resize (tmpNumEnt);
6255 indInTmp.resize (tmpNumEnt);
6258 ! unpackRow (valInTmp.getRawPtr (), indInTmp.getRawPtr (), tmpNumEnt,
6259 valBeg, indBeg, numEnt, lclRow, combineMode);
6262 firstBadOffset = offset;
6263 firstBadExpectedNumBytes = expectedNumBytes;
6264 firstBadNumBytes = numBytes;
6265 firstBadNumEnt = numEnt;
6272 if (wrongNumBytes || outOfBounds || unpackErr) {
6273 std::ostringstream os;
6274 os <<
" importLIDs[i]: " << importLIDs[firstBadIndex]
6275 <<
", bufSize: " << bufSize
6276 <<
", offset: " << firstBadOffset
6277 <<
", numBytes: " << firstBadNumBytes
6278 <<
", expectedNumBytes: " << firstBadExpectedNumBytes
6279 <<
", numEnt: " << firstBadNumEnt;
6280 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6281 wrongNumBytes, std::logic_error,
"At index i = " << firstBadIndex
6282 <<
", expectedNumBytes > numBytes." << os.str ());
6283 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6284 outOfBounds, std::logic_error,
"First invalid offset into 'imports' " 6285 "unpack buffer at index i = " << firstBadIndex <<
"." << os.str ());
6286 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6287 unpackErr, std::logic_error,
"First error in unpackRow() at index i = " 6288 << firstBadIndex <<
"." << os.str ());
6292 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6293 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
6296 const bool force)
const 6298 using Teuchos::null;
6302 TEUCHOS_TEST_FOR_EXCEPTION(
6303 ! this->hasColMap (), std::runtime_error,
"Tpetra::CrsMatrix::getColumn" 6304 "MapMultiVector: You may only call this method if the matrix has a " 6305 "column Map. If the matrix does not yet have a column Map, you should " 6306 "first call fillComplete (with domain and range Map if necessary).");
6310 TEUCHOS_TEST_FOR_EXCEPTION(
6311 ! this->getGraph ()->isFillComplete (), std::runtime_error,
"Tpetra::" 6312 "CrsMatrix::getColumnMapMultiVector: You may only call this method if " 6313 "this matrix's graph is fill complete.");
6316 RCP<const import_type> importer = this->getGraph ()->getImporter ();
6317 RCP<const map_type> colMap = this->getColMap ();
6330 if (! importer.is_null () || force) {
6331 if (importMV_.is_null () || importMV_->getNumVectors () != numVecs) {
6332 X_colMap = rcp (
new MV (colMap, numVecs));
6335 importMV_ = X_colMap;
6338 X_colMap = importMV_;
6349 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6350 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
6353 const bool force)
const 6355 using Teuchos::null;
6361 TEUCHOS_TEST_FOR_EXCEPTION(
6362 ! this->getGraph ()->isFillComplete (), std::runtime_error,
"Tpetra::" 6363 "CrsMatrix::getRowMapMultiVector: You may only call this method if this " 6364 "matrix's graph is fill complete.");
6367 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
6371 RCP<const map_type> rowMap = this->getRowMap ();
6383 if (! exporter.is_null () || force) {
6384 if (exportMV_.is_null () || exportMV_->getNumVectors () != numVecs) {
6385 Y_rowMap = rcp (
new MV (rowMap, numVecs));
6386 exportMV_ = Y_rowMap;
6389 Y_rowMap = exportMV_;
6395 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6400 TEUCHOS_TEST_FOR_EXCEPTION(
6401 myGraph_.is_null (), std::logic_error,
"Tpetra::CrsMatrix::" 6402 "removeEmptyProcessesInPlace: This method does not work when the matrix " 6403 "was created with a constant graph (that is, when it was created using " 6404 "the version of its constructor that takes an RCP<const CrsGraph>). " 6405 "This is because the matrix is not allowed to modify the graph in that " 6406 "case, but removing empty processes requires modifying the graph.");
6407 myGraph_->removeEmptyProcessesInPlace (newMap);
6411 this->map_ = this->getRowMap ();
6415 staticGraph_ = Teuchos::rcp_const_cast<
const Graph> (myGraph_);
6418 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6419 Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
6424 const Teuchos::RCP<const map_type>& domainMap,
6425 const Teuchos::RCP<const map_type>& rangeMap,
6426 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6428 using Teuchos::Array;
6429 using Teuchos::ArrayRCP;
6430 using Teuchos::ArrayView;
6431 using Teuchos::ParameterList;
6434 using Teuchos::rcp_implicit_cast;
6435 using Teuchos::sublist;
6436 typedef LocalOrdinal LO;
6437 typedef GlobalOrdinal GO;
6441 const crs_matrix_type& B = *
this;
6442 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
6443 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
6450 RCP<const map_type> A_rangeMap = A.
getRangeMap ();
6451 RCP<const map_type> B_domainMap = B.getDomainMap ();
6452 RCP<const map_type> B_rangeMap = B.getRangeMap ();
6454 RCP<const map_type> theDomainMap = domainMap;
6455 RCP<const map_type> theRangeMap = rangeMap;
6457 if (domainMap.is_null ()) {
6458 if (B_domainMap.is_null ()) {
6459 TEUCHOS_TEST_FOR_EXCEPTION(
6460 A_domainMap.is_null (), std::invalid_argument,
6461 "Tpetra::CrsMatrix::add: If neither A nor B have a domain Map, " 6462 "then you must supply a nonnull domain Map to this method.");
6463 theDomainMap = A_domainMap;
6465 theDomainMap = B_domainMap;
6468 if (rangeMap.is_null ()) {
6469 if (B_rangeMap.is_null ()) {
6470 TEUCHOS_TEST_FOR_EXCEPTION(
6471 A_rangeMap.is_null (), std::invalid_argument,
6472 "Tpetra::CrsMatrix::add: If neither A nor B have a range Map, " 6473 "then you must supply a nonnull range Map to this method.");
6474 theRangeMap = A_rangeMap;
6476 theRangeMap = B_rangeMap;
6480 #ifdef HAVE_TPETRA_DEBUG 6484 if (! A_domainMap.is_null () && ! A_rangeMap.is_null ()) {
6485 if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
6486 TEUCHOS_TEST_FOR_EXCEPTION(
6487 ! B_domainMap->isSameAs (*A_domainMap), std::invalid_argument,
6488 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a domain Map " 6489 "which is the same as (isSameAs) this RowMatrix's domain Map.");
6490 TEUCHOS_TEST_FOR_EXCEPTION(
6491 ! B_rangeMap->isSameAs (*A_rangeMap), std::invalid_argument,
6492 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a range Map " 6493 "which is the same as (isSameAs) this RowMatrix's range Map.");
6494 TEUCHOS_TEST_FOR_EXCEPTION(
6495 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
6496 std::invalid_argument,
6497 "Tpetra::CrsMatrix::add: The input domain Map must be the same as " 6498 "(isSameAs) this RowMatrix's domain Map.");
6499 TEUCHOS_TEST_FOR_EXCEPTION(
6500 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
6501 std::invalid_argument,
6502 "Tpetra::CrsMatrix::add: The input range Map must be the same as " 6503 "(isSameAs) this RowMatrix's range Map.");
6506 else if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
6507 TEUCHOS_TEST_FOR_EXCEPTION(
6508 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
6509 std::invalid_argument,
6510 "Tpetra::CrsMatrix::add: The input domain Map must be the same as " 6511 "(isSameAs) this RowMatrix's domain Map.");
6512 TEUCHOS_TEST_FOR_EXCEPTION(
6513 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
6514 std::invalid_argument,
6515 "Tpetra::CrsMatrix::add: The input range Map must be the same as " 6516 "(isSameAs) this RowMatrix's range Map.");
6519 TEUCHOS_TEST_FOR_EXCEPTION(
6520 domainMap.is_null () || rangeMap.is_null (), std::invalid_argument,
6521 "Tpetra::CrsMatrix::add: If neither A nor B have a domain and range " 6522 "Map, then you must supply a nonnull domain and range Map to this " 6525 #endif // HAVE_TPETRA_DEBUG 6530 bool callFillComplete =
true;
6531 RCP<ParameterList> constructorSublist;
6532 RCP<ParameterList> fillCompleteSublist;
6533 if (! params.is_null ()) {
6534 callFillComplete = params->get (
"Call fillComplete", callFillComplete);
6535 constructorSublist = sublist (params,
"Constructor parameters");
6536 fillCompleteSublist = sublist (params,
"fillComplete parameters");
6539 RCP<const map_type> A_rowMap = A.
getRowMap ();
6540 RCP<const map_type> B_rowMap = B.getRowMap ();
6541 RCP<const map_type> C_rowMap = B_rowMap;
6542 RCP<crs_matrix_type> C;
6549 if (A_rowMap->isSameAs (*B_rowMap)) {
6550 const LO localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
6551 ArrayRCP<size_t> C_maxNumEntriesPerRow (localNumRows, 0);
6554 if (alpha !=
ZERO) {
6555 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
6557 C_maxNumEntriesPerRow[localRow] += A_numEntries;
6562 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
6563 const size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
6564 C_maxNumEntriesPerRow[localRow] += B_numEntries;
6568 if (constructorSublist.is_null ()) {
6569 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
6572 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
6583 if (constructorSublist.is_null ()) {
6587 constructorSublist));
6591 #ifdef HAVE_TPETRA_DEBUG 6592 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
6593 "Tpetra::RowMatrix::add: C should not be null at this point. " 6594 "Please report this bug to the Tpetra developers.");
6595 #endif // HAVE_TPETRA_DEBUG 6602 if (alpha !=
ZERO) {
6603 const LO A_localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
6604 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
6606 const GO globalRow = A_rowMap->getGlobalElement (localRow);
6607 if (A_numEntries > static_cast<size_t> (ind.size ())) {
6608 ind.resize (A_numEntries);
6609 val.resize (A_numEntries);
6611 ArrayView<GO> indView = ind (0, A_numEntries);
6612 ArrayView<Scalar> valView = val (0, A_numEntries);
6616 for (
size_t k = 0; k < A_numEntries; ++k) {
6617 valView[k] *= alpha;
6620 C->insertGlobalValues (globalRow, indView, valView);
6625 const LO B_localNumRows =
static_cast<LO
> (B_rowMap->getNodeNumElements ());
6626 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
6627 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
6628 const GO globalRow = B_rowMap->getGlobalElement (localRow);
6629 if (B_numEntries > static_cast<size_t> (ind.size ())) {
6630 ind.resize (B_numEntries);
6631 val.resize (B_numEntries);
6633 ArrayView<GO> indView = ind (0, B_numEntries);
6634 ArrayView<Scalar> valView = val (0, B_numEntries);
6635 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
6638 for (
size_t k = 0; k < B_numEntries; ++k) {
6642 C->insertGlobalValues (globalRow, indView, valView);
6646 if (callFillComplete) {
6647 if (fillCompleteSublist.is_null ()) {
6648 C->fillComplete (theDomainMap, theRangeMap);
6650 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
6653 return rcp_implicit_cast<row_matrix_type> (C);
6656 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6660 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
6661 const Teuchos::RCP<const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node> > & domainTransfer,
6662 const Teuchos::RCP<const map_type>& domainMap,
6663 const Teuchos::RCP<const map_type>& rangeMap,
6664 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6667 using Teuchos::ArrayRCP;
6668 using Teuchos::ArrayView;
6669 using Teuchos::Comm;
6670 using Teuchos::ParameterList;
6672 typedef LocalOrdinal LO;
6673 typedef GlobalOrdinal GO;
6674 typedef node_type NT;
6678 #ifdef HAVE_TPETRA_MMM_TIMINGS 6680 if(!params.is_null())
6681 label = params->get(
"Timer Label",label);
6682 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
6683 using Teuchos::TimeMonitor;
6684 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-1"))));
6692 const import_type* xferAsImport =
dynamic_cast<const import_type*
> (&rowTransfer);
6693 const export_type* xferAsExport =
dynamic_cast<const export_type*
> (&rowTransfer);
6694 TEUCHOS_TEST_FOR_EXCEPTION(
6695 xferAsImport == NULL && xferAsExport == NULL, std::invalid_argument,
6696 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' input " 6697 "argument must be either an Import or an Export, and its template " 6698 "parameters must match the corresponding template parameters of the " 6706 Teuchos::RCP<const import_type> xferDomainAsImport = Teuchos::rcp_dynamic_cast<
const import_type> (domainTransfer);
6707 Teuchos::RCP<const export_type> xferDomainAsExport = Teuchos::rcp_dynamic_cast<
const export_type> (domainTransfer);
6709 if(! domainTransfer.is_null()) {
6710 TEUCHOS_TEST_FOR_EXCEPTION(
6711 (xferDomainAsImport.is_null() && xferDomainAsExport.is_null()), std::invalid_argument,
6712 "Tpetra::CrsMatrix::transferAndFillComplete: The 'domainTransfer' input " 6713 "argument must be either an Import or an Export, and its template " 6714 "parameters must match the corresponding template parameters of the " 6717 TEUCHOS_TEST_FOR_EXCEPTION(
6718 ( xferAsImport != NULL || ! xferDomainAsImport.is_null() ) &&
6719 (( xferAsImport != NULL && xferDomainAsImport.is_null() ) ||
6720 ( xferAsImport == NULL && ! xferDomainAsImport.is_null() )), std::invalid_argument,
6721 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input " 6722 "arguments must be of the same type (either Import or Export).");
6724 TEUCHOS_TEST_FOR_EXCEPTION(
6725 ( xferAsExport != NULL || ! xferDomainAsExport.is_null() ) &&
6726 (( xferAsExport != NULL && xferDomainAsExport.is_null() ) ||
6727 ( xferAsExport == NULL && ! xferDomainAsExport.is_null() )), std::invalid_argument,
6728 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input " 6729 "arguments must be of the same type (either Import or Export).");
6735 const bool communication_needed = rowTransfer.getSourceMap ()->isDistributed ();
6741 bool reverseMode =
false;
6742 bool restrictComm =
false;
6743 RCP<ParameterList> matrixparams;
6744 if (! params.is_null ()) {
6745 reverseMode = params->get (
"Reverse Mode", reverseMode);
6746 restrictComm = params->get (
"Restrict Communicator", restrictComm);
6747 matrixparams = sublist (params,
"CrsMatrix");
6752 RCP<const map_type> MyRowMap = reverseMode ?
6753 rowTransfer.getSourceMap () : rowTransfer.getTargetMap ();
6754 RCP<const map_type> MyColMap;
6755 RCP<const map_type> MyDomainMap = ! domainMap.is_null () ?
6756 domainMap : getDomainMap ();
6757 RCP<const map_type> MyRangeMap = ! rangeMap.is_null () ?
6758 rangeMap : getRangeMap ();
6759 RCP<const map_type> BaseRowMap = MyRowMap;
6760 RCP<const map_type> BaseDomainMap = MyDomainMap;
6768 if (! destMat.is_null ()) {
6779 const bool NewFlag = ! destMat->getGraph ()->isLocallyIndexed () &&
6780 ! destMat->getGraph ()->isGloballyIndexed ();
6781 TEUCHOS_TEST_FOR_EXCEPTION(
6782 ! NewFlag, std::invalid_argument,
"Tpetra::CrsMatrix::" 6783 "transferAndFillComplete: The input argument 'destMat' is only allowed " 6784 "to be nonnull, if its graph is empty (neither locally nor globally " 6793 TEUCHOS_TEST_FOR_EXCEPTION(
6794 ! destMat->getRowMap ()->isSameAs (*MyRowMap), std::invalid_argument,
6795 "Tpetra::CrsMatrix::transferAndFillComplete: The (row) Map of the " 6796 "input argument 'destMat' is not the same as the (row) Map specified " 6797 "by the input argument 'rowTransfer'.");
6798 TEUCHOS_TEST_FOR_EXCEPTION(
6799 ! destMat->checkSizes (*
this), std::invalid_argument,
6800 "Tpetra::CrsMatrix::transferAndFillComplete: You provided a nonnull " 6801 "destination matrix, but checkSizes() indicates that it is not a legal " 6802 "legal target for redistribution from the source matrix (*this). This " 6803 "may mean that they do not have the same dimensions.");
6817 TEUCHOS_TEST_FOR_EXCEPTION(
6818 ! (reverseMode || getRowMap ()->isSameAs (*rowTransfer.getSourceMap ())),
6819 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: " 6820 "rowTransfer->getSourceMap() must match this->getRowMap() in forward mode.");
6821 TEUCHOS_TEST_FOR_EXCEPTION(
6822 ! (! reverseMode || getRowMap ()->isSameAs (*rowTransfer.getTargetMap ())),
6823 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: " 6824 "rowTransfer->getTargetMap() must match this->getRowMap() in reverse mode.");
6827 TEUCHOS_TEST_FOR_EXCEPTION(
6828 ! xferDomainAsImport.is_null() && ! xferDomainAsImport->getTargetMap()->isSameAs(*domainMap),
6829 std::invalid_argument,
6830 "Tpetra::CrsMatrix::transferAndFillComplete: The target map of the 'domainTransfer' input " 6831 "argument must be the same as the rebalanced domain map 'domainMap'");
6833 TEUCHOS_TEST_FOR_EXCEPTION(
6834 ! xferDomainAsExport.is_null() && ! xferDomainAsExport->getSourceMap()->isSameAs(*domainMap),
6835 std::invalid_argument,
6836 "Tpetra::CrsMatrix::transferAndFillComplete: The source map of the 'domainTransfer' input " 6837 "argument must be the same as the rebalanced domain map 'domainMap'");
6850 const size_t NumSameIDs = rowTransfer.getNumSameIDs();
6851 ArrayView<const LO> ExportLIDs = reverseMode ?
6852 rowTransfer.getRemoteLIDs () : rowTransfer.getExportLIDs ();
6853 ArrayView<const LO> RemoteLIDs = reverseMode ?
6854 rowTransfer.getExportLIDs () : rowTransfer.getRemoteLIDs ();
6855 ArrayView<const LO> PermuteToLIDs = reverseMode ?
6856 rowTransfer.getPermuteFromLIDs () : rowTransfer.getPermuteToLIDs ();
6857 ArrayView<const LO> PermuteFromLIDs = reverseMode ?
6858 rowTransfer.getPermuteToLIDs () : rowTransfer.getPermuteFromLIDs ();
6859 Distributor& Distor = rowTransfer.getDistributor ();
6862 Teuchos::Array<int> SourcePids;
6863 Teuchos::Array<int> TargetPids;
6864 int MyPID = getComm ()->getRank ();
6867 RCP<const map_type> ReducedRowMap, ReducedColMap,
6868 ReducedDomainMap, ReducedRangeMap;
6869 RCP<const Comm<int> > ReducedComm;
6873 if (destMat.is_null ()) {
6874 destMat = rcp (
new this_type (MyRowMap, 0,
StaticProfile, matrixparams));
6881 ReducedRowMap = MyRowMap->removeEmptyProcesses ();
6882 ReducedComm = ReducedRowMap.is_null () ?
6884 ReducedRowMap->getComm ();
6885 destMat->removeEmptyProcessesInPlace (ReducedRowMap);
6887 ReducedDomainMap = MyRowMap.getRawPtr () == MyDomainMap.getRawPtr () ?
6889 MyDomainMap->replaceCommWithSubset (ReducedComm);
6890 ReducedRangeMap = MyRowMap.getRawPtr () == MyRangeMap.getRawPtr () ?
6892 MyRangeMap->replaceCommWithSubset (ReducedComm);
6895 MyRowMap = ReducedRowMap;
6896 MyDomainMap = ReducedDomainMap;
6897 MyRangeMap = ReducedRangeMap;
6900 if (! ReducedComm.is_null ()) {
6901 MyPID = ReducedComm->getRank ();
6908 ReducedComm = MyRowMap->getComm ();
6914 #ifdef HAVE_TPETRA_MMM_TIMINGS 6915 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ImportSetup"))));
6918 RCP<const import_type> MyImporter = getGraph ()->getImporter ();
6921 bool bSameDomainMap = BaseDomainMap->isSameAs (*getDomainMap ());
6923 if (! restrictComm && ! MyImporter.is_null () && bSameDomainMap ) {
6930 Import_Util::getPids (*MyImporter, SourcePids,
false);
6932 else if (restrictComm && ! MyImporter.is_null () && bSameDomainMap) {
6935 IntVectorType SourceDomain_pids(getDomainMap (),
true);
6936 IntVectorType SourceCol_pids(getColMap());
6938 SourceDomain_pids.putScalar(MyPID);
6940 SourceCol_pids.doImport (SourceDomain_pids, *MyImporter,
INSERT);
6941 SourcePids.resize (getColMap ()->getNodeNumElements ());
6942 SourceCol_pids.get1dCopy (SourcePids ());
6944 else if (MyImporter.is_null () && bSameDomainMap) {
6946 SourcePids.resize (getColMap ()->getNodeNumElements ());
6947 SourcePids.assign (getColMap ()->getNodeNumElements (), MyPID);
6949 else if ( ! MyImporter.is_null () &&
6950 ! domainTransfer.is_null () ) {
6957 IntVectorType TargetDomain_pids (domainMap);
6958 TargetDomain_pids.putScalar (MyPID);
6961 IntVectorType SourceDomain_pids (getDomainMap ());
6964 IntVectorType SourceCol_pids (getColMap ());
6966 if (! reverseMode && ! xferDomainAsImport.is_null() ) {
6967 SourceDomain_pids.doExport (TargetDomain_pids, *xferDomainAsImport,
INSERT);
6969 else if (reverseMode && ! xferDomainAsExport.is_null() ) {
6970 SourceDomain_pids.doExport (TargetDomain_pids, *xferDomainAsExport,
INSERT);
6972 else if (! reverseMode && ! xferDomainAsExport.is_null() ) {
6973 SourceDomain_pids.doImport (TargetDomain_pids, *xferDomainAsExport,
INSERT);
6975 else if (reverseMode && ! xferDomainAsImport.is_null() ) {
6976 SourceDomain_pids.doImport (TargetDomain_pids, *xferDomainAsImport,
INSERT);
6979 TEUCHOS_TEST_FOR_EXCEPTION(
6980 true, std::logic_error,
"Tpetra::CrsMatrix::" 6981 "transferAndFillComplete: Should never get here! " 6982 "Please report this bug to a Tpetra developer.");
6984 SourceCol_pids.doImport (SourceDomain_pids, *MyImporter,
INSERT);
6985 SourcePids.resize (getColMap ()->getNodeNumElements ());
6986 SourceCol_pids.get1dCopy (SourcePids ());
6988 else if (BaseDomainMap->isSameAs (*BaseRowMap) &&
6989 getDomainMap ()->isSameAs (*getRowMap ())) {
6991 IntVectorType TargetRow_pids (domainMap);
6992 IntVectorType SourceRow_pids (getRowMap ());
6993 IntVectorType SourceCol_pids (getColMap ());
6995 TargetRow_pids.putScalar (MyPID);
6996 if (! reverseMode && xferAsImport != NULL) {
6997 SourceRow_pids.doExport (TargetRow_pids, *xferAsImport,
INSERT);
6999 else if (reverseMode && xferAsExport != NULL) {
7000 SourceRow_pids.doExport (TargetRow_pids, *xferAsExport,
INSERT);
7002 else if (! reverseMode && xferAsExport != NULL) {
7003 SourceRow_pids.doImport (TargetRow_pids, *xferAsExport,
INSERT);
7005 else if (reverseMode && xferAsImport != NULL) {
7006 SourceRow_pids.doImport (TargetRow_pids, *xferAsImport,
INSERT);
7009 TEUCHOS_TEST_FOR_EXCEPTION(
7010 true, std::logic_error,
"Tpetra::CrsMatrix::" 7011 "transferAndFillComplete: Should never get here! " 7012 "Please report this bug to a Tpetra developer.");
7014 SourceCol_pids.doImport (SourceRow_pids, *MyImporter,
INSERT);
7015 SourcePids.resize (getColMap ()->getNodeNumElements ());
7016 SourceCol_pids.get1dCopy (SourcePids ());
7019 TEUCHOS_TEST_FOR_EXCEPTION(
7020 true, std::invalid_argument,
"Tpetra::CrsMatrix::" 7021 "transferAndFillComplete: This method only allows either domainMap == " 7022 "getDomainMap (), or (domainMap == rowTransfer.getTargetMap () and " 7023 "getDomainMap () == getRowMap ()).");
7025 #ifdef HAVE_TPETRA_MMM_TIMINGS 7026 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-2"))));
7030 size_t constantNumPackets = destMat->constantNumberOfPackets ();
7031 if (constantNumPackets == 0) {
7034 execution_space::fence ();
7035 destMat->numExportPacketsPerLID_ =
7036 decltype (destMat->numExportPacketsPerLID_) (
"numExportPacketsPerLID",
7037 ExportLIDs.size ());
7038 execution_space::fence ();
7039 destMat->numImportPacketsPerLID_ =
7040 decltype (destMat->numImportPacketsPerLID_) (
"numImportPacketsPerLID",
7041 RemoteLIDs.size ());
7042 execution_space::fence ();
7049 const size_t rbufLen = RemoteLIDs.size() * constantNumPackets;
7050 destMat->reallocImportsIfNeeded (rbufLen);
7054 #ifdef HAVE_TPETRA_DEBUG 7056 using Teuchos::outArg;
7057 using Teuchos::REDUCE_MAX;
7058 using Teuchos::reduceAll;
7061 RCP<const Teuchos::Comm<int> > comm = this->getComm ();
7062 const int myRank = comm->getRank ();
7063 const int numProcs = comm->getSize ();
7065 std::ostringstream os;
7069 destMat->numExportPacketsPerLID_.template modify<Kokkos::HostSpace> ();
7070 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
7072 Import_Util::packAndPrepareWithOwningPIDs (*
this, ExportLIDs,
7074 numExportPacketsPerLID,
7075 constantNumPackets, Distor,
7078 catch (std::exception& e) {
7079 os <<
"Proc " << myRank <<
": " << e.what ();
7083 if (! comm.is_null ()) {
7084 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
7088 cerr <<
"packAndPrepareWithOwningPIDs threw an exception: " << endl;
7090 std::ostringstream err;
7091 for (
int r = 0; r < numProcs; ++r) {
7092 if (r == myRank && lclErr != 0) {
7093 cerr << os.str () << endl;
7100 TEUCHOS_TEST_FOR_EXCEPTION(
7101 true, std::logic_error,
"packAndPrepareWithOwningPIDs threw an " 7109 destMat->numExportPacketsPerLID_.template modify<Kokkos::HostSpace> ();
7110 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
7112 Import_Util::packAndPrepareWithOwningPIDs (*
this, ExportLIDs,
7114 numExportPacketsPerLID,
7115 constantNumPackets, Distor,
7118 #endif // HAVE_TPETRA_DEBUG 7121 #ifdef HAVE_TPETRA_MMM_TIMINGS 7122 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Transfer"))));
7125 if (communication_needed) {
7127 if (constantNumPackets == 0) {
7131 destMat->numExportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
7132 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
7134 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
7135 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
7137 Distor.doReversePostsAndWaits (numExportPacketsPerLID, 1,
7138 numImportPacketsPerLID);
7139 size_t totalImportPackets = 0;
7140 for (
Array_size_type i = 0; i < numImportPacketsPerLID.size (); ++i) {
7141 totalImportPackets += numImportPacketsPerLID[i];
7146 destMat->reallocImportsIfNeeded (totalImportPackets);
7147 destMat->imports_.template modify<Kokkos::HostSpace> ();
7148 Teuchos::ArrayView<char> hostImports =
7152 destMat->exports_.template sync<Kokkos::HostSpace> ();
7153 Teuchos::ArrayView<const char> hostExports =
7155 Distor.doReversePostsAndWaits (hostExports,
7156 numExportPacketsPerLID,
7158 numImportPacketsPerLID);
7161 destMat->imports_.template modify<Kokkos::HostSpace> ();
7162 Teuchos::ArrayView<char> hostImports =
7166 destMat->exports_.template sync<Kokkos::HostSpace> ();
7167 Teuchos::ArrayView<const char> hostExports =
7169 Distor.doReversePostsAndWaits (hostExports,
7175 if (constantNumPackets == 0) {
7179 destMat->numExportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
7180 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
7182 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
7183 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
7185 Distor.doPostsAndWaits (numExportPacketsPerLID, 1,
7186 numImportPacketsPerLID);
7187 size_t totalImportPackets = 0;
7188 for (
Array_size_type i = 0; i < numImportPacketsPerLID.size (); ++i) {
7189 totalImportPackets += numImportPacketsPerLID[i];
7194 destMat->reallocImportsIfNeeded (totalImportPackets);
7195 destMat->imports_.template modify<Kokkos::HostSpace> ();
7196 Teuchos::ArrayView<char> hostImports =
7200 destMat->exports_.template sync<Kokkos::HostSpace> ();
7201 Teuchos::ArrayView<const char> hostExports =
7203 Distor.doPostsAndWaits (hostExports,
7204 numExportPacketsPerLID,
7206 numImportPacketsPerLID);
7209 destMat->imports_.template modify<Kokkos::HostSpace> ();
7210 Teuchos::ArrayView<char> hostImports =
7214 destMat->exports_.template sync<Kokkos::HostSpace> ();
7215 Teuchos::ArrayView<const char> hostExports =
7217 Distor.doPostsAndWaits (hostExports,
7228 #ifdef HAVE_TPETRA_MMM_TIMINGS 7229 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-1"))));
7233 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
7234 Teuchos::ArrayView<const size_t> numImportPacketsPerLID =
7236 destMat->imports_.template sync<Kokkos::HostSpace> ();
7237 Teuchos::ArrayView<const char> hostImports =
7240 Import_Util::unpackAndCombineWithOwningPIDsCount (*
this, RemoteLIDs,
7242 numImportPacketsPerLID,
7248 size_t N = BaseRowMap->getNodeNumElements ();
7251 ArrayRCP<size_t> CSR_rowptr(N+1);
7252 ArrayRCP<GO> CSR_colind_GID;
7253 ArrayRCP<LO> CSR_colind_LID;
7254 ArrayRCP<Scalar> CSR_vals;
7255 CSR_colind_GID.resize (mynnz);
7256 CSR_vals.resize (mynnz);
7260 if (
typeid (LO) ==
typeid (GO)) {
7261 CSR_colind_LID = Teuchos::arcp_reinterpret_cast<LO> (CSR_colind_GID);
7264 CSR_colind_LID.resize (mynnz);
7272 Import_Util::unpackAndCombineIntoCrsArrays (*
this, RemoteLIDs, hostImports,
7273 numImportPacketsPerLID,
7274 constantNumPackets, Distor,
INSERT, NumSameIDs,
7275 PermuteToLIDs, PermuteFromLIDs, N, mynnz, MyPID,
7276 CSR_rowptr (), CSR_colind_GID (), CSR_vals (),
7277 SourcePids (), TargetPids);
7282 #ifdef HAVE_TPETRA_MMM_TIMINGS 7283 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-2"))));
7288 Teuchos::Array<int> RemotePids;
7289 Import_Util::lowCommunicationMakeColMapAndReindex (CSR_rowptr (),
7293 TargetPids, RemotePids,
7300 ReducedColMap = (MyRowMap.getRawPtr () == MyColMap.getRawPtr ()) ?
7302 MyColMap->replaceCommWithSubset (ReducedComm);
7303 MyColMap = ReducedColMap;
7307 destMat->replaceColMap (MyColMap);
7314 if (ReducedComm.is_null ()) {
7321 #ifdef HAVE_TPETRA_MMM_TIMINGS 7322 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-3"))));
7324 Import_Util::sortCrsEntries (CSR_rowptr (),
7327 if ((! reverseMode && xferAsImport != NULL) ||
7328 (reverseMode && xferAsExport != NULL)) {
7329 Import_Util::sortCrsEntries (CSR_rowptr (),
7333 else if ((! reverseMode && xferAsExport != NULL) ||
7334 (reverseMode && xferAsImport != NULL)) {
7335 Import_Util::sortAndMergeCrsEntries (CSR_rowptr (),
7338 if (CSR_rowptr[N] != mynnz) {
7339 CSR_colind_LID.resize (CSR_rowptr[N]);
7340 CSR_vals.resize (CSR_rowptr[N]);
7344 TEUCHOS_TEST_FOR_EXCEPTION(
7345 true, std::logic_error,
"Tpetra::CrsMatrix::" 7346 "transferAndFillComplete: Should never get here! " 7347 "Please report this bug to a Tpetra developer.");
7358 destMat->setAllValues (CSR_rowptr, CSR_colind_LID, CSR_vals);
7364 Teuchos::ParameterList esfc_params;
7365 #ifdef HAVE_TPETRA_MMM_TIMINGS 7366 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC CreateImporter"))));
7368 RCP<import_type> MyImport = rcp (
new import_type (MyDomainMap, MyColMap, RemotePids));
7369 #ifdef HAVE_TPETRA_MMM_TIMINGS 7370 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ESFC"))));
7372 esfc_params.set(
"Timer Label",prefix + std::string(
"TAFC"));
7375 destMat->expertStaticFillComplete (MyDomainMap, MyRangeMap, MyImport,Teuchos::null,rcp(&esfc_params,
false));
7378 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
7383 const Teuchos::RCP<const map_type>& domainMap,
7384 const Teuchos::RCP<const map_type>& rangeMap,
7385 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 7387 transferAndFillComplete (destMatrix, importer, Teuchos::null, domainMap, rangeMap, params);
7390 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
7396 const Teuchos::RCP<const map_type>& domainMap,
7397 const Teuchos::RCP<const map_type>& rangeMap,
7398 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 7400 transferAndFillComplete (destMatrix, rowImporter, Teuchos::rcpFromRef(domainImporter), domainMap, rangeMap, params);
7403 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
7408 const Teuchos::RCP<const map_type>& domainMap,
7409 const Teuchos::RCP<const map_type>& rangeMap,
7410 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 7412 transferAndFillComplete (destMatrix, exporter, Teuchos::null, domainMap, rangeMap, params);
7415 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
7421 const Teuchos::RCP<const map_type>& domainMap,
7422 const Teuchos::RCP<const map_type>& rangeMap,
7423 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 7425 transferAndFillComplete (destMatrix, rowExporter, Teuchos::rcpFromRef(domainExporter), domainMap, rangeMap, params);
7436 #define TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR,LO,GO,NODE) \ 7438 template class CrsMatrix< SCALAR , LO , GO , NODE >; \ 7439 template Teuchos::RCP< CrsMatrix< SCALAR , LO , GO , NODE > > \ 7440 CrsMatrix< SCALAR , LO , GO , NODE >::convert< SCALAR > () const; 7442 #define TPETRA_CRSMATRIX_CONVERT_INSTANT(SO,SI,LO,GO,NODE) \ 7444 template Teuchos::RCP< CrsMatrix< SO , LO , GO , NODE > > \ 7445 CrsMatrix< SI , LO , GO , NODE >::convert< SO > () const; 7447 #define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 7449 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 7450 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 7451 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7452 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7453 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& importer, \ 7454 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7455 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7456 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 7457 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7458 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7459 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 7460 const Teuchos::RCP<Teuchos::ParameterList>& params); 7462 #define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \ 7464 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 7465 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 7466 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7467 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7468 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowImporter, \ 7469 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7470 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7471 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainImporter, \ 7472 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7473 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7474 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 7475 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7476 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7477 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 7478 const Teuchos::RCP<Teuchos::ParameterList>& params); 7481 #define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 7483 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 7484 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 7485 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7486 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7487 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& exporter, \ 7488 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7489 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7490 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 7491 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7492 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7493 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 7494 const Teuchos::RCP<Teuchos::ParameterList>& params); 7496 #define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \ 7498 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 7499 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 7500 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7501 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7502 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowExporter, \ 7503 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7504 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7505 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainExporter, \ 7506 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7507 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7508 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 7509 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 7510 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 7511 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 7512 const Teuchos::RCP<Teuchos::ParameterList>& params); 7515 #define TPETRA_CRSMATRIX_INSTANT(SCALAR, LO, GO ,NODE) \ 7516 TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \ 7517 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 7518 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 7519 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \ 7520 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) 7522 #endif // TPETRA_CRSMATRIX_DEF_HPP void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
ProfileType getProfileType() const
Returns true if the matrix was allocated with static data structures.
LocalOrdinal replaceGlobalValues(const GlobalOrdinal globalRow, const typename UnmanagedView< GlobalIndicesViewType >::type &inputInds, const typename UnmanagedView< ImplScalarViewType >::type &inputVals) const
Replace one or more entries' values, using global indices.
void getGlobalRowView(GlobalOrdinal GlobalRow, Teuchos::ArrayView< const GlobalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const
Get a constant, nonpersisting view of a row of this matrix, using global row and column indices...
Kokkos::CrsMatrix< impl_scalar_type, LocalOrdinal, execution_space, void, typename local_graph_type::size_type > local_matrix_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
void reindexColumns(crs_graph_type *const graph, const Teuchos::RCP< const map_type > &newColMap, const Teuchos::RCP< const import_type > &newImport=Teuchos::null, const bool sortEachRow=true)
Reindex the column indices in place, and replace the column Map. Optionally, replace the Import objec...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
KOKKOS_INLINE_FUNCTION bool outOfBounds(const IntegerType x, const IntegerType exclusiveUpperBound)
Is x out of bounds? That is, is x less than zero, or greater than or equal to the given exclusive upp...
Functor for the the ABSMAX CombineMode of Import and Export operations.
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
void copyOffsets(const OutputViewType &dst, const InputViewType &src)
Copy row offsets (in a sparse graph or matrix) from src to dst. The offsets may have different types...
virtual bool isLocallyIndexed() const =0
Whether matrix indices are locally indexed.
std::string description() const
A one-line description of this object.
mag_type getFrobeniusNorm() const
Compute and return the Frobenius norm of the matrix.
LocalOrdinal local_ordinal_type
This class' second template parameter; the type of local indices.
size_t getNodeNumEntries() const
The local number of entries in this matrix.
void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to alpha.
LocalOrdinal sumIntoLocalValues(const LocalOrdinal localRow, const typename UnmanagedView< LocalIndicesViewType >::type &inputInds, const typename UnmanagedView< ImplScalarViewType >::type &inputVals, const bool atomic=useAtomicUpdatesByDefault) const
Sum into one or more sparse matrix entries, using local row and column indices.
LocalOrdinal replaceLocalValues(const LocalOrdinal localRow, const typename UnmanagedView< LocalIndicesViewType >::type &inputInds, const typename UnmanagedView< ImplScalarViewType >::type &inputVals) const
Replace one or more entries' values, using local row and column indices.
Teuchos::RCP< const RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const
This matrix's graph, as a RowGraph.
global_size_t getGlobalNumDiags() const
Returns the number of global diagonal entries, based on global row/column index comparisons.
bool isFillActive() const
Whether the matrix is not fill complete.
void sortEntries()
Sort the entries of each row by their column indices.
bool isNodeGlobalElement(GlobalOrdinal globalIndex) const
Whether the given global index is owned by this Map on the calling process.
Teuchos::RCP< CrsMatrix< T, LocalOrdinal, GlobalOrdinal, Node, classic > > convert() const
Return another CrsMatrix with the same entries, but converted to a different Scalar type T...
virtual void copyAndPermute(const SrcDistObject &source, size_t numSameIDs, const Teuchos::ArrayView< const LocalOrdinal > &permuteToLIDs, const Teuchos::ArrayView< const LocalOrdinal > &permuteFromLIDs)
Perform copies and permutations that are local to this process.
size_t getNodeNumDiags() const
Returns the number of local diagonal entries, based on global row/column index comparisons.
void resumeFill(const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Resume operations that may change the values or structure of the matrix.
Teuchos::ArrayView< const impl_scalar_type > getView(RowInfo rowinfo) const
Constant view of all entries (including extra space) in the given row.
global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
void gaussSeidelCopy(MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps, const bool zeroInitialGuess) const
Version of gaussSeidel(), with fewer requirements on X.
void rightScale(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &x)
bool isLocallyComplete() const
Do all source Map indices on the calling process exist on at least one process (not necessarily this ...
Teuchos::ArrayView< impl_scalar_type > getViewNonConst(const RowInfo &rowinfo) const
Nonconst view of all entries (including extra space) in the given row.
size_t getNodeNumCols() const
The number of columns connected to the locally owned rows of this matrix.
Teuchos::RCP< const map_type > getRowMap() const
The Map that describes the row distribution in this matrix.
bool hasColMap() const
Indicates whether the matrix has a well-defined column map.
One or more distributed dense vectors.
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
The current number of entries on the calling process in the specified local row.
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries on this node in the specified global row.
Teuchos::RCP< node_type > getNode() const
The Kokkos Node instance.
GlobalOrdinal global_ordinal_type
This class' third template parameter; the type of global indices.
void mergeRedundantEntries()
Merge entries in each row with the same column indices.
Declare and define Tpetra::Details::copyOffsets, an implementation detail of Tpetra (in particular...
void deep_copy(MultiVector< DS, DL, DG, DN, dstClassic > &dst, const MultiVector< SS, SL, SG, SN, srcClassic > &src)
Copy the contents of the MultiVector src into dst.
Teuchos::RCP< MV > getColumnMapMultiVector(const MV &X_domainMap, const bool force=false) const
Create a (or fetch a cached) column Map MultiVector.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Compute a sparse matrix-MultiVector multiply.
void replaceColMap(const Teuchos::RCP< const map_type > &newColMap)
Replace the matrix's column Map with the given Map.
virtual bool checkSizes(const SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
local_map_type getLocalMap() const
Get the local Map for Kokkos kernels.
Node node_type
This class' fourth template parameter; the Kokkos device type.
bool isLowerTriangular() const
Indicates whether the matrix is lower triangular.
Node::device_type device_type
The Kokkos device type.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
Allocation information for a locally owned row in a CrsGraph or CrsMatrix.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
virtual bool supportsRowViews() const
Return true if getLocalRowView() and getGlobalRowView() are valid for this object.
void getLocalDiagCopy(Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &diag) const
Get a copy of the diagonal entries of the matrix.
Teuchos_Ordinal Array_size_type
Size type for Teuchos Array objects.
virtual void removeEmptyProcessesInPlace(const Teuchos::RCP< const map_type > &newMap)
Remove processes owning zero rows from the Maps and their communicator.
device_type::execution_space execution_space
The Kokkos execution space.
Implementation details of Tpetra.
void insertGlobalValues(const GlobalOrdinal globalRow, const Teuchos::ArrayView< const GlobalOrdinal > &cols, const Teuchos::ArrayView< const Scalar > &vals)
Insert one or more entries into the matrix, using global column indices.
void reduce()
Sum values of a locally replicated multivector across all processes.
void fillLocalMatrix(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Fill data into the local matrix.
size_t global_size_t
Global size_t object.
LocalOrdinal sumIntoGlobalValues(const GlobalOrdinal globalRow, const Teuchos::ArrayView< const GlobalOrdinal > &cols, const Teuchos::ArrayView< const Scalar > &vals, const bool atomic=useAtomicUpdatesByDefault)
Sum into one or more sparse matrix entries, using global indices.
void merge2(IT1 &indResultOut, IT2 &valResultOut, IT1 indBeg, IT1 indEnd, IT2 valBeg, IT2 valEnd)
Merge values in place, additively, with the same index.
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView< GlobalOrdinal > &Indices, const Teuchos::ArrayView< Scalar > &Values, size_t &NumEntries) const =0
Get a copy of the given global row's entries.
Kokkos::StaticCrsGraph< LocalOrdinal, Kokkos::LayoutLeft, execution_space > local_graph_type
The type of the part of the sparse graph on each MPI process.
Traits class for "invalid" (flag) values of integer types that Tpetra uses as local ordinals or globa...
Insert new values that don't currently exist.
Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > createOneToOne(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &M)
Creates a one-to-one version of the given Map where each GID is owned by only one process...
void exportAndFillComplete(Teuchos::RCP< CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > &destMatrix, const export_type &exporter, const Teuchos::RCP< const map_type > &domainMap=Teuchos::null, const Teuchos::RCP< const map_type > &rangeMap=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null) const
Export from this to the given destination matrix, and make the result fill complete.
global_size_t getGlobalNumCols() const
The number of global columns in the matrix.
bool isConstantStride() const
Whether this multivector has constant stride between columns.
void applyTranspose(const MV &X_in, MV &Y_in, const Teuchos::ETransp mode, Scalar alpha, Scalar beta) const
Special case of apply() for mode != Teuchos::NO_TRANS.
void getLocalDiagOffsets(Teuchos::ArrayRCP< size_t > &offsets) const
Get offsets of the diagonal entries in the matrix.
global_size_t getGlobalNumRows() const
Number of global elements in the row map of this matrix.
bool isFillComplete() const
Whether the matrix is fill complete.
ESweepDirection
Sweep direction for Gauss-Seidel or Successive Over-Relaxation (SOR).
void allocateValues(ELocalGlobal lg, GraphAllocationStatus gas)
Allocate values (and optionally indices) using the Node.
Declare and define the function Tpetra::Details::computeOffsetsFromCounts, an implementation detail o...
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
void scale(const Scalar &alpha)
Scale the matrix's values: this := alpha*this.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
#define TPETRA_ABUSE_WARNING(throw_exception_test, Exception, msg)
Handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WA...
void unpackAndCombine(const Teuchos::ArrayView< const LocalOrdinal > &importLIDs, const Teuchos::ArrayView< const char > &imports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t constantNumPackets, Distributor &distor, CombineMode combineMode)
Unpack the imported column indices and values, and combine into matrix.
Sets up and executes a communication plan for a Tpetra DistObject.
bool isStorageOptimized() const
Returns true if storage has been optimized.
CombineMode
Rule for combining data in an Import or Export.
Sum new values into existing values.
void setAllValues(const typename local_matrix_type::row_map_type &ptr, const typename local_graph_type::entries_type::non_const_type &ind, const typename local_matrix_type::values_type &val)
Set the local matrix using three (compressed sparse row) arrays.
LocalOrdinal getLocalRowViewRaw(const LocalOrdinal lclRow, LocalOrdinal &numEnt, const LocalOrdinal *&lclColInds, const Scalar *&vals) const
Get a constant, nonpersisting, locally indexed view of the given row of the matrix, using "raw" pointers instead of Teuchos::ArrayView.
bool isUpperTriangular() const
Indicates whether the matrix is upper triangular.
OffsetsViewType::non_const_value_type computeOffsetsFromCounts(const OffsetsViewType &ptr, const CountsViewType &counts)
Compute offsets from counts.
Utility functions for packing and unpacking sparse matrix entries.
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
bool isStaticGraph() const
Indicates that the graph is static, so that new entries cannot be added to this matrix.
virtual ~CrsMatrix()
Destructor.
void getLocalRowView(LocalOrdinal LocalRow, Teuchos::ArrayView< const LocalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const
Get a constant, nonpersisting view of a row of this matrix, using local row and column indices...
Replace old value with maximum of magnitudes of old and new values.
GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const
The global index corresponding to the given local index.
Abstract base class for objects that can be the source of an Import or Export operation.
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
void replaceDomainMapAndImporter(const Teuchos::RCP< const map_type > &newDomainMap, Teuchos::RCP< const import_type > &newImporter)
Replace the current domain Map and Import with the given objects.
static LocalMapType::local_ordinal_type getDiagCopyWithoutOffsets(const DiagType &D, const LocalMapType &rowMap, const LocalMapType &colMap, const CrsMatrixType &A)
Given a locally indexed, local sparse matrix, and corresponding local row and column Maps...
LocalOrdinal getViewRaw(impl_scalar_type *&vals, LocalOrdinal &numEnt, const RowInfo &rowinfo) const
Nonconst pointer to all entries (including extra space) in the given row.
void fillComplete(const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Tell the matrix that you are done changing its structure or values, and that you are ready to do comp...
LO getLocalDiagCopyWithoutOffsetsNotFillComplete(::Tpetra::Vector< SC, LO, GO, NT > &diag, const ::Tpetra::RowMatrix< SC, LO, GO, NT > &A, const bool debug=false)
Given a locally indexed, global sparse matrix, extract the matrix's diagonal entries into a Tpetra::V...
Replace existing values with new values.
#define TPETRA_EFFICIENCY_WARNING(throw_exception_test, Exception, msg)
Print or throw an efficency warning.
void reindexColumns(const Teuchos::RCP< const map_type > &newColMap, const Teuchos::RCP< const import_type > &newImport=Teuchos::null, const bool sortIndicesInEachRow=true)
Reindex the column indices in place, and replace the column Map. Optionally, replace the Import objec...
void computeGlobalConstants()
Compute matrix properties that require collectives.
bool hasTransposeApply() const
Whether apply() allows applying the transpose or conjugate transpose.
Replace old values with zero.
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
void reorderedGaussSeidel(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Teuchos::ArrayView< LocalOrdinal > &rowIndices, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps) const
Reordered "Hybrid" Jacobi + (Gauss-Seidel or SOR) on .
void importAndFillComplete(Teuchos::RCP< CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > &destMatrix, const import_type &importer, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null) const
Import from this to the given destination matrix, and make the result fill complete.
void getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView< GlobalOrdinal > &Indices, const Teuchos::ArrayView< Scalar > &Values, size_t &NumEntries) const
Fill given arrays with a deep copy of the locally owned entries of the matrix in a given row...
void getLocalRowCopy(LocalOrdinal localRow, const Teuchos::ArrayView< LocalOrdinal > &colInds, const Teuchos::ArrayView< Scalar > &vals, size_t &numEntries) const
Fill given arrays with a deep copy of the locally owned entries of the matrix in a given row...
Kokkos::Details::ArithTraits< impl_scalar_type >::mag_type mag_type
Type of a norm result.
size_t getNodeNumRows() const
The number of matrix rows owned by the calling process.
void leftScale(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &x)
void sort2(const IT1 &first1, const IT1 &last1, const IT2 &first2)
Sort the first array, and apply the resulting permutation to the second array.
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
void checkInternalState() const
Check that this object's state is sane; throw if it's not.
void reorderedGaussSeidelCopy(MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Teuchos::ArrayView< LocalOrdinal > &rowIndices, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps, const bool zeroInitialGuess) const
Version of reorderedGaussSeidel(), with fewer requirements on X.
Describes a parallel distribution of objects over processes.
Teuchos::ArrayView< typename DualViewType::t_dev::value_type > getArrayViewFromDualView(const DualViewType &x)
Get a Teuchos::ArrayView which views the host Kokkos::View of the input 1-D Kokkos::DualView.
void applyNonTranspose(const MV &X_in, MV &Y_in, Scalar alpha, Scalar beta) const
Special case of apply() for mode == Teuchos::NO_TRANS.
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects...
size_t getNumVectors() const
Number of columns in the multivector.
Teuchos::RCP< Node > getNode() const
Get this Map's Node object.
A read-only, row-oriented interface to a sparse matrix.
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
Scalar operator()(const Scalar &x, const Scalar &y)
Return the maximum of the magnitudes (absolute values) of x and y.
A distributed dense vector.
bool isGloballyIndexed() const
Whether the matrix is globally indexed on the calling process.
void gaussSeidel(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps) const
"Hybrid" Jacobi + (Gauss-Seidel or SOR) on .
void insertLocalValues(const LocalOrdinal localRow, const Teuchos::ArrayView< const LocalOrdinal > &cols, const Teuchos::ArrayView< const Scalar > &vals)
Insert one or more entries into the matrix, using local column indices.
Teuchos::RCP< const map_type > getColMap() const
The Map that describes the column distribution in this matrix.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this matrix.
void doExport(const SrcDistObject &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)
Export data into this object using an Export object ("forward mode").
void globalAssemble()
Communicate nonlocal contributions to other processes.
Kokkos::Details::ArithTraits< Scalar >::val_type impl_scalar_type
The type used internally in place of Scalar.
Teuchos::RCP< MV > getRowMapMultiVector(const MV &Y_rangeMap, const bool force=false) const
Create a (or fetch a cached) row Map MultiVector.
bool isDistributed() const
Whether this is a globally distributed object.
Teuchos::RCP< const map_type > getRangeMap() const
The range Map of this matrix.
void expertStaticFillComplete(const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< const import_type > &importer=Teuchos::null, const Teuchos::RCP< const export_type > &exporter=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Perform a fillComplete on a matrix that already has data.
bool isNodeLocalElement(LocalOrdinal localIndex) const
Whether the given local index is valid for this Map on the calling process.
virtual Teuchos::RCP< RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > add(const Scalar &alpha, const RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms) const
Implementation of RowMatrix::add: return alpha*A + beta*this.
Teuchos::RCP< const crs_graph_type > getCrsGraph() const
This matrix's graph, as a CrsGraph.
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
void clearGlobalConstants()
Clear matrix properties that require collectives.
LocalOrdinal getViewRawConst(const impl_scalar_type *&vals, LocalOrdinal &numEnt, const RowInfo &rowinfo) const
Const pointer to all entries (including extra space) in the given row.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
The Map that describes the distribution of rows over processes.
virtual void pack(const Teuchos::ArrayView< const LocalOrdinal > &exportLIDs, Teuchos::Array< char > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t &constantNumPackets, Distributor &distor) const
Pack this object's data for an Import or Export.
bool isLocallyIndexed() const
Whether the matrix is locally indexed on the calling process.
local_matrix_type lclMatrix_
The local sparse matrix.
void fillLocalGraphAndMatrix(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Fill data into the local graph and matrix.
GlobalOrdinal getIndexBase() const
The index base for global indices for this matrix.