42 #ifndef TPETRA_CRSMATRIX_DEF_HPP
43 #define TPETRA_CRSMATRIX_DEF_HPP
53 #include "Tpetra_RowMatrix.hpp"
61 #include "Tpetra_Details_getDiagCopyWithoutOffsets.hpp"
62 #include "Tpetra_Details_gathervPrint.hpp"
63 #include "Tpetra_Details_leftScaleLocalCrsMatrix.hpp"
65 #include "Tpetra_Details_rightScaleLocalCrsMatrix.hpp"
68 #include "Teuchos_SerialDenseMatrix.hpp"
69 #include "KokkosSparse_getDiagCopy.hpp"
70 #include "Tpetra_Details_copyConvert.hpp"
72 #include "Tpetra_Details_packCrsMatrix.hpp"
73 #include "Tpetra_Details_unpackCrsMatrixAndCombine.hpp"
83 template<
class T,
class BinaryFunction>
84 T atomic_binary_function_update (
volatile T*
const dest,
98 T newVal = f (assume, inputVal);
99 oldVal = Kokkos::atomic_compare_exchange (dest, assume, newVal);
100 }
while (assume != oldVal);
120 template<
class Scalar>
124 typedef Teuchos::ScalarTraits<Scalar> STS;
125 return std::max (STS::magnitude (x), STS::magnitude (y));
135 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
138 size_t maxNumEntriesPerRow,
140 const Teuchos::RCP<Teuchos::ParameterList>& params) :
145 fillComplete_ (false),
146 frobNorm_ (-STM::one ())
148 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, size_t, "
149 "ProfileType[, RCP<ParameterList>]): ";
150 Teuchos::RCP<crs_graph_type> graph;
152 graph = Teuchos::rcp (
new crs_graph_type (rowMap, maxNumEntriesPerRow,
155 catch (std::exception& e) {
156 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
157 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
158 "size_t, ProfileType[, RCP<ParameterList>]) threw an exception: "
165 staticGraph_ = myGraph_;
167 checkInternalState ();
170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
173 const Teuchos::ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
175 const Teuchos::RCP<Teuchos::ParameterList>& params) :
180 fillComplete_ (false),
181 frobNorm_ (-STM::one ())
183 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, "
184 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]): ";
185 Teuchos::RCP<crs_graph_type> graph;
187 graph = Teuchos::rcp (
new crs_graph_type (rowMap, NumEntriesPerRowToAlloc,
190 catch (std::exception &e) {
191 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
192 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
193 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]) threw "
194 "an exception: " << e.what ());
200 staticGraph_ = graph;
202 checkInternalState ();
205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
208 const Teuchos::RCP<const map_type>& colMap,
209 size_t maxNumEntriesPerRow,
211 const Teuchos::RCP<Teuchos::ParameterList>& params) :
216 fillComplete_ (false),
217 frobNorm_ (-STM::one ())
219 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, RCP<const Map>, "
220 "size_t, ProfileType[, RCP<ParameterList>]): ";
222 #ifdef HAVE_TPETRA_DEBUG
224 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
225 (! staticGraph_.is_null (), std::logic_error,
226 "staticGraph_ is not null at the beginning of the constructor. "
227 "Please report this bug to the Tpetra developers.");
228 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
229 (! myGraph_.is_null (), std::logic_error,
230 "myGraph_ is not null at the beginning of the constructor. "
231 "Please report this bug to the Tpetra developers.");
232 #endif // HAVE_TPETRA_DEBUG
234 Teuchos::RCP<crs_graph_type> graph;
240 catch (std::exception &e) {
241 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
242 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
243 "RCP<const Map>, size_t, ProfileType[, RCP<ParameterList>]) threw an "
244 "exception: " << e.what ());
250 staticGraph_ = myGraph_;
252 checkInternalState ();
255 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
258 const Teuchos::RCP<const map_type>& colMap,
259 const Teuchos::ArrayRCP<const size_t>& numEntPerRow,
261 const Teuchos::RCP<Teuchos::ParameterList>& params) :
266 fillComplete_ (false),
267 frobNorm_ (-STM::one ())
269 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, RCP<const Map>, "
270 "ArrayRCP<const size_t>, ProfileType[, RCP<ParameterList>]): ";
271 Teuchos::RCP<crs_graph_type> graph;
273 graph = Teuchos::rcp (
new crs_graph_type (rowMap, colMap, numEntPerRow,
276 catch (std::exception &e) {
277 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
278 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
279 "RCP<const Map>, ArrayRCP<const size_t>, ProfileType[, "
280 "RCP<ParameterList>]) threw an exception: " << e.what ());
286 staticGraph_ = graph;
288 checkInternalState ();
291 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
293 CrsMatrix (
const Teuchos::RCP<const crs_graph_type>& graph,
294 const Teuchos::RCP<Teuchos::ParameterList>& ) :
296 staticGraph_ (graph),
298 fillComplete_ (false),
299 frobNorm_ (-STM::one ())
301 typedef typename local_matrix_type::values_type values_type;
302 const char tfecfFuncName[] =
"CrsMatrix(RCP<const CrsGraph>[, "
303 "RCP<ParameterList>]): ";
304 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
305 (graph.is_null (), std::runtime_error,
"Input graph is null.");
306 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
307 (! graph->isFillComplete (), std::runtime_error,
"Input graph is not "
308 "fill complete. You must call fillComplete on the graph before using "
309 "it to construct a CrsMatrix. Note that calling resumeFill on the "
310 "graph makes it not fill complete, even if you had previously called "
311 "fillComplete. In that case, you must call fillComplete on the graph "
320 const size_t numCols = graph->getColMap ()->getNodeNumElements ();
321 auto lclGraph = graph->getLocalGraph ();
322 const size_t numEnt = lclGraph.entries.extent (0);
323 values_type val (
"Tpetra::CrsMatrix::val", numEnt);
326 numCols, val, lclGraph);
330 this->k_values1D_ = this->lclMatrix_.values;
332 checkInternalState ();
335 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
338 const Teuchos::RCP<const map_type>& colMap,
339 const typename local_matrix_type::row_map_type& rowPointers,
340 const typename local_graph_type::entries_type::non_const_type& columnIndices,
341 const typename local_matrix_type::values_type& values,
342 const Teuchos::RCP<Teuchos::ParameterList>& params) :
345 fillComplete_ (false),
346 frobNorm_ (-STM::one ())
349 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
350 "RCP<const Map>, ptr, ind, val[, params]): ";
351 const char suffix[] =
". Please report this bug to the Tpetra developers.";
357 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
358 (values.extent (0) != columnIndices.extent (0),
359 std::invalid_argument,
"Input arrays don't have matching dimensions. "
360 "values.extent(0) = " << values.extent (0) <<
" != "
361 "columnIndices.extent(0) = " << columnIndices.extent (0) <<
".");
362 #ifdef HAVE_TPETRA_DEBUG
363 if (rowPointers.extent (0) != 0) {
364 const size_t numEnt =
365 ::Tpetra::Details::getEntryOnHost (rowPointers, rowPointers.extent (0) - 1);
366 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
367 (numEnt != static_cast<size_t> (columnIndices.extent (0)) ||
368 numEnt != static_cast<size_t> (values.extent (0)),
369 std::invalid_argument,
"Last entry of rowPointers says that the matrix"
370 " has " << numEnt <<
" entr" << (numEnt != 1 ?
"ies" :
"y") <<
", but "
371 "the dimensions of columnIndices and values don't match this. "
372 "columnIndices.extent(0) = " << columnIndices.extent (0) <<
373 " and values.extent(0) = " << values.extent (0) <<
".");
375 #endif // HAVE_TPETRA_DEBUG
377 RCP<crs_graph_type> graph;
379 graph = Teuchos::rcp (
new crs_graph_type (rowMap, colMap, rowPointers,
380 columnIndices, params));
382 catch (std::exception& e) {
383 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
384 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
385 "RCP<const Map>, ptr, ind[, params]) threw an exception: "
393 auto lclGraph = graph->getLocalGraph ();
394 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
395 (lclGraph.row_map.extent (0) != rowPointers.extent (0) ||
396 lclGraph.entries.extent (0) != columnIndices.extent (0),
397 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, "
398 "ind[, params]) did not set the local graph correctly." << suffix);
399 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
400 (lclGraph.entries.extent (0) != values.extent (0),
401 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, ind[, "
402 "params]) did not set the local graph correctly. "
403 "lclGraph.entries.extent(0) = " << lclGraph.entries.extent (0)
404 <<
" != values.extent(0) = " << values.extent (0) << suffix);
410 staticGraph_ = graph;
419 const size_t numCols = graph->getColMap ()->getNodeNumElements ();
421 numCols, values, lclGraph);
422 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
423 (lclMatrix_.values.extent (0) != values.extent (0),
424 std::logic_error,
"Local matrix's constructor did not set the values "
425 "correctly. lclMatrix_.values.extent(0) = " <<
426 lclMatrix_.values.extent (0) <<
" != values.extent(0) = " <<
427 values.extent (0) << suffix);
432 this->k_values1D_ = this->lclMatrix_.values;
434 checkInternalState ();
437 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
440 const Teuchos::RCP<const map_type>& colMap,
441 const Teuchos::ArrayRCP<size_t>& ptr,
442 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
443 const Teuchos::ArrayRCP<Scalar>& val,
444 const Teuchos::RCP<Teuchos::ParameterList>& params) :
447 fillComplete_ (false),
448 frobNorm_ (-STM::one ())
450 using Kokkos::Compat::getKokkosViewDeepCopy;
451 using Teuchos::av_reinterpret_cast;
453 typedef typename local_matrix_type::values_type values_type;
455 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
456 "RCP<const Map>, ptr, ind, val[, params]): ";
458 RCP<crs_graph_type> graph;
463 catch (std::exception& e) {
464 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
465 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
466 "RCP<const Map>, ArrayRCP<size_t>, ArrayRCP<LocalOrdinal>[, "
467 "RCP<ParameterList>]) threw an exception: " << e.what ());
473 staticGraph_ = graph;
486 auto lclGraph = staticGraph_->getLocalGraph ();
487 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
488 (static_cast<size_t> (lclGraph.row_map.extent (0)) != static_cast<size_t> (ptr.size ()) ||
489 static_cast<size_t> (lclGraph.entries.extent (0)) != static_cast<size_t> (ind.size ()),
490 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, "
491 "ind[, params]) did not set the local graph correctly. Please "
492 "report this bug to the Tpetra developers.");
494 const size_t numCols = staticGraph_->getColMap ()->getNodeNumElements ();
495 values_type valIn = getKokkosViewDeepCopy<device_type> (av_reinterpret_cast<IST> (val ()));
497 numCols, valIn, lclGraph);
501 this->k_values1D_ = this->lclMatrix_.values;
503 checkInternalState ();
506 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
509 const Teuchos::RCP<const map_type>& colMap,
511 const Teuchos::RCP<Teuchos::ParameterList>& params) :
513 lclMatrix_ (lclMatrix),
514 k_values1D_ (lclMatrix.values),
516 fillComplete_ (true),
517 frobNorm_ (-STM::one ())
519 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
520 "RCP<const Map>, local_matrix_type[, RCP<ParameterList>]): ";
521 Teuchos::RCP<crs_graph_type> graph;
524 lclMatrix.graph, params));
526 catch (std::exception& e) {
527 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
528 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
529 "RCP<const Map>, local_graph_type[, RCP<ParameterList>]) threw an "
530 "exception: " << e.what ());
532 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
533 (!graph->isFillComplete (), std::logic_error,
"CrsGraph constructor (RCP"
534 "<const Map>, RCP<const Map>, local_graph_type[, RCP<ParameterList>]) "
535 "did not produce a fill-complete graph. Please report this bug to the "
536 "Tpetra developers.");
541 staticGraph_ = graph;
543 const bool callComputeGlobalConstants = params.get () ==
nullptr ||
544 params->get (
"compute global constants",
true);
545 if (callComputeGlobalConstants) {
546 this->computeGlobalConstants ();
550 #ifdef HAVE_TPETRA_DEBUG
551 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isFillActive (), std::logic_error,
552 "We're at the end of fillComplete(), but isFillActive() is true. "
553 "Please report this bug to the Tpetra developers.");
554 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillComplete (), std::logic_error,
555 "We're at the end of fillComplete(), but isFillComplete() is false. "
556 "Please report this bug to the Tpetra developers.");
557 #endif // HAVE_TPETRA_DEBUG
558 checkInternalState ();
561 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
564 const Teuchos::RCP<const map_type>& rowMap,
565 const Teuchos::RCP<const map_type>& colMap,
566 const Teuchos::RCP<const map_type>& domainMap,
567 const Teuchos::RCP<const map_type>& rangeMap,
568 const Teuchos::RCP<Teuchos::ParameterList>& params) :
570 lclMatrix_ (lclMatrix),
571 k_values1D_ (lclMatrix.values),
573 fillComplete_ (true),
574 frobNorm_ (-STM::one ())
576 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
577 "RCP<const Map>, RCP<const Map>, RCP<const Map>, local_matrix_type[, "
578 "RCP<ParameterList>]): ";
579 Teuchos::RCP<crs_graph_type> graph;
581 graph = Teuchos::rcp (
new crs_graph_type (lclMatrix.graph, rowMap, colMap,
582 domainMap, rangeMap, params));
584 catch (std::exception& e) {
585 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
586 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
587 "RCP<const Map>, RCP<const Map>, RCP<const Map>, local_graph_type[, "
588 "RCP<ParameterList>]) threw an exception: " << e.what ());
590 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
591 (!graph->isFillComplete (), std::logic_error,
"CrsGraph constructor (RCP"
592 "<const Map>, RCP<const Map>, RCP<const Map>, RCP<const Map>, local_graph_type[, "
593 "RCP<ParameterList>]) did not produce a fill-complete graph. Please report this "
594 "bug to the Tpetra developers.");
599 staticGraph_ = graph;
601 const bool callComputeGlobalConstants = params.get () ==
nullptr ||
602 params->get (
"compute global constants",
true);
603 if (callComputeGlobalConstants) {
604 this->computeGlobalConstants ();
608 #ifdef HAVE_TPETRA_DEBUG
609 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isFillActive (), std::logic_error,
610 "We're at the end of fillComplete(), but isFillActive() is true. "
611 "Please report this bug to the Tpetra developers.");
612 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillComplete (), std::logic_error,
613 "We're at the end of fillComplete(), but isFillComplete() is false. "
614 "Please report this bug to the Tpetra developers.");
615 #endif // HAVE_TPETRA_DEBUG
616 checkInternalState ();
619 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
624 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
625 Teuchos::RCP<const Teuchos::Comm<int> >
628 return getCrsGraphRef ().
getComm ();
631 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
635 return getCrsGraphRef ().
getNode ();
638 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
645 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
649 return fillComplete_;
652 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
656 return ! fillComplete_;
659 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
666 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
673 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
680 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
687 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
694 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
701 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
708 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
715 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
722 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
729 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
735 return dynamic_cast<const HDM&> (G).getGlobalNumDiagsImpl ();
738 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
742 return this->getGlobalNumDiagsImpl ();
745 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
751 return dynamic_cast<const HDM&> (G).getNodeNumDiagsImpl ();
754 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
758 return this->getNodeNumDiagsImpl ();
761 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
768 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
775 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
782 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
789 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
796 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
797 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
803 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
804 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
810 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
811 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
817 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
818 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
824 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
825 Teuchos::RCP<const RowGraph<LocalOrdinal, GlobalOrdinal, Node> >
828 if (staticGraph_ != Teuchos::null) {
834 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
835 Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> >
838 if (staticGraph_ != Teuchos::null) {
844 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
848 if (! this->staticGraph_.is_null ()) {
849 return * (this->staticGraph_);
852 #ifdef HAVE_TPETRA_DEBUG
853 const char tfecfFuncName[] =
"getCrsGraphRef: ";
854 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
855 (this->myGraph_.is_null (), std::logic_error,
856 "Both staticGraph_ and myGraph_ are null. "
857 "Please report this bug to the Tpetra developers.");
858 #endif // HAVE_TPETRA_DEBUG
859 return * (this->myGraph_);
863 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
865 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
866 isLowerTriangularImpl ()
const {
869 return dynamic_cast<const HDM&> (G).isLowerTriangularImpl ();
872 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
876 return this->isLowerTriangularImpl ();
879 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
885 return dynamic_cast<const HDM&> (G).isUpperTriangularImpl ();
888 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
892 return this->isUpperTriangularImpl ();
895 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
899 return myGraph_.is_null ();
902 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
909 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
916 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
917 Teuchos::ArrayRCP<Teuchos::Array<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::impl_scalar_type> >
922 using Teuchos::Array;
923 using Teuchos::ArrayRCP;
924 typedef impl_scalar_type IST;
925 typedef LocalOrdinal LO;
926 const char tfecfFuncName[] =
"allocateValues2D: ";
928 const crs_graph_type& graph = this->getCrsGraphRef ();
929 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
930 (! graph.indicesAreAllocated (), std::runtime_error,
931 "Graph indices must be allocated before values.");
932 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
934 "Graph indices must be allocated in a dynamic profile.");
936 const LO lclNumRows = graph.getNodeNumRows ();
937 Teuchos::ArrayRCP<Teuchos::Array<IST> > values2D (lclNumRows);
938 if (! graph.lclInds2D_.is_null ()) {
939 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
940 values2D[lclRow].resize (graph.lclInds2D_[lclRow].size ());
943 else if (! graph.gblInds2D_.is_null ()) {
944 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
945 values2D[lclRow].resize (graph.gblInds2D_[lclRow].size ());
951 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
953 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
954 allocateValues (ELocalGlobal lg, GraphAllocationStatus gas)
956 using ::Tpetra::Details::ProfilingRegion;
957 const char tfecfFuncName[] =
"allocateValues: ";
958 ProfilingRegion regionAllocateValues (
"Tpetra::CrsMatrix::allocateValues");
960 #ifdef HAVE_TPETRA_DEBUG
961 const char suffix[] =
" Please report this bug to the Tpetra developers.";
963 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
964 (this->staticGraph_.is_null (), std::logic_error,
965 "staticGraph_ is null." << suffix);
970 if ((gas == GraphAlreadyAllocated) != this->staticGraph_->indicesAreAllocated ()) {
971 const char err1[] =
"The caller has asserted that the graph is ";
972 const char err2[] =
"already allocated, but the static graph says "
973 "that its indices are ";
974 const char err3[] =
"already allocated. Please report this bug to "
975 "the Tpetra developers.";
976 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
977 (gas == GraphAlreadyAllocated && ! this->staticGraph_->indicesAreAllocated (),
978 std::logic_error, err1 << err2 <<
"not " << err3);
979 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
980 (gas != GraphAlreadyAllocated && this->staticGraph_->indicesAreAllocated (),
981 std::logic_error, err1 <<
"not " << err2 << err3);
989 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
990 (! this->staticGraph_->indicesAreAllocated () &&
991 this->myGraph_.is_null (), std::logic_error,
992 "The static graph says that its indices are not allocated, "
993 "but the graph is not owned by the matrix." << suffix);
994 #endif // HAVE_TPETRA_DEBUG
996 if (gas == GraphNotYetAllocated) {
997 #ifdef HAVE_TPETRA_DEBUG
998 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
999 (this->myGraph_.is_null (), std::logic_error,
1000 "gas = GraphNotYetAllocated, but myGraph_ is null." << suffix);
1001 #endif // HAVE_TPETRA_DEBUG
1003 this->myGraph_->allocateIndices (lg);
1005 catch (std::exception& e) {
1006 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1007 (
true, std::runtime_error,
"CrsGraph::allocateIndices "
1008 "threw an exception: " << e.what ());
1011 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1012 (
true, std::runtime_error,
"CrsGraph::allocateIndices "
1013 "threw an exception not a subclass of std::exception.");
1025 #ifdef HAVE_TPETRA_DEBUG
1026 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1027 (this->staticGraph_.is_null (), std::logic_error,
1028 "this->getProfileType() == StaticProfile, but staticGraph_ is null."
1030 #endif // HAVE_TPETRA_DEBUG
1032 const size_t lclNumRows = this->staticGraph_->getNodeNumRows ();
1033 typename Graph::local_graph_type::row_map_type k_ptrs =
1034 this->staticGraph_->k_rowPtrs_;
1035 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1036 (k_ptrs.extent (0) != lclNumRows+1, std::logic_error,
1037 "With StaticProfile, row offsets array has length "
1038 << k_ptrs.extent (0) <<
" != (lclNumRows+1) = "
1039 << (lclNumRows+1) <<
".");
1041 const size_t lclTotalNumEntries =
1042 ::Tpetra::Details::getEntryOnHost (k_ptrs, lclNumRows);
1045 typedef typename local_matrix_type::values_type values_type;
1047 values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1055 this->values2D_ = this->allocateValues2D ();
1059 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1062 getAllValues (Teuchos::ArrayRCP<const size_t>& rowPointers,
1063 Teuchos::ArrayRCP<const LocalOrdinal>& columnIndices,
1064 Teuchos::ArrayRCP<const Scalar>& values)
const
1067 const char tfecfFuncName[] =
"getAllValues: ";
1068 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1069 columnIndices.size () != values.size (), std::runtime_error,
1070 "Requires that columnIndices and values are the same size.");
1072 RCP<const crs_graph_type> relevantGraph = getCrsGraph ();
1073 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1074 relevantGraph.is_null (), std::runtime_error,
1075 "Requires that getCrsGraph() is not null.");
1077 rowPointers = relevantGraph->getNodeRowPtrs ();
1079 catch (std::exception &e) {
1080 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1081 true, std::runtime_error,
1082 "Caught exception while calling graph->getNodeRowPtrs(): "
1086 columnIndices = relevantGraph->getNodePackedIndices ();
1088 catch (std::exception &e) {
1089 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1090 true, std::runtime_error,
1091 "Caught exception while calling graph->getNodePackedIndices(): "
1094 Teuchos::ArrayRCP<const impl_scalar_type> vals =
1095 Kokkos::Compat::persistingView (k_values1D_);
1096 values = Teuchos::arcp_reinterpret_cast<const Scalar> (vals);
1099 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1101 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1102 fillLocalGraphAndMatrix (
const Teuchos::RCP<Teuchos::ParameterList>& params)
1105 using ::Tpetra::Details::ProfilingRegion;
1106 using Kokkos::create_mirror_view;
1107 using Teuchos::arcp_const_cast;
1108 using Teuchos::Array;
1109 using Teuchos::ArrayRCP;
1110 using Teuchos::null;
1113 typedef typename local_matrix_type::row_map_type row_map_type;
1114 typedef typename Graph::local_graph_type::entries_type::non_const_type lclinds_1d_type;
1115 typedef typename local_matrix_type::values_type values_type;
1116 ProfilingRegion regionFLGAM (
"Tpetra::CrsGraph::fillLocalGraphAndMatrix");
1118 #ifdef HAVE_TPETRA_DEBUG
1119 const char tfecfFuncName[] =
"fillLocalGraphAndMatrix (called from "
1120 "fillComplete or expertStaticFillComplete): ";
1121 #endif // HAVE_TPETRA_DEBUG
1123 #ifdef HAVE_TPETRA_DEBUG
1126 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1127 (myGraph_.is_null (), std::logic_error,
"The nonconst graph (myGraph_) "
1128 "is null. This means that the matrix has a const (a.k.a. \"static\") "
1129 "graph. fillComplete or expertStaticFillComplete should never call "
1130 "fillLocalGraphAndMatrix in that case. "
1131 "Please report this bug to the Tpetra developers.");
1132 #endif // HAVE_TPETRA_DEBUG
1134 const size_t lclNumRows = this->getNodeNumRows ();
1142 typename row_map_type::non_const_type k_ptrs;
1143 row_map_type k_ptrs_const;
1144 lclinds_1d_type k_inds;
1150 lclinds_1d_type k_lclInds1D_ = myGraph_->k_lclInds1D_;
1152 typedef decltype (myGraph_->k_numRowEntries_) row_entries_type;
1173 typename row_entries_type::const_type numRowEnt_h =
1174 myGraph_->k_numRowEntries_;
1175 #ifdef HAVE_TPETRA_DEBUG
1176 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1177 (static_cast<size_t> (numRowEnt_h.extent (0)) != lclNumRows,
1178 std::logic_error,
"(DynamicProfile branch) numRowEnt_h has the "
1179 "wrong length. numRowEnt_h.extent(0) = "
1180 << numRowEnt_h.extent (0) <<
" != getNodeNumRows() = "
1181 << lclNumRows <<
".");
1182 #endif // HAVE_TPETRA_DEBUG
1187 k_ptrs =
typename row_map_type::non_const_type (
"Tpetra::CrsGraph::ptr",
1189 typename row_map_type::non_const_type::HostMirror h_ptrs =
1190 create_mirror_view (k_ptrs);
1198 const size_t lclTotalNumEntries =
1200 #ifdef HAVE_TPETRA_DEBUG
1201 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1202 (static_cast<size_t> (h_ptrs.extent (0)) != lclNumRows + 1,
1203 std::logic_error,
"(DynamicProfile branch) After packing h_ptrs, "
1204 "h_ptrs.extent(0) = " << h_ptrs.extent (0) <<
" != "
1205 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
1207 const size_t h_ptrs_lastEnt = h_ptrs(lclNumRows);
1208 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1209 (h_ptrs_lastEnt != lclTotalNumEntries, std::logic_error,
1210 "(DynamicProfile branch) After packing h_ptrs, h_ptrs(lclNumRows="
1211 << lclNumRows <<
") = " << h_ptrs_lastEnt <<
" != total number "
1212 "of entries on the calling process = " << lclTotalNumEntries <<
".");
1214 #endif // HAVE_TPETRA_DEBUG
1217 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
1218 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1221 typename lclinds_1d_type::HostMirror h_inds = create_mirror_view (k_inds);
1222 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
1225 ArrayRCP<Array<LocalOrdinal> > lclInds2D = myGraph_->lclInds2D_;
1226 for (
size_t row = 0; row < lclNumRows; ++row) {
1227 const size_t numEnt = numRowEnt_h(row);
1228 std::copy (lclInds2D[row].begin(),
1229 lclInds2D[row].begin() + numEnt,
1230 h_inds.data() + h_ptrs(row));
1231 std::copy (values2D_[row].begin(),
1232 values2D_[row].begin() + numEnt,
1233 h_vals.data() + h_ptrs(row));
1242 k_ptrs_const = k_ptrs;
1244 #ifdef HAVE_TPETRA_DEBUG
1246 if (k_ptrs.extent (0) != 0) {
1247 const size_t numOffsets = static_cast<size_t> (k_ptrs.extent (0));
1248 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1249 (numOffsets != lclNumRows + 1, std::logic_error,
"(DynamicProfile "
1250 "branch) After copying into k_ptrs, k_ptrs.extent(0) = " <<
1251 numOffsets <<
" != (lclNumRows+1) = " << (lclNumRows+1) <<
".");
1253 const auto valToCheck = ::Tpetra::Details::getEntryOnHost (k_ptrs, numOffsets-1);
1254 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1255 (static_cast<size_t> (valToCheck) != k_vals.extent (0),
1256 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs("
1257 << (numOffsets-1) <<
") = " << valToCheck <<
" != "
1258 "k_vals.extent(0) = " << k_vals.extent (0) <<
".");
1259 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1260 (static_cast<size_t> (valToCheck) != k_inds.extent (0),
1261 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs("
1262 << (numOffsets-1) <<
") = " << valToCheck <<
" != "
1263 "k_inds.extent(0) = " << k_inds.extent (0) <<
".");
1265 #endif // HAVE_TPETRA_DEBUG
1274 typename Graph::local_graph_type::row_map_type curRowOffsets =
1275 myGraph_->k_rowPtrs_;
1277 #ifdef HAVE_TPETRA_DEBUG
1278 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1279 (curRowOffsets.extent (0) == 0, std::logic_error,
1280 "(StaticProfile branch) curRowOffsets.extent(0) == 0.");
1281 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1282 (curRowOffsets.extent (0) != lclNumRows + 1, std::logic_error,
1283 "(StaticProfile branch) curRowOffsets.extent(0) = "
1284 << curRowOffsets.extent (0) <<
" != lclNumRows + 1 = "
1285 << (lclNumRows + 1) <<
".")
1287 const size_t numOffsets = curRowOffsets.extent (0);
1288 const auto valToCheck =
1289 ::Tpetra::Details::getEntryOnHost (curRowOffsets, numOffsets - 1);
1290 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1292 myGraph_->k_lclInds1D_.extent (0) != valToCheck,
1293 std::logic_error,
"(StaticProfile branch) numOffsets = " <<
1294 numOffsets <<
" != 0 and myGraph_->k_lclInds1D_.extent(0) = "
1295 << myGraph_->k_lclInds1D_.extent (0) <<
" != curRowOffsets("
1296 << numOffsets <<
") = " << valToCheck <<
".");
1298 #endif // HAVE_TPETRA_DEBUG
1300 if (myGraph_->getNodeNumEntries () != myGraph_->getNodeAllocationSize ()) {
1307 #ifdef HAVE_TPETRA_DEBUG
1308 if (curRowOffsets.extent (0) != 0) {
1309 const size_t numOffsets =
1310 static_cast<size_t> (curRowOffsets.extent (0));
1311 const auto valToCheck =
1312 ::Tpetra::Details::getEntryOnHost (curRowOffsets, numOffsets-1);
1313 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1314 (static_cast<size_t> (valToCheck) !=
1315 static_cast<size_t> (k_values1D_.extent (0)),
1316 std::logic_error,
"(StaticProfile unpacked branch) Before "
1317 "allocating or packing, curRowOffsets(" << (numOffsets-1) <<
") = "
1318 << valToCheck <<
" != k_values1D_.extent(0)"
1319 " = " << k_values1D_.extent (0) <<
".");
1320 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1321 (static_cast<size_t> (valToCheck) !=
1322 static_cast<size_t> (myGraph_->k_lclInds1D_.extent (0)),
1323 std::logic_error,
"(StaticProfile unpacked branch) Before "
1324 "allocating or packing, curRowOffsets(" << (numOffsets-1) <<
") = "
1326 <<
" != myGraph_->k_lclInds1D_.extent(0) = "
1327 << myGraph_->k_lclInds1D_.extent (0) <<
".");
1329 #endif // HAVE_TPETRA_DEBUG
1337 size_t lclTotalNumEntries = 0;
1339 typename row_map_type::non_const_type::HostMirror h_ptrs;
1344 typename row_map_type::non_const_type
1345 packedRowOffsets (
"Tpetra::CrsGraph::ptr", lclNumRows + 1);
1346 typename row_entries_type::const_type numRowEnt_h =
1347 myGraph_->k_numRowEntries_;
1350 lclTotalNumEntries =
1354 k_ptrs = packedRowOffsets;
1355 k_ptrs_const = k_ptrs;
1358 #ifdef HAVE_TPETRA_DEBUG
1359 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1360 (static_cast<size_t> (k_ptrs.extent (0)) != lclNumRows + 1,
1362 "(StaticProfile unpacked branch) After packing k_ptrs, "
1363 "k_ptrs.extent(0) = " << k_ptrs.extent (0) <<
" != "
1364 "lclNumRows+1 = " << (lclNumRows+1) <<
".");
1366 const auto valToCheck = ::Tpetra::Details::getEntryOnHost (k_ptrs, lclNumRows);
1367 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1368 (valToCheck != lclTotalNumEntries, std::logic_error,
1369 "(StaticProfile unpacked branch) After filling k_ptrs, "
1370 "k_ptrs(lclNumRows=" << lclNumRows <<
") = " << valToCheck
1371 <<
" != total number of entries on the calling process = "
1372 << lclTotalNumEntries <<
".");
1374 #endif // HAVE_TPETRA_DEBUG
1377 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
1378 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1390 typedef pack_functor<
typename Graph::local_graph_type::entries_type::non_const_type,
1391 typename Graph::local_graph_type::row_map_type>
1393 inds_packer_type indsPacker (k_inds, myGraph_->k_lclInds1D_,
1394 k_ptrs, curRowOffsets);
1396 typedef Kokkos::RangePolicy<exec_space, LocalOrdinal> range_type;
1397 Kokkos::parallel_for (range_type (0, lclNumRows), indsPacker);
1401 typedef pack_functor<values_type, row_map_type> vals_packer_type;
1402 vals_packer_type valsPacker (k_vals, this->k_values1D_,
1403 k_ptrs, curRowOffsets);
1404 Kokkos::parallel_for (range_type (0, lclNumRows), valsPacker);
1406 #ifdef HAVE_TPETRA_DEBUG
1407 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1408 (k_ptrs.extent (0) == 0, std::logic_error,
1409 "(StaticProfile \"Optimize Storage\" = "
1410 "true branch) After packing, k_ptrs.extent(0) = 0. This "
1411 "probably means that k_rowPtrs_ was never allocated.");
1412 if (k_ptrs.extent (0) != 0) {
1413 const size_t numOffsets = static_cast<size_t> (k_ptrs.extent (0));
1414 const auto valToCheck = ::Tpetra::Details::getEntryOnHost (k_ptrs, numOffsets - 1);
1415 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1416 (static_cast<size_t> (valToCheck) != k_vals.extent (0),
1418 "(StaticProfile \"Optimize Storage\"=true branch) After packing, "
1419 "k_ptrs(" << (numOffsets-1) <<
") = " << valToCheck <<
1420 " != k_vals.extent(0) = " << k_vals.extent (0) <<
".");
1421 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1422 (static_cast<size_t> (valToCheck) != k_inds.extent (0),
1424 "(StaticProfile \"Optimize Storage\"=true branch) After packing, "
1425 "k_ptrs(" << (numOffsets-1) <<
") = " << valToCheck <<
1426 " != k_inds.extent(0) = " << k_inds.extent (0) <<
".");
1428 #endif // HAVE_TPETRA_DEBUG
1431 k_ptrs_const = myGraph_->k_rowPtrs_;
1432 k_inds = myGraph_->k_lclInds1D_;
1433 k_vals = this->k_values1D_;
1435 #ifdef HAVE_TPETRA_DEBUG
1436 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1437 (k_ptrs_const.extent (0) == 0, std::logic_error,
1438 "(StaticProfile \"Optimize Storage\"=false branch) "
1439 "k_ptrs_const.extent(0) = 0. This probably means that "
1440 "k_rowPtrs_ was never allocated.");
1441 if (k_ptrs_const.extent (0) != 0) {
1442 const size_t numOffsets = static_cast<size_t> (k_ptrs_const.extent (0));
1443 const auto valToCheck = ::Tpetra::Details::getEntryOnHost (k_ptrs_const, numOffsets - 1);
1444 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1445 (static_cast<size_t> (valToCheck) != k_vals.extent (0),
1447 "(StaticProfile \"Optimize Storage\"=false branch) "
1448 "k_ptrs_const(" << (numOffsets-1) <<
") = " << valToCheck
1449 <<
" != k_vals.extent(0) = " << k_vals.extent (0) <<
".");
1450 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1451 (static_cast<size_t> (valToCheck) != k_inds.extent (0),
1453 "(StaticProfile \"Optimize Storage\" = false branch) "
1454 "k_ptrs_const(" << (numOffsets-1) <<
") = " << valToCheck
1455 <<
" != k_inds.extent(0) = " << k_inds.extent (0) <<
".");
1457 #endif // HAVE_TPETRA_DEBUG
1461 #ifdef HAVE_TPETRA_DEBUG
1463 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1464 (static_cast<size_t> (k_ptrs_const.extent (0)) != lclNumRows + 1,
1465 std::logic_error,
"After packing, k_ptrs_const.extent(0) = " <<
1466 k_ptrs_const.extent (0) <<
" != lclNumRows+1 = " << (lclNumRows+1)
1468 if (k_ptrs_const.extent (0) != 0) {
1469 const size_t numOffsets = static_cast<size_t> (k_ptrs_const.extent (0));
1470 const size_t k_ptrs_const_numOffsetsMinus1 =
1471 ::Tpetra::Details::getEntryOnHost (k_ptrs_const, numOffsets - 1);
1472 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1473 (k_ptrs_const_numOffsetsMinus1 != k_vals.extent (0),
1474 std::logic_error,
"After packing, k_ptrs_const(" << (numOffsets-1) <<
1475 ") = " << k_ptrs_const_numOffsetsMinus1 <<
" != k_vals.extent(0)"
1476 " = " << k_vals.extent (0) <<
".");
1477 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1478 (k_ptrs_const_numOffsetsMinus1 != k_inds.extent (0),
1479 std::logic_error,
"After packing, k_ptrs_const(" << (numOffsets-1) <<
1480 ") = " << k_ptrs_const_numOffsetsMinus1 <<
" != k_inds.extent(0)"
1481 " = " << k_inds.extent (0) <<
".");
1483 #endif // HAVE_TPETRA_DEBUG
1489 const bool defaultOptStorage =
1490 ! isStaticGraph () || staticGraph_->isStorageOptimized ();
1491 const bool requestOptimizedStorage =
1492 (! params.is_null () && params->get (
"Optimize Storage", defaultOptStorage)) ||
1493 (params.is_null () && defaultOptStorage);
1500 if (requestOptimizedStorage) {
1506 myGraph_->lclInds2D_ =
null;
1507 myGraph_->k_numRowEntries_ = row_entries_type ();
1510 this->values2D_ =
null;
1513 myGraph_->k_rowPtrs_ = k_ptrs_const;
1514 myGraph_->k_lclInds1D_ = k_inds;
1515 this->k_values1D_ = k_vals;
1519 myGraph_->storageStatus_ = ::Tpetra::Details::STORAGE_1D_PACKED;
1520 this->storageStatus_ = ::Tpetra::Details::STORAGE_1D_PACKED;
1531 myGraph_->lclGraph_ =
1536 getNodeNumCols (), k_vals,
1537 myGraph_->lclGraph_);
1540 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1545 using ::Tpetra::Details::ProfilingRegion;
1546 using Kokkos::create_mirror_view;
1547 using Teuchos::ArrayRCP;
1548 using Teuchos::Array;
1549 using Teuchos::null;
1552 typedef LocalOrdinal LO;
1553 typedef typename Graph::local_graph_type::row_map_type row_map_type;
1554 typedef typename row_map_type::non_const_type non_const_row_map_type;
1555 typedef typename local_matrix_type::values_type values_type;
1556 #ifdef HAVE_TPETRA_DEBUG
1557 const char tfecfFuncName[] =
"fillLocalMatrix (called from fillComplete): ";
1558 #endif // HAVE_TPETRA_DEBUG
1559 ProfilingRegion regionFLM (
"Tpetra::CrsMatrix::fillLocalMatrix");
1561 const size_t lclNumRows = getNodeNumRows();
1562 const map_type& rowMap = * (getRowMap ());
1563 RCP<node_type> node = rowMap.
getNode ();
1574 ArrayRCP<Array<LO> > lclInds2D = staticGraph_->lclInds2D_;
1575 size_t nodeNumEntries = staticGraph_->getNodeNumEntries ();
1576 size_t nodeNumAllocated = staticGraph_->getNodeAllocationSize ();
1577 row_map_type k_rowPtrs_ = staticGraph_->lclGraph_.row_map;
1579 row_map_type k_ptrs;
1585 bool requestOptimizedStorage =
true;
1586 const bool default_OptimizeStorage =
1587 ! isStaticGraph () || staticGraph_->isStorageOptimized ();
1588 if (! params.is_null () && ! params->get (
"Optimize Storage", default_OptimizeStorage)) {
1589 requestOptimizedStorage =
false;
1596 if (! staticGraph_->isStorageOptimized () && requestOptimizedStorage) {
1598 "You requested optimized storage by setting the"
1599 "\"Optimize Storage\" flag to \"true\" in the parameter list, or by virtue"
1600 "of default behavior. However, the associated CrsGraph was filled separately"
1601 "and requested not to optimize storage. Therefore, the CrsMatrix cannot"
1602 "optimize storage.");
1603 requestOptimizedStorage =
false;
1606 typedef decltype (staticGraph_->k_numRowEntries_) row_entries_type;
1631 size_t lclTotalNumEntries = 0;
1633 typename non_const_row_map_type::HostMirror h_ptrs;
1635 typename row_entries_type::const_type numRowEnt_h =
1636 staticGraph_->k_numRowEntries_;
1638 non_const_row_map_type packedRowOffsets (
"Tpetra::CrsGraph::ptr",
1642 h_ptrs = create_mirror_view (packedRowOffsets);
1646 k_ptrs = packedRowOffsets;
1649 #ifdef HAVE_TPETRA_DEBUG
1650 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1651 (static_cast<size_t> (k_ptrs.extent (0)) != lclNumRows + 1,
1652 std::logic_error,
"In DynamicProfile branch, after packing k_ptrs, "
1653 "k_ptrs.extent(0) = " << k_ptrs.extent (0) <<
" != "
1654 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
1655 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1656 (static_cast<size_t> (h_ptrs.extent (0)) != lclNumRows + 1,
1657 std::logic_error,
"In DynamicProfile branch, after packing h_ptrs, "
1658 "h_ptrs.extent(0) = " << h_ptrs.extent (0) <<
" != "
1659 "(lclNumRows+1) = " << (lclNumRows+1) <<
".");
1661 const auto valToCheck = ::Tpetra::Details::getEntryOnHost (k_ptrs, lclNumRows);
1662 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1663 (static_cast<size_t> (valToCheck) != lclTotalNumEntries,
1664 std::logic_error,
"(DynamicProfile branch) After packing k_ptrs, "
1665 "k_ptrs(lclNumRows = " << lclNumRows <<
") = " << valToCheck
1666 <<
" != total number of entries on the calling process = "
1667 << lclTotalNumEntries <<
".");
1669 #endif // HAVE_TPETRA_DEBUG
1672 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1674 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
1676 for (
size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) {
1677 const size_t numEnt = numRowEnt_h(lclRow);
1678 std::copy (values2D_[lclRow].begin(),
1679 values2D_[lclRow].begin() + numEnt,
1680 h_vals.data() + h_ptrs(lclRow));
1685 #ifdef HAVE_TPETRA_DEBUG
1687 if (k_ptrs.extent (0) != 0) {
1688 const size_t numOffsets = static_cast<size_t> (k_ptrs.extent (0));
1689 const auto valToCheck =
1690 ::Tpetra::Details::getEntryOnHost (k_ptrs, numOffsets - 1);
1691 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1692 (static_cast<size_t> (valToCheck) != k_vals.extent (0),
1693 std::logic_error,
"(DynamicProfile branch) After packing, k_ptrs("
1694 << (numOffsets-1) <<
") = " << valToCheck <<
" != "
1695 "k_vals.extent(0) = " << k_vals.extent (0) <<
".");
1697 #endif // HAVE_TPETRA_DEBUG
1717 if (nodeNumEntries != nodeNumAllocated) {
1720 non_const_row_map_type tmpk_ptrs (
"Tpetra::CrsGraph::ptr",
1725 size_t lclTotalNumEntries = 0;
1728 typename row_entries_type::const_type numRowEnt_d =
1729 staticGraph_->k_numRowEntries_;
1737 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1740 typedef pack_functor<values_type, row_map_type> packer_type;
1741 packer_type valsPacker (k_vals, k_values1D_, tmpk_ptrs, k_rowPtrs_);
1744 typedef Kokkos::RangePolicy<exec_space, LocalOrdinal> range_type;
1745 Kokkos::parallel_for (range_type (0, lclNumRows), valsPacker);
1748 k_vals = k_values1D_;
1753 if (requestOptimizedStorage) {
1757 k_values1D_ = k_vals;
1758 this->storageStatus_ = ::Tpetra::Details::STORAGE_1D_PACKED;
1766 getColMap ()->getNodeNumElements (),
1768 staticGraph_->getLocalGraph ());
1771 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1776 const typename crs_graph_type::SLocalGlobalViews& newInds,
1777 const Teuchos::ArrayView<impl_scalar_type>& oldRowVals,
1778 const Teuchos::ArrayView<const impl_scalar_type>& newRowVals,
1779 const ELocalGlobal lg,
1780 const ELocalGlobal I)
1782 const size_t oldNumEnt = rowInfo.numEntries;
1783 const size_t numInserted = graph.insertIndices (rowInfo, newInds, lg, I);
1789 if (numInserted > 0) {
1790 const size_t startOffset = oldNumEnt;
1791 memcpy (&oldRowVals[startOffset], &newRowVals[0],
1792 numInserted *
sizeof (impl_scalar_type));
1796 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1798 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1799 insertLocalValues (
const LocalOrdinal lclRow,
1800 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1801 const Teuchos::ArrayView<const Scalar>& values)
1805 const char tfecfFuncName[] =
"insertLocalValues: ";
1807 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1808 (! this->isFillActive (), std::runtime_error,
1809 "Fill is not active. After calling fillComplete, you must call "
1810 "resumeFill before you may insert entries into the matrix again.");
1811 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1812 (this->isStaticGraph (), std::runtime_error,
1813 "Cannot insert indices with static graph; use replaceLocalValues() "
1817 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1818 (graph.
colMap_.is_null (), std::runtime_error,
1819 "Cannot insert local indices without a column map.");
1820 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1822 std::runtime_error,
"Graph indices are global; use "
1823 "insertGlobalValues().");
1824 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1825 (values.size () != indices.size (), std::runtime_error,
1826 "values.size() = " << values.size ()
1827 <<
" != indices.size() = " << indices.size () <<
".");
1828 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1829 ! graph.
rowMap_->isNodeLocalElement (lclRow), std::runtime_error,
1830 "Local row index " << lclRow <<
" does not belong to this process.");
1832 if (! graph.indicesAreAllocated ()) {
1833 this->allocateValues (LocalIndices, GraphNotYetAllocated);
1836 const size_t numEntriesToAdd = static_cast<size_t> (indices.size ());
1837 #ifdef HAVE_TPETRA_DEBUG
1842 using Teuchos::toString;
1845 Teuchos::Array<LocalOrdinal> badColInds;
1846 bool allInColMap =
true;
1847 for (
size_t k = 0; k < numEntriesToAdd; ++k) {
1849 allInColMap =
false;
1850 badColInds.push_back (indices[k]);
1853 if (! allInColMap) {
1854 std::ostringstream os;
1855 os <<
"You attempted to insert entries in owned row " << lclRow
1856 <<
", at the following column indices: " << toString (indices)
1858 os <<
"Of those, the following indices are not in the column Map on "
1859 "this process: " << toString (badColInds) <<
"." << endl <<
"Since "
1860 "the matrix has a column Map already, it is invalid to insert "
1861 "entries at those locations.";
1862 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1863 (
true, std::invalid_argument, os.str ());
1866 #endif // HAVE_TPETRA_DEBUG
1869 const size_t curNumEnt = rowInfo.numEntries;
1870 const size_t newNumEnt = curNumEnt + numEntriesToAdd;
1871 if (newNumEnt > rowInfo.allocSize) {
1872 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1873 (this->getProfileType () ==
StaticProfile, std::runtime_error,
1874 "New indices exceed statically allocated graph structure.");
1876 Teuchos::Array<IST>& curVals = this->values2D_[lclRow];
1880 graph.
lclInds2D_[rowInfo.localRow].resize (newNumEnt);
1881 curVals.resize (newNumEnt);
1882 rowInfo.allocSize = newNumEnt;
1884 typename crs_graph_type::SLocalGlobalViews indsView;
1885 indsView.linds = indices;
1887 Teuchos::ArrayView<IST> valsView = this->getViewNonConst (rowInfo);
1888 Teuchos::ArrayView<const IST> valsIn =
1889 Teuchos::av_reinterpret_cast<const IST> (values);
1890 this->insertIndicesAndValues (graph, rowInfo, indsView, valsView,
1891 valsIn, LocalIndices, LocalIndices);
1892 #ifdef HAVE_TPETRA_DEBUG
1894 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1895 (chkNewNumEnt != newNumEnt, std::logic_error,
1896 "The row should have " << newNumEnt <<
" entries after insert, but "
1897 "instead has " << chkNewNumEnt <<
". Please report this bug to "
1898 "the Tpetra developers.");
1899 #endif // HAVE_TPETRA_DEBUG
1902 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1906 const LocalOrdinal numEnt,
1907 const Scalar vals[],
1908 const LocalOrdinal cols[])
1910 Teuchos::ArrayView<const LocalOrdinal> colsT (cols, numEnt);
1911 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
1912 this->insertLocalValues (localRow, colsT, valsT);
1915 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1920 const GlobalOrdinal gblColInds[],
1921 const impl_scalar_type vals[],
1922 const size_t numInputEnt)
1924 typedef impl_scalar_type IST;
1925 typedef GlobalOrdinal GO;
1926 const char tfecfFuncName[] =
"insertGlobalValuesImpl: ";
1928 #ifdef HAVE_TPETRA_DEBUG
1930 #endif // HAVE_TPETRA_DEBUG
1932 if (! graph.indicesAreAllocated ()) {
1933 this->allocateValues (GlobalIndices, GraphNotYetAllocated);
1938 rowInfo = graph.getRowInfo (rowInfo.localRow);
1941 const size_t curNumEnt = rowInfo.numEntries;
1942 const size_t newNumEnt = curNumEnt + numInputEnt;
1943 if (newNumEnt > rowInfo.allocSize) {
1944 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1945 (this->getProfileType () == StaticProfile &&
1946 newNumEnt > rowInfo.allocSize, std::runtime_error,
1947 "New indices exceed statically allocated graph structure. "
1948 "curNumEnt (" << curNumEnt <<
") + numInputEnt ("
1949 << numInputEnt <<
") > allocSize (" << rowInfo.allocSize
1953 Teuchos::Array<IST>& curVals = this->values2D_[rowInfo.localRow];
1956 graph.gblInds2D_[rowInfo.localRow].resize (newNumEnt);
1957 curVals.resize (newNumEnt);
1958 rowInfo.allocSize = newNumEnt;
1961 using Teuchos::ArrayView;
1962 typename crs_graph_type::SLocalGlobalViews inputIndsAV;
1963 inputIndsAV.ginds = ArrayView<const GO> (gblColInds, numInputEnt);
1964 ArrayView<IST> curValsAV = this->getViewNonConst (rowInfo);
1965 ArrayView<const IST> inputValsAV (vals, numInputEnt);
1967 const ELocalGlobal curIndexingStatus =
1968 this->isGloballyIndexed () ? GlobalIndices : LocalIndices;
1974 this->insertIndicesAndValues (graph, rowInfo, inputIndsAV, curValsAV,
1975 inputValsAV, GlobalIndices,
1977 #ifdef HAVE_TPETRA_DEBUG
1978 const size_t chkNewNumEnt =
1979 graph.getNumEntriesInLocalRow (rowInfo.localRow);
1980 if (chkNewNumEnt != newNumEnt) {
1981 std::ostringstream os;
1982 os << std::endl <<
"newNumEnt = " << newNumEnt
1983 <<
" != graph.getNumEntriesInLocalRow(" << rowInfo.localRow
1984 <<
") = " << chkNewNumEnt <<
"." << std::endl
1985 <<
"\torigNumEnt: " << origNumEnt << std::endl
1986 <<
"\tnumInputEnt: " << numInputEnt << std::endl
1987 <<
"\tgblColInds: [";
1988 for (
size_t k = 0; k < numInputEnt; ++k) {
1989 os << gblColInds[k];
1990 if (k +
size_t (1) < numInputEnt) {
1994 os <<
"]" << std::endl
1996 for (
size_t k = 0; k < numInputEnt; ++k) {
1998 if (k +
size_t (1) < numInputEnt) {
2002 os <<
"]" << std::endl;
2004 if (this->supportsRowViews ()) {
2005 Teuchos::ArrayView<const Scalar> vals2;
2006 if (this->isGloballyIndexed ()) {
2007 Teuchos::ArrayView<const GlobalOrdinal> gblColInds2;
2008 const GlobalOrdinal gblRow =
2009 graph.rowMap_->getGlobalElement (rowInfo.localRow);
2010 if (gblRow == Tpetra::Details::OrdinalTraits<GlobalOrdinal>::invalid ()) {
2011 os <<
"Local row index " << rowInfo.localRow <<
" is invalid!" << std::endl;
2014 bool getViewThrew =
false;
2016 this->getGlobalRowView (gblRow, gblColInds2, vals2);
2018 catch (std::exception& e) {
2019 getViewThrew =
true;
2020 os <<
"getGlobalRowView threw exception:" << std::endl
2021 << e.what () << std::endl;
2023 if (! getViewThrew) {
2024 os <<
"\tNew global column indices: "
2025 << Teuchos::toString (gblColInds2) << std::endl
2026 <<
"\tNew values: " << Teuchos::toString (vals2) << std::endl;
2030 else if (this->isLocallyIndexed ()) {
2031 Teuchos::ArrayView<const LocalOrdinal> lclColInds2;
2032 this->getLocalRowView (rowInfo.localRow, lclColInds2, vals2);
2033 os <<
"\tNew local column indices: " << Teuchos::toString (lclColInds2)
2035 os <<
"\tNew values: " << Teuchos::toString (vals2) << std::endl;
2039 os <<
"Please report this bug to the Tpetra developers.";
2040 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2041 (
true, std::logic_error, os.str ());
2043 #endif // HAVE_TPETRA_DEBUG
2046 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2048 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2049 insertGlobalValues (
const GlobalOrdinal gblRow,
2050 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2051 const Teuchos::ArrayView<const Scalar>& values)
2053 using Teuchos::toString;
2056 typedef LocalOrdinal LO;
2057 typedef GlobalOrdinal GO;
2058 typedef Tpetra::Details::OrdinalTraits<LO> OTLO;
2059 typedef typename Teuchos::ArrayView<const GO>::size_type size_type;
2060 const char tfecfFuncName[] =
"insertGlobalValues: ";
2062 #ifdef HAVE_TPETRA_DEBUG
2063 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2064 (values.size () != indices.size (), std::runtime_error,
2065 "values.size() = " << values.size () <<
" != indices.size() = "
2066 << indices.size () <<
".");
2067 #endif // HAVE_TPETRA_DEBUG
2071 const map_type& rowMap = * (this->getCrsGraphRef ().rowMap_);
2074 if (lclRow == OTLO::invalid ()) {
2081 this->insertNonownedGlobalValues (gblRow, indices, values);
2084 if (this->isStaticGraph ()) {
2086 const int myRank = rowMap.getComm ()->getRank ();
2087 const int numProcs = rowMap.getComm ()->getSize ();
2088 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2089 (
true, std::runtime_error,
2090 "The matrix was constructed with a constant (\"static\") graph, "
2091 "yet the given global row index " << gblRow <<
" is in the row "
2092 "Map on the calling process (with rank " << myRank <<
", of " <<
2093 numProcs <<
" process(es)). In this case, you may not insert "
2094 "new entries into rows owned by the calling process.");
2098 const IST*
const inputVals =
2099 reinterpret_cast<const IST*> (values.getRawPtr ());
2100 const GO*
const inputGblColInds = indices.getRawPtr ();
2101 const size_t numInputEnt = indices.size ();
2110 if (! graph.
colMap_.is_null ()) {
2116 #ifdef HAVE_TPETRA_DEBUG
2117 Teuchos::Array<GO> badColInds;
2118 #endif // HAVE_TPETRA_DEBUG
2119 const size_type numEntriesToInsert = indices.size ();
2120 bool allInColMap =
true;
2121 for (size_type k = 0; k < numEntriesToInsert; ++k) {
2123 allInColMap =
false;
2124 #ifdef HAVE_TPETRA_DEBUG
2125 badColInds.push_back (indices[k]);
2128 #endif // HAVE_TPETRA_DEBUG
2131 if (! allInColMap) {
2132 std::ostringstream os;
2133 os <<
"You attempted to insert entries in owned row " << gblRow
2134 <<
", at the following column indices: " << toString (indices)
2136 #ifdef HAVE_TPETRA_DEBUG
2137 os <<
"Of those, the following indices are not in the column Map "
2138 "on this process: " << toString (badColInds) <<
"." << endl
2139 <<
"Since the matrix has a column Map already, it is invalid "
2140 "to insert entries at those locations.";
2142 os <<
"At least one of those indices is not in the column Map "
2143 "on this process." << endl <<
"It is invalid to insert into "
2144 "columns not in the column Map on the process that owns the "
2146 #endif // HAVE_TPETRA_DEBUG
2147 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2148 (
true, std::invalid_argument, os.str ());
2152 this->insertGlobalValuesImpl (graph, rowInfo, inputGblColInds,
2153 inputVals, numInputEnt);
2158 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2162 const LocalOrdinal numEnt,
2163 const Scalar vals[],
2164 const GlobalOrdinal inds[])
2166 Teuchos::ArrayView<const GlobalOrdinal> indsT (inds, numEnt);
2167 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
2168 this->insertGlobalValues (globalRow, indsT, valsT);
2172 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2176 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2177 const Teuchos::ArrayView<const Scalar>& values)
2179 typedef impl_scalar_type IST;
2180 typedef LocalOrdinal LO;
2181 typedef GlobalOrdinal GO;
2182 typedef Tpetra::Details::OrdinalTraits<LO> OTLO;
2183 const char tfecfFuncName[] =
"insertGlobalValuesFiltered: ";
2185 #ifdef HAVE_TPETRA_DEBUG
2186 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2187 (values.size () != indices.size (), std::runtime_error,
2188 "values.size() = " << values.size () <<
" != indices.size() = "
2189 << indices.size () <<
".");
2190 #endif // HAVE_TPETRA_DEBUG
2194 const map_type& rowMap = * (this->getCrsGraphRef ().rowMap_);
2195 const LO lclRow = rowMap.getLocalElement (gblRow);
2197 if (lclRow == OTLO::invalid ()) {
2204 this->insertNonownedGlobalValues (gblRow, indices, values);
2207 if (this->isStaticGraph ()) {
2209 const int myRank = rowMap.getComm ()->getRank ();
2210 const int numProcs = rowMap.getComm ()->getSize ();
2211 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2212 (
true, std::runtime_error,
2213 "The matrix was constructed with a constant (\"static\") graph, "
2214 "yet the given global row index " << gblRow <<
" is in the row "
2215 "Map on the calling process (with rank " << myRank <<
", of " <<
2216 numProcs <<
" process(es)). In this case, you may not insert "
2217 "new entries into rows owned by the calling process.");
2220 crs_graph_type& graph = * (this->myGraph_);
2221 const IST*
const inputVals =
2222 reinterpret_cast<const IST*> (values.getRawPtr ());
2223 const GO*
const inputGblColInds = indices.getRawPtr ();
2224 const size_t numInputEnt = indices.size ();
2225 RowInfo rowInfo = graph.getRowInfo (lclRow);
2227 if (! graph.colMap_.is_null ()) {
2228 const map_type& colMap = * (graph.colMap_);
2229 size_t curOffset = 0;
2230 while (curOffset < numInputEnt) {
2234 size_t endOffset = curOffset;
2235 for ( ; endOffset < numInputEnt &&
2236 colMap.getLocalElement (inputGblColInds[endOffset]) != OTLO::invalid ();
2242 const LO numIndInSeq = (endOffset - curOffset);
2243 if (numIndInSeq != 0) {
2244 this->insertGlobalValuesImpl (graph, rowInfo,
2245 inputGblColInds + curOffset,
2246 inputVals + curOffset,
2252 #ifdef HAVE_TPETRA_DEBUG
2253 const bool invariant = endOffset == numInputEnt ||
2254 colMap.getLocalElement (inputGblColInds[endOffset]) == OTLO::invalid ();
2255 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2256 (! invariant, std::logic_error, std::endl <<
"Invariant failed!");
2257 #endif // HAVE_TPETRA_DEBUG
2258 curOffset = endOffset + 1;
2262 this->insertGlobalValuesImpl (graph, rowInfo, inputGblColInds,
2263 inputVals, numInputEnt);
2268 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2270 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2271 replaceLocalValuesImpl (impl_scalar_type rowVals[],
2272 const crs_graph_type& graph,
2273 const RowInfo& rowInfo,
2274 const LocalOrdinal inds[],
2275 const impl_scalar_type newVals[],
2276 const LocalOrdinal numElts)
const
2278 typedef LocalOrdinal LO;
2279 typedef GlobalOrdinal GO;
2281 const bool sorted = graph.isSorted ();
2290 if (graph.isLocallyIndexed ()) {
2293 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2295 for (LO j = 0; j < numElts; ++j) {
2296 const LO lclColInd = inds[j];
2297 const size_t offset =
2298 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2299 lclColInd, hint, sorted);
2300 if (offset != rowInfo.numEntries) {
2301 rowVals[offset] = newVals[j];
2307 else if (graph.isGloballyIndexed ()) {
2308 if (graph.colMap_.is_null ()) {
2309 return Teuchos::OrdinalTraits<LO>::invalid ();
2311 const map_type colMap = * (graph.colMap_);
2315 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
2317 for (LO j = 0; j < numElts; ++j) {
2318 const GO gblColInd = colMap.getGlobalElement (inds[j]);
2319 if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
2320 const size_t offset =
2321 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2322 gblColInd, hint, sorted);
2323 if (offset != rowInfo.numEntries) {
2324 rowVals[offset] = newVals[j];
2343 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2345 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2346 replaceLocalValues (
const LocalOrdinal localRow,
2347 const Teuchos::ArrayView<const LocalOrdinal>& lclCols,
2348 const Teuchos::ArrayView<const Scalar>& vals)
const
2350 typedef LocalOrdinal LO;
2352 const LO numInputEnt = static_cast<LO> (lclCols.size ());
2353 if (static_cast<LO> (vals.size ()) != numInputEnt) {
2354 return Teuchos::OrdinalTraits<LO>::invalid ();
2356 const LO*
const inputInds = lclCols.getRawPtr ();
2357 const Scalar*
const inputVals = vals.getRawPtr ();
2358 return this->replaceLocalValues (localRow, numInputEnt,
2359 inputVals, inputInds);
2362 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2366 const LocalOrdinal numEnt,
2367 const Scalar inputVals[],
2368 const LocalOrdinal inputCols[])
const
2371 typedef LocalOrdinal LO;
2373 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2375 return Teuchos::OrdinalTraits<LO>::invalid ();
2380 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2383 return static_cast<LO> (0);
2385 auto curRowVals = this->getRowViewNonConst (rowInfo);
2386 const IST*
const inVals = reinterpret_cast<const IST*> (inputVals);
2387 return this->replaceLocalValuesImpl (curRowVals.data (), graph, rowInfo,
2388 inputCols, inVals, numEnt);
2391 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2395 const crs_graph_type& graph,
2397 const GlobalOrdinal inds[],
2398 const impl_scalar_type newVals[],
2399 const LocalOrdinal numElts)
const
2401 typedef LocalOrdinal LO;
2402 typedef GlobalOrdinal GO;
2404 const bool sorted = graph.isSorted ();
2413 if (graph.isLocallyIndexed ()) {
2417 if (graph.colMap_.is_null ()) {
2423 const map_type& colMap = * (graph.colMap_);
2427 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2428 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
2429 for (LO j = 0; j < numElts; ++j) {
2430 const LO lclColInd = colMap.getLocalElement (inds[j]);
2431 if (lclColInd != LINV) {
2432 const size_t offset =
2433 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2434 lclColInd, hint, sorted);
2435 if (offset != rowInfo.numEntries) {
2436 rowVals[offset] = newVals[j];
2443 else if (graph.isGloballyIndexed ()) {
2446 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
2448 for (LO j = 0; j < numElts; ++j) {
2449 const GO gblColInd = inds[j];
2450 const size_t offset =
2451 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2452 gblColInd, hint, sorted);
2453 if (offset != rowInfo.numEntries) {
2454 rowVals[offset] = newVals[j];
2467 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2469 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2470 replaceGlobalValues (
const GlobalOrdinal globalRow,
2471 const Teuchos::ArrayView<const GlobalOrdinal>& inputGblColInds,
2472 const Teuchos::ArrayView<const Scalar>& inputVals)
const
2474 typedef LocalOrdinal LO;
2476 const LO numInputEnt = static_cast<LO> (inputGblColInds.size ());
2477 if (static_cast<LO> (inputVals.size ()) != numInputEnt) {
2478 return Teuchos::OrdinalTraits<LO>::invalid ();
2480 return this->replaceGlobalValues (globalRow, numInputEnt,
2481 inputVals.getRawPtr (),
2482 inputGblColInds.getRawPtr ());
2485 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2489 const LocalOrdinal numEnt,
2490 const Scalar inputVals[],
2491 const GlobalOrdinal inputGblColInds[])
const
2494 typedef LocalOrdinal LO;
2496 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2498 return Teuchos::OrdinalTraits<LO>::invalid ();
2503 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2506 return static_cast<LO> (0);
2509 auto curRowVals = this->getRowViewNonConst (rowInfo);
2510 const IST*
const inVals = reinterpret_cast<const IST*> (inputVals);
2511 return this->replaceGlobalValuesImpl (curRowVals.data (), graph, rowInfo,
2512 inputGblColInds, inVals, numEnt);
2515 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2519 const crs_graph_type& graph,
2521 const GlobalOrdinal inds[],
2522 const impl_scalar_type newVals[],
2523 const LocalOrdinal numElts,
2524 const bool atomic)
const
2526 typedef LocalOrdinal LO;
2527 typedef GlobalOrdinal GO;
2529 const bool sorted = graph.isSorted ();
2538 if (graph.isLocallyIndexed ()) {
2542 if (graph.colMap_.is_null ()) {
2548 const map_type& colMap = * (graph.colMap_);
2552 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2553 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
2555 for (LO j = 0; j < numElts; ++j) {
2556 const LO lclColInd = colMap.getLocalElement (inds[j]);
2557 if (lclColInd != LINV) {
2558 const size_t offset =
2559 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2560 lclColInd, hint, sorted);
2561 if (offset != rowInfo.numEntries) {
2563 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2566 rowVals[offset] += newVals[j];
2574 else if (graph.isGloballyIndexed ()) {
2577 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
2579 for (LO j = 0; j < numElts; ++j) {
2580 const GO gblColInd = inds[j];
2581 const size_t offset =
2582 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2583 gblColInd, hint, sorted);
2584 if (offset != rowInfo.numEntries) {
2586 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2589 rowVals[offset] += newVals[j];
2603 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2605 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2606 sumIntoGlobalValues (
const GlobalOrdinal gblRow,
2607 const Teuchos::ArrayView<const GlobalOrdinal>& inputGblColInds,
2608 const Teuchos::ArrayView<const Scalar>& inputVals,
2611 typedef LocalOrdinal LO;
2613 const LO numInputEnt = static_cast<LO> (inputGblColInds.size ());
2614 if (static_cast<LO> (inputVals.size ()) != numInputEnt) {
2615 return Teuchos::OrdinalTraits<LO>::invalid ();
2617 return this->sumIntoGlobalValues (gblRow, numInputEnt,
2618 inputVals.getRawPtr (),
2619 inputGblColInds.getRawPtr (),
2623 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2627 const LocalOrdinal numInputEnt,
2628 const Scalar inputVals[],
2629 const GlobalOrdinal inputGblColInds[],
2633 typedef LocalOrdinal LO;
2634 typedef GlobalOrdinal GO;
2636 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2638 return Teuchos::OrdinalTraits<LO>::invalid ();
2643 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2648 using Teuchos::ArrayView;
2649 ArrayView<const GO> inputGblColInds_av (numInputEnt == 0 ? NULL :
2650 inputGblColInds, numInputEnt);
2651 ArrayView<const Scalar> inputVals_av (numInputEnt == 0 ? NULL :
2652 inputVals, numInputEnt);
2657 this->insertNonownedGlobalValues (gblRow, inputGblColInds_av,
2668 auto curRowVals = this->getRowViewNonConst (rowInfo);
2669 const IST*
const inVals = reinterpret_cast<const IST*> (inputVals);
2670 return this->sumIntoGlobalValuesImpl (curRowVals.data (), graph, rowInfo,
2671 inputGblColInds, inVals,
2672 numInputEnt, atomic);
2676 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2680 const LocalOrdinal numInputEnt,
2681 const impl_scalar_type inputVals[],
2682 const LocalOrdinal inputCols[],
2683 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2684 const bool atomic)
const
2686 using Tpetra::Details::OrdinalTraits;
2687 typedef LocalOrdinal LO;
2689 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2691 return Teuchos::OrdinalTraits<LO>::invalid ();
2693 const crs_graph_type& graph = * (this->staticGraph_);
2694 const RowInfo rowInfo = graph.getRowInfo (lclRow);
2696 if (rowInfo.localRow == OrdinalTraits<size_t>::invalid ()) {
2699 return static_cast<LO> (0);
2701 auto curRowVals = this->getRowViewNonConst (rowInfo);
2702 return this->transformLocalValues (curRowVals.data (), graph,
2703 rowInfo, inputCols, inputVals,
2704 numInputEnt, f, atomic);
2707 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2709 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2710 transformGlobalValues (
const GlobalOrdinal gblRow,
2711 const LocalOrdinal numInputEnt,
2712 const impl_scalar_type inputVals[],
2713 const GlobalOrdinal inputCols[],
2714 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2715 const bool atomic)
const
2717 using Tpetra::Details::OrdinalTraits;
2718 typedef LocalOrdinal LO;
2720 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2722 return OrdinalTraits<LO>::invalid ();
2724 const crs_graph_type& graph = * (this->staticGraph_);
2725 const RowInfo rowInfo = graph.getRowInfoFromGlobalRowIndex (gblRow);
2727 if (rowInfo.localRow == OrdinalTraits<size_t>::invalid ()) {
2730 return static_cast<LO> (0);
2732 auto curRowVals = this->getRowViewNonConst (rowInfo);
2733 return this->transformGlobalValues (curRowVals.data (), graph,
2734 rowInfo, inputCols, inputVals,
2735 numInputEnt, f, atomic);
2738 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2740 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2741 transformLocalValues (impl_scalar_type rowVals[],
2742 const crs_graph_type& graph,
2743 const RowInfo& rowInfo,
2744 const LocalOrdinal inds[],
2745 const impl_scalar_type newVals[],
2746 const LocalOrdinal numElts,
2747 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2748 const bool atomic)
const
2750 typedef impl_scalar_type ST;
2751 typedef LocalOrdinal LO;
2752 typedef GlobalOrdinal GO;
2759 const bool sorted = graph.isSorted ();
2764 if (graph.isLocallyIndexed ()) {
2767 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2769 for (LO j = 0; j < numElts; ++j) {
2770 const LO lclColInd = inds[j];
2771 const size_t offset =
2772 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2773 lclColInd, hint, sorted);
2774 if (offset != rowInfo.numEntries) {
2783 volatile ST*
const dest = &rowVals[offset];
2784 (void) atomic_binary_function_update (dest, newVals[j], f);
2788 rowVals[offset] = f (rowVals[offset], newVals[j]);
2795 else if (graph.isGloballyIndexed ()) {
2799 if (graph.colMap_.is_null ()) {
2806 const map_type& colMap = * (graph.colMap_);
2809 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
2811 const GO GINV = Teuchos::OrdinalTraits<GO>::invalid ();
2812 for (LO j = 0; j < numElts; ++j) {
2813 const GO gblColInd = colMap.getGlobalElement (inds[j]);
2814 if (gblColInd != GINV) {
2815 const size_t offset =
2816 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2817 gblColInd, hint, sorted);
2818 if (offset != rowInfo.numEntries) {
2827 volatile ST*
const dest = &rowVals[offset];
2828 (void) atomic_binary_function_update (dest, newVals[j], f);
2832 rowVals[offset] = f (rowVals[offset], newVals[j]);
2847 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2849 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2850 transformGlobalValues (impl_scalar_type rowVals[],
2851 const crs_graph_type& graph,
2852 const RowInfo& rowInfo,
2853 const GlobalOrdinal inds[],
2854 const impl_scalar_type newVals[],
2855 const LocalOrdinal numElts,
2856 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2857 const bool atomic)
const
2859 typedef impl_scalar_type ST;
2860 typedef LocalOrdinal LO;
2861 typedef GlobalOrdinal GO;
2868 const bool sorted = graph.isSorted ();
2873 if (graph.isGloballyIndexed ()) {
2876 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
2878 for (LO j = 0; j < numElts; ++j) {
2879 const GO gblColInd = inds[j];
2880 const size_t offset =
2881 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2882 gblColInd, hint, sorted);
2883 if (offset != rowInfo.numEntries) {
2892 volatile ST*
const dest = &rowVals[offset];
2893 (void) atomic_binary_function_update (dest, newVals[j], f);
2897 rowVals[offset] = f (rowVals[offset], newVals[j]);
2904 else if (graph.isLocallyIndexed ()) {
2908 if (graph.colMap_.is_null ()) {
2914 const map_type& colMap = * (graph.colMap_);
2917 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2919 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
2920 for (LO j = 0; j < numElts; ++j) {
2921 const LO lclColInd = colMap.getLocalElement (inds[j]);
2922 if (lclColInd != LINV) {
2923 const size_t offset =
2924 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2925 lclColInd, hint, sorted);
2926 if (offset != rowInfo.numEntries) {
2935 volatile ST*
const dest = &rowVals[offset];
2936 (void) atomic_binary_function_update (dest, newVals[j], f);
2940 rowVals[offset] = f (rowVals[offset], newVals[j]);
2955 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2957 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2958 sumIntoLocalValuesImpl (impl_scalar_type rowVals[],
2959 const crs_graph_type& graph,
2960 const RowInfo& rowInfo,
2961 const LocalOrdinal inds[],
2962 const impl_scalar_type newVals[],
2963 const LocalOrdinal numElts,
2964 const bool atomic)
const
2966 typedef LocalOrdinal LO;
2967 typedef GlobalOrdinal GO;
2969 const bool sorted = graph.isSorted ();
2978 if (graph.isLocallyIndexed ()) {
2981 auto colInds = graph.getLocalKokkosRowView (rowInfo);
2983 for (LO j = 0; j < numElts; ++j) {
2984 const LO lclColInd = inds[j];
2985 const size_t offset =
2986 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2987 lclColInd, hint, sorted);
2988 if (offset != rowInfo.numEntries) {
2990 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2993 rowVals[offset] += newVals[j];
3000 else if (graph.isGloballyIndexed ()) {
3001 if (graph.colMap_.is_null ()) {
3002 return Teuchos::OrdinalTraits<LO>::invalid ();
3004 const map_type colMap = * (graph.colMap_);
3008 auto colInds = graph.getGlobalKokkosRowView (rowInfo);
3010 for (LO j = 0; j < numElts; ++j) {
3011 const GO gblColInd = colMap.getGlobalElement (inds[j]);
3012 if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
3013 const size_t offset =
3014 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
3015 gblColInd, hint, sorted);
3016 if (offset != rowInfo.numEntries) {
3018 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
3021 rowVals[offset] += newVals[j];
3041 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3043 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
3044 sumIntoLocalValues (
const LocalOrdinal localRow,
3045 const Teuchos::ArrayView<const LocalOrdinal>& indices,
3046 const Teuchos::ArrayView<const Scalar>& values,
3047 const bool atomic)
const
3049 typedef LocalOrdinal LO;
3051 const LO numInputEnt = static_cast<LO> (indices.size ());
3052 if (static_cast<LO> (values.size ()) != numInputEnt) {
3053 return Teuchos::OrdinalTraits<LO>::invalid ();
3055 const LO*
const inputInds = indices.getRawPtr ();
3056 const Scalar*
const inputVals = values.getRawPtr ();
3057 return this->sumIntoLocalValues (localRow, numInputEnt,
3058 inputVals, inputInds, atomic);
3061 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3065 const LocalOrdinal numEnt,
3066 const Scalar vals[],
3067 const LocalOrdinal cols[],
3068 const bool atomic)
const
3071 typedef LocalOrdinal LO;
3073 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
3075 return Teuchos::OrdinalTraits<LO>::invalid ();
3080 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
3083 return static_cast<LO> (0);
3085 auto curRowVals = this->getRowViewNonConst (rowInfo);
3086 const IST*
const inputVals = reinterpret_cast<const IST*> (vals);
3087 return this->sumIntoLocalValuesImpl (curRowVals.data (), graph, rowInfo,
3088 cols, inputVals, numEnt, atomic);
3091 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3092 Teuchos::ArrayView<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::impl_scalar_type>
3096 using Kokkos::MemoryUnmanaged;
3098 using Teuchos::ArrayView;
3100 typedef std::pair<size_t, size_t> range_type;
3102 if (k_values1D_.extent (0) != 0 && rowinfo.allocSize > 0) {
3103 #ifdef HAVE_TPETRA_DEBUG
3104 TEUCHOS_TEST_FOR_EXCEPTION(
3105 rowinfo.offset1D + rowinfo.allocSize > k_values1D_.extent (0),
3106 std::range_error,
"Tpetra::CrsMatrix::getView: Invalid access "
3107 "to 1-D storage of values." << std::endl <<
"rowinfo.offset1D (" <<
3108 rowinfo.offset1D <<
") + rowinfo.allocSize (" << rowinfo.allocSize <<
3109 ") > k_values1D_.extent(0) (" << k_values1D_.extent (0) <<
").");
3110 #endif // HAVE_TPETRA_DEBUG
3111 range_type range (rowinfo.offset1D, rowinfo.offset1D + rowinfo.allocSize);
3112 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
3119 subview_type sv = Kokkos::subview (subview_type (k_values1D_), range);
3120 const ST*
const sv_raw = (rowinfo.allocSize == 0) ? NULL : sv.data ();
3121 return ArrayView<const ST> (sv_raw, rowinfo.allocSize);
3123 else if (values2D_ != Teuchos::null) {
3124 return values2D_[rowinfo.localRow] ();
3127 return ArrayView<impl_scalar_type> ();
3132 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3136 LocalOrdinal& numEnt,
3139 if (k_values1D_.extent (0) != 0 && rowinfo.allocSize > 0) {
3140 #ifdef HAVE_TPETRA_DEBUG
3141 if (rowinfo.offset1D + rowinfo.allocSize > k_values1D_.extent (0)) {
3144 return Teuchos::OrdinalTraits<LocalOrdinal>::invalid ();
3146 #endif // HAVE_TPETRA_DEBUG
3147 vals = k_values1D_.data () + rowinfo.offset1D;
3148 numEnt = rowinfo.allocSize;
3150 else if (! values2D_.is_null ()) {
3151 #ifdef HAVE_TPETRA_DEBUG
3152 if (rowinfo.localRow >= static_cast<size_t> (values2D_.size ())) {
3155 return Teuchos::OrdinalTraits<LocalOrdinal>::invalid ();
3157 #endif // HAVE_TPETRA_DEBUG
3160 const auto& curRow = values2D_[rowinfo.localRow];
3161 vals = curRow.getRawPtr ();
3162 numEnt = curRow.size ();
3169 return static_cast<LocalOrdinal> (0);
3172 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3176 LocalOrdinal& numEnt,
3180 const LocalOrdinal err = this->getViewRawConst (valsConst, numEnt, rowinfo);
3181 vals = const_cast<impl_scalar_type*> (valsConst);
3185 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3186 Kokkos::View<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::impl_scalar_type*,
3188 Kokkos::MemoryUnmanaged>
3192 using Kokkos::MemoryUnmanaged;
3194 typedef impl_scalar_type ST;
3195 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
3196 typedef std::pair<size_t, size_t> range_type;
3198 if (k_values1D_.extent (0) != 0 && rowInfo.allocSize > 0) {
3199 #ifdef HAVE_TPETRA_DEBUG
3200 TEUCHOS_TEST_FOR_EXCEPTION
3201 (rowInfo.offset1D + rowInfo.allocSize > this->k_values1D_.extent (0),
3202 std::range_error,
"Tpetra::CrsMatrix::getRowView: Invalid access "
3203 "to 1-D storage of values. rowInfo.offset1D ("
3204 << rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize
3205 <<
") > this->k_values1D_.extent(0) ("
3206 << this->k_values1D_.extent (0) <<
").");
3207 #endif // HAVE_TPETRA_DEBUG
3208 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
3215 return Kokkos::subview (subview_type (this->k_values1D_), range);
3217 else if (this->values2D_ != Teuchos::null) {
3221 auto& rowView = this->values2D_[rowInfo.localRow];
3222 return subview_type (rowView.getRawPtr (), rowView.size ());
3225 return subview_type ();
3229 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3230 Kokkos::View<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::impl_scalar_type*,
3231 typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::execution_space,
3232 Kokkos::MemoryUnmanaged>
3233 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
3234 getRowViewNonConst (
const RowInfo& rowInfo)
const
3236 using Kokkos::MemoryUnmanaged;
3238 typedef impl_scalar_type ST;
3239 typedef View<ST*, execution_space, MemoryUnmanaged> subview_type;
3240 typedef std::pair<size_t, size_t> range_type;
3242 if (k_values1D_.extent (0) != 0 && rowInfo.allocSize > 0) {
3243 #ifdef HAVE_TPETRA_DEBUG
3244 TEUCHOS_TEST_FOR_EXCEPTION
3245 (rowInfo.offset1D + rowInfo.allocSize > this->k_values1D_.extent (0),
3246 std::range_error,
"Tpetra::CrsMatrix::getRowViewNonConst: Invalid "
3247 "access to 1-D storage of values. rowInfo.offset1D ("
3248 << rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize
3249 <<
") > this->k_values1D_.extent(0) ("
3250 << this->k_values1D_.extent (0) <<
").");
3251 #endif // HAVE_TPETRA_DEBUG
3252 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
3259 return Kokkos::subview (subview_type (this->k_values1D_), range);
3261 else if (this->values2D_ != Teuchos::null) {
3265 auto& rowView = this->values2D_[rowInfo.localRow];
3266 return subview_type (rowView.getRawPtr (), rowView.size ());
3269 return subview_type ();
3273 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3274 Teuchos::ArrayView<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::impl_scalar_type>
3275 CrsMatrix<Scalar, LocalOrdinal,GlobalOrdinal, Node>::
3278 return Teuchos::av_const_cast<impl_scalar_type> (this->getView (rowinfo));
3281 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3285 const Teuchos::ArrayView<LocalOrdinal>& indices,
3286 const Teuchos::ArrayView<Scalar>& values,
3287 size_t& numEntries)
const
3289 using Teuchos::ArrayView;
3290 using Teuchos::av_reinterpret_cast;
3291 const char tfecfFuncName[] =
"getLocalRowCopy: ";
3293 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3294 (! this->hasColMap (), std::runtime_error,
3295 "The matrix does not have a column Map yet. This means we don't have "
3296 "local indices for columns yet, so it doesn't make sense to call this "
3297 "method. If the matrix doesn't have a column Map yet, you should call "
3298 "fillComplete on it first.");
3300 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
3301 const size_t theNumEntries = rowinfo.numEntries;
3302 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3303 (static_cast<size_t> (indices.size ()) < theNumEntries ||
3304 static_cast<size_t> (values.size ()) < theNumEntries,
3305 std::runtime_error,
"Row with local index " << localRow <<
" has " <<
3306 theNumEntries <<
" entry/ies, but indices.size() = " <<
3307 indices.size () <<
" and values.size() = " << values.size () <<
".");
3308 numEntries = theNumEntries;
3310 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
3311 if (staticGraph_->isLocallyIndexed ()) {
3312 const LocalOrdinal* curLclInds;
3314 LocalOrdinal numSpots;
3319 #ifdef HAVE_TPETRA_DEBUG
3321 staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
3322 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3323 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3324 "staticGraph_->getLocalViewRawConst returned nonzero error code "
3326 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3327 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
3328 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
3330 const LocalOrdinal numSpotsBefore = numSpots;
3331 err = getViewRawConst (curVals, numSpots, rowinfo);
3332 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3333 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3334 "getViewRaw returned nonzero error code " << err <<
".");
3335 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3336 (numSpotsBefore != numSpots, std::logic_error,
3337 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = "
3338 << numSpots <<
".");
3340 (void) staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
3341 (void) getViewRawConst (curVals, numSpots, rowinfo);
3342 #endif // HAVE_TPETRA_DEBUG
3344 for (
size_t j = 0; j < theNumEntries; ++j) {
3345 values[j] = curVals[j];
3346 indices[j] = curLclInds[j];
3349 else if (staticGraph_->isGloballyIndexed ()) {
3351 const map_type& colMap = * (staticGraph_->colMap_);
3352 const GlobalOrdinal* curGblInds;
3354 LocalOrdinal numSpots;
3359 #ifdef HAVE_TPETRA_DEBUG
3361 staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
3362 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3363 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3364 "staticGraph_->getGlobalViewRawConst returned nonzero error code "
3366 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3367 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
3368 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
3370 const LocalOrdinal numSpotsBefore = numSpots;
3371 err = getViewRawConst (curVals, numSpots, rowinfo);
3372 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3373 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3374 "getViewRawConst returned nonzero error code " << err <<
".");
3375 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3376 (numSpotsBefore != numSpots, std::logic_error,
3377 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = "
3378 << numSpots <<
".");
3380 (void) staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
3381 (void) getViewRawConst (curVals, numSpots, rowinfo);
3382 #endif //HAVE_TPETRA_DEBUG
3384 for (
size_t j = 0; j < theNumEntries; ++j) {
3385 values[j] = curVals[j];
3392 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3396 const Teuchos::ArrayView<GlobalOrdinal>& indices,
3397 const Teuchos::ArrayView<Scalar>& values,
3398 size_t& numEntries)
const
3400 using Teuchos::ArrayView;
3401 using Teuchos::av_reinterpret_cast;
3402 const char tfecfFuncName[] =
"getGlobalRowCopy: ";
3405 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
3406 const size_t theNumEntries = rowinfo.numEntries;
3407 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3408 static_cast<size_t> (indices.size ()) < theNumEntries ||
3409 static_cast<size_t> (values.size ()) < theNumEntries,
3410 std::runtime_error,
"Row with global index " << globalRow <<
" has "
3411 << theNumEntries <<
" entry/ies, but indices.size() = " <<
3412 indices.size () <<
" and values.size() = " << values.size () <<
".");
3413 numEntries = theNumEntries;
3415 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
3416 if (staticGraph_->isLocallyIndexed ()) {
3417 const map_type& colMap = * (staticGraph_->colMap_);
3418 const LocalOrdinal* curLclInds;
3420 LocalOrdinal numSpots;
3425 #ifdef HAVE_TPETRA_DEBUG
3427 staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
3428 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3429 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3430 "staticGraph_->getLocalViewRawConst returned nonzero error code "
3432 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3433 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
3434 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
3436 const LocalOrdinal numSpotsBefore = numSpots;
3437 err = getViewRawConst (curVals, numSpots, rowinfo);
3438 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3439 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3440 "getViewRaw returned nonzero error code " << err <<
".");
3441 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3442 (numSpotsBefore != numSpots, std::logic_error,
3443 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = "
3444 << numSpots <<
".");
3446 (void) staticGraph_->getLocalViewRawConst (curLclInds, numSpots, rowinfo);
3447 (void) getViewRawConst (curVals, numSpots, rowinfo);
3448 #endif //HAVE_TPETRA_DEBUG
3450 for (
size_t j = 0; j < theNumEntries; ++j) {
3451 values[j] = curVals[j];
3455 else if (staticGraph_->isGloballyIndexed ()) {
3456 const GlobalOrdinal* curGblInds;
3458 LocalOrdinal numSpots;
3463 #ifdef HAVE_TPETRA_DEBUG
3465 staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
3466 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3467 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3468 "staticGraph_->getGlobalViewRawConst returned nonzero error code "
3470 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3471 (static_cast<size_t> (numSpots) < theNumEntries, std::logic_error,
3472 "numSpots = " << numSpots <<
" < theNumEntries = " << theNumEntries
3474 const LocalOrdinal numSpotsBefore = numSpots;
3475 err = getViewRawConst (curVals, numSpots, rowinfo);
3476 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3477 (err != static_cast<LocalOrdinal> (0), std::logic_error,
3478 "getViewRawConst returned nonzero error code " << err <<
".");
3479 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3480 (numSpotsBefore != numSpots, std::logic_error,
3481 "numSpotsBefore = " << numSpotsBefore <<
" != numSpots = "
3482 << numSpots <<
".");
3484 (void) staticGraph_->getGlobalViewRawConst (curGblInds, numSpots, rowinfo);
3485 (void) getViewRawConst (curVals, numSpots, rowinfo);
3486 #endif //HAVE_TPETRA_DEBUG
3488 for (
size_t j = 0; j < theNumEntries; ++j) {
3489 values[j] = curVals[j];
3490 indices[j] = curGblInds[j];
3496 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3500 Teuchos::ArrayView<const LocalOrdinal>& indices,
3501 Teuchos::ArrayView<const Scalar>& values)
const
3503 using Teuchos::ArrayView;
3504 using Teuchos::av_reinterpret_cast;
3505 typedef LocalOrdinal LO;
3506 const char tfecfFuncName[] =
"getLocalRowView: ";
3508 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3509 isGloballyIndexed (), std::runtime_error,
"The matrix currently stores "
3510 "its indices as global indices, so you cannot get a view with local "
3511 "column indices. If the matrix has a column Map, you may call "
3512 "getLocalRowCopy() to get local column indices; otherwise, you may get "
3513 "a view with global column indices by calling getGlobalRowCopy().");
3514 indices = Teuchos::null;
3515 values = Teuchos::null;
3516 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
3517 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
3518 rowinfo.numEntries > 0) {
3519 ArrayView<const LO> indTmp = staticGraph_->getLocalView (rowinfo);
3520 ArrayView<const Scalar> valTmp =
3521 av_reinterpret_cast<const Scalar> (this->getView (rowinfo));
3522 indices = indTmp (0, rowinfo.numEntries);
3523 values = valTmp (0, rowinfo.numEntries);
3526 #ifdef HAVE_TPETRA_DEBUG
3527 const char suffix[] =
". This should never happen. Please report this "
3528 "bug to the Tpetra developers.";
3529 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3530 (static_cast<size_t> (indices.size ()) !=
3531 static_cast<size_t> (values.size ()), std::logic_error,
3532 "At the end of this method, for local row " << localRow <<
", "
3533 "indices.size() = " << indices.size () <<
" != values.size () = "
3534 << values.size () << suffix);
3535 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3536 (static_cast<size_t> (indices.size ()) !=
3537 static_cast<size_t> (rowinfo.numEntries), std::logic_error,
3538 "At the end of this method, for local row " << localRow <<
", "
3539 "indices.size() = " << indices.size () <<
" != rowinfo.numEntries = "
3540 << rowinfo.numEntries << suffix);
3541 const size_t expectedNumEntries = getNumEntriesInLocalRow (localRow);
3542 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3543 (rowinfo.numEntries != expectedNumEntries, std::logic_error,
"At the end "
3544 "of this method, for local row " << localRow <<
", rowinfo.numEntries = "
3545 << rowinfo.numEntries <<
" != getNumEntriesInLocalRow(localRow) = " <<
3546 expectedNumEntries << suffix);
3547 #endif // HAVE_TPETRA_DEBUG
3550 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3554 LocalOrdinal& numEnt,
3556 const LocalOrdinal*& ind)
const
3558 typedef LocalOrdinal LO;
3565 if (staticGraph_.is_null () || staticGraph_->isGloballyIndexed ()) {
3566 return Tpetra::Details::OrdinalTraits<LO>::invalid ();
3569 const RowInfo rowInfo = staticGraph_->getRowInfo (lclRow);
3570 if (rowInfo.localRow == Tpetra::Details::OrdinalTraits<size_t>::invalid ()) {
3575 return static_cast<LO> (1);
3578 numEnt = static_cast<LO> (rowInfo.numEntries);
3579 auto lclColInds = staticGraph_->getLocalKokkosRowView (rowInfo);
3580 ind = lclColInds.data ();
3581 const LO err = this->getViewRawConst (val, numEnt, rowInfo);
3587 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3591 LocalOrdinal& numEnt,
3592 const LocalOrdinal*& lclColInds,
3593 const Scalar*& vals)
const
3596 const LocalOrdinal errCode =
3597 this->getLocalRowView (lclRow, numEnt, vals_ist, lclColInds);
3598 vals = reinterpret_cast<const Scalar*> (vals_ist);
3602 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3606 Teuchos::ArrayView<const GlobalOrdinal>& indices,
3607 Teuchos::ArrayView<const Scalar>& values)
const
3609 using Teuchos::ArrayView;
3610 using Teuchos::av_reinterpret_cast;
3611 typedef GlobalOrdinal GO;
3612 const char tfecfFuncName[] =
"getGlobalRowView: ";
3614 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3615 isLocallyIndexed (), std::runtime_error,
3616 "The matrix is locally indexed, so we cannot return a view of the row "
3617 "with global column indices. Use getGlobalRowCopy() instead.");
3618 indices = Teuchos::null;
3619 values = Teuchos::null;
3621 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
3622 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
3623 rowinfo.numEntries > 0) {
3624 ArrayView<const GO> indTmp = staticGraph_->getGlobalView (rowinfo);
3625 ArrayView<const Scalar> valTmp =
3626 av_reinterpret_cast<const Scalar> (this->getView (rowinfo));
3627 #ifdef HAVE_TPETRA_DEBUG
3628 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3629 (static_cast<size_t> (indTmp.size ()) < rowinfo.numEntries ||
3630 static_cast<size_t> (valTmp.size ()) < rowinfo.numEntries,
3631 std::logic_error, std::endl <<
"rowinfo.numEntries not accurate. "
3632 << std::endl <<
"indTmp.size() = " << indTmp.size ()
3633 <<
", valTmp.size() = " << valTmp.size ()
3634 <<
", rowinfo.numEntries = " << rowinfo.numEntries <<
".");
3635 #endif // HAVE_TPETRA_DEBUG
3636 indices = indTmp (0, rowinfo.numEntries);
3637 values = valTmp (0, rowinfo.numEntries);
3640 #ifdef HAVE_TPETRA_DEBUG
3641 const char suffix[] =
". This should never happen. Please report this "
3642 "bug to the Tpetra developers.";
3643 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3644 (static_cast<size_t> (indices.size ()) !=
3645 static_cast<size_t> (values.size ()), std::logic_error,
3646 "At the end of this method, for global row " << globalRow <<
", "
3647 "indices.size() = " << indices.size () <<
" != values.size () = "
3648 << values.size () << suffix);
3649 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3650 (static_cast<size_t> (indices.size ()) !=
3651 static_cast<size_t> (rowinfo.numEntries), std::logic_error,
3652 "At the end of this method, for global row " << globalRow <<
", "
3653 "indices.size() = " << indices.size () <<
" != rowinfo.numEntries = "
3654 << rowinfo.numEntries << suffix);
3655 const size_t expectedNumEntries = getNumEntriesInGlobalRow (globalRow);
3656 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3657 (rowinfo.numEntries != expectedNumEntries, std::logic_error,
"At the end "
3658 "of this method, for global row " << globalRow <<
", rowinfo.numEntries "
3659 "= " << rowinfo.numEntries <<
" != getNumEntriesInGlobalRow(globalRow) ="
3660 " " << expectedNumEntries << suffix);
3661 #endif // HAVE_TPETRA_DEBUG
3664 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3669 typedef LocalOrdinal LO;
3670 typedef typename Teuchos::Array<Scalar>::size_type size_type;
3671 const char tfecfFuncName[] =
"scale: ";
3674 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3675 ! isFillActive (), std::runtime_error,
3676 "Fill must be active before you may call this method. "
3677 "Please call resumeFill() to make fill active.");
3679 const size_t nlrs = staticGraph_->getNodeNumRows ();
3680 const size_t numEntries = staticGraph_->getNodeNumEntries ();
3681 if (! staticGraph_->indicesAreAllocated () ||
3682 nlrs == 0 || numEntries == 0) {
3687 const LO lclNumRows = lclMatrix_.numRows ();
3688 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
3689 auto row_i = lclMatrix_.row (lclRow);
3690 for (LO k = 0; k < row_i.length; ++k) {
3692 row_i.value (k) *= theAlpha;
3697 for (
size_t row = 0; row < nlrs; ++row) {
3698 const size_type numEnt = getNumEntriesInLocalRow (row);
3699 Teuchos::ArrayView<impl_scalar_type> rowVals = values2D_[row] ();
3700 for (size_type k = 0; k < numEnt; ++k) {
3701 rowVals[k] *= theAlpha;
3708 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3713 const char tfecfFuncName[] =
"setAllToScalar: ";
3715 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3716 ! isFillActive (), std::runtime_error,
3717 "Fill must be active before you may call this method. "
3718 "Please call resumeFill() to make fill active.");
3724 const size_t nlrs = staticGraph_->getNodeNumRows();
3725 const size_t numEntries = staticGraph_->getNodeNumEntries();
3726 if (! staticGraph_->indicesAreAllocated () || numEntries == 0) {
3730 const ProfileType profType = staticGraph_->getProfileType ();
3738 for (
size_t row = 0; row < nlrs; ++row) {
3739 std::fill (values2D_[row].begin (), values2D_[row].end (), theAlpha);
3745 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3749 const typename local_graph_type::entries_type::non_const_type& columnIndices,
3750 const typename local_matrix_type::values_type& values)
3752 const char tfecfFuncName[] =
"setAllValues: ";
3753 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3754 (columnIndices.size () != values.size (), std::invalid_argument,
3755 "columnIndices.size() = " << columnIndices.size () <<
" != values.size()"
3756 " = " << values.size () <<
".");
3757 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3758 (myGraph_.is_null (), std::runtime_error,
"myGraph_ must not be null.");
3761 myGraph_->setAllIndices (rowPointers, columnIndices);
3763 catch (std::exception &e) {
3764 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3765 (
true, std::runtime_error,
"myGraph_->setAllIndices() threw an "
3766 "exception: " << e.what ());
3772 auto lclGraph = myGraph_->getLocalGraph ();
3773 const size_t numEnt = lclGraph.entries.extent (0);
3774 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3775 (lclGraph.row_map.extent (0) != rowPointers.extent (0) ||
3776 numEnt != static_cast<size_t> (columnIndices.extent (0)),
3777 std::logic_error,
"myGraph_->setAllIndices() did not correctly create "
3778 "local graph. Please report this bug to the Tpetra developers.");
3780 const size_t numCols = myGraph_->getColMap ()->getNodeNumElements ();
3782 numCols, values, lclGraph);
3786 this->k_values1D_ = this->lclMatrix_.values;
3790 this->storageStatus_ = ::Tpetra::Details::STORAGE_1D_PACKED;
3792 checkInternalState ();
3795 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3799 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
3800 const Teuchos::ArrayRCP<Scalar>& val)
3802 using Kokkos::Compat::getKokkosViewDeepCopy;
3803 using Teuchos::ArrayRCP;
3804 using Teuchos::av_reinterpret_cast;
3807 typedef typename local_matrix_type::row_map_type row_map_type;
3809 const char tfecfFuncName[] =
"setAllValues(ArrayRCP<size_t>, ArrayRCP<LO>, ArrayRCP<Scalar>): ";
3815 typename row_map_type::non_const_type ptrNative (
"ptr", ptr.size ());
3816 Kokkos::View<
const size_t*,
3817 typename row_map_type::array_layout,
3819 Kokkos::MemoryUnmanaged> ptrSizeT (ptr.getRawPtr (), ptr.size ());
3822 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3823 (ptrNative.extent (0) != ptrSizeT.extent (0),
3824 std::logic_error,
"ptrNative.extent(0) = " <<
3825 ptrNative.extent (0) <<
" != ptrSizeT.extent(0) = "
3826 << ptrSizeT.extent (0) <<
". Please report this bug to the "
3827 "Tpetra developers.");
3829 auto indIn = getKokkosViewDeepCopy<DT> (ind ());
3830 auto valIn = getKokkosViewDeepCopy<DT> (av_reinterpret_cast<IST> (val ()));
3831 this->setAllValues (ptrNative, indIn, valIn);
3834 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3839 const char tfecfFuncName[] =
"getLocalDiagOffsets: ";
3840 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3841 (staticGraph_.is_null (), std::runtime_error,
"The matrix has no graph.");
3849 if (static_cast<size_t> (offsets.size ()) < lclNumRows) {
3850 offsets.resize (lclNumRows);
3856 typedef typename device_type::memory_space memory_space;
3857 if (std::is_same<memory_space, Kokkos::HostSpace>::value) {
3862 Kokkos::MemoryUnmanaged> output_type;
3863 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3864 staticGraph_->getLocalDiagOffsets (offsetsOut);
3867 Kokkos::View<size_t*, device_type> offsetsTmp (
"diagOffsets", lclNumRows);
3868 staticGraph_->getLocalDiagOffsets (offsetsTmp);
3869 typedef Kokkos::View<
size_t*, Kokkos::HostSpace,
3870 Kokkos::MemoryUnmanaged> output_type;
3871 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3876 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3881 using Teuchos::ArrayRCP;
3882 using Teuchos::ArrayView;
3883 using Teuchos::av_reinterpret_cast;
3884 const char tfecfFuncName[] =
"getLocalDiagCopy (1-arg): ";
3888 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3889 staticGraph_.is_null (), std::runtime_error,
3890 "This method requires that the matrix have a graph.");
3891 auto rowMapPtr = this->getRowMap ();
3892 if (rowMapPtr.is_null () || rowMapPtr->getComm ().is_null ()) {
3898 auto colMapPtr = this->getColMap ();
3899 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3900 (! this->hasColMap () || colMapPtr.is_null (), std::runtime_error,
3901 "This method requires that the matrix have a column Map.");
3902 const map_type& rowMap = * rowMapPtr;
3903 const map_type& colMap = * colMapPtr;
3904 const LO myNumRows = static_cast<LO> (this->getNodeNumRows ());
3906 #ifdef HAVE_TPETRA_DEBUG
3909 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3910 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3911 "The input Vector's Map must be compatible with the CrsMatrix's row "
3912 "Map. You may check this by using Map's isCompatible method: "
3913 "diag.getMap ()->isCompatible (A.getRowMap ());");
3914 #endif // HAVE_TPETRA_DEBUG
3916 if (this->isFillComplete ()) {
3917 diag.template modify<device_type> ();
3918 const auto D_lcl = diag.template getLocalView<device_type> ();
3920 const auto D_lcl_1d =
3921 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3923 const auto lclRowMap = rowMap.getLocalMap ();
3925 const auto lclMatrix = this->lclMatrix_;
3928 lclColMap, lclMatrix);
3936 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3941 Kokkos::MemoryUnmanaged>& offsets)
const
3943 typedef LocalOrdinal LO;
3945 #ifdef HAVE_TPETRA_DEBUG
3946 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3947 const map_type& rowMap = * (this->getRowMap ());
3950 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3951 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3952 "The input Vector's Map must be compatible with (in the sense of Map::"
3953 "isCompatible) the CrsMatrix's row Map.");
3954 #endif // HAVE_TPETRA_DEBUG
3963 diag.template modify<device_type> ();
3964 auto D_lcl = diag.template getLocalView<device_type> ();
3965 const LO myNumRows = static_cast<LO> (this->getNodeNumRows ());
3968 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3970 KokkosSparse::getDiagCopy (D_lcl_1d, offsets, this->lclMatrix_);
3973 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3977 const Teuchos::ArrayView<const size_t>& offsets)
const
3979 typedef LocalOrdinal LO;
3982 typedef typename vec_type::dual_view_type dual_view_type;
3983 typedef typename dual_view_type::host_mirror_space::execution_space host_execution_space;
3985 #ifdef HAVE_TPETRA_DEBUG
3986 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3987 const map_type& rowMap = * (this->getRowMap ());
3990 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3991 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3992 "The input Vector's Map must be compatible with (in the sense of Map::"
3993 "isCompatible) the CrsMatrix's row Map.");
3994 #endif // HAVE_TPETRA_DEBUG
4000 diag_dv.modified_device () = 0;
4005 diag.template modify<host_execution_space> ();
4006 auto lclVecHost = diag.template getLocalView<host_execution_space> ();
4008 auto lclVecHost1d = Kokkos::subview (lclVecHost, Kokkos::ALL (), 0);
4010 Kokkos::View<
const size_t*, Kokkos::HostSpace,
4011 Kokkos::MemoryTraits<Kokkos::Unmanaged> >
4012 h_offsets (offsets.getRawPtr (), offsets.size ());
4014 const LO myNumRows = static_cast<LO> (this->getNodeNumRows ());
4015 typedef Kokkos::RangePolicy<host_execution_space, LO> policy_type;
4016 const size_t INV = Tpetra::Details::OrdinalTraits<size_t>::invalid ();
4018 Kokkos::parallel_for (policy_type (0, myNumRows), [&] (
const LO& lclRow) {
4019 lclVecHost1d(lclRow) = STS::zero ();
4020 if (h_offsets[lclRow] != INV) {
4021 auto curRow = lclMatrix_.rowConst (lclRow);
4022 lclVecHost1d(lclRow) = static_cast<IST> (curRow.value(h_offsets[lclRow]));
4025 diag.template sync<execution_space> ();
4029 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4034 using ::Tpetra::Details::ProfilingRegion;
4035 using Teuchos::ArrayRCP;
4036 using Teuchos::ArrayView;
4037 using Teuchos::null;
4040 using Teuchos::rcpFromRef;
4043 const char tfecfFuncName[] =
"leftScale: ";
4045 ProfilingRegion region (
"Tpetra::CrsMatrix::leftScale");
4047 RCP<const vec_type> xp;
4048 if (this->getRangeMap ()->isSameAs (* (x.
getMap ()))) {
4051 auto exporter = this->getCrsGraphRef ().getExporter ();
4052 if (exporter.get () !=
nullptr) {
4053 RCP<vec_type> tempVec (
new vec_type (this->getRowMap ()));
4054 tempVec->doImport (x, *exporter,
REPLACE);
4058 xp = rcpFromRef (x);
4061 else if (this->getRowMap ()->isSameAs (* (x.
getMap ()))) {
4062 xp = rcpFromRef (x);
4065 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4066 (
true, std::invalid_argument,
"x's Map must be the same as "
4067 "either the row Map or the range Map of the CrsMatrix.");
4075 const LO lclNumRows =
4076 static_cast<LO> (this->getRowMap ()->getNodeNumElements ());
4077 const bool validLocalMatrix = this->lclMatrix_.numRows () == lclNumRows;
4079 if (validLocalMatrix) {
4080 using dev_memory_space =
typename device_type::memory_space;
4081 if (xp->template need_sync<dev_memory_space> ()) {
4082 using Teuchos::rcp_const_cast;
4083 rcp_const_cast<vec_type> (xp)->template sync<dev_memory_space> ();
4085 auto x_lcl = xp->template getLocalView<dev_memory_space> ();
4086 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
4087 ::Tpetra::Details::leftScaleLocalCrsMatrix (this->lclMatrix_, x_lcl_1d,
false,
false);
4090 execution_space::fence ();
4092 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
4093 ArrayView<impl_scalar_type> rowValues = Teuchos::null;
4094 for (LocalOrdinal i = 0; i < lclNumRows; ++i) {
4095 const RowInfo rowinfo = this->staticGraph_->getRowInfo (i);
4096 rowValues = this->getViewNonConst (rowinfo);
4097 const impl_scalar_type scaleValue = static_cast<impl_scalar_type> (vectorVals[i]);
4098 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
4099 rowValues[j] *= scaleValue;
4102 execution_space::fence ();
4106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4111 using ::Tpetra::Details::ProfilingRegion;
4112 using Teuchos::ArrayRCP;
4113 using Teuchos::ArrayView;
4114 using Teuchos::null;
4117 using Teuchos::rcpFromRef;
4120 const char tfecfFuncName[] =
"rightScale: ";
4122 ProfilingRegion region (
"Tpetra::CrsMatrix::rightScale");
4124 RCP<const vec_type> xp;
4125 if (this->getDomainMap ()->isSameAs (* (x.
getMap ()))) {
4128 auto importer = this->getCrsGraphRef ().getImporter ();
4129 if (importer.get () !=
nullptr) {
4130 RCP<vec_type> tempVec (
new vec_type (this->getColMap ()));
4131 tempVec->doImport (x, *importer,
REPLACE);
4135 xp = rcpFromRef (x);
4138 else if (this->getColMap ()->isSameAs (* (x.
getMap ()))) {
4139 xp = rcpFromRef (x);
4141 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4142 (
true, std::runtime_error,
"x's Map must be the same as "
4143 "either the domain Map or the column Map of the CrsMatrix.");
4151 const LO lclNumRows =
4152 static_cast<LO> (this->getRowMap ()->getNodeNumElements ());
4153 const bool validLocalMatrix = this->lclMatrix_.numRows () == lclNumRows;
4155 if (validLocalMatrix) {
4156 using dev_memory_space =
typename device_type::memory_space;
4157 if (xp->template need_sync<dev_memory_space> ()) {
4158 using Teuchos::rcp_const_cast;
4159 rcp_const_cast<vec_type> (xp)->template sync<dev_memory_space> ();
4161 auto x_lcl = xp->template getLocalView<dev_memory_space> ();
4162 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
4163 ::Tpetra::Details::rightScaleLocalCrsMatrix (this->lclMatrix_, x_lcl_1d,
false,
false);
4166 execution_space::fence ();
4168 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
4169 ArrayView<impl_scalar_type> rowValues =
null;
4170 for (LO i = 0; i < lclNumRows; ++i) {
4171 const RowInfo rowinfo = this->staticGraph_->getRowInfo (i);
4172 rowValues = this->getViewNonConst (rowinfo);
4173 ArrayView<const LO> colInds;
4174 this->getCrsGraphRef ().getLocalRowView (i, colInds);
4175 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
4176 rowValues[j] *= static_cast<impl_scalar_type> (vectorVals[colInds[j]]);
4179 execution_space::fence ();
4183 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4188 using Teuchos::ArrayView;
4189 using Teuchos::outArg;
4190 using Teuchos::REDUCE_SUM;
4191 using Teuchos::reduceAll;
4192 typedef typename Teuchos::ArrayRCP<const impl_scalar_type>::size_type size_type;
4199 mag_type frobNorm = frobNorm_;
4200 if (frobNorm == -STM::one ()) {
4201 mag_type mySum = STM::zero ();
4202 if (getNodeNumEntries() > 0) {
4203 if (isStorageOptimized ()) {
4206 const size_type numEntries =
4207 static_cast<size_type> (getNodeNumEntries ());
4208 for (size_type k = 0; k < numEntries; ++k) {
4214 const mag_type val_abs = STS::abs (val);
4215 mySum += val_abs * val_abs;
4219 const LocalOrdinal numRows =
4220 static_cast<LocalOrdinal> (this->getNodeNumRows ());
4221 for (LocalOrdinal r = 0; r < numRows; ++r) {
4222 const RowInfo rowInfo = myGraph_->getRowInfo (r);
4223 const size_type numEntries =
4224 static_cast<size_type> (rowInfo.numEntries);
4225 ArrayView<const impl_scalar_type> A_r =
4226 this->getView (rowInfo).view (0, numEntries);
4227 for (size_type k = 0; k < numEntries; ++k) {
4229 const mag_type val_abs = STS::abs (val);
4230 mySum += val_abs * val_abs;
4235 mag_type totalSum = STM::zero ();
4236 reduceAll<int, mag_type> (* (getComm ()), REDUCE_SUM,
4237 mySum, outArg (totalSum));
4238 frobNorm = STM::sqrt (totalSum);
4240 if (isFillComplete ()) {
4244 frobNorm_ = frobNorm;
4249 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4254 const char tfecfFuncName[] =
"replaceColMap: ";
4258 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4259 myGraph_.is_null (), std::runtime_error,
4260 "This method does not work if the matrix has a const graph. The whole "
4261 "idea of a const graph is that you are not allowed to change it, but "
4262 "this method necessarily must modify the graph, since the graph owns "
4263 "the matrix's column Map.");
4267 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4271 const Teuchos::RCP<const map_type>& newColMap,
4272 const Teuchos::RCP<const import_type>& newImport,
4273 const bool sortEachRow)
4275 const char tfecfFuncName[] =
"reindexColumns: ";
4276 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4277 graph == NULL && myGraph_.is_null (), std::invalid_argument,
4278 "The input graph is NULL, but the matrix does not own its graph.");
4281 const bool sortGraph =
false;
4283 if (sortEachRow && theGraph.isLocallyIndexed () && ! theGraph.isSorted ()) {
4284 const LocalOrdinal lclNumRows =
4285 static_cast<LocalOrdinal> (theGraph.getNodeNumRows ());
4286 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
4287 const RowInfo rowInfo = theGraph.getRowInfo (row);
4288 auto lclColInds = theGraph.getLocalKokkosRowViewNonConst (rowInfo);
4289 auto vals = this->getRowViewNonConst (rowInfo);
4292 sort2 (lclColInds.data (),
4293 lclColInds.data () + rowInfo.numEntries,
4296 theGraph.indicesAreSorted_ =
true;
4300 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4304 Teuchos::RCP<const import_type>& newImporter)
4306 const char tfecfFuncName[] =
"replaceDomainMapAndImporter: ";
4307 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4308 myGraph_.is_null (), std::runtime_error,
4309 "This method does not work if the matrix has a const graph. The whole "
4310 "idea of a const graph is that you are not allowed to change it, but this"
4311 " method necessarily must modify the graph, since the graph owns the "
4312 "matrix's domain Map and Import objects.");
4316 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4320 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
4321 const Teuchos::ArrayView<const Scalar>& values)
4323 using Teuchos::Array;
4324 typedef GlobalOrdinal GO;
4325 typedef typename Array<GO>::size_type size_type;
4327 const size_type numToInsert = indices.size ();
4330 std::pair<Array<GO>, Array<Scalar> >& curRow = nonlocals_[globalRow];
4331 Array<GO>& curRowInds = curRow.first;
4332 Array<Scalar>& curRowVals = curRow.second;
4333 const size_type newCapacity = curRowInds.size () + numToInsert;
4334 curRowInds.reserve (newCapacity);
4335 curRowVals.reserve (newCapacity);
4336 for (size_type k = 0; k < numToInsert; ++k) {
4337 curRowInds.push_back (indices[k]);
4338 curRowVals.push_back (values[k]);
4342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4344 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
4347 using ::Tpetra::Details::ProfilingRegion;
4348 using Teuchos::Comm;
4349 using Teuchos::outArg;
4352 using Teuchos::REDUCE_MAX;
4353 using Teuchos::REDUCE_MIN;
4354 using Teuchos::reduceAll;
4357 typedef GlobalOrdinal GO;
4358 typedef typename Teuchos::Array<GO>::size_type size_type;
4359 const char tfecfFuncName[] =
"globalAssemble: ";
4360 ProfilingRegion regionGlobalAssemble (
"Tpetra::CrsMatrix::globalAssemble");
4362 RCP<const Comm<int> > comm = getComm ();
4364 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4365 (! isFillActive (), std::runtime_error,
"Fill must be active before "
4366 "you may call this method.");
4368 const size_t myNumNonlocalRows = nonlocals_.size ();
4375 const int iHaveNonlocalRows = (myNumNonlocalRows == 0) ? 0 : 1;
4376 int someoneHasNonlocalRows = 0;
4377 reduceAll<int, int> (*comm, REDUCE_MAX, iHaveNonlocalRows,
4378 outArg (someoneHasNonlocalRows));
4379 if (someoneHasNonlocalRows == 0) {
4393 RCP<const map_type> nonlocalRowMap;
4395 Teuchos::ArrayRCP<size_t> numEntPerNonlocalRow (myNumNonlocalRows);
4397 Teuchos::Array<GO> myNonlocalGblRows (myNumNonlocalRows);
4398 size_type curPos = 0;
4399 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
4400 ++mapIter, ++curPos) {
4401 myNonlocalGblRows[curPos] = mapIter->first;
4404 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4405 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4412 sort2 (gblCols.begin (), gblCols.end (), vals.begin ());
4413 typename Teuchos::Array<GO>::iterator gblCols_newEnd;
4414 typename Teuchos::Array<Scalar>::iterator vals_newEnd;
4415 merge2 (gblCols_newEnd, vals_newEnd,
4416 gblCols.begin (), gblCols.end (),
4417 vals.begin (), vals.end ());
4418 gblCols.erase (gblCols_newEnd, gblCols.end ());
4419 vals.erase (vals_newEnd, vals.end ());
4420 numEntPerNonlocalRow[curPos] = gblCols.size ();
4431 GO myMinNonlocalGblRow = std::numeric_limits<GO>::max ();
4433 auto iter = std::min_element (myNonlocalGblRows.begin (),
4434 myNonlocalGblRows.end ());
4435 if (iter != myNonlocalGblRows.end ()) {
4436 myMinNonlocalGblRow = *iter;
4439 GO gblMinNonlocalGblRow = 0;
4440 reduceAll<int, GO> (*comm, REDUCE_MIN, myMinNonlocalGblRow,
4441 outArg (gblMinNonlocalGblRow));
4442 const GO indexBase = gblMinNonlocalGblRow;
4443 const global_size_t INV = Teuchos::OrdinalTraits<global_size_t>::invalid ();
4444 nonlocalRowMap = rcp (
new map_type (INV, myNonlocalGblRows (), indexBase, comm));
4452 RCP<crs_matrix_type> nonlocalMatrix =
4453 rcp (
new crs_matrix_type (nonlocalRowMap, numEntPerNonlocalRow,
4456 size_type curPos = 0;
4457 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
4458 ++mapIter, ++curPos) {
4459 const GO gblRow = mapIter->first;
4461 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4462 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4464 nonlocalMatrix->insertGlobalValues (gblRow, gblCols (), vals ());
4476 auto origRowMap = this->getRowMap ();
4477 const bool origRowMapIsOneToOne = origRowMap->isOneToOne ();
4479 int isLocallyComplete = 1;
4481 if (origRowMapIsOneToOne) {
4482 export_type exportToOrig (nonlocalRowMap, origRowMap);
4484 isLocallyComplete = 0;
4486 this->doExport (*nonlocalMatrix, exportToOrig,
Tpetra::ADD);
4495 export_type exportToOneToOne (nonlocalRowMap, oneToOneRowMap);
4497 isLocallyComplete = 0;
4504 crs_matrix_type oneToOneMatrix (oneToOneRowMap, 0);
4506 oneToOneMatrix.doExport (*nonlocalMatrix, exportToOneToOne,
Tpetra::ADD);
4510 nonlocalMatrix = Teuchos::null;
4513 import_type importToOrig (oneToOneRowMap, origRowMap);
4514 this->doImport (oneToOneMatrix, importToOrig,
Tpetra::ADD);
4521 decltype (nonlocals_) newNonlocals;
4522 std::swap (nonlocals_, newNonlocals);
4531 int isGloballyComplete = 0;
4532 reduceAll<int, int> (*comm, REDUCE_MIN, isLocallyComplete,
4533 outArg (isGloballyComplete));
4534 TEUCHOS_TEST_FOR_EXCEPTION
4535 (isGloballyComplete != 1, std::runtime_error,
"On at least one process, "
4536 "you called insertGlobalValues with a global row index which is not in "
4537 "the matrix's row Map on any process in its communicator.");
4540 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4545 if (! isStaticGraph ()) {
4548 clearGlobalConstants ();
4549 fillComplete_ =
false;
4552 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4566 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4573 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4585 frobNorm_ = -STM::one ();
4588 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4593 const char tfecfFuncName[] =
"fillComplete(params): ";
4595 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4596 (this->getCrsGraph ().is_null (), std::logic_error,
4597 "getCrsGraph() returns null. This should not happen at this point. "
4598 "Please report this bug to the Tpetra developers.");
4608 Teuchos::RCP<const map_type> rangeMap = graph.
getRowMap ();
4609 Teuchos::RCP<const map_type> domainMap = rangeMap;
4610 this->fillComplete (domainMap, rangeMap, params);
4614 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4618 const Teuchos::RCP<const map_type>& rangeMap,
4619 const Teuchos::RCP<Teuchos::ParameterList>& params)
4621 using ::Tpetra::Details::ProfilingRegion;
4622 using Teuchos::ArrayRCP;
4625 const char tfecfFuncName[] =
"fillComplete: ";
4626 ProfilingRegion regionFillComplete (
"Tpetra::CrsMatrix::fillComplete");
4628 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4629 (! this->isFillActive () || this->isFillComplete (), std::runtime_error,
4630 "Matrix fill state must be active (isFillActive() "
4631 "must be true) before you may call fillComplete().");
4632 const int numProcs = this->getComm ()->getSize ();
4640 bool assertNoNonlocalInserts =
false;
4643 bool sortGhosts =
true;
4645 if (! params.is_null ()) {
4646 assertNoNonlocalInserts = params->get (
"No Nonlocal Changes",
4647 assertNoNonlocalInserts);
4648 if (params->isParameter (
"sort column map ghost gids")) {
4649 sortGhosts = params->get (
"sort column map ghost gids", sortGhosts);
4651 else if (params->isParameter (
"Sort column Map ghost GIDs")) {
4652 sortGhosts = params->get (
"Sort column Map ghost GIDs", sortGhosts);
4657 const bool needGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
4659 if (! this->myGraph_.is_null ()) {
4660 this->myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
4663 if (! this->getCrsGraphRef ().indicesAreAllocated ()) {
4664 if (this->hasColMap ()) {
4666 this->allocateValues (LocalIndices, GraphNotYetAllocated);
4669 this->allocateValues (GlobalIndices, GraphNotYetAllocated);
4674 if (needGlobalAssemble) {
4675 this->globalAssemble ();
4678 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4679 (numProcs == 1 && nonlocals_.size() > 0,
4680 std::runtime_error,
"Cannot have nonlocal entries on a serial run. "
4681 "An invalid entry (i.e., with row index not in the row Map) must have "
4682 "been submitted to the CrsMatrix.");
4685 if (this->isStaticGraph ()) {
4692 #ifdef HAVE_TPETRA_DEBUG
4710 const bool domainMapsMatch =
4711 this->staticGraph_->getDomainMap ()->isSameAs (*domainMap);
4712 const bool rangeMapsMatch =
4713 this->staticGraph_->getRangeMap ()->isSameAs (*rangeMap);
4715 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4716 (! domainMapsMatch, std::runtime_error,
4717 "The CrsMatrix's domain Map does not match the graph's domain Map. "
4718 "The graph cannot be changed because it was given to the CrsMatrix "
4719 "constructor as const. You can fix this by passing in the graph's "
4720 "domain Map and range Map to the matrix's fillComplete call.");
4722 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4723 (! rangeMapsMatch, std::runtime_error,
4724 "The CrsMatrix's range Map does not match the graph's range Map. "
4725 "The graph cannot be changed because it was given to the CrsMatrix "
4726 "constructor as const. You can fix this by passing in the graph's "
4727 "domain Map and range Map to the matrix's fillComplete call.");
4728 #endif // HAVE_TPETRA_DEBUG
4732 this->fillLocalMatrix (params);
4739 this->myGraph_->setDomainRangeMaps (domainMap, rangeMap);
4742 Teuchos::Array<int> remotePIDs (0);
4743 const bool mustBuildColMap = ! this->hasColMap ();
4744 if (mustBuildColMap) {
4745 this->myGraph_->makeColMap (remotePIDs);
4750 const std::pair<size_t, std::string> makeIndicesLocalResult =
4751 this->myGraph_->makeIndicesLocal ();
4756 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4757 (makeIndicesLocalResult.first != 0, std::runtime_error,
4758 makeIndicesLocalResult.second);
4760 const bool sorted = this->myGraph_->isSorted ();
4761 const bool merged = this->myGraph_->isMerged ();
4762 this->sortAndMergeIndicesAndValues (sorted, merged);
4767 this->myGraph_->makeImportExport (remotePIDs, mustBuildColMap);
4771 this->fillLocalGraphAndMatrix (params);
4773 const bool callGraphComputeGlobalConstants = params.get () ==
nullptr ||
4774 params->get (
"compute global constants",
true);
4775 const bool computeLocalTriangularConstants = params.get () ==
nullptr ||
4776 params->get (
"compute local triangular constants",
true);
4777 if (callGraphComputeGlobalConstants) {
4778 this->myGraph_->computeGlobalConstants (computeLocalTriangularConstants);
4781 this->myGraph_->computeLocalConstants (computeLocalTriangularConstants);
4783 this->myGraph_->fillComplete_ =
true;
4784 this->myGraph_->checkInternalState ();
4787 const bool callComputeGlobalConstants = params.get () ==
nullptr ||
4788 params->get (
"compute global constants",
true);
4789 if (callComputeGlobalConstants) {
4790 this->computeGlobalConstants ();
4795 this->fillComplete_ =
true;
4796 this->checkInternalState ();
4799 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4803 const Teuchos::RCP<const map_type> & rangeMap,
4804 const Teuchos::RCP<const import_type>& importer,
4805 const Teuchos::RCP<const export_type>& exporter,
4806 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
4808 #ifdef HAVE_TPETRA_MMM_TIMINGS
4810 if(!params.is_null())
4811 label = params->get(
"Timer Label",label);
4812 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
4813 using Teuchos::TimeMonitor;
4814 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-Graph"))));
4818 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
4819 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
4820 std::runtime_error,
"Matrix fill state must be active (isFillActive() "
4821 "must be true) before calling fillComplete().");
4822 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4823 myGraph_.is_null (), std::logic_error,
"myGraph_ is null. This is not allowed.");
4827 myGraph_->expertStaticFillComplete (domainMap, rangeMap, importer, exporter,params);
4829 const bool callComputeGlobalConstants = params.get () ==
nullptr ||
4830 params->get (
"compute global constants",
true);
4831 if (callComputeGlobalConstants) {
4832 #ifdef HAVE_TPETRA_MMM_TIMINGS
4833 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cGC"))));
4835 this->computeGlobalConstants ();
4838 #ifdef HAVE_TPETRA_MMM_TIMINGS
4839 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-fLGAM"))));
4843 fillLocalGraphAndMatrix (params);
4848 fillComplete_ =
true;
4851 #ifdef HAVE_TPETRA_DEBUG
4852 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isFillActive(), std::logic_error,
4853 ": We're at the end of fillComplete(), but isFillActive() is true. "
4854 "Please report this bug to the Tpetra developers.");
4855 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillComplete(), std::logic_error,
4856 ": We're at the end of fillComplete(), but isFillActive() is true. "
4857 "Please report this bug to the Tpetra developers.");
4858 #endif // HAVE_TPETRA_DEBUG
4860 #ifdef HAVE_TPETRA_MMM_TIMINGS
4861 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cIS"))));
4864 checkInternalState();
4867 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4873 #ifdef HAVE_TPETRA_DEBUG
4874 const char tfecfFuncName[] =
"mergeRowIndicesAndValues: ";
4875 #endif // HAVE_TPETRA_DEBUG
4877 auto rowValues = this->getRowViewNonConst (rowInfo);
4878 typedef typename std::decay<decltype (rowValues[0]) >::type value_type;
4879 value_type* rowValueIter = rowValues.data ();
4880 auto inds_view = graph.getLocalKokkosRowViewNonConst (rowInfo);
4883 LocalOrdinal* beg = inds_view.data ();
4884 LocalOrdinal* end = inds_view.data () + rowInfo.numEntries;
4886 #ifdef HAVE_TPETRA_DEBUG
4887 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4888 (rowInfo.allocSize != static_cast<size_t> (inds_view.extent (0)) ||
4889 rowInfo.allocSize != static_cast<size_t> (rowValues.extent (0)),
4890 std::runtime_error,
"rowInfo.allocSize = " << rowInfo.allocSize
4891 <<
" != inds_view.extent(0) = " << inds_view.extent (0)
4892 <<
" || rowInfo.allocSize = " << rowInfo.allocSize
4893 <<
" != rowValues.extent(0) = " << rowValues.extent (0) <<
".");
4894 #endif // HAVE_TPETRA_DEBUG
4896 LocalOrdinal* newend = beg;
4898 LocalOrdinal* cur = beg + 1;
4899 value_type* vcur = rowValueIter + 1;
4900 value_type* vend = rowValueIter;
4902 while (cur != end) {
4903 if (*cur != *newend) {
4920 const size_t mergedEntries = newend - beg;
4922 const size_t numDups = rowInfo.numEntries - mergedEntries;
4926 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4931 using ::Tpetra::Details::ProfilingRegion;
4932 typedef LocalOrdinal LO;
4933 typedef typename Kokkos::View<LO*, device_type>::HostMirror::execution_space
4934 host_execution_space;
4935 typedef Kokkos::RangePolicy<host_execution_space, LO> range_type;
4937 const char tfecfFuncName[] =
"sortAndMergeIndicesAndValues: ";
4938 ProfilingRegion regionSAM (
"Tpetra::CrsMatrix::sortAndMergeIndicesAndValues");
4940 if (! sorted || ! merged) {
4941 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4942 (this->isStaticGraph (), std::runtime_error,
"Cannot sort or merge with "
4943 "\"static\" (const) graph, since the matrix does not own the graph.");
4944 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4945 (this->myGraph_.is_null (), std::logic_error,
"myGraph_ is null, but "
4946 "this matrix claims ! isStaticGraph(). "
4947 "Please report this bug to the Tpetra developers.");
4948 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4949 (this->isStorageOptimized (), std::logic_error,
"It is invalid to call "
4950 "this method if the graph's storage has already been optimized. "
4951 "Please report this bug to the Tpetra developers.");
4954 const LO lclNumRows = static_cast<LO> (this->getNodeNumRows ());
4955 size_t totalNumDups = 0;
4957 Kokkos::parallel_reduce (range_type (0, lclNumRows),
4958 [
this, &graph, sorted, merged] (
const LO& lclRow,
size_t& numDups) {
4961 auto lclColInds = graph.getLocalKokkosRowViewNonConst (rowInfo);
4962 auto vals = this->getRowViewNonConst (rowInfo);
4965 sort2 (lclColInds.data (),
4966 lclColInds.data () + rowInfo.numEntries,
4970 numDups += this->mergeRowIndicesAndValues (graph, rowInfo);
4982 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4991 using Teuchos::null;
4994 using Teuchos::rcp_const_cast;
4995 using Teuchos::rcpFromRef;
4996 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4997 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
5003 if (alpha ==
ZERO) {
5006 }
else if (beta != ONE) {
5020 RCP<const import_type> importer = this->getGraph ()->getImporter ();
5021 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
5027 const bool Y_is_overwritten = (beta ==
ZERO);
5030 const bool Y_is_replicated =
5031 (! Y_in.
isDistributed () && this->getComm ()->getSize () != 1);
5039 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
5046 RCP<const MV> X_colMap;
5047 if (importer.is_null ()) {
5055 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in,
true);
5057 X_colMap = rcp_const_cast<const MV> (X_colMapNonConst);
5062 X_colMap = rcpFromRef (X_in);
5066 ProfilingRegion regionImport (
"Tpetra::CrsMatrix::apply: Import");
5072 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in);
5075 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
5076 X_colMap = rcp_const_cast<const MV> (X_colMapNonConst);
5083 RCP<MV> Y_rowMap = getRowMapMultiVector (Y_in);
5090 if (! exporter.is_null ()) {
5091 this->localApply (*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha,
ZERO);
5093 ProfilingRegion regionExport (
"Tpetra::CrsMatrix::apply: Export");
5099 if (Y_is_overwritten) {
5125 Y_rowMap = getRowMapMultiVector (Y_in,
true);
5132 this->localApply (*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
5136 this->localApply (*X_colMap, Y_in, Teuchos::NO_TRANS, alpha, beta);
5144 if (Y_is_replicated) {
5145 ProfilingRegion regionReduce (
"Tpetra::CrsMatrix::apply: Reduce Y");
5150 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5155 const Teuchos::ETransp mode,
5160 using Teuchos::null;
5163 using Teuchos::rcp_const_cast;
5164 using Teuchos::rcpFromRef;
5165 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
5168 if (alpha ==
ZERO) {
5191 RCP<const import_type> importer = this->getGraph ()->getImporter ();
5192 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
5198 const bool Y_is_overwritten = (beta ==
ZERO);
5199 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
5205 X = rcp (
new MV (X_in, Teuchos::Copy));
5207 X = rcpFromRef (X_in);
5211 if (importer != Teuchos::null) {
5212 if (importMV_ != Teuchos::null && importMV_->getNumVectors() != numVectors) {
5215 if (importMV_ ==
null) {
5216 importMV_ = rcp (
new MV (this->getColMap (), numVectors));
5219 if (exporter != Teuchos::null) {
5220 if (exportMV_ != Teuchos::null && exportMV_->getNumVectors() != numVectors) {
5223 if (exportMV_ ==
null) {
5224 exportMV_ = rcp (
new MV (this->getRowMap (), numVectors));
5230 if (! exporter.is_null ()) {
5231 ProfilingRegion regionImport (
"Tpetra::CrsMatrix::apply (transpose): Import");
5232 exportMV_->doImport (X_in, *exporter,
INSERT);
5239 if (importer != Teuchos::null) {
5240 ProfilingRegion regionExport (
"Tpetra::CrsMatrix::apply (transpose): Export");
5247 importMV_->putScalar (
ZERO);
5249 this->localApply (*X, *importMV_, mode, alpha,
ZERO);
5250 if (Y_is_overwritten) {
5267 MV Y (Y_in, Teuchos::Copy);
5268 this->localApply (*X, Y, mode, alpha, beta);
5271 this->localApply (*X, Y_in, mode, alpha, beta);
5278 if (Y_is_replicated) {
5279 ProfilingRegion regionReduce (
"Tpetra::CrsMatrix::apply (transpose): Reduce Y");
5284 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5289 const Teuchos::ETransp mode,
5290 const Scalar& alpha,
5291 const Scalar& beta)
const
5294 ProfilingRegion regionLocalApply (
"Tpetra::CrsMatrix::localApply");
5295 this->
template localMultiply<Scalar, Scalar> (X, Y, mode, alpha, beta);
5298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5303 Teuchos::ETransp mode,
5308 const char fnName[] =
"Tpetra::CrsMatrix::apply";
5310 TEUCHOS_TEST_FOR_EXCEPTION
5311 (! isFillComplete (), std::runtime_error,
5312 fnName <<
": Cannot call apply() until fillComplete() "
5313 "has been called.");
5315 if (mode == Teuchos::NO_TRANS) {
5316 ProfilingRegion regionNonTranspose (fnName);
5317 this->applyNonTranspose (X, Y, alpha, beta);
5320 ProfilingRegion regionTranspose (
"Tpetra::CrsMatrix::apply (transpose)");
5327 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
5331 this->applyTranspose (X, Y, mode, alpha, beta);
5335 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5341 const Scalar& dampingFactor,
5343 const int numSweeps)
const
5345 reorderedGaussSeidel (B, X, D, Teuchos::null, dampingFactor, direction, numSweeps);
5348 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5354 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
5355 const Scalar& dampingFactor,
5357 const int numSweeps)
const
5359 using Teuchos::null;
5362 using Teuchos::rcp_const_cast;
5363 using Teuchos::rcpFromRef;
5366 TEUCHOS_TEST_FOR_EXCEPTION(
5367 isFillComplete() ==
false, std::runtime_error,
5368 "Tpetra::CrsMatrix::gaussSeidel: cannot call this method until "
5369 "fillComplete() has been called.");
5370 TEUCHOS_TEST_FOR_EXCEPTION(
5372 std::invalid_argument,
5373 "Tpetra::CrsMatrix::gaussSeidel: The number of sweeps must be , "
5374 "nonnegative but you provided numSweeps = " << numSweeps <<
" < 0.");
5379 if (direction == Forward) {
5380 localDirection = Forward;
5382 else if (direction == Backward) {
5383 localDirection = Backward;
5385 else if (direction == Symmetric) {
5387 localDirection = Forward;
5390 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
5391 "Tpetra::CrsMatrix::gaussSeidel: The 'direction' enum does not have "
5392 "any of its valid values: Forward, Backward, or Symmetric.");
5395 if (numSweeps == 0) {
5402 RCP<const import_type> importer = this->getGraph()->getImporter();
5403 RCP<const export_type> exporter = this->getGraph()->getExporter();
5404 TEUCHOS_TEST_FOR_EXCEPTION(
5405 ! exporter.is_null (), std::runtime_error,
5406 "Tpetra's gaussSeidel implementation requires that the row, domain, "
5407 "and range Maps be the same. This cannot be the case, because the "
5408 "matrix has a nontrivial Export object.");
5410 RCP<const map_type> domainMap = this->getDomainMap ();
5411 RCP<const map_type> rangeMap = this->getRangeMap ();
5412 RCP<const map_type> rowMap = this->getGraph ()->getRowMap ();
5413 RCP<const map_type> colMap = this->getGraph ()->getColMap ();
5415 #ifdef HAVE_TEUCHOS_DEBUG
5420 TEUCHOS_TEST_FOR_EXCEPTION(
5421 ! X.
getMap ()->isSameAs (*domainMap),
5423 "Tpetra::CrsMatrix::gaussSeidel requires that the input "
5424 "multivector X be in the domain Map of the matrix.");
5425 TEUCHOS_TEST_FOR_EXCEPTION(
5426 ! B.
getMap ()->isSameAs (*rangeMap),
5428 "Tpetra::CrsMatrix::gaussSeidel requires that the input "
5429 "B be in the range Map of the matrix.");
5430 TEUCHOS_TEST_FOR_EXCEPTION(
5431 ! D.
getMap ()->isSameAs (*rowMap),
5433 "Tpetra::CrsMatrix::gaussSeidel requires that the input "
5434 "D be in the row Map of the matrix.");
5435 TEUCHOS_TEST_FOR_EXCEPTION(
5436 ! rowMap->isSameAs (*rangeMap),
5438 "Tpetra::CrsMatrix::gaussSeidel requires that the row Map and the "
5439 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
5440 TEUCHOS_TEST_FOR_EXCEPTION(
5441 ! domainMap->isSameAs (*rangeMap),
5443 "Tpetra::CrsMatrix::gaussSeidel requires that the domain Map and "
5444 "the range Map of the matrix be the same.");
5450 #endif // HAVE_TEUCHOS_DEBUG
5460 B_in = rcpFromRef (B);
5467 RCP<MV> B_in_nonconst = getRowMapMultiVector (B,
true);
5469 B_in = rcp_const_cast<const MV> (B_in_nonconst);
5474 "gaussSeidel: The current implementation of the Gauss-Seidel kernel "
5475 "requires that X and B both have constant stride. Since B does not "
5476 "have constant stride, we had to make a copy. This is a limitation of "
5477 "the current implementation and not your fault, but we still report it "
5478 "as an efficiency warning for your information.");
5487 RCP<MV> X_domainMap;
5489 bool copiedInput =
false;
5491 if (importer.is_null ()) {
5493 X_domainMap = rcpFromRef (X);
5494 X_colMap = X_domainMap;
5495 copiedInput =
false;
5501 X_colMap = getColumnMapMultiVector (X,
true);
5502 X_domainMap = X_colMap;
5507 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the "
5508 "Gauss-Seidel kernel requires that X and B both have constant "
5509 "stride. Since X does not have constant stride, we had to make a "
5510 "copy. This is a limitation of the current implementation and not "
5511 "your fault, but we still report it as an efficiency warning for "
5512 "your information.");
5517 X_domainMap = rcpFromRef (X);
5521 X_colMap = X_domainMap->offsetViewNonConst (colMap, 0);
5531 X_colMap->doImport (X, *importer,
INSERT);
5532 copiedInput =
false;
5539 X_colMap = getColumnMapMultiVector (X,
true);
5540 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
5541 X_colMap->doImport (X, *importer,
INSERT);
5545 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the "
5546 "Gauss-Seidel kernel requires that X and B both have constant stride. "
5547 "Since X does not have constant stride, we had to make a copy. "
5548 "This is a limitation of the current implementation and not your fault, "
5549 "but we still report it as an efficiency warning for your information.");
5553 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
5554 if (! importer.is_null () && sweep > 0) {
5556 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
5560 if (direction != Symmetric) {
5561 if (rowIndices.is_null ()) {
5562 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5567 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5574 const bool doImportBetweenDirections =
false;
5575 if (rowIndices.is_null ()) {
5576 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5583 if (doImportBetweenDirections) {
5585 if (! importer.is_null ()) {
5586 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
5589 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5594 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5598 if (doImportBetweenDirections) {
5600 if (! importer.is_null ()) {
5601 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
5604 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5617 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5623 const Scalar& dampingFactor,
5625 const int numSweeps,
5626 const bool zeroInitialGuess)
const
5628 reorderedGaussSeidelCopy (X, B, D, Teuchos::null, dampingFactor, direction,
5629 numSweeps, zeroInitialGuess);
5632 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5638 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
5639 const Scalar& dampingFactor,
5641 const int numSweeps,
5642 const bool zeroInitialGuess)
const
5644 using Teuchos::null;
5647 using Teuchos::rcpFromRef;
5648 using Teuchos::rcp_const_cast;
5650 const char prefix[] =
"Tpetra::CrsMatrix::(reordered)gaussSeidelCopy: ";
5651 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
5653 TEUCHOS_TEST_FOR_EXCEPTION(
5654 ! isFillComplete (), std::runtime_error,
5655 prefix <<
"The matrix is not fill complete.");
5656 TEUCHOS_TEST_FOR_EXCEPTION(
5657 numSweeps < 0, std::invalid_argument,
5658 prefix <<
"The number of sweeps must be nonnegative, "
5659 "but you provided numSweeps = " << numSweeps <<
" < 0.");
5664 if (direction == Forward) {
5665 localDirection = Forward;
5667 else if (direction == Backward) {
5668 localDirection = Backward;
5670 else if (direction == Symmetric) {
5672 localDirection = Forward;
5675 TEUCHOS_TEST_FOR_EXCEPTION(
5676 true, std::invalid_argument,
5677 prefix <<
"The 'direction' enum does not have any of its valid "
5678 "values: Forward, Backward, or Symmetric.");
5681 if (numSweeps == 0) {
5685 RCP<const import_type> importer = this->getGraph ()->getImporter ();
5686 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
5687 TEUCHOS_TEST_FOR_EXCEPTION(
5688 ! exporter.is_null (), std::runtime_error,
5689 "This method's implementation currently requires that the matrix's row, "
5690 "domain, and range Maps be the same. This cannot be the case, because "
5691 "the matrix has a nontrivial Export object.");
5693 RCP<const map_type> domainMap = this->getDomainMap ();
5694 RCP<const map_type> rangeMap = this->getRangeMap ();
5695 RCP<const map_type> rowMap = this->getGraph ()->getRowMap ();
5696 RCP<const map_type> colMap = this->getGraph ()->getColMap ();
5698 #ifdef HAVE_TEUCHOS_DEBUG
5703 TEUCHOS_TEST_FOR_EXCEPTION(
5704 ! X.
getMap ()->isSameAs (*domainMap), std::runtime_error,
5705 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
5706 "multivector X be in the domain Map of the matrix.");
5707 TEUCHOS_TEST_FOR_EXCEPTION(
5708 ! B.
getMap ()->isSameAs (*rangeMap), std::runtime_error,
5709 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
5710 "B be in the range Map of the matrix.");
5711 TEUCHOS_TEST_FOR_EXCEPTION(
5712 ! D.
getMap ()->isSameAs (*rowMap), std::runtime_error,
5713 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input "
5714 "D be in the row Map of the matrix.");
5715 TEUCHOS_TEST_FOR_EXCEPTION(
5716 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
5717 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the "
5718 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
5719 TEUCHOS_TEST_FOR_EXCEPTION(
5720 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
5721 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and "
5722 "the range Map of the matrix be the same.");
5728 #endif // HAVE_TEUCHOS_DEBUG
5739 RCP<MV> X_domainMap;
5740 bool copyBackOutput =
false;
5741 if (importer.is_null ()) {
5743 X_colMap = rcpFromRef (X);
5744 X_domainMap = rcpFromRef (X);
5750 if (zeroInitialGuess) {
5751 X_colMap->putScalar (
ZERO);
5759 X_colMap = getColumnMapMultiVector (X,
true);
5763 X_domainMap = X_colMap;
5764 if (! zeroInitialGuess) {
5767 }
catch (std::exception& e) {
5768 std::ostringstream os;
5769 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: "
5770 "deep_copy(*X_domainMap, X) threw an exception: "
5771 << e.what () <<
".";
5772 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
5775 copyBackOutput =
true;
5779 "gaussSeidelCopy: The current implementation of the Gauss-Seidel "
5780 "kernel requires that X and B both have constant stride. Since X "
5781 "does not have constant stride, we had to make a copy. This is a "
5782 "limitation of the current implementation and not your fault, but we "
5783 "still report it as an efficiency warning for your information.");
5787 X_colMap = getColumnMapMultiVector (X);
5788 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
5790 #ifdef HAVE_TPETRA_DEBUG
5791 auto X_colMap_host_view =
5792 X_colMap->template getLocalView<Kokkos::HostSpace> ();
5793 auto X_domainMap_host_view =
5794 X_domainMap->template getLocalView<Kokkos::HostSpace> ();
5796 if (X_colMap->getLocalLength () != 0 && X_domainMap->getLocalLength ()) {
5797 TEUCHOS_TEST_FOR_EXCEPTION
5798 (X_colMap_host_view.data () != X_domainMap_host_view.data (),
5799 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: Pointer to "
5800 "start of column Map view of X is not equal to pointer to start of "
5801 "(domain Map view of) X. This may mean that Tpetra::MultiVector::"
5802 "offsetViewNonConst is broken. "
5803 "Please report this bug to the Tpetra developers.");
5806 TEUCHOS_TEST_FOR_EXCEPTION(
5807 X_colMap_host_view.extent (0) < X_domainMap_host_view.extent (0) ||
5808 X_colMap->getLocalLength () < X_domainMap->getLocalLength (),
5809 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: "
5810 "X_colMap has fewer local rows than X_domainMap. "
5811 "X_colMap_host_view.extent(0) = " << X_colMap_host_view.extent (0)
5812 <<
", X_domainMap_host_view.extent(0) = "
5813 << X_domainMap_host_view.extent (0)
5814 <<
", X_colMap->getLocalLength() = " << X_colMap->getLocalLength ()
5815 <<
", and X_domainMap->getLocalLength() = "
5816 << X_domainMap->getLocalLength ()
5817 <<
". This means that Tpetra::MultiVector::offsetViewNonConst "
5818 "is broken. Please report this bug to the Tpetra developers.");
5820 TEUCHOS_TEST_FOR_EXCEPTION(
5821 X_colMap->getNumVectors () != X_domainMap->getNumVectors (),
5822 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: "
5823 "X_colMap has a different number of columns than X_domainMap. "
5824 "X_colMap->getNumVectors() = " << X_colMap->getNumVectors ()
5825 <<
" != X_domainMap->getNumVectors() = "
5826 << X_domainMap->getNumVectors ()
5827 <<
". This means that Tpetra::MultiVector::offsetViewNonConst "
5828 "is broken. Please report this bug to the Tpetra developers.");
5829 #endif // HAVE_TPETRA_DEBUG
5831 if (zeroInitialGuess) {
5833 X_colMap->putScalar (
ZERO);
5843 X_colMap->doImport (X, *importer,
INSERT);
5845 copyBackOutput =
true;
5853 B_in = rcpFromRef (B);
5859 RCP<MV> B_in_nonconst = getRowMapMultiVector (B,
true);
5862 }
catch (std::exception& e) {
5863 std::ostringstream os;
5864 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: "
5865 "deep_copy(*B_in_nonconst, B) threw an exception: "
5866 << e.what () <<
".";
5867 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
5869 B_in = rcp_const_cast<const MV> (B_in_nonconst);
5874 "gaussSeidelCopy: The current implementation requires that B have "
5875 "constant stride. Since B does not have constant stride, we had to "
5876 "copy it into a separate constant-stride multivector. This is a "
5877 "limitation of the current implementation and not your fault, but we "
5878 "still report it as an efficiency warning for your information.");
5881 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
5882 if (! importer.is_null () && sweep > 0) {
5885 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
5889 if (direction != Symmetric) {
5890 if (rowIndices.is_null ()) {
5891 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5896 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5903 if (rowIndices.is_null ()) {
5904 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5912 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
5917 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5921 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
5930 if (copyBackOutput) {
5933 }
catch (std::exception& e) {
5934 TEUCHOS_TEST_FOR_EXCEPTION(
5935 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) "
5936 "threw an exception: " << e.what ());
5941 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5943 Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node> >
5949 const char tfecfFuncName[] =
"convert: ";
5951 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5952 (! this->isFillComplete (), std::runtime_error,
"This matrix (the source "
5953 "of the conversion) is not fill complete. You must first call "
5954 "fillComplete() (possibly with the domain and range Map) without an "
5955 "intervening call to resumeFill(), before you may call this method.");
5956 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5957 (! this->isStaticGraph (), std::logic_error,
"This matrix (the source "
5958 "of the conversion) claims to be fill complete, but does not have a "
5959 "static (i.e., constant) graph. Please report this bug to the Tpetra "
5962 RCP<output_matrix_type> newMatrix
5963 (
new output_matrix_type (this->getCrsGraph ()));
5966 ::Tpetra::Details::copyConvert (newMatrix->lclMatrix_.values,
5967 this->lclMatrix_.values);
5971 newMatrix->fillComplete (this->getDomainMap (), this->getRangeMap ());
5977 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5982 #ifdef HAVE_TPETRA_DEBUG
5983 const char tfecfFuncName[] =
"checkInternalState: ";
5984 const char err[] =
"Internal state is not consistent. "
5985 "Please report this bug to the Tpetra developers.";
5989 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5990 staticGraph_.is_null (),
5991 std::logic_error, err);
5995 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5996 ! myGraph_.is_null () && myGraph_ != staticGraph_,
5997 std::logic_error, err);
5999 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6000 isFillComplete () && ! staticGraph_->isFillComplete (),
6001 std::logic_error, err <<
" Specifically, the matrix is fill complete, "
6002 "but its graph is NOT fill complete.");
6004 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6005 isStorageOptimized () && ! values2D_.is_null (),
6006 std::logic_error, err);
6008 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6009 getProfileType() ==
StaticProfile && values2D_ != Teuchos::null,
6010 std::logic_error, err);
6012 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6013 getProfileType() ==
DynamicProfile && k_values1D_.extent (0) > 0,
6014 std::logic_error, err);
6017 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6018 staticGraph_->indicesAreAllocated () &&
6019 staticGraph_->getNodeAllocationSize() > 0 &&
6020 staticGraph_->getNodeNumRows() > 0
6021 && values2D_.is_null () &&
6022 k_values1D_.extent (0) == 0,
6023 std::logic_error, err);
6025 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
6026 k_values1D_.extent (0) > 0 && values2D_ != Teuchos::null,
6027 std::logic_error, err <<
" Specifically, k_values1D_ is allocated (has "
6028 "size " << k_values1D_.extent (0) <<
" > 0) and values2D_ is also "
6029 "allocated. CrsMatrix is not suppose to have both a 1-D and a 2-D "
6030 "allocation at the same time.");
6034 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6039 std::ostringstream os;
6041 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
6042 if (this->getObjectLabel () !=
"") {
6043 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
6045 if (isFillComplete ()) {
6046 os <<
"isFillComplete: true"
6047 <<
", global dimensions: [" << getGlobalNumRows () <<
", "
6048 << getGlobalNumCols () <<
"]"
6049 <<
", global number of entries: " << getGlobalNumEntries ()
6053 os <<
"isFillComplete: false"
6054 <<
", global dimensions: [" << getGlobalNumRows () <<
", "
6055 << getGlobalNumCols () <<
"]}";
6060 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6064 const Teuchos::EVerbosityLevel verbLevel)
const
6068 using Teuchos::ArrayView;
6069 using Teuchos::Comm;
6071 using Teuchos::TypeNameTraits;
6072 using Teuchos::VERB_DEFAULT;
6073 using Teuchos::VERB_NONE;
6074 using Teuchos::VERB_LOW;
6075 using Teuchos::VERB_MEDIUM;
6076 using Teuchos::VERB_HIGH;
6077 using Teuchos::VERB_EXTREME;
6079 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
6081 if (vl == VERB_NONE) {
6085 Teuchos::OSTab tab0 (out);
6087 RCP<const Comm<int> > comm = this->getComm();
6088 const int myRank = comm->getRank();
6089 const int numProcs = comm->getSize();
6091 for (
size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
6094 width = std::max<size_t> (width, static_cast<size_t> (11)) + 2;
6104 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
6106 Teuchos::OSTab tab1 (out);
6109 if (this->getObjectLabel () !=
"") {
6110 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
6113 out <<
"Template parameters:" << endl;
6114 Teuchos::OSTab tab2 (out);
6115 out <<
"Scalar: " << TypeNameTraits<Scalar>::name () << endl
6116 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name () << endl
6117 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name () << endl
6118 <<
"Node: " << TypeNameTraits<Node>::name () << endl;
6120 if (isFillComplete()) {
6121 out <<
"isFillComplete: true" << endl
6122 <<
"Global dimensions: [" << getGlobalNumRows () <<
", "
6123 << getGlobalNumCols () <<
"]" << endl
6124 <<
"Global number of entries: " << getGlobalNumEntries () << endl
6125 << endl <<
"Global max number of entries in a row: "
6126 << getGlobalMaxNumRowEntries () << endl;
6129 out <<
"isFillComplete: false" << endl
6130 <<
"Global dimensions: [" << getGlobalNumRows () <<
", "
6131 << getGlobalNumCols () <<
"]" << endl;
6135 if (vl < VERB_MEDIUM) {
6141 out << endl <<
"Row Map:" << endl;
6143 if (getRowMap ().is_null ()) {
6145 out <<
"null" << endl;
6152 getRowMap ()->describe (out, vl);
6157 out <<
"Column Map: ";
6159 if (getColMap ().is_null ()) {
6161 out <<
"null" << endl;
6163 }
else if (getColMap () == getRowMap ()) {
6165 out <<
"same as row Map" << endl;
6171 getColMap ()->describe (out, vl);
6176 out <<
"Domain Map: ";
6178 if (getDomainMap ().is_null ()) {
6180 out <<
"null" << endl;
6182 }
else if (getDomainMap () == getRowMap ()) {
6184 out <<
"same as row Map" << endl;
6186 }
else if (getDomainMap () == getColMap ()) {
6188 out <<
"same as column Map" << endl;
6194 getDomainMap ()->describe (out, vl);
6199 out <<
"Range Map: ";
6201 if (getRangeMap ().is_null ()) {
6203 out <<
"null" << endl;
6205 }
else if (getRangeMap () == getDomainMap ()) {
6207 out <<
"same as domain Map" << endl;
6209 }
else if (getRangeMap () == getRowMap ()) {
6211 out <<
"same as row Map" << endl;
6217 getRangeMap ()->describe (out, vl);
6221 for (
int curRank = 0; curRank < numProcs; ++curRank) {
6222 if (myRank == curRank) {
6223 out <<
"Process rank: " << curRank << endl;
6224 Teuchos::OSTab tab2 (out);
6225 if (! staticGraph_->indicesAreAllocated ()) {
6226 out <<
"Graph indices not allocated" << endl;
6229 out <<
"Number of allocated entries: "
6230 << staticGraph_->getNodeAllocationSize () << endl;
6232 out <<
"Number of entries: " << getNodeNumEntries () << endl
6233 <<
"Max number of entries per row: " << getNodeMaxNumRowEntries ()
6242 if (vl < VERB_HIGH) {
6247 for (
int curRank = 0; curRank < numProcs; ++curRank) {
6248 if (myRank == curRank) {
6249 out << std::setw(width) <<
"Proc Rank"
6250 << std::setw(width) <<
"Global Row"
6251 << std::setw(width) <<
"Num Entries";
6252 if (vl == VERB_EXTREME) {
6253 out << std::setw(width) <<
"(Index,Value)";
6256 for (
size_t r = 0; r < getNodeNumRows (); ++r) {
6257 const size_t nE = getNumEntriesInLocalRow(r);
6258 GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
6259 out << std::setw(width) << myRank
6260 << std::setw(width) << gid
6261 << std::setw(width) << nE;
6262 if (vl == VERB_EXTREME) {
6263 if (isGloballyIndexed()) {
6264 ArrayView<const GlobalOrdinal> rowinds;
6265 ArrayView<const Scalar> rowvals;
6266 getGlobalRowView (gid, rowinds, rowvals);
6267 for (
size_t j = 0; j < nE; ++j) {
6268 out <<
" (" << rowinds[j]
6269 <<
", " << rowvals[j]
6273 else if (isLocallyIndexed()) {
6274 ArrayView<const LocalOrdinal> rowinds;
6275 ArrayView<const Scalar> rowvals;
6276 getLocalRowView (r, rowinds, rowvals);
6277 for (
size_t j=0; j < nE; ++j) {
6278 out <<
" (" << getColMap()->getGlobalElement(rowinds[j])
6279 <<
", " << rowvals[j]
6295 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6308 const row_matrix_type* srcRowMat =
6309 dynamic_cast<const row_matrix_type*> (&source);
6310 return (srcRowMat != NULL);
6313 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6321 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6325 const size_t numSameIDs,
6326 const LocalOrdinal permuteToLIDs[],
6327 const LocalOrdinal permuteFromLIDs[],
6328 const size_t numPermutes)
6331 using Teuchos::Array;
6332 using Teuchos::ArrayView;
6333 typedef LocalOrdinal LO;
6334 typedef GlobalOrdinal GO;
6335 #ifdef HAVE_TPETRA_DEBUG
6337 const char tfecfFuncName[] =
"copyAndPermuteImpl: ";
6338 #endif // HAVE_TPETRA_DEBUG
6339 ProfilingRegion regionCAP (
"Tpetra::CrsMatrix::copyAndPermuteImpl");
6346 const map_type& srcRowMap = * (srcMat.
getRowMap ());
6348 Array<Scalar> rowVals;
6349 const LO numSameIDs_as_LID = static_cast<LO> (numSameIDs);
6350 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
6354 const GO sourceGID = srcRowMap.getGlobalElement (sourceLID);
6355 const GO targetGID = sourceGID;
6358 ArrayView<const GO> rowIndsConstView;
6359 ArrayView<const Scalar> rowValsConstView;
6361 if (sourceIsLocallyIndexed) {
6363 if (rowLength > static_cast<size_t> (rowInds.size())) {
6364 rowInds.resize (rowLength);
6365 rowVals.resize (rowLength);
6369 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
6370 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
6375 size_t checkRowLength = 0;
6376 srcMat.
getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
6378 #ifdef HAVE_TPETRA_DEBUG
6379 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
6380 std::logic_error,
"For global row index " << sourceGID <<
", the source"
6381 " matrix's getNumEntriesInGlobalRow() method returns a row length of "
6382 << rowLength <<
", but the getGlobalRowCopy() method reports that "
6383 "the row length is " << checkRowLength <<
". Please report this bug "
6384 "to the Tpetra developers.");
6385 #endif // HAVE_TPETRA_DEBUG
6387 rowIndsConstView = rowIndsView.view (0, rowLength);
6388 rowValsConstView = rowValsView.view (0, rowLength);
6395 if (this->isStaticGraph ()) {
6398 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
REPLACE);
6404 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
INSERT);
6411 const map_type& tgtRowMap = * (this->getRowMap ());
6412 for (
size_t p = 0; p < numPermutes; ++p) {
6413 const GO sourceGID = srcRowMap.getGlobalElement (permuteFromLIDs[p]);
6414 const GO targetGID = tgtRowMap.getGlobalElement (permuteToLIDs[p]);
6417 ArrayView<const GO> rowIndsConstView;
6418 ArrayView<const Scalar> rowValsConstView;
6420 if (sourceIsLocallyIndexed) {
6422 if (rowLength > static_cast<size_t> (rowInds.size ())) {
6423 rowInds.resize (rowLength);
6424 rowVals.resize (rowLength);
6428 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
6429 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
6434 size_t checkRowLength = 0;
6435 srcMat.
getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
6437 #ifdef HAVE_TPETRA_DEBUG
6438 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
6439 std::logic_error,
"For the source matrix's global row index "
6440 << sourceGID <<
", the source matrix's getNumEntriesInGlobalRow() "
6441 "method returns a row length of " << rowLength <<
", but the "
6442 "getGlobalRowCopy() method reports that the row length is "
6443 << checkRowLength <<
". Please report this bug to the Tpetra "
6445 #endif // HAVE_TPETRA_DEBUG
6447 rowIndsConstView = rowIndsView.view (0, rowLength);
6448 rowValsConstView = rowValsView.view (0, rowLength);
6455 if (isStaticGraph()) {
6456 this->combineGlobalValues (targetGID, rowIndsConstView,
6460 this->combineGlobalValues (targetGID, rowIndsConstView,
6461 rowValsConstView,
INSERT);
6466 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6468 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6469 copyAndPermuteNew (
const SrcDistObject& srcObj,
6470 const size_t numSameIDs,
6471 const Kokkos::DualView<const local_ordinal_type*, device_type>& permuteToLIDs,
6472 const Kokkos::DualView<const local_ordinal_type*, device_type>& permuteFromLIDs)
6478 typedef Kokkos::HostSpace host_mem_space;
6479 typedef typename device_type::memory_space dev_mem_space;
6481 const char tfecfFuncName[] =
"copyAndPermuteNew: ";
6482 ProfilingRegion regionCAP (
"Tpetra::CrsMatrix::copyAndPermuteNew");
6488 std::unique_ptr<std::string> prefix;
6491 auto map = this->getMap ();
6492 if (! map.is_null ()) {
6493 auto comm = map->getComm ();
6494 if (! comm.is_null ()) {
6495 myRank = comm->getRank ();
6500 prefix = [myRank] () {
6501 std::ostringstream pfxStrm;
6502 pfxStrm <<
"(Proc " << myRank <<
") ";
6503 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
6505 std::ostringstream os;
6506 os << *prefix <<
"Tpetra::CrsMatrix::copyAndPermuteNew: " << endl
6511 std::cerr << os.str ();
6514 const auto numPermute = permuteToLIDs.extent (0);
6515 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6516 (numPermute != permuteFromLIDs.extent (0),
6517 std::invalid_argument,
"permuteToLIDs.extent(0) = "
6518 << numPermute <<
"!= permuteFromLIDs.extent(0) = "
6519 << permuteFromLIDs.extent (0) <<
".");
6525 const bool permuteToLIDs_sync_back =
6526 permuteToLIDs.modified_device () >= permuteToLIDs.modified_host ();
6528 permuteToLIDs_nc.template sync<host_mem_space> ();
6529 auto permuteToLIDs_h = permuteToLIDs.template view<host_mem_space> ();
6531 const bool permuteFromLIDs_sync_back =
6532 permuteFromLIDs.modified_device () >= permuteFromLIDs.modified_host ();
6534 permuteFromLIDs_nc.template sync<host_mem_space> ();
6535 auto permuteFromLIDs_h = permuteFromLIDs.template view<host_mem_space> ();
6538 std::ostringstream os;
6539 os << *prefix <<
"permuteToLIDs_sync_back: "
6540 << (permuteToLIDs_sync_back ?
"true" :
"false") <<
", "
6541 <<
"permuteFromLIDs_sync_back: "
6542 << (permuteFromLIDs_sync_back ?
"true" :
"false") << endl;
6543 std::cerr << os.str ();
6548 typedef ::Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> RMT;
6549 const RMT& srcMat = dynamic_cast<const RMT&> (srcObj);
6551 this->copyAndPermuteImpl (srcMat, numSameIDs, permuteToLIDs_h.data (),
6552 permuteFromLIDs_h.data (), numPermute);
6554 if (permuteToLIDs_sync_back) {
6555 permuteToLIDs_nc.template sync<dev_mem_space> ();
6557 if (permuteFromLIDs_sync_back) {
6558 permuteFromLIDs_nc.template sync<dev_mem_space> ();
6562 std::ostringstream os;
6563 os << *prefix <<
"copyAndPermuteNew: after:" << endl
6568 std::cerr << os.str ();
6572 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6574 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6575 packAndPrepareNew (
const SrcDistObject& source,
6576 const Kokkos::DualView<const local_ordinal_type*, device_type>& exportLIDs,
6577 Kokkos::DualView<char*, buffer_device_type>& exports,
6578 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
6579 size_t& constantNumPackets,
6580 Distributor& distor)
6584 using Teuchos::outArg;
6585 using Teuchos::REDUCE_MAX;
6586 using Teuchos::reduceAll;
6588 typedef LocalOrdinal LO;
6589 typedef GlobalOrdinal GO;
6590 const char tfecfFuncName[] =
"packAndPrepareNew: ";
6591 ProfilingRegion regionPAP (
"Tpetra::CrsMatrix::packAndPrepareNew");
6597 Teuchos::RCP<const Teuchos::Comm<int> > pComm = this->getComm ();
6598 if (pComm.is_null ()) {
6601 const Teuchos::Comm<int>& comm = *pComm;
6602 const int myRank = comm.getSize ();
6604 std::unique_ptr<std::string> prefix;
6607 prefix = [myRank] () {
6608 std::ostringstream pfxStrm;
6609 pfxStrm <<
"(Proc " << myRank <<
") ";
6610 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
6613 std::ostringstream os;
6614 os << *prefix <<
"Tpetra::CrsMatrix::packAndPrepareNew: " << endl
6624 std::cerr << os.str ();
6647 std::ostringstream msg;
6650 typedef CrsMatrix<Scalar, LO, GO, Node> crs_matrix_type;
6651 const crs_matrix_type* srcCrsMat =
6652 dynamic_cast<const crs_matrix_type*> (&source);
6653 if (srcCrsMat != NULL) {
6655 std::ostringstream os;
6656 os << *prefix <<
"Source matrix same (CrsMatrix) type as target; "
6657 "calling packNew" << endl;
6658 std::cerr << os.str ();
6661 srcCrsMat->packNew (exportLIDs, exports, numPacketsPerLID,
6662 constantNumPackets, distor);
6664 catch (std::exception& e) {
6666 msg <<
"Proc " << myRank <<
": " << e.what () << std::endl;
6671 using Kokkos::HostSpace;
6672 using Kokkos::subview;
6673 typedef Kokkos::DualView<char*, buffer_device_type> exports_type;
6674 typedef Kokkos::pair<size_t, size_t> range_type;
6677 std::ostringstream os;
6678 os << *prefix <<
"Source matrix NOT same (CrsMatrix) type as target"
6680 std::cerr << os.str ();
6683 typedef RowMatrix<Scalar, LO, GO, Node> row_matrix_type;
6684 const row_matrix_type* srcRowMat =
6685 dynamic_cast<const row_matrix_type*> (&source);
6686 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6687 (srcRowMat == NULL, std::invalid_argument,
6688 "The source object of the Import or Export operation is neither a "
6689 "CrsMatrix (with the same template parameters as the target object), "
6690 "nor a RowMatrix (with the same first four template parameters as the "
6703 exportLIDs_nc.template sync<HostSpace> ();
6705 auto exportLIDs_h = exportLIDs.template view<HostSpace> ();
6706 Teuchos::ArrayView<const LO> exportLIDs_av (exportLIDs_h.data (),
6707 exportLIDs_h.size ());
6711 Teuchos::Array<char> exports_a;
6717 auto numPacketsPerLID_nc = numPacketsPerLID;
6718 numPacketsPerLID_nc.modified_device() = 0;
6719 numPacketsPerLID_nc.modified_host() = 1;
6721 auto numPacketsPerLID_h = numPacketsPerLID.template view<HostSpace> ();
6722 Teuchos::ArrayView<size_t> numPacketsPerLID_av (numPacketsPerLID_h.data (),
6723 numPacketsPerLID_h.size ());
6728 srcRowMat->pack (exportLIDs_av, exports_a, numPacketsPerLID_av,
6729 constantNumPackets, distor);
6731 catch (std::exception& e) {
6733 msg <<
"Proc " << myRank <<
": " << e.what () << std::endl;
6737 const size_t newAllocSize = static_cast<size_t> (exports_a.size ());
6738 if (static_cast<size_t> (exports.extent (0)) < newAllocSize) {
6739 const std::string oldLabel = exports.d_view.label ();
6740 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6741 exports = exports_type (newLabel, newAllocSize);
6745 exports.modified_device() = 0;
6746 exports.modified_host() = 1;
6748 auto exports_h = exports.template view<HostSpace> ();
6749 auto exports_h_sub = subview (exports_h, range_type (0, newAllocSize));
6753 typedef typename exports_type::t_host::execution_space HES;
6754 typedef Kokkos::Device<HES, HostSpace> host_device_type;
6755 Kokkos::View<const char*, host_device_type>
6756 exports_a_kv (exports_a.getRawPtr (), newAllocSize);
6762 reduceAll<int, int> (comm, REDUCE_MAX, lclBad, outArg (gblBad));
6765 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6766 (
true, std::logic_error,
"packNew() or pack() threw an exception on "
6767 "one or more participating processes.");
6771 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6772 (lclBad != 0, std::logic_error,
"packNew threw an exception on one "
6773 "or more participating processes. Here is this process' error "
6774 "message: " << msg.str ());
6778 std::ostringstream os;
6779 os << *prefix <<
"packAndPrepareNew: Done!" << endl
6789 std::cerr << os.str ();
6793 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6795 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6796 packRow (
char exports[],
6797 const size_t offset,
6798 const size_t numEnt,
6799 const GlobalOrdinal gidsIn[],
6800 const impl_scalar_type valsIn[],
6801 const size_t numBytesPerValue)
const
6804 using Kokkos::subview;
6806 typedef LocalOrdinal LO;
6807 typedef GlobalOrdinal GO;
6808 typedef impl_scalar_type ST;
6809 typedef typename View<int*, device_type>::HostMirror::execution_space HES;
6817 const LO numEntLO = static_cast<size_t> (numEnt);
6819 const size_t numEntBeg = offset;
6820 const size_t numEntLen = PackTraits<LO, HES>::packValueCount (numEntLO);
6821 const size_t gidsBeg = numEntBeg + numEntLen;
6822 const size_t gidsLen = numEnt * PackTraits<GO, HES>::packValueCount (gid);
6823 const size_t valsBeg = gidsBeg + gidsLen;
6824 const size_t valsLen = numEnt * numBytesPerValue;
6826 char*
const numEntOut = exports + numEntBeg;
6827 char*
const gidsOut = exports + gidsBeg;
6828 char*
const valsOut = exports + valsBeg;
6830 size_t numBytesOut = 0;
6832 numBytesOut += PackTraits<LO, HES>::packValue (numEntOut, numEntLO);
6835 Kokkos::pair<int, size_t> p;
6836 p = PackTraits<GO, HES>::packArray (gidsOut, gidsIn, numEnt);
6837 errorCode += p.first;
6838 numBytesOut += p.second;
6840 p = PackTraits<ST, HES>::packArray (valsOut, valsIn, numEnt);
6841 errorCode += p.first;
6842 numBytesOut += p.second;
6845 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6846 TEUCHOS_TEST_FOR_EXCEPTION
6847 (numBytesOut != expectedNumBytes, std::logic_error,
"packRow: "
6848 "numBytesOut = " << numBytesOut <<
" != expectedNumBytes = "
6849 << expectedNumBytes <<
".");
6850 TEUCHOS_TEST_FOR_EXCEPTION
6851 (errorCode != 0, std::runtime_error,
"packRow: "
6852 "PackTraits::packArray returned a nonzero error code");
6857 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6859 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6860 unpackRow (GlobalOrdinal gidsOut[],
6861 impl_scalar_type valsOut[],
6862 const char imports[],
6863 const size_t offset,
6864 const size_t numBytes,
6865 const size_t numEnt,
6866 const size_t numBytesPerValue)
6869 using Kokkos::subview;
6871 typedef LocalOrdinal LO;
6872 typedef GlobalOrdinal GO;
6873 typedef impl_scalar_type ST;
6874 typedef typename View<int*, device_type>::HostMirror::execution_space HES;
6876 if (numBytes == 0) {
6879 const int myRank = this->getMap ()->getComm ()->getRank ();
6880 TEUCHOS_TEST_FOR_EXCEPTION
6881 (
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6882 "unpackRow: The number of bytes to unpack numBytes=0, but the "
6883 "number of entries to unpack (as reported by numPacketsPerLID) "
6884 "for this row numEnt=" << numEnt <<
" != 0.");
6889 if (numEnt == 0 && numBytes != 0) {
6890 const int myRank = this->getMap ()->getComm ()->getRank ();
6891 TEUCHOS_TEST_FOR_EXCEPTION
6892 (
true, std::logic_error,
"(Proc " << myRank <<
") CrsMatrix::"
6893 "unpackRow: The number of entries to unpack (as reported by "
6894 "numPacketsPerLID) numEnt=0, but the number of bytes to unpack "
6895 "numBytes=" << numBytes <<
" != 0.");
6901 const size_t numEntBeg = offset;
6902 const size_t numEntLen = PackTraits<LO, HES>::packValueCount (lid);
6903 const size_t gidsBeg = numEntBeg + numEntLen;
6904 const size_t gidsLen = numEnt * PackTraits<GO, HES>::packValueCount (gid);
6905 const size_t valsBeg = gidsBeg + gidsLen;
6906 const size_t valsLen = numEnt * numBytesPerValue;
6908 const char*
const numEntIn = imports + numEntBeg;
6909 const char*
const gidsIn = imports + gidsBeg;
6910 const char*
const valsIn = imports + valsBeg;
6912 size_t numBytesOut = 0;
6915 numBytesOut += PackTraits<LO, HES>::unpackValue (numEntOut, numEntIn);
6916 if (static_cast<size_t> (numEntOut) != numEnt ||
6917 numEntOut == static_cast<LO> (0)) {
6918 const int myRank = this->getMap ()->getComm ()->getRank ();
6919 std::ostringstream os;
6920 os <<
"(Proc " << myRank <<
") CrsMatrix::unpackRow: ";
6921 bool firstErrorCondition =
false;
6922 if (static_cast<size_t> (numEntOut) != numEnt) {
6923 os <<
"Number of entries from numPacketsPerLID numEnt=" << numEnt
6924 <<
" does not equal number of entries unpacked from imports "
6925 "buffer numEntOut=" << numEntOut <<
".";
6926 firstErrorCondition =
true;
6928 if (numEntOut == static_cast<LO> (0)) {
6929 if (firstErrorCondition) {
6932 os <<
"Number of entries unpacked from imports buffer numEntOut=0, "
6933 "but number of bytes to unpack for this row numBytes=" << numBytes
6934 <<
" != 0. This should never happen, since packRow should only "
6935 "ever pack rows with a nonzero number of entries. In this case, "
6936 "the number of entries from numPacketsPerLID is numEnt=" << numEnt
6939 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error, os.str ());
6943 Kokkos::pair<int, size_t> p;
6944 p = PackTraits<GO, HES>::unpackArray (gidsOut, gidsIn, numEnt);
6945 errorCode += p.first;
6946 numBytesOut += p.second;
6948 p = PackTraits<ST, HES>::unpackArray (valsOut, valsIn, numEnt);
6949 errorCode += p.first;
6950 numBytesOut += p.second;
6953 TEUCHOS_TEST_FOR_EXCEPTION
6954 (numBytesOut != numBytes, std::logic_error,
"unpackRow: numBytesOut = "
6955 << numBytesOut <<
" != numBytes = " << numBytes <<
".");
6957 const size_t expectedNumBytes = numEntLen + gidsLen + valsLen;
6958 TEUCHOS_TEST_FOR_EXCEPTION
6959 (numBytesOut != expectedNumBytes, std::logic_error,
"unpackRow: "
6960 "numBytesOut = " << numBytesOut <<
" != expectedNumBytes = "
6961 << expectedNumBytes <<
".");
6963 TEUCHOS_TEST_FOR_EXCEPTION
6964 (errorCode != 0, std::runtime_error,
"unpackRow: "
6965 "PackTraits::unpackArray returned a nonzero error code");
6970 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6972 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6973 allocatePackSpaceNew (Kokkos::DualView<char*, buffer_device_type>& exports,
6974 size_t& totalNumEntries,
6975 const Kokkos::DualView<const local_ordinal_type*, device_type>& exportLIDs)
const
6979 typedef impl_scalar_type IST;
6980 typedef LocalOrdinal LO;
6981 typedef GlobalOrdinal GO;
6988 std::unique_ptr<std::string> prefix;
6991 auto map = this->getMap ();
6992 if (! map.is_null ()) {
6993 auto comm = map->getComm ();
6994 if (! comm.is_null ()) {
6995 myRank = comm->getRank ();
6999 prefix = [myRank] () {
7000 std::ostringstream pfxStrm;
7001 pfxStrm <<
"(Proc " << myRank <<
") ";
7002 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
7005 std::ostringstream os;
7006 os << *prefix <<
"Tpetra::CrsMatrix::allocatePackSpaceNew: Before:"
7014 std::cerr << os.str ();
7019 const LO numExportLIDs = static_cast<LO> (exportLIDs.extent (0));
7025 Kokkos::DualView<local_ordinal_type*, device_type> exportLIDs_nc =
7027 exportLIDs_nc.template sync<Kokkos::HostSpace> ();
7029 auto exportLIDs_h = exportLIDs.template view<Kokkos::HostSpace> ();
7032 totalNumEntries = 0;
7033 for (LO i = 0; i < numExportLIDs; ++i) {
7034 const LO lclRow = exportLIDs_h[i];
7035 size_t curNumEntries = this->getNumEntriesInLocalRow (lclRow);
7038 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
7041 totalNumEntries += curNumEntries;
7052 const size_t allocSize =
7053 static_cast<size_t> (numExportLIDs) *
sizeof (LO) +
7054 totalNumEntries * (
sizeof (IST) +
sizeof (GO));
7055 if (static_cast<size_t> (exports.extent (0)) < allocSize) {
7056 typedef Kokkos::DualView<char*, buffer_device_type> exports_type;
7058 const std::string oldLabel = exports.d_view.label ();
7059 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
7060 exports = exports_type (newLabel, allocSize);
7064 std::ostringstream os;
7065 os << *prefix <<
"Tpetra::CrsMatrix::allocatePackSpaceNew: After:"
7073 std::cerr << os.str ();
7077 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7079 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7080 packNew (
const Kokkos::DualView<const local_ordinal_type*, device_type>& exportLIDs,
7081 Kokkos::DualView<char*, buffer_device_type>& exports,
7082 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
7083 size_t& constantNumPackets,
7087 if (this->isStaticGraph ()) {
7088 using ::Tpetra::Details::packCrsMatrixNew;
7089 packCrsMatrixNew (*
this, exports, numPacketsPerLID, exportLIDs,
7090 constantNumPackets, dist);
7093 this->packNonStaticNew (exportLIDs, exports, numPacketsPerLID,
7094 constantNumPackets, dist);
7098 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7101 packNonStaticNew (
const Kokkos::DualView<const local_ordinal_type*, device_type>& exportLIDs,
7102 Kokkos::DualView<char*, buffer_device_type>& exports,
7103 const Kokkos::DualView<size_t*, buffer_device_type>& numPacketsPerLID,
7104 size_t& constantNumPackets,
7112 typedef LocalOrdinal LO;
7113 typedef GlobalOrdinal GO;
7114 typedef impl_scalar_type ST;
7115 typedef typename View<int*, device_type>::HostMirror::execution_space HES;
7116 const char tfecfFuncName[] =
"packNonStaticNew: ";
7122 std::unique_ptr<std::string> prefix;
7125 auto map = this->getMap ();
7126 if (! map.is_null ()) {
7127 auto comm = map->getComm ();
7128 if (! comm.is_null ()) {
7129 myRank = comm->getRank ();
7133 prefix = [myRank] () {
7134 std::ostringstream pfxStrm;
7135 pfxStrm <<
"(Proc " << myRank <<
") ";
7136 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
7139 std::ostringstream os;
7140 os << *prefix <<
"Tpetra::CrsMatrix::packNonStaticNew:" << endl;
7141 std::cerr << os.str ();
7144 const size_t numExportLIDs = static_cast<size_t> (exportLIDs.extent (0));
7145 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7146 (numExportLIDs != static_cast<size_t> (numPacketsPerLID.extent (0)),
7147 std::invalid_argument,
"exportLIDs.size() = " << numExportLIDs
7148 <<
" != numPacketsPerLID.size() = " << numPacketsPerLID.extent (0)
7154 constantNumPackets = 0;
7159 size_t totalNumEntries = 0;
7160 this->allocatePackSpaceNew (exports, totalNumEntries, exportLIDs);
7161 const size_t bufSize = static_cast<size_t> (exports.extent (0));
7164 exports.modified_device() = 0;
7165 exports.modified_host() = 1;
7166 auto exports_h = exports.template view<Kokkos::HostSpace> ();
7168 std::ostringstream os;
7169 os << *prefix <<
"After marking exports as modified on host, "
7171 std::cerr << os.str ();
7175 auto exportLIDs_h = exportLIDs.template view<Kokkos::HostSpace> ();
7178 numPacketsPerLID.modified_device() = 0;
7179 numPacketsPerLID.modified_host() = 1;
7180 auto numPacketsPerLID_h = numPacketsPerLID.template view<Kokkos::HostSpace> ();
7186 for (
size_t i = 0; i < numExportLIDs; ++i) {
7187 const LO lclRow = exportLIDs_h[i];
7190 numEnt = this->getNumEntriesInLocalRow (lclRow);
7197 numPacketsPerLID_h[i] = 0;
7202 View<GO*, HES> gidsIn_k;
7205 gidsIn_k = PackTraits<GO, HES>::allocateArray(gid, numEnt,
"gids");
7208 Teuchos::ArrayView<const Scalar> valsIn;
7209 if (this->isLocallyIndexed ()) {
7213 Teuchos::ArrayView<const LO> lidsIn;
7214 this->getLocalRowView (lclRow, lidsIn, valsIn);
7215 const map_type& colMap = * (this->getColMap ());
7216 for (
size_t k = 0; k < numEnt; ++k) {
7217 gidsIn_k[k] = colMap.getGlobalElement (lidsIn[k]);
7220 else if (this->isGloballyIndexed ()) {
7226 Teuchos::ArrayView<const GO> gblIndView;;
7227 const map_type& rowMap = * (this->getRowMap ());
7228 const GO gblRow = rowMap.getGlobalElement (lclRow);
7229 this->getGlobalRowView (gblRow, gblIndView, valsIn);
7230 for (
size_t k = 0; k < numEnt; ++k) {
7231 gidsIn_k[k] = gblIndView[k];
7238 typename HES::device_type outputDevice;
7241 reinterpret_cast<const ST*> (valsIn.getRawPtr ()),
7244 const size_t numBytesPerValue =
7245 PackTraits<ST,HES>::packValueCount (valsIn[0]);
7246 const size_t numBytes =
7247 this->packRow (exports_h.data (), offset, numEnt, gidsIn_k.data (),
7248 valsIn_k.data (), numBytesPerValue);
7249 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7250 (offset > bufSize || offset + numBytes > bufSize, std::logic_error,
7251 "First invalid offset into 'exports' pack buffer at index i = " << i
7252 <<
". exportLIDs_h[i]: " << exportLIDs_h[i] <<
", bufSize: " <<
7253 bufSize <<
", offset: " << offset <<
", numBytes: " << numBytes <<
7258 numPacketsPerLID_h[i] = numBytes;
7263 std::ostringstream os;
7264 os << *prefix <<
"Tpetra::CrsMatrix::packNonStaticNew: After:" << endl
7271 std::cerr << os.str ();
7275 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7277 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7278 combineGlobalValuesRaw (
const LocalOrdinal lclRow,
7279 const LocalOrdinal numEnt,
7280 const impl_scalar_type vals[],
7281 const GlobalOrdinal cols[],
7284 typedef GlobalOrdinal GO;
7289 const GO gblRow = this->myGraph_->rowMap_->getGlobalElement (lclRow);
7290 Teuchos::ArrayView<const GO> cols_av (numEnt == 0 ? NULL : cols, numEnt);
7291 Teuchos::ArrayView<const Scalar> vals_av (numEnt == 0 ? NULL : reinterpret_cast<const Scalar*> (vals), numEnt);
7296 this->combineGlobalValues (gblRow, cols_av, vals_av, combineMode);
7300 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7302 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7303 combineGlobalValues (
const GlobalOrdinal globalRowIndex,
7304 const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
7305 const Teuchos::ArrayView<const Scalar>& values,
7308 const char tfecfFuncName[] =
"combineGlobalValues: ";
7310 if (isStaticGraph ()) {
7314 if (combineMode ==
ADD) {
7315 sumIntoGlobalValues (globalRowIndex, columnIndices, values);
7317 else if (combineMode ==
REPLACE) {
7318 replaceGlobalValues (globalRowIndex, columnIndices, values);
7320 else if (combineMode ==
ABSMAX) {
7321 using ::Tpetra::Details::AbsMax;
7323 this->
template transformGlobalValues<AbsMax<Scalar> > (globalRowIndex,
7327 else if (combineMode ==
INSERT) {
7328 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
7329 isStaticGraph () && combineMode ==
INSERT, std::invalid_argument,
7330 "INSERT combine mode is not allowed if the matrix has a static graph "
7331 "(i.e., was constructed with the CrsMatrix constructor that takes a "
7332 "const CrsGraph pointer).");
7335 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
7336 true, std::logic_error,
"Invalid combine mode; should never get "
7337 "here! Please report this bug to the Tpetra developers.");
7341 if (combineMode ==
ADD || combineMode ==
INSERT) {
7349 this->insertGlobalValuesFiltered (globalRowIndex, columnIndices,
7352 catch (std::exception& e) {
7353 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7354 (
true, std::runtime_error, std::endl
7355 <<
"insertGlobalValuesFiltered(" << globalRowIndex <<
", "
7356 << std::endl << Teuchos::toString (columnIndices) <<
", "
7357 << std::endl << Teuchos::toString (values)
7358 <<
") threw an exception: " << std::endl << e.what ());
7370 else if (combineMode ==
ABSMAX) {
7371 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
7372 ! isStaticGraph () && combineMode ==
ABSMAX, std::logic_error,
7373 "ABSMAX combine mode when the matrix has a dynamic graph is not yet "
7376 else if (combineMode ==
REPLACE) {
7377 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
7378 ! isStaticGraph () && combineMode ==
REPLACE, std::logic_error,
7379 "REPLACE combine mode when the matrix has a dynamic graph is not yet "
7383 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
7384 true, std::logic_error,
"Should never get here! Please report this "
7385 "bug to the Tpetra developers.");
7390 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7392 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7393 unpackAndCombineNew (
const Kokkos::DualView<const local_ordinal_type*, device_type>& importLIDs,
7394 const Kokkos::DualView<const char*, buffer_device_type>& imports,
7395 const Kokkos::DualView<const size_t*, buffer_device_type>& numPacketsPerLID,
7396 const size_t constantNumPackets,
7403 const char tfecfFuncName[] =
"unpackAndCombineNew: ";
7404 ProfilingRegion regionUAC (
"Tpetra::CrsMatrix::unpackAndCombineNew");
7408 constexpr
int numValidModes = 5;
7411 const char* validModeNames[numValidModes] =
7412 {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT",
"ZERO"};
7414 std::unique_ptr<std::string> prefix;
7417 auto map = this->getMap ();
7418 if (! map.is_null ()) {
7419 auto comm = map->getComm ();
7420 if (! comm.is_null ()) {
7421 myRank = comm->getRank ();
7424 prefix = [myRank] () {
7425 std::ostringstream pfxStrm;
7426 pfxStrm <<
"(Proc " << myRank <<
") ";
7427 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
7430 std::ostringstream os;
7431 os << *prefix <<
"Tpetra::CrsMatrix::unpackAndCombineNew: " << endl
7441 << *prefix <<
" constantNumPackets: " << constantNumPackets
7445 std::cerr << os.str ();
7449 if (std::find (validModes, validModes+numValidModes, combineMode) ==
7450 validModes+numValidModes) {
7451 std::ostringstream os;
7452 os <<
"Invalid combine mode. Valid modes are {";
7453 for (
int k = 0; k < numValidModes; ++k) {
7454 os << validModeNames[k];
7455 if (k < numValidModes - 1) {
7460 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7461 (
true, std::invalid_argument, os.str ());
7465 if (combineMode ==
ZERO) {
7470 using Teuchos::reduceAll;
7471 std::unique_ptr<std::ostringstream> msg (
new std::ostringstream ());
7474 this->unpackAndCombineNewImpl (importLIDs, imports, numPacketsPerLID,
7475 constantNumPackets, distor, combineMode);
7476 }
catch (std::exception& e) {
7481 const Teuchos::Comm<int>& comm = * (this->getComm ());
7482 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
7483 lclBad, Teuchos::outArg (gblBad));
7489 std::ostringstream os;
7490 os <<
"(Proc " << comm.getRank () <<
") " << msg->str () << endl;
7491 msg = std::unique_ptr<std::ostringstream> (
new std::ostringstream ());
7492 ::Tpetra::Details::gathervPrint (*msg, os.str (), comm);
7493 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7494 (
true, std::logic_error, std::endl <<
"unpackAndCombineNewImpl() "
7495 "threw an exception on one or more participating processes: "
7496 << endl << msg->str ());
7500 this->unpackAndCombineNewImpl (importLIDs, imports, numPacketsPerLID,
7501 constantNumPackets, distor, combineMode);
7505 std::ostringstream os;
7506 os << *prefix <<
"unpackAndCombineNew: Done!" << endl
7516 std::cerr << os.str ();
7520 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7524 const Kokkos::DualView<const char*, buffer_device_type>& imports,
7525 const Kokkos::DualView<const size_t*, buffer_device_type>& numPacketsPerLID,
7526 const size_t constantNumPackets,
7533 if (this->isStaticGraph ()) {
7534 using ::Tpetra::Details::unpackCrsMatrixAndCombineNew;
7535 unpackCrsMatrixAndCombineNew (*
this, imports, numPacketsPerLID,
7536 importLIDs, constantNumPackets,
7537 distor, combineMode, atomic);
7540 this->unpackAndCombineNewImplNonStatic (importLIDs, imports,
7543 distor, combineMode);
7547 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7549 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7550 unpackAndCombineNewImplNonStatic (
const Kokkos::DualView<const LocalOrdinal*, device_type>& importLIDs,
7551 const Kokkos::DualView<const char*, buffer_device_type>& imports,
7552 const Kokkos::DualView<const size_t*, buffer_device_type>& numPacketsPerLID,
7553 const size_t constantNumPackets,
7554 Distributor& distor,
7558 using Kokkos::subview;
7559 using Kokkos::MemoryUnmanaged;
7564 typedef LocalOrdinal LO;
7565 typedef GlobalOrdinal GO;
7566 typedef impl_scalar_type ST;
7567 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
7568 typedef typename View<int*, device_type>::HostMirror::execution_space HES;
7569 typedef std::pair<typename View<int*, HES>::size_type,
7570 typename View<int*, HES>::size_type> pair_type;
7571 typedef View<GO*, HES, MemoryUnmanaged> gids_out_type;
7572 typedef View<ST*, HES, MemoryUnmanaged> vals_out_type;
7573 const char tfecfFuncName[] =
"unpackAndCombineNewImplNonStatic: ";
7579 std::unique_ptr<std::string> prefix;
7582 auto map = this->getMap ();
7583 if (! map.is_null ()) {
7584 auto comm = map->getComm ();
7585 if (! comm.is_null ()) {
7586 myRank = comm->getRank ();
7590 prefix = [myRank] () {
7591 std::ostringstream pfxStrm;
7592 pfxStrm <<
"(Proc " << myRank <<
") ";
7593 return std::unique_ptr<std::string> (
new std::string (pfxStrm.str ()));
7596 std::ostringstream os;
7597 os << *prefix <<
"Tpetra::CrsMatrix::unpackAndCombineNewImplNonStatic:"
7599 std::cerr << os.str ();
7602 const size_type numImportLIDs = importLIDs.extent (0);
7603 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7604 (numImportLIDs != static_cast<size_type> (numPacketsPerLID.extent (0)),
7605 std::invalid_argument,
"importLIDs.size() = " << numImportLIDs
7606 <<
" != numPacketsPerLID.size() = " << numPacketsPerLID.extent (0)
7608 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7610 combineMode !=
ABSMAX && combineMode !=
ZERO, std::invalid_argument,
7611 "Invalid CombineMode value " << combineMode <<
". Valid "
7612 <<
"values include ADD, INSERT, REPLACE, ABSMAX, and ZERO.");
7613 if (combineMode ==
ZERO || numImportLIDs == 0) {
7620 imports_nc.template sync<Kokkos::HostSpace> ();
7622 auto imports_h = imports.template view<Kokkos::HostSpace> ();
7627 numPacketsPerLID_nc.template sync<Kokkos::HostSpace> ();
7629 auto numPacketsPerLID_h = numPacketsPerLID.template view<Kokkos::HostSpace> ();
7634 importLIDs_nc.template sync<Kokkos::HostSpace> ();
7636 auto importLIDs_h = importLIDs.template view<Kokkos::HostSpace> ();
7638 size_t numBytesPerValue;
7649 numBytesPerValue = PackTraits<ST, HES>::packValueCount (val);
7654 size_t maxRowNumEnt = 0;
7655 for (size_type i = 0; i < numImportLIDs; ++i) {
7656 const size_t numBytes = numPacketsPerLID_h[i];
7657 if (numBytes == 0) {
7661 #ifdef HAVE_TPETRA_DEBUG
7662 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7663 (offset + numBytes > static_cast<size_t> (imports_h.extent (0)),
7664 std::logic_error,
"At local row index importLIDs_h[i=" << i <<
"]="
7665 << importLIDs_h[i] <<
", offset (=" << offset <<
") + numBytes (="
7666 << numBytes <<
") > imports_h.extent(0)="
7667 << imports_h.extent (0) <<
".");
7668 #endif // HAVE_TPETRA_DEBUG
7672 #ifdef HAVE_TPETRA_DEBUG
7673 const size_t theNumBytes = PackTraits<LO, HES>::packValueCount (numEntLO);
7674 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7675 (theNumBytes > numBytes, std::logic_error,
"theNumBytes = "
7676 << theNumBytes <<
" > numBytes = " << numBytes <<
".");
7677 #endif // HAVE_TPETRA_DEBUG
7679 const char*
const inBuf = imports_h.data () + offset;
7680 const size_t actualNumBytes =
7681 PackTraits<LO, HES>::unpackValue (numEntLO, inBuf);
7683 #ifdef HAVE_TPETRA_DEBUG
7684 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7685 (actualNumBytes > numBytes, std::logic_error,
"At i = " << i
7686 <<
", actualNumBytes=" << actualNumBytes
7687 <<
" > numBytes=" << numBytes <<
".");
7688 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7689 (numEntLO == 0, std::logic_error,
"At local row index importLIDs_h[i="
7690 << i <<
"]=" << importLIDs_h[i] <<
", the number of entries read "
7691 "from the packed data is numEntLO=" << numEntLO <<
", but numBytes="
7692 << numBytes <<
" != 0.");
7694 (void) actualNumBytes;
7695 #endif // HAVE_TPETRA_DEBUG
7697 maxRowNumEnt = std::max (static_cast<size_t> (numEntLO), maxRowNumEnt);
7705 View<GO*, HES> gblColInds;
7706 View<LO*, HES> lclColInds;
7707 View<ST*, HES> vals;
7720 gblColInds = PackTraits<GO, HES>::allocateArray (gid, maxRowNumEnt,
"gids");
7721 lclColInds = PackTraits<LO, HES>::allocateArray (lid, maxRowNumEnt,
"lids");
7722 vals = PackTraits<ST, HES>::allocateArray (val, maxRowNumEnt,
"vals");
7726 for (size_type i = 0; i < numImportLIDs; ++i) {
7727 const size_t numBytes = numPacketsPerLID_h[i];
7728 if (numBytes == 0) {
7732 const char*
const inBuf = imports_h.data () + offset;
7733 const size_t actualNumBytes = PackTraits<LO, HES>::unpackValue (numEntLO, inBuf);
7734 (void) actualNumBytes;
7736 const size_t numEnt = static_cast<size_t>(numEntLO);;
7737 const LO lclRow = importLIDs_h[i];
7739 gids_out_type gidsOut = subview (gblColInds, pair_type (0, numEnt));
7740 vals_out_type valsOut = subview (vals, pair_type (0, numEnt));
7742 const size_t numBytesOut =
7743 unpackRow (gidsOut.data (), valsOut.data (), imports_h.data (),
7744 offset, numBytes, numEnt, numBytesPerValue);
7745 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
7746 (numBytes != numBytesOut, std::logic_error,
"At i = " << i <<
", "
7747 <<
"numBytes = " << numBytes <<
" != numBytesOut = " << numBytesOut
7750 const ST*
const valsRaw = const_cast<const ST*> (valsOut.data ());
7751 const GO*
const gidsRaw = const_cast<const GO*> (gidsOut.data ());
7752 this->combineGlobalValuesRaw (lclRow, numEnt, valsRaw, gidsRaw, combineMode);
7759 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7760 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
7761 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
7762 getColumnMapMultiVector (
const MV& X_domainMap,
7763 const bool force)
const
7765 using Teuchos::null;
7769 TEUCHOS_TEST_FOR_EXCEPTION(
7770 ! this->hasColMap (), std::runtime_error,
"Tpetra::CrsMatrix::getColumn"
7771 "MapMultiVector: You may only call this method if the matrix has a "
7772 "column Map. If the matrix does not yet have a column Map, you should "
7773 "first call fillComplete (with domain and range Map if necessary).");
7777 TEUCHOS_TEST_FOR_EXCEPTION(
7778 ! this->getGraph ()->isFillComplete (), std::runtime_error,
"Tpetra::"
7779 "CrsMatrix::getColumnMapMultiVector: You may only call this method if "
7780 "this matrix's graph is fill complete.");
7783 RCP<const import_type> importer = this->getGraph ()->getImporter ();
7784 RCP<const map_type> colMap = this->getColMap ();
7797 if (! importer.is_null () || force) {
7798 if (importMV_.is_null () || importMV_->getNumVectors () != numVecs) {
7799 X_colMap = rcp (
new MV (colMap, numVecs));
7802 importMV_ = X_colMap;
7805 X_colMap = importMV_;
7816 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7817 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
7820 const bool force)
const
7822 using Teuchos::null;
7828 TEUCHOS_TEST_FOR_EXCEPTION(
7829 ! this->getGraph ()->isFillComplete (), std::runtime_error,
"Tpetra::"
7830 "CrsMatrix::getRowMapMultiVector: You may only call this method if this "
7831 "matrix's graph is fill complete.");
7834 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
7838 RCP<const map_type> rowMap = this->getRowMap ();
7850 if (! exporter.is_null () || force) {
7851 if (exportMV_.is_null () || exportMV_->getNumVectors () != numVecs) {
7852 Y_rowMap = rcp (
new MV (rowMap, numVecs));
7853 exportMV_ = Y_rowMap;
7856 Y_rowMap = exportMV_;
7862 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7867 TEUCHOS_TEST_FOR_EXCEPTION(
7868 myGraph_.is_null (), std::logic_error,
"Tpetra::CrsMatrix::"
7869 "removeEmptyProcessesInPlace: This method does not work when the matrix "
7870 "was created with a constant graph (that is, when it was created using "
7871 "the version of its constructor that takes an RCP<const CrsGraph>). "
7872 "This is because the matrix is not allowed to modify the graph in that "
7873 "case, but removing empty processes requires modifying the graph.");
7874 myGraph_->removeEmptyProcessesInPlace (newMap);
7878 this->map_ = this->getRowMap ();
7882 staticGraph_ = Teuchos::rcp_const_cast<const Graph> (myGraph_);
7885 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
7886 Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
7891 const Teuchos::RCP<const map_type>& domainMap,
7892 const Teuchos::RCP<const map_type>& rangeMap,
7893 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
7895 using Teuchos::Array;
7896 using Teuchos::ArrayRCP;
7897 using Teuchos::ArrayView;
7898 using Teuchos::ParameterList;
7901 using Teuchos::rcp_implicit_cast;
7902 using Teuchos::sublist;
7903 typedef LocalOrdinal LO;
7904 typedef GlobalOrdinal GO;
7908 const crs_matrix_type& B = *
this;
7909 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
7910 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
7917 RCP<const map_type> A_rangeMap = A.
getRangeMap ();
7918 RCP<const map_type> B_domainMap = B.getDomainMap ();
7919 RCP<const map_type> B_rangeMap = B.getRangeMap ();
7921 RCP<const map_type> theDomainMap = domainMap;
7922 RCP<const map_type> theRangeMap = rangeMap;
7924 if (domainMap.is_null ()) {
7925 if (B_domainMap.is_null ()) {
7926 TEUCHOS_TEST_FOR_EXCEPTION(
7927 A_domainMap.is_null (), std::invalid_argument,
7928 "Tpetra::CrsMatrix::add: If neither A nor B have a domain Map, "
7929 "then you must supply a nonnull domain Map to this method.");
7930 theDomainMap = A_domainMap;
7932 theDomainMap = B_domainMap;
7935 if (rangeMap.is_null ()) {
7936 if (B_rangeMap.is_null ()) {
7937 TEUCHOS_TEST_FOR_EXCEPTION(
7938 A_rangeMap.is_null (), std::invalid_argument,
7939 "Tpetra::CrsMatrix::add: If neither A nor B have a range Map, "
7940 "then you must supply a nonnull range Map to this method.");
7941 theRangeMap = A_rangeMap;
7943 theRangeMap = B_rangeMap;
7947 #ifdef HAVE_TPETRA_DEBUG
7951 if (! A_domainMap.is_null () && ! A_rangeMap.is_null ()) {
7952 if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
7953 TEUCHOS_TEST_FOR_EXCEPTION(
7954 ! B_domainMap->isSameAs (*A_domainMap), std::invalid_argument,
7955 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a domain Map "
7956 "which is the same as (isSameAs) this RowMatrix's domain Map.");
7957 TEUCHOS_TEST_FOR_EXCEPTION(
7958 ! B_rangeMap->isSameAs (*A_rangeMap), std::invalid_argument,
7959 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a range Map "
7960 "which is the same as (isSameAs) this RowMatrix's range Map.");
7961 TEUCHOS_TEST_FOR_EXCEPTION(
7962 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
7963 std::invalid_argument,
7964 "Tpetra::CrsMatrix::add: The input domain Map must be the same as "
7965 "(isSameAs) this RowMatrix's domain Map.");
7966 TEUCHOS_TEST_FOR_EXCEPTION(
7967 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
7968 std::invalid_argument,
7969 "Tpetra::CrsMatrix::add: The input range Map must be the same as "
7970 "(isSameAs) this RowMatrix's range Map.");
7973 else if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
7974 TEUCHOS_TEST_FOR_EXCEPTION(
7975 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
7976 std::invalid_argument,
7977 "Tpetra::CrsMatrix::add: The input domain Map must be the same as "
7978 "(isSameAs) this RowMatrix's domain Map.");
7979 TEUCHOS_TEST_FOR_EXCEPTION(
7980 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
7981 std::invalid_argument,
7982 "Tpetra::CrsMatrix::add: The input range Map must be the same as "
7983 "(isSameAs) this RowMatrix's range Map.");
7986 TEUCHOS_TEST_FOR_EXCEPTION(
7987 domainMap.is_null () || rangeMap.is_null (), std::invalid_argument,
7988 "Tpetra::CrsMatrix::add: If neither A nor B have a domain and range "
7989 "Map, then you must supply a nonnull domain and range Map to this "
7992 #endif // HAVE_TPETRA_DEBUG
7997 bool callFillComplete =
true;
7998 RCP<ParameterList> constructorSublist;
7999 RCP<ParameterList> fillCompleteSublist;
8000 if (! params.is_null ()) {
8001 callFillComplete = params->get (
"Call fillComplete", callFillComplete);
8002 constructorSublist = sublist (params,
"Constructor parameters");
8003 fillCompleteSublist = sublist (params,
"fillComplete parameters");
8006 RCP<const map_type> A_rowMap = A.
getRowMap ();
8007 RCP<const map_type> B_rowMap = B.getRowMap ();
8008 RCP<const map_type> C_rowMap = B_rowMap;
8009 RCP<crs_matrix_type> C;
8016 if (A_rowMap->isSameAs (*B_rowMap)) {
8017 const LO localNumRows = static_cast<LO> (A_rowMap->getNodeNumElements ());
8018 ArrayRCP<size_t> C_maxNumEntriesPerRow (localNumRows, 0);
8021 if (alpha !=
ZERO) {
8022 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
8024 C_maxNumEntriesPerRow[localRow] += A_numEntries;
8029 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
8030 const size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
8031 C_maxNumEntriesPerRow[localRow] += B_numEntries;
8035 if (constructorSublist.is_null ()) {
8036 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
8039 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
8050 if (constructorSublist.is_null ()) {
8054 constructorSublist));
8058 #ifdef HAVE_TPETRA_DEBUG
8059 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
8060 "Tpetra::RowMatrix::add: C should not be null at this point. "
8061 "Please report this bug to the Tpetra developers.");
8062 #endif // HAVE_TPETRA_DEBUG
8069 if (alpha !=
ZERO) {
8070 const LO A_localNumRows = static_cast<LO> (A_rowMap->getNodeNumElements ());
8071 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
8073 const GO globalRow = A_rowMap->getGlobalElement (localRow);
8074 if (A_numEntries > static_cast<size_t> (ind.size ())) {
8075 ind.resize (A_numEntries);
8076 val.resize (A_numEntries);
8078 ArrayView<GO> indView = ind (0, A_numEntries);
8079 ArrayView<Scalar> valView = val (0, A_numEntries);
8083 for (
size_t k = 0; k < A_numEntries; ++k) {
8084 valView[k] *= alpha;
8087 C->insertGlobalValues (globalRow, indView, valView);
8092 const LO B_localNumRows = static_cast<LO> (B_rowMap->getNodeNumElements ());
8093 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
8094 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
8095 const GO globalRow = B_rowMap->getGlobalElement (localRow);
8096 if (B_numEntries > static_cast<size_t> (ind.size ())) {
8097 ind.resize (B_numEntries);
8098 val.resize (B_numEntries);
8100 ArrayView<GO> indView = ind (0, B_numEntries);
8101 ArrayView<Scalar> valView = val (0, B_numEntries);
8102 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
8105 for (
size_t k = 0; k < B_numEntries; ++k) {
8109 C->insertGlobalValues (globalRow, indView, valView);
8113 if (callFillComplete) {
8114 if (fillCompleteSublist.is_null ()) {
8115 C->fillComplete (theDomainMap, theRangeMap);
8117 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
8120 return rcp_implicit_cast<row_matrix_type> (C);
8123 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8127 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
8128 const Teuchos::RCP<const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node> > & domainTransfer,
8129 const Teuchos::RCP<const map_type>& domainMap,
8130 const Teuchos::RCP<const map_type>& rangeMap,
8131 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
8137 using Teuchos::ArrayRCP;
8138 using Teuchos::ArrayView;
8139 using Teuchos::Comm;
8140 using Teuchos::ParameterList;
8142 typedef LocalOrdinal LO;
8143 typedef GlobalOrdinal GO;
8144 typedef node_type NT;
8148 #ifdef HAVE_TPETRA_MMM_TIMINGS
8150 if(!params.is_null())
8151 label = params->get(
"Timer Label",label);
8152 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
8153 using Teuchos::TimeMonitor;
8154 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-1"))));
8162 const import_type* xferAsImport = dynamic_cast<const import_type*> (&rowTransfer);
8163 const export_type* xferAsExport = dynamic_cast<const export_type*> (&rowTransfer);
8164 TEUCHOS_TEST_FOR_EXCEPTION(
8165 xferAsImport == NULL && xferAsExport == NULL, std::invalid_argument,
8166 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' input "
8167 "argument must be either an Import or an Export, and its template "
8168 "parameters must match the corresponding template parameters of the "
8176 Teuchos::RCP<const import_type> xferDomainAsImport = Teuchos::rcp_dynamic_cast<const import_type> (domainTransfer);
8177 Teuchos::RCP<const export_type> xferDomainAsExport = Teuchos::rcp_dynamic_cast<const export_type> (domainTransfer);
8179 if(! domainTransfer.is_null()) {
8180 TEUCHOS_TEST_FOR_EXCEPTION(
8181 (xferDomainAsImport.is_null() && xferDomainAsExport.is_null()), std::invalid_argument,
8182 "Tpetra::CrsMatrix::transferAndFillComplete: The 'domainTransfer' input "
8183 "argument must be either an Import or an Export, and its template "
8184 "parameters must match the corresponding template parameters of the "
8187 TEUCHOS_TEST_FOR_EXCEPTION(
8188 ( xferAsImport != NULL || ! xferDomainAsImport.is_null() ) &&
8189 (( xferAsImport != NULL && xferDomainAsImport.is_null() ) ||
8190 ( xferAsImport == NULL && ! xferDomainAsImport.is_null() )), std::invalid_argument,
8191 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
8192 "arguments must be of the same type (either Import or Export).");
8194 TEUCHOS_TEST_FOR_EXCEPTION(
8195 ( xferAsExport != NULL || ! xferDomainAsExport.is_null() ) &&
8196 (( xferAsExport != NULL && xferDomainAsExport.is_null() ) ||
8197 ( xferAsExport == NULL && ! xferDomainAsExport.is_null() )), std::invalid_argument,
8198 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' and 'domainTransfer' input "
8199 "arguments must be of the same type (either Import or Export).");
8205 const bool communication_needed = rowTransfer.getSourceMap ()->isDistributed ();
8211 bool reverseMode =
false;
8212 bool restrictComm =
false;
8213 RCP<ParameterList> matrixparams;
8214 if (! params.is_null ()) {
8215 reverseMode = params->get (
"Reverse Mode", reverseMode);
8216 restrictComm = params->get (
"Restrict Communicator", restrictComm);
8217 matrixparams = sublist (params,
"CrsMatrix");
8222 RCP<const map_type> MyRowMap = reverseMode ?
8223 rowTransfer.getSourceMap () : rowTransfer.getTargetMap ();
8224 RCP<const map_type> MyColMap;
8225 RCP<const map_type> MyDomainMap = ! domainMap.is_null () ?
8226 domainMap : getDomainMap ();
8227 RCP<const map_type> MyRangeMap = ! rangeMap.is_null () ?
8228 rangeMap : getRangeMap ();
8229 RCP<const map_type> BaseRowMap = MyRowMap;
8230 RCP<const map_type> BaseDomainMap = MyDomainMap;
8238 if (! destMat.is_null ()) {
8249 const bool NewFlag = ! destMat->getGraph ()->isLocallyIndexed () &&
8250 ! destMat->getGraph ()->isGloballyIndexed ();
8251 TEUCHOS_TEST_FOR_EXCEPTION(
8252 ! NewFlag, std::invalid_argument,
"Tpetra::CrsMatrix::"
8253 "transferAndFillComplete: The input argument 'destMat' is only allowed "
8254 "to be nonnull, if its graph is empty (neither locally nor globally "
8263 TEUCHOS_TEST_FOR_EXCEPTION(
8264 ! destMat->getRowMap ()->isSameAs (*MyRowMap), std::invalid_argument,
8265 "Tpetra::CrsMatrix::transferAndFillComplete: The (row) Map of the "
8266 "input argument 'destMat' is not the same as the (row) Map specified "
8267 "by the input argument 'rowTransfer'.");
8268 TEUCHOS_TEST_FOR_EXCEPTION(
8269 ! destMat->checkSizes (*
this), std::invalid_argument,
8270 "Tpetra::CrsMatrix::transferAndFillComplete: You provided a nonnull "
8271 "destination matrix, but checkSizes() indicates that it is not a legal "
8272 "legal target for redistribution from the source matrix (*this). This "
8273 "may mean that they do not have the same dimensions.");
8287 TEUCHOS_TEST_FOR_EXCEPTION(
8288 ! (reverseMode || getRowMap ()->isSameAs (*rowTransfer.getSourceMap ())),
8289 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: "
8290 "rowTransfer->getSourceMap() must match this->getRowMap() in forward mode.");
8291 TEUCHOS_TEST_FOR_EXCEPTION(
8292 ! (! reverseMode || getRowMap ()->isSameAs (*rowTransfer.getTargetMap ())),
8293 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: "
8294 "rowTransfer->getTargetMap() must match this->getRowMap() in reverse mode.");
8297 TEUCHOS_TEST_FOR_EXCEPTION(
8298 ! xferDomainAsImport.is_null() && ! xferDomainAsImport->getTargetMap()->isSameAs(*domainMap),
8299 std::invalid_argument,
8300 "Tpetra::CrsMatrix::transferAndFillComplete: The target map of the 'domainTransfer' input "
8301 "argument must be the same as the rebalanced domain map 'domainMap'");
8303 TEUCHOS_TEST_FOR_EXCEPTION(
8304 ! xferDomainAsExport.is_null() && ! xferDomainAsExport->getSourceMap()->isSameAs(*domainMap),
8305 std::invalid_argument,
8306 "Tpetra::CrsMatrix::transferAndFillComplete: The source map of the 'domainTransfer' input "
8307 "argument must be the same as the rebalanced domain map 'domainMap'");
8320 const size_t NumSameIDs = rowTransfer.getNumSameIDs();
8321 ArrayView<const LO> ExportLIDs = reverseMode ?
8322 rowTransfer.getRemoteLIDs () : rowTransfer.getExportLIDs ();
8323 ArrayView<const LO> RemoteLIDs = reverseMode ?
8324 rowTransfer.getExportLIDs () : rowTransfer.getRemoteLIDs ();
8325 ArrayView<const LO> PermuteToLIDs = reverseMode ?
8326 rowTransfer.getPermuteFromLIDs () : rowTransfer.getPermuteToLIDs ();
8327 ArrayView<const LO> PermuteFromLIDs = reverseMode ?
8328 rowTransfer.getPermuteToLIDs () : rowTransfer.getPermuteFromLIDs ();
8329 Distributor& Distor = rowTransfer.getDistributor ();
8332 Teuchos::Array<int> SourcePids;
8333 Teuchos::Array<int> TargetPids;
8334 int MyPID = getComm ()->getRank ();
8337 RCP<const map_type> ReducedRowMap, ReducedColMap,
8338 ReducedDomainMap, ReducedRangeMap;
8339 RCP<const Comm<int> > ReducedComm;
8343 if (destMat.is_null ()) {
8344 destMat = rcp (
new this_type (MyRowMap, 0, StaticProfile, matrixparams));
8351 ReducedRowMap = MyRowMap->removeEmptyProcesses ();
8352 ReducedComm = ReducedRowMap.is_null () ?
8354 ReducedRowMap->getComm ();
8355 destMat->removeEmptyProcessesInPlace (ReducedRowMap);
8357 ReducedDomainMap = MyRowMap.getRawPtr () == MyDomainMap.getRawPtr () ?
8359 MyDomainMap->replaceCommWithSubset (ReducedComm);
8360 ReducedRangeMap = MyRowMap.getRawPtr () == MyRangeMap.getRawPtr () ?
8362 MyRangeMap->replaceCommWithSubset (ReducedComm);
8365 MyRowMap = ReducedRowMap;
8366 MyDomainMap = ReducedDomainMap;
8367 MyRangeMap = ReducedRangeMap;
8370 if (! ReducedComm.is_null ()) {
8371 MyPID = ReducedComm->getRank ();
8378 ReducedComm = MyRowMap->getComm ();
8384 #ifdef HAVE_TPETRA_MMM_TIMINGS
8385 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ImportSetup"))));
8388 RCP<const import_type> MyImporter = getGraph ()->getImporter ();
8391 bool bSameDomainMap = BaseDomainMap->isSameAs (*getDomainMap ());
8393 if (! restrictComm && ! MyImporter.is_null () && bSameDomainMap ) {
8400 Import_Util::getPids (*MyImporter, SourcePids,
false);
8402 else if (restrictComm && ! MyImporter.is_null () && bSameDomainMap) {
8405 IntVectorType SourceDomain_pids(getDomainMap (),
true);
8406 IntVectorType SourceCol_pids(getColMap());
8408 SourceDomain_pids.putScalar(MyPID);
8410 SourceCol_pids.doImport (SourceDomain_pids, *MyImporter,
INSERT);
8411 SourcePids.resize (getColMap ()->getNodeNumElements ());
8412 SourceCol_pids.get1dCopy (SourcePids ());
8414 else if (MyImporter.is_null () && bSameDomainMap) {
8416 SourcePids.resize (getColMap ()->getNodeNumElements ());
8417 SourcePids.assign (getColMap ()->getNodeNumElements (), MyPID);
8419 else if ( ! MyImporter.is_null () &&
8420 ! domainTransfer.is_null () ) {
8427 IntVectorType TargetDomain_pids (domainMap);
8428 TargetDomain_pids.putScalar (MyPID);
8431 IntVectorType SourceDomain_pids (getDomainMap ());
8434 IntVectorType SourceCol_pids (getColMap ());
8436 if (! reverseMode && ! xferDomainAsImport.is_null() ) {
8437 SourceDomain_pids.doExport (TargetDomain_pids, *xferDomainAsImport,
INSERT);
8439 else if (reverseMode && ! xferDomainAsExport.is_null() ) {
8440 SourceDomain_pids.doExport (TargetDomain_pids, *xferDomainAsExport,
INSERT);
8442 else if (! reverseMode && ! xferDomainAsExport.is_null() ) {
8443 SourceDomain_pids.doImport (TargetDomain_pids, *xferDomainAsExport,
INSERT);
8445 else if (reverseMode && ! xferDomainAsImport.is_null() ) {
8446 SourceDomain_pids.doImport (TargetDomain_pids, *xferDomainAsImport,
INSERT);
8449 TEUCHOS_TEST_FOR_EXCEPTION(
8450 true, std::logic_error,
"Tpetra::CrsMatrix::"
8451 "transferAndFillComplete: Should never get here! "
8452 "Please report this bug to a Tpetra developer.");
8454 SourceCol_pids.doImport (SourceDomain_pids, *MyImporter,
INSERT);
8455 SourcePids.resize (getColMap ()->getNodeNumElements ());
8456 SourceCol_pids.get1dCopy (SourcePids ());
8458 else if (BaseDomainMap->isSameAs (*BaseRowMap) &&
8459 getDomainMap ()->isSameAs (*getRowMap ())) {
8461 IntVectorType TargetRow_pids (domainMap);
8462 IntVectorType SourceRow_pids (getRowMap ());
8463 IntVectorType SourceCol_pids (getColMap ());
8465 TargetRow_pids.putScalar (MyPID);
8466 if (! reverseMode && xferAsImport != NULL) {
8467 SourceRow_pids.doExport (TargetRow_pids, *xferAsImport,
INSERT);
8469 else if (reverseMode && xferAsExport != NULL) {
8470 SourceRow_pids.doExport (TargetRow_pids, *xferAsExport,
INSERT);
8472 else if (! reverseMode && xferAsExport != NULL) {
8473 SourceRow_pids.doImport (TargetRow_pids, *xferAsExport,
INSERT);
8475 else if (reverseMode && xferAsImport != NULL) {
8476 SourceRow_pids.doImport (TargetRow_pids, *xferAsImport,
INSERT);
8479 TEUCHOS_TEST_FOR_EXCEPTION(
8480 true, std::logic_error,
"Tpetra::CrsMatrix::"
8481 "transferAndFillComplete: Should never get here! "
8482 "Please report this bug to a Tpetra developer.");
8484 SourceCol_pids.doImport (SourceRow_pids, *MyImporter,
INSERT);
8485 SourcePids.resize (getColMap ()->getNodeNumElements ());
8486 SourceCol_pids.get1dCopy (SourcePids ());
8489 TEUCHOS_TEST_FOR_EXCEPTION(
8490 true, std::invalid_argument,
"Tpetra::CrsMatrix::"
8491 "transferAndFillComplete: This method only allows either domainMap == "
8492 "getDomainMap (), or (domainMap == rowTransfer.getTargetMap () and "
8493 "getDomainMap () == getRowMap ()).");
8495 #ifdef HAVE_TPETRA_MMM_TIMINGS
8496 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-2"))));
8500 size_t constantNumPackets = destMat->constantNumberOfPackets ();
8501 if (constantNumPackets == 0) {
8502 destMat->reallocArraysForNumPacketsPerLid (ExportLIDs.size (),
8503 RemoteLIDs.size ());
8510 const size_t rbufLen = RemoteLIDs.size() * constantNumPackets;
8511 destMat->reallocImportsIfNeeded (rbufLen);
8515 #ifdef HAVE_TPETRA_DEBUG
8517 using Teuchos::outArg;
8518 using Teuchos::REDUCE_MAX;
8519 using Teuchos::reduceAll;
8522 RCP<const Teuchos::Comm<int> > comm = this->getComm ();
8523 const int myRank = comm->getRank ();
8524 const int numProcs = comm->getSize ();
8526 std::ostringstream os;
8530 destMat->numExportPacketsPerLID_.template modify<Kokkos::HostSpace> ();
8531 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
8534 numExportPacketsPerLID, ExportLIDs,
8535 SourcePids, constantNumPackets, Distor);
8537 catch (std::exception& e) {
8538 os <<
"Proc " << myRank <<
": " << e.what ();
8542 if (! comm.is_null ()) {
8543 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
8547 cerr <<
"packCrsMatrixWithOwningPIDs threw an exception: " << endl;
8549 std::ostringstream err;
8550 for (
int r = 0; r < numProcs; ++r) {
8551 if (r == myRank && lclErr != 0) {
8552 cerr << os.str () << endl;
8559 TEUCHOS_TEST_FOR_EXCEPTION(
8560 true, std::logic_error,
"packCrsMatrixWithOwningPIDs threw an "
8568 destMat->numExportPacketsPerLID_.template modify<Kokkos::HostSpace> ();
8569 Teuchos::ArrayView<size_t> numExportPacketsPerLID =
8572 numExportPacketsPerLID, ExportLIDs,
8573 SourcePids, constantNumPackets, Distor);
8575 #endif // HAVE_TPETRA_DEBUG
8578 #ifdef HAVE_TPETRA_MMM_TIMINGS
8579 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Transfer"))));
8582 if (communication_needed) {
8584 if (constantNumPackets == 0) {
8588 destMat->numExportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
8589 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8591 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
8592 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8594 Distor.doReversePostsAndWaits (numExportPacketsPerLID, 1,
8595 numImportPacketsPerLID);
8596 size_t totalImportPackets = 0;
8597 for (Array_size_type i = 0; i < numImportPacketsPerLID.size (); ++i) {
8598 totalImportPackets += numImportPacketsPerLID[i];
8603 destMat->reallocImportsIfNeeded (totalImportPackets);
8604 destMat->imports_.template modify<Kokkos::HostSpace> ();
8605 Teuchos::ArrayView<char> hostImports =
8609 destMat->exports_.template sync<Kokkos::HostSpace> ();
8610 Teuchos::ArrayView<const char> hostExports =
8612 Distor.doReversePostsAndWaits (hostExports,
8613 numExportPacketsPerLID,
8615 numImportPacketsPerLID);
8618 destMat->imports_.template modify<Kokkos::HostSpace> ();
8619 Teuchos::ArrayView<char> hostImports =
8623 destMat->exports_.template sync<Kokkos::HostSpace> ();
8624 Teuchos::ArrayView<const char> hostExports =
8626 Distor.doReversePostsAndWaits (hostExports,
8632 if (constantNumPackets == 0) {
8636 destMat->numExportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
8637 Teuchos::ArrayView<const size_t> numExportPacketsPerLID =
8639 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
8640 Teuchos::ArrayView<size_t> numImportPacketsPerLID =
8642 Distor.doPostsAndWaits (numExportPacketsPerLID, 1,
8643 numImportPacketsPerLID);
8644 size_t totalImportPackets = 0;
8645 for (Array_size_type i = 0; i < numImportPacketsPerLID.size (); ++i) {
8646 totalImportPackets += numImportPacketsPerLID[i];
8651 destMat->reallocImportsIfNeeded (totalImportPackets);
8652 destMat->imports_.template modify<Kokkos::HostSpace> ();
8653 Teuchos::ArrayView<char> hostImports =
8657 destMat->exports_.template sync<Kokkos::HostSpace> ();
8658 Teuchos::ArrayView<const char> hostExports =
8660 Distor.doPostsAndWaits (hostExports,
8661 numExportPacketsPerLID,
8663 numImportPacketsPerLID);
8666 destMat->imports_.template modify<Kokkos::HostSpace> ();
8667 Teuchos::ArrayView<char> hostImports =
8671 destMat->exports_.template sync<Kokkos::HostSpace> ();
8672 Teuchos::ArrayView<const char> hostExports =
8674 Distor.doPostsAndWaits (hostExports,
8685 #ifdef HAVE_TPETRA_MMM_TIMINGS
8686 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-1"))));
8690 destMat->numImportPacketsPerLID_.template sync<Kokkos::HostSpace> ();
8691 Teuchos::ArrayView<const size_t> numImportPacketsPerLID =
8693 destMat->imports_.template sync<Kokkos::HostSpace> ();
8694 Teuchos::ArrayView<const char> hostImports =
8698 numImportPacketsPerLID,
8699 constantNumPackets, Distor,
INSERT,
8700 NumSameIDs, PermuteToLIDs, PermuteFromLIDs);
8701 size_t N = BaseRowMap->getNodeNumElements ();
8704 ArrayRCP<size_t> CSR_rowptr(N+1);
8705 ArrayRCP<GO> CSR_colind_GID;
8706 ArrayRCP<LO> CSR_colind_LID;
8707 ArrayRCP<Scalar> CSR_vals;
8708 CSR_colind_GID.resize (mynnz);
8709 CSR_vals.resize (mynnz);
8713 if (
typeid (LO) ==
typeid (GO)) {
8714 CSR_colind_LID = Teuchos::arcp_reinterpret_cast<LO> (CSR_colind_GID);
8717 CSR_colind_LID.resize (mynnz);
8726 numImportPacketsPerLID, constantNumPackets,
8727 Distor,
INSERT, NumSameIDs, PermuteToLIDs,
8728 PermuteFromLIDs, N, mynnz, MyPID,
8729 CSR_rowptr (), CSR_colind_GID (),
8730 Teuchos::av_reinterpret_cast<impl_scalar_type> (CSR_vals ()),
8731 SourcePids (), TargetPids);
8736 #ifdef HAVE_TPETRA_MMM_TIMINGS
8737 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-2"))));
8742 Teuchos::Array<int> RemotePids;
8743 Import_Util::lowCommunicationMakeColMapAndReindex (CSR_rowptr (),
8747 TargetPids, RemotePids,
8754 ReducedColMap = (MyRowMap.getRawPtr () == MyColMap.getRawPtr ()) ?
8756 MyColMap->replaceCommWithSubset (ReducedComm);
8757 MyColMap = ReducedColMap;
8761 destMat->replaceColMap (MyColMap);
8768 if (ReducedComm.is_null ()) {
8775 #ifdef HAVE_TPETRA_MMM_TIMINGS
8776 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-3"))));
8778 if ((! reverseMode && xferAsImport != NULL) ||
8779 (reverseMode && xferAsExport != NULL)) {
8780 Import_Util::sortCrsEntries (CSR_rowptr (),
8784 else if ((! reverseMode && xferAsExport != NULL) ||
8785 (reverseMode && xferAsImport != NULL)) {
8786 Import_Util::sortAndMergeCrsEntries (CSR_rowptr (),
8789 if (CSR_rowptr[N] != mynnz) {
8790 CSR_colind_LID.resize (CSR_rowptr[N]);
8791 CSR_vals.resize (CSR_rowptr[N]);
8795 TEUCHOS_TEST_FOR_EXCEPTION(
8796 true, std::logic_error,
"Tpetra::CrsMatrix::"
8797 "transferAndFillComplete: Should never get here! "
8798 "Please report this bug to a Tpetra developer.");
8809 destMat->setAllValues (CSR_rowptr, CSR_colind_LID, CSR_vals);
8815 Teuchos::ParameterList esfc_params;
8816 #ifdef HAVE_TPETRA_MMM_TIMINGS
8817 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC CreateImporter"))));
8819 RCP<import_type> MyImport = rcp (
new import_type (MyDomainMap, MyColMap, RemotePids));
8820 #ifdef HAVE_TPETRA_MMM_TIMINGS
8821 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ESFC"))));
8823 esfc_params.set(
"Timer Label",prefix + std::string(
"TAFC"));
8825 if(!params.is_null())
8826 esfc_params.set(
"compute global constants",params->get(
"compute global constants",
true));
8828 destMat->expertStaticFillComplete (MyDomainMap, MyRangeMap, MyImport,Teuchos::null,rcp(&esfc_params,
false));
8831 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8833 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
8836 const Teuchos::RCP<const map_type>& domainMap,
8837 const Teuchos::RCP<const map_type>& rangeMap,
8838 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
8840 transferAndFillComplete (destMatrix, importer, Teuchos::null, domainMap, rangeMap, params);
8843 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8849 const Teuchos::RCP<const map_type>& domainMap,
8850 const Teuchos::RCP<const map_type>& rangeMap,
8851 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
8853 transferAndFillComplete (destMatrix, rowImporter, Teuchos::rcpFromRef(domainImporter), domainMap, rangeMap, params);
8856 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8861 const Teuchos::RCP<const map_type>& domainMap,
8862 const Teuchos::RCP<const map_type>& rangeMap,
8863 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
8865 transferAndFillComplete (destMatrix, exporter, Teuchos::null, domainMap, rangeMap, params);
8868 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
8874 const Teuchos::RCP<const map_type>& domainMap,
8875 const Teuchos::RCP<const map_type>& rangeMap,
8876 const Teuchos::RCP<Teuchos::ParameterList>& params)
const
8878 transferAndFillComplete (destMatrix, rowExporter, Teuchos::rcpFromRef(domainExporter), domainMap, rangeMap, params);
8890 #define TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR,LO,GO,NODE) \
8892 namespace Classes { \
8893 template class CrsMatrix< SCALAR , LO , GO , NODE >; \
8895 template Teuchos::RCP< CrsMatrix< SCALAR , LO , GO , NODE > > \
8896 CrsMatrix< SCALAR , LO , GO , NODE >::convert< SCALAR > () const;
8898 #define TPETRA_CRSMATRIX_CONVERT_INSTANT(SO,SI,LO,GO,NODE) \
8900 template Teuchos::RCP< CrsMatrix< SO , LO , GO , NODE > > \
8901 CrsMatrix< SI , LO , GO , NODE >::convert< SO > () const;
8903 #define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
8905 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \
8906 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \
8907 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8908 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8909 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& importer, \
8910 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8911 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8912 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \
8913 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8914 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8915 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \
8916 const Teuchos::RCP<Teuchos::ParameterList>& params);
8918 #define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
8920 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \
8921 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \
8922 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8923 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8924 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowImporter, \
8925 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8926 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8927 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainImporter, \
8928 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8929 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8930 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \
8931 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8932 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8933 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \
8934 const Teuchos::RCP<Teuchos::ParameterList>& params);
8937 #define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
8939 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \
8940 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \
8941 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8942 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8943 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& exporter, \
8944 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8945 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8946 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \
8947 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8948 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8949 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \
8950 const Teuchos::RCP<Teuchos::ParameterList>& params);
8952 #define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
8954 Teuchos::RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \
8955 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \
8956 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8957 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8958 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& rowExporter, \
8959 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8960 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8961 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& domainExporter, \
8962 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8963 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8964 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \
8965 const Teuchos::RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \
8966 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \
8967 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \
8968 const Teuchos::RCP<Teuchos::ParameterList>& params);
8971 #define TPETRA_CRSMATRIX_INSTANT(SCALAR, LO, GO ,NODE) \
8972 TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \
8973 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
8974 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \
8975 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE) \
8976 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT_TWO(SCALAR, LO, GO, NODE)
8978 #endif // TPETRA_CRSMATRIX_DEF_HPP