10 #ifndef TPETRA_CRSMATRIX_DEF_HPP
11 #define TPETRA_CRSMATRIX_DEF_HPP
23 #include "Tpetra_RowMatrix.hpp"
24 #include "Tpetra_LocalCrsMatrixOperator.hpp"
32 #include "Tpetra_Details_getDiagCopyWithoutOffsets.hpp"
40 #include "Tpetra_Details_packCrsMatrix.hpp"
41 #include "Tpetra_Details_unpackCrsMatrixAndCombine.hpp"
43 #include "Teuchos_FancyOStream.hpp"
44 #include "Teuchos_RCP.hpp"
45 #include "Teuchos_DataAccess.hpp"
46 #include "Teuchos_SerialDenseMatrix.hpp"
47 #include "KokkosBlas1_scal.hpp"
48 #include "KokkosSparse_getDiagCopy.hpp"
49 #include "KokkosSparse_spmv.hpp"
61 template<
class T,
class BinaryFunction>
62 T atomic_binary_function_update (
volatile T*
const dest,
76 T newVal = f (assume, inputVal);
77 oldVal = Kokkos::atomic_compare_exchange (dest, assume, newVal);
78 }
while (assume != oldVal);
98 template<
class Scalar>
102 typedef Teuchos::ScalarTraits<Scalar> STS;
103 return std::max (STS::magnitude (x), STS::magnitude (y));
112 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
114 CrsMatrix (
const Teuchos::RCP<const map_type>& rowMap,
115 size_t maxNumEntriesPerRow,
116 const Teuchos::RCP<Teuchos::ParameterList>& params) :
119 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, size_t "
120 "[, RCP<ParameterList>]): ";
121 Teuchos::RCP<crs_graph_type> graph;
123 graph = Teuchos::rcp (
new crs_graph_type (rowMap, maxNumEntriesPerRow,
126 catch (std::exception& e) {
127 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
128 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
129 "size_t [, RCP<ParameterList>]) threw an exception: "
136 staticGraph_ = myGraph_;
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 const Teuchos::ArrayView<const size_t>& numEntPerRowToAlloc,
145 const Teuchos::RCP<Teuchos::ParameterList>& params) :
148 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, "
149 "ArrayView<const size_t>[, RCP<ParameterList>]): ";
150 Teuchos::RCP<crs_graph_type> graph;
156 catch (std::exception& e) {
157 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
158 (
true, std::runtime_error,
"CrsGraph constructor "
159 "(RCP<const Map>, ArrayView<const size_t>"
160 "[, RCP<ParameterList>]) threw an exception: "
167 staticGraph_ = graph;
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
175 const Teuchos::RCP<const map_type>& colMap,
176 const size_t maxNumEntPerRow,
177 const Teuchos::RCP<Teuchos::ParameterList>& params) :
180 const char tfecfFuncName[] =
"CrsMatrix(RCP<const Map>, "
181 "RCP<const Map>, size_t[, RCP<ParameterList>]): ";
182 const char suffix[] =
183 " Please report this bug to the Tpetra developers.";
186 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
187 (! staticGraph_.is_null (), std::logic_error,
188 "staticGraph_ is not null at the beginning of the constructor."
190 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
191 (! myGraph_.is_null (), std::logic_error,
192 "myGraph_ is not null at the beginning of the constructor."
194 Teuchos::RCP<crs_graph_type> graph;
200 catch (std::exception& e) {
201 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
202 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
203 "RCP<const Map>, size_t[, RCP<ParameterList>]) threw an "
204 "exception: " << e.what ());
210 staticGraph_ = myGraph_;
215 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
218 const Teuchos::RCP<const map_type>& colMap,
219 const Teuchos::ArrayView<const size_t>& numEntPerRowToAlloc,
220 const Teuchos::RCP<Teuchos::ParameterList>& params) :
223 const char tfecfFuncName[] =
224 "CrsMatrix(RCP<const Map>, RCP<const Map>, "
225 "ArrayView<const size_t>[, RCP<ParameterList>]): ";
226 Teuchos::RCP<crs_graph_type> graph;
232 catch (std::exception& e) {
233 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
234 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
235 "RCP<const Map>, ArrayView<const size_t>[, "
236 "RCP<ParameterList>]) threw an exception: " << e.what ());
242 staticGraph_ = graph;
248 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
250 CrsMatrix (
const Teuchos::RCP<const crs_graph_type>& graph,
251 const Teuchos::RCP<Teuchos::ParameterList>& ) :
253 staticGraph_ (graph),
254 storageStatus_ (Details::STORAGE_1D_PACKED)
257 typedef typename local_matrix_device_type::values_type values_type;
258 const char tfecfFuncName[] =
"CrsMatrix(RCP<const CrsGraph>[, "
259 "RCP<ParameterList>]): ";
262 std::unique_ptr<std::string> prefix;
264 prefix = this->createPrefix(
"CrsMatrix",
"CrsMatrix(graph,params)");
265 std::ostringstream os;
266 os << *prefix <<
"Start" << endl;
267 std::cerr << os.str ();
270 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
271 (graph.is_null (), std::runtime_error,
"Input graph is null.");
272 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
273 (! graph->isFillComplete (), std::runtime_error,
"Input graph "
274 "is not fill complete. You must call fillComplete on the "
275 "graph before using it to construct a CrsMatrix. Note that "
276 "calling resumeFill on the graph makes it not fill complete, "
277 "even if you had previously called fillComplete. In that "
278 "case, you must call fillComplete on the graph again.");
286 const size_t numEnt = graph->lclIndsPacked_wdv.extent (0);
288 std::ostringstream os;
289 os << *prefix <<
"Allocate values: " << numEnt << endl;
290 std::cerr << os.str ();
293 values_type val (
"Tpetra::CrsMatrix::values", numEnt);
295 valuesUnpacked_wdv = valuesPacked_wdv;
300 std::ostringstream os;
301 os << *prefix <<
"Done" << endl;
302 std::cerr << os.str ();
306 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
309 const Teuchos::RCP<const crs_graph_type>& graph,
310 const Teuchos::RCP<Teuchos::ParameterList>& params) :
312 staticGraph_ (graph),
313 storageStatus_ (matrix.storageStatus_)
315 const char tfecfFuncName[] =
"CrsMatrix(RCP<const CrsGraph>, "
316 "local_matrix_device_type::values_type, "
317 "[,RCP<ParameterList>]): ";
318 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
319 (graph.is_null (), std::runtime_error,
"Input graph is null.");
320 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
321 (! graph->isFillComplete (), std::runtime_error,
"Input graph "
322 "is not fill complete. You must call fillComplete on the "
323 "graph before using it to construct a CrsMatrix. Note that "
324 "calling resumeFill on the graph makes it not fill complete, "
325 "even if you had previously called fillComplete. In that "
326 "case, you must call fillComplete on the graph again.");
328 size_t numValuesPacked = graph->lclIndsPacked_wdv.extent(0);
329 valuesPacked_wdv =
values_wdv_type(matrix.valuesPacked_wdv, 0, numValuesPacked);
331 size_t numValuesUnpacked = graph->lclIndsUnpacked_wdv.extent(0);
332 valuesUnpacked_wdv =
values_wdv_type(matrix.valuesUnpacked_wdv, 0, numValuesUnpacked);
338 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
340 CrsMatrix (
const Teuchos::RCP<const crs_graph_type>& graph,
341 const typename local_matrix_device_type::values_type& values,
342 const Teuchos::RCP<Teuchos::ParameterList>& ) :
344 staticGraph_ (graph),
345 storageStatus_ (Details::STORAGE_1D_PACKED)
347 const char tfecfFuncName[] =
"CrsMatrix(RCP<const CrsGraph>, "
348 "local_matrix_device_type::values_type, "
349 "[,RCP<ParameterList>]): ";
350 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
351 (graph.is_null (), std::runtime_error,
"Input graph is null.");
352 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
353 (! graph->isFillComplete (), std::runtime_error,
"Input graph "
354 "is not fill complete. You must call fillComplete on the "
355 "graph before using it to construct a CrsMatrix. Note that "
356 "calling resumeFill on the graph makes it not fill complete, "
357 "even if you had previously called fillComplete. In that "
358 "case, you must call fillComplete on the graph again.");
367 valuesUnpacked_wdv = valuesPacked_wdv;
378 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
381 const Teuchos::RCP<const map_type>& colMap,
382 const typename local_graph_device_type::row_map_type& rowPointers,
383 const typename local_graph_device_type::entries_type::non_const_type& columnIndices,
384 const typename local_matrix_device_type::values_type& values,
385 const Teuchos::RCP<Teuchos::ParameterList>& params) :
387 storageStatus_ (Details::STORAGE_1D_PACKED)
389 using Details::getEntryOnHost;
392 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
393 "RCP<const Map>, ptr, ind, val[, params]): ";
394 const char suffix[] =
395 ". Please report this bug to the Tpetra developers.";
399 std::unique_ptr<std::string> prefix;
401 prefix = this->createPrefix(
402 "CrsMatrix",
"CrsMatrix(rowMap,colMap,ptr,ind,val[,params])");
403 std::ostringstream os;
404 os << *prefix <<
"Start" << endl;
405 std::cerr << os.str ();
412 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
413 (values.extent(0) != columnIndices.extent(0),
414 std::invalid_argument,
"values.extent(0)=" << values.extent(0)
415 <<
" != columnIndices.extent(0) = " << columnIndices.extent(0)
417 if (debug && rowPointers.extent(0) != 0) {
418 const size_t numEnt =
419 getEntryOnHost(rowPointers, rowPointers.extent(0) - 1);
420 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
421 (numEnt !=
size_t(columnIndices.extent(0)) ||
422 numEnt !=
size_t(values.extent(0)),
423 std::invalid_argument,
"Last entry of rowPointers says that "
424 "the matrix has " << numEnt <<
" entr"
425 << (numEnt != 1 ?
"ies" :
"y") <<
", but the dimensions of "
426 "columnIndices and values don't match this. "
427 "columnIndices.extent(0)=" << columnIndices.extent (0)
428 <<
" and values.extent(0)=" << values.extent (0) <<
".");
431 RCP<crs_graph_type> graph;
433 graph = Teuchos::rcp (
new crs_graph_type (rowMap, colMap, rowPointers,
434 columnIndices, params));
436 catch (std::exception& e) {
437 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
438 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
439 "RCP<const Map>, ptr, ind[, params]) threw an exception: "
447 auto lclGraph = graph->getLocalGraphDevice ();
448 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
449 (lclGraph.row_map.extent (0) != rowPointers.extent (0) ||
450 lclGraph.entries.extent (0) != columnIndices.extent (0),
451 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, "
452 "ind[, params]) did not set the local graph correctly." << suffix);
453 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
454 (lclGraph.entries.extent (0) != values.extent (0),
455 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, ptr, ind[, "
456 "params]) did not set the local graph correctly. "
457 "lclGraph.entries.extent(0) = " << lclGraph.entries.extent (0)
458 <<
" != values.extent(0) = " << values.extent (0) << suffix);
464 staticGraph_ = graph;
474 valuesUnpacked_wdv = valuesPacked_wdv;
483 std::ostringstream os;
484 os << *prefix <<
"Done" << endl;
485 std::cerr << os.str();
489 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
492 const Teuchos::RCP<const map_type>& colMap,
493 const Teuchos::ArrayRCP<size_t>& ptr,
494 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
495 const Teuchos::ArrayRCP<Scalar>& val,
496 const Teuchos::RCP<Teuchos::ParameterList>& params) :
498 storageStatus_ (Details::STORAGE_1D_PACKED)
500 using Kokkos::Compat::getKokkosViewDeepCopy;
501 using Teuchos::av_reinterpret_cast;
503 using values_type =
typename local_matrix_device_type::values_type;
505 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
506 "RCP<const Map>, ptr, ind, val[, params]): ";
508 RCP<crs_graph_type> graph;
513 catch (std::exception& e) {
514 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
515 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
516 "RCP<const Map>, ArrayRCP<size_t>, ArrayRCP<LocalOrdinal>[, "
517 "RCP<ParameterList>]) threw an exception: " << e.what ());
523 staticGraph_ = graph;
536 auto lclGraph = staticGraph_->getLocalGraphDevice ();
537 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
538 (
size_t (lclGraph.row_map.extent (0)) !=
size_t (ptr.size ()) ||
539 size_t (lclGraph.entries.extent (0)) !=
size_t (ind.size ()),
540 std::logic_error,
"CrsGraph's constructor (rowMap, colMap, "
541 "ptr, ind[, params]) did not set the local graph correctly. "
542 "Please report this bug to the Tpetra developers.");
545 getKokkosViewDeepCopy<device_type> (av_reinterpret_cast<IST> (val ()));
547 valuesUnpacked_wdv = valuesPacked_wdv;
557 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
560 const Teuchos::RCP<const map_type>& colMap,
562 const Teuchos::RCP<Teuchos::ParameterList>& params) :
564 storageStatus_ (Details::STORAGE_1D_PACKED),
567 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
568 "RCP<const Map>, local_matrix_device_type[, RCP<ParameterList>]): ";
569 const char suffix[] =
570 " Please report this bug to the Tpetra developers.";
572 Teuchos::RCP<crs_graph_type> graph;
575 lclMatrix.graph, params));
577 catch (std::exception& e) {
578 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
579 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
580 "RCP<const Map>, local_graph_device_type[, RCP<ParameterList>]) threw an "
581 "exception: " << e.what ());
583 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
584 (!graph->isFillComplete (), std::logic_error,
"CrsGraph constructor (RCP"
585 "<const Map>, RCP<const Map>, local_graph_device_type[, RCP<ParameterList>]) "
586 "did not produce a fill-complete graph. Please report this bug to the "
587 "Tpetra developers.");
592 staticGraph_ = graph;
595 valuesUnpacked_wdv = valuesPacked_wdv;
597 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
599 "At the end of a CrsMatrix constructor that should produce "
600 "a fillComplete matrix, isFillActive() is true." << suffix);
601 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
603 "CrsMatrix constructor that should produce a fillComplete "
604 "matrix, isFillComplete() is false." << suffix);
608 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
611 const Teuchos::RCP<const map_type>& rowMap,
612 const Teuchos::RCP<const map_type>& colMap,
613 const Teuchos::RCP<const map_type>& domainMap,
614 const Teuchos::RCP<const map_type>& rangeMap,
615 const Teuchos::RCP<Teuchos::ParameterList>& params) :
617 storageStatus_ (Details::STORAGE_1D_PACKED),
620 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(RCP<const Map>, "
621 "RCP<const Map>, RCP<const Map>, RCP<const Map>, "
622 "local_matrix_device_type[, RCP<ParameterList>]): ";
623 const char suffix[] =
624 " Please report this bug to the Tpetra developers.";
626 Teuchos::RCP<crs_graph_type> graph;
628 graph = Teuchos::rcp (
new crs_graph_type (lclMatrix.graph, rowMap, colMap,
629 domainMap, rangeMap, params));
631 catch (std::exception& e) {
632 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
633 (
true, std::runtime_error,
"CrsGraph constructor (RCP<const Map>, "
634 "RCP<const Map>, RCP<const Map>, RCP<const Map>, local_graph_device_type[, "
635 "RCP<ParameterList>]) threw an exception: " << e.what ());
637 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
638 (! graph->isFillComplete (), std::logic_error,
"CrsGraph "
639 "constructor (RCP<const Map>, RCP<const Map>, RCP<const Map>, "
640 "RCP<const Map>, local_graph_device_type[, RCP<ParameterList>]) did "
641 "not produce a fillComplete graph." << suffix);
646 staticGraph_ = graph;
649 valuesUnpacked_wdv = valuesPacked_wdv;
651 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
653 "At the end of a CrsMatrix constructor that should produce "
654 "a fillComplete matrix, isFillActive() is true." << suffix);
655 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
657 "CrsMatrix constructor that should produce a fillComplete "
658 "matrix, isFillComplete() is false." << suffix);
662 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
665 const Teuchos::RCP<const map_type>& rowMap,
666 const Teuchos::RCP<const map_type>& colMap,
667 const Teuchos::RCP<const map_type>& domainMap,
668 const Teuchos::RCP<const map_type>& rangeMap,
669 const Teuchos::RCP<const import_type>& importer,
670 const Teuchos::RCP<const export_type>& exporter,
671 const Teuchos::RCP<Teuchos::ParameterList>& params) :
673 storageStatus_ (Details::STORAGE_1D_PACKED),
677 const char tfecfFuncName[] =
"Tpetra::CrsMatrix"
678 "(lclMat,Map,Map,Map,Map,Import,Export,params): ";
679 const char suffix[] =
680 " Please report this bug to the Tpetra developers.";
682 Teuchos::RCP<crs_graph_type> graph;
685 domainMap, rangeMap, importer,
688 catch (std::exception& e) {
689 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
690 (
true, std::runtime_error,
"CrsGraph constructor "
691 "(local_graph_device_type, Map, Map, Map, Map, Import, Export, "
692 "params) threw: " << e.what ());
694 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
695 (!graph->isFillComplete (), std::logic_error,
"CrsGraph "
696 "constructor (local_graph_device_type, Map, Map, Map, Map, Import, "
697 "Export, params) did not produce a fill-complete graph. "
698 "Please report this bug to the Tpetra developers.");
703 staticGraph_ = graph;
706 valuesUnpacked_wdv = valuesPacked_wdv;
708 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
710 "At the end of a CrsMatrix constructor that should produce "
711 "a fillComplete matrix, isFillActive() is true." << suffix);
712 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
714 "CrsMatrix constructor that should produce a fillComplete "
715 "matrix, isFillComplete() is false." << suffix);
719 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
722 const Teuchos::DataAccess copyOrView):
724 staticGraph_ (source.getCrsGraph()),
725 storageStatus_ (source.storageStatus_)
727 const char tfecfFuncName[] =
"Tpetra::CrsMatrix("
728 "const CrsMatrix&, const Teuchos::DataAccess): ";
729 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
731 "Source graph must be fillComplete().");
733 if (copyOrView == Teuchos::Copy) {
734 using values_type =
typename local_matrix_device_type::values_type;
736 using Kokkos::view_alloc;
737 using Kokkos::WithoutInitializing;
738 values_type newvals (view_alloc (
"val", WithoutInitializing),
743 valuesUnpacked_wdv = valuesPacked_wdv;
746 else if (copyOrView == Teuchos::View) {
752 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
753 (
true, std::invalid_argument,
"Second argument 'copyOrView' "
754 "has an invalid value " << copyOrView <<
". Valid values "
755 "include Teuchos::Copy = " << Teuchos::Copy <<
" and "
756 "Teuchos::View = " << Teuchos::View <<
".");
761 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
766 std::swap(crs_matrix.
importMV_, this->importMV_);
767 std::swap(crs_matrix.
exportMV_, this->exportMV_);
768 std::swap(crs_matrix.staticGraph_, this->staticGraph_);
769 std::swap(crs_matrix.myGraph_, this->myGraph_);
770 std::swap(crs_matrix.valuesPacked_wdv, this->valuesPacked_wdv);
771 std::swap(crs_matrix.valuesUnpacked_wdv, this->valuesUnpacked_wdv);
774 std::swap(crs_matrix.
nonlocals_, this->nonlocals_);
777 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
778 Teuchos::RCP<const Teuchos::Comm<int> >
781 return getCrsGraphRef ().getComm ();
784 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
788 return fillComplete_;
791 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
795 return ! fillComplete_;
798 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
802 return this->getCrsGraphRef ().isStorageOptimized ();
805 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
809 return getCrsGraphRef ().isLocallyIndexed ();
812 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
816 return getCrsGraphRef ().isGloballyIndexed ();
819 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
823 return getCrsGraphRef ().hasColMap ();
826 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
830 return getCrsGraphRef ().getGlobalNumEntries ();
833 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
837 return getCrsGraphRef ().getLocalNumEntries ();
840 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
844 return getCrsGraphRef ().getGlobalNumRows ();
847 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
851 return getCrsGraphRef ().getGlobalNumCols ();
854 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
858 return getCrsGraphRef ().getLocalNumRows ();
862 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
866 return getCrsGraphRef ().getLocalNumCols ();
870 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
874 return getCrsGraphRef ().getNumEntriesInGlobalRow (globalRow);
877 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
881 return getCrsGraphRef ().getNumEntriesInLocalRow (localRow);
884 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
888 return getCrsGraphRef ().getGlobalMaxNumRowEntries ();
891 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
895 return getCrsGraphRef ().getLocalMaxNumRowEntries ();
898 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
902 return getRowMap ()->getIndexBase ();
905 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
906 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
909 return getCrsGraphRef ().getRowMap ();
912 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
913 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
916 return getCrsGraphRef ().getColMap ();
919 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
920 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
923 return getCrsGraphRef ().getDomainMap ();
926 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
927 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
930 return getCrsGraphRef ().getRangeMap ();
933 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
934 Teuchos::RCP<const RowGraph<LocalOrdinal, GlobalOrdinal, Node> >
937 if (staticGraph_ != Teuchos::null) {
943 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
944 Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node> >
947 if (staticGraph_ != Teuchos::null) {
953 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
958 #ifdef HAVE_TPETRA_DEBUG
959 constexpr
bool debug =
true;
961 constexpr
bool debug =
false;
962 #endif // HAVE_TPETRA_DEBUG
964 if (! this->staticGraph_.is_null ()) {
965 return * (this->staticGraph_);
969 const char tfecfFuncName[] =
"getCrsGraphRef: ";
970 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
971 (this->myGraph_.is_null (), std::logic_error,
972 "Both staticGraph_ and myGraph_ are null. "
973 "Please report this bug to the Tpetra developers.");
975 return * (this->myGraph_);
979 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
980 typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_device_type
984 auto numCols = staticGraph_->getColMap()->getLocalNumElements();
987 valuesPacked_wdv.getDeviceView(Access::ReadWrite),
988 staticGraph_->getLocalGraphDevice());
991 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
992 typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_matrix_host_type
996 auto numCols = staticGraph_->
getColMap()->getLocalNumElements();
997 return local_matrix_host_type(
"Tpetra::CrsMatrix::lclMatrixHost", numCols,
998 valuesPacked_wdv.getHostView(Access::ReadWrite),
999 staticGraph_->getLocalGraphHost());
1002 #if KOKKOSKERNELS_VERSION < 40299
1004 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1005 std::shared_ptr<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_multiply_op_type>
1009 auto localMatrix = getLocalMatrixDevice();
1010 #if defined(KOKKOSKERNELS_ENABLE_TPL_CUSPARSE) || defined(KOKKOSKERNELS_ENABLE_TPL_ROCSPARSE) || defined(KOKKOSKERNELS_ENABLE_TPL_MKL)
1011 if(this->getLocalNumEntries() <=
size_t(Teuchos::OrdinalTraits<LocalOrdinal>::max()))
1013 if(this->ordinalRowptrs.data() ==
nullptr)
1015 auto originalRowptrs = localMatrix.graph.row_map;
1018 this->ordinalRowptrs = ordinal_rowptrs_type(
1019 Kokkos::ViewAllocateWithoutInitializing(
"CrsMatrix::ordinalRowptrs"), originalRowptrs.extent(0));
1020 auto ordinalRowptrs_ = this->ordinalRowptrs;
1021 Kokkos::parallel_for(
"CrsMatrix::getLocalMultiplyOperator::convertRowptrs",
1022 Kokkos::RangePolicy<execution_space>(0, originalRowptrs.extent(0)),
1023 KOKKOS_LAMBDA(LocalOrdinal i)
1025 ordinalRowptrs_(i) = originalRowptrs(i);
1029 return std::make_shared<local_multiply_op_type>(
1030 std::make_shared<local_matrix_device_type>(localMatrix), this->ordinalRowptrs);
1034 return std::make_shared<local_multiply_op_type>(
1035 std::make_shared<local_matrix_device_type>(localMatrix));
1039 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1043 return myGraph_.is_null ();
1046 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1053 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1060 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1069 const char tfecfFuncName[] =
"allocateValues: ";
1070 const char suffix[] =
1071 " Please report this bug to the Tpetra developers.";
1072 ProfilingRegion region(
"Tpetra::CrsMatrix::allocateValues");
1074 std::unique_ptr<std::string> prefix;
1076 prefix = this->createPrefix(
"CrsMatrix",
"allocateValues");
1077 std::ostringstream os;
1078 os << *prefix <<
"lg: "
1079 << (lg == LocalIndices ?
"Local" :
"Global") <<
"Indices"
1081 << (gas == GraphAlreadyAllocated ?
"Already" :
"NotYet")
1082 <<
"Allocated" << endl;
1083 std::cerr << os.str();
1086 const bool debug = Behavior::debug(
"CrsMatrix");
1088 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1089 (this->staticGraph_.is_null (), std::logic_error,
1090 "staticGraph_ is null." << suffix);
1095 if ((gas == GraphAlreadyAllocated) !=
1096 staticGraph_->indicesAreAllocated ()) {
1097 const char err1[] =
"The caller has asserted that the graph "
1099 const char err2[] =
"already allocated, but the static graph "
1100 "says that its indices are ";
1101 const char err3[] =
"already allocated. ";
1102 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1103 (gas == GraphAlreadyAllocated &&
1104 ! staticGraph_->indicesAreAllocated (), std::logic_error,
1105 err1 << err2 <<
"not " << err3 << suffix);
1106 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1107 (gas != GraphAlreadyAllocated &&
1108 staticGraph_->indicesAreAllocated (), std::logic_error,
1109 err1 <<
"not " << err2 << err3 << suffix);
1117 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1118 (! this->staticGraph_->indicesAreAllocated () &&
1119 this->myGraph_.is_null (), std::logic_error,
1120 "The static graph says that its indices are not allocated, "
1121 "but the graph is not owned by the matrix." << suffix);
1124 if (gas == GraphNotYetAllocated) {
1126 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1127 (this->myGraph_.is_null (), std::logic_error,
1128 "gas = GraphNotYetAllocated, but myGraph_ is null." << suffix);
1131 this->myGraph_->allocateIndices (lg, verbose);
1133 catch (std::exception& e) {
1134 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1135 (
true, std::runtime_error,
"CrsGraph::allocateIndices "
1136 "threw an exception: " << e.what ());
1139 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1140 (
true, std::runtime_error,
"CrsGraph::allocateIndices "
1141 "threw an exception not a subclass of std::exception.");
1146 const size_t lclTotalNumEntries = this->staticGraph_->getLocalAllocationSize();
1148 const size_t lclNumRows = this->staticGraph_->getLocalNumRows ();
1149 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1150 (this->staticGraph_->getRowPtrsUnpackedHost()(lclNumRows) != lclTotalNumEntries, std::logic_error,
1151 "length of staticGraph's lclIndsUnpacked does not match final entry of rowPtrsUnapcked_host." << suffix);
1155 using values_type =
typename local_matrix_device_type::values_type;
1157 std::ostringstream os;
1158 os << *prefix <<
"Allocate values_wdv: Pre "
1159 << valuesUnpacked_wdv.extent(0) <<
", post "
1160 << lclTotalNumEntries << endl;
1161 std::cerr << os.str();
1165 values_type(
"Tpetra::CrsMatrix::values",
1166 lclTotalNumEntries));
1170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1176 using ::Tpetra::Details::getEntryOnHost;
1177 using Teuchos::arcp_const_cast;
1178 using Teuchos::Array;
1179 using Teuchos::ArrayRCP;
1180 using Teuchos::null;
1184 using row_map_type =
typename local_graph_device_type::row_map_type;
1185 using lclinds_1d_type =
typename Graph::local_graph_device_type::entries_type::non_const_type;
1186 using values_type =
typename local_matrix_device_type::values_type;
1188 (
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix");
1190 const char tfecfFuncName[] =
"fillLocalGraphAndMatrix (called from "
1191 "fillComplete or expertStaticFillComplete): ";
1192 const char suffix[] =
1193 " Please report this bug to the Tpetra developers.";
1197 std::unique_ptr<std::string> prefix;
1199 prefix = this->createPrefix(
"CrsMatrix",
"fillLocalGraphAndMatrix");
1200 std::ostringstream os;
1201 os << *prefix << endl;
1202 std::cerr << os.str ();
1208 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1209 (myGraph_.is_null (), std::logic_error,
"The nonconst graph "
1210 "(myGraph_) is null. This means that the matrix has a "
1211 "const (a.k.a. \"static\") graph. fillComplete or "
1212 "expertStaticFillComplete should never call "
1213 "fillLocalGraphAndMatrix in that case." << suffix);
1216 const size_t lclNumRows = this->getLocalNumRows ();
1231 typename Graph::local_graph_device_type::row_map_type curRowOffsets =
1232 myGraph_->rowPtrsUnpacked_dev_;
1235 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1236 (curRowOffsets.extent (0) == 0, std::logic_error,
1237 "curRowOffsets.extent(0) == 0.");
1238 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1239 (curRowOffsets.extent (0) != lclNumRows + 1, std::logic_error,
1240 "curRowOffsets.extent(0) = "
1241 << curRowOffsets.extent (0) <<
" != lclNumRows + 1 = "
1242 << (lclNumRows + 1) <<
".");
1243 const size_t numOffsets = curRowOffsets.extent (0);
1244 const auto valToCheck = myGraph_->getRowPtrsUnpackedHost()(numOffsets - 1);
1245 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1247 myGraph_->lclIndsUnpacked_wdv.extent (0) != valToCheck,
1248 std::logic_error,
"numOffsets = " <<
1249 numOffsets <<
" != 0 and myGraph_->lclIndsUnpacked_wdv.extent(0) = "
1250 << myGraph_->lclIndsUnpacked_wdv.extent (0) <<
" != curRowOffsets("
1251 << numOffsets <<
") = " << valToCheck <<
".");
1254 if (myGraph_->getLocalNumEntries() !=
1255 myGraph_->getLocalAllocationSize()) {
1259 typename row_map_type::non_const_type k_ptrs;
1260 row_map_type k_ptrs_const;
1261 lclinds_1d_type k_inds;
1265 std::ostringstream os;
1266 const auto numEnt = myGraph_->getLocalNumEntries();
1267 const auto allocSize = myGraph_->getLocalAllocationSize();
1268 os << *prefix <<
"Unpacked 1-D storage: numEnt=" << numEnt
1269 <<
", allocSize=" << allocSize << endl;
1270 std::cerr << os.str ();
1278 if (debug && curRowOffsets.extent (0) != 0) {
1279 const size_t numOffsets =
1280 static_cast<size_t> (curRowOffsets.extent (0));
1281 const auto valToCheck = myGraph_->getRowPtrsUnpackedHost()(numOffsets - 1);
1282 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1283 (static_cast<size_t> (valToCheck) !=
1284 static_cast<size_t> (valuesUnpacked_wdv.extent (0)),
1285 std::logic_error,
"(unpacked branch) Before "
1286 "allocating or packing, curRowOffsets(" << (numOffsets-1)
1287 <<
") = " << valToCheck <<
" != valuesUnpacked_wdv.extent(0)"
1288 " = " << valuesUnpacked_wdv.extent (0) <<
".");
1289 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1290 (static_cast<size_t> (valToCheck) !=
1291 static_cast<size_t> (myGraph_->lclIndsUnpacked_wdv.extent (0)),
1292 std::logic_error,
"(unpacked branch) Before "
1293 "allocating or packing, curRowOffsets(" << (numOffsets-1)
1294 <<
") = " << valToCheck
1295 <<
" != myGraph_->lclIndsUnpacked_wdv.extent(0) = "
1296 << myGraph_->lclIndsUnpacked_wdv.extent (0) <<
".");
1304 size_t lclTotalNumEntries = 0;
1310 std::ostringstream os;
1311 os << *prefix <<
"Allocate packed row offsets: "
1312 << (lclNumRows+1) << endl;
1313 std::cerr << os.str ();
1315 typename row_map_type::non_const_type
1316 packedRowOffsets (
"Tpetra::CrsGraph::ptr", lclNumRows + 1);
1317 typename row_entries_type::const_type numRowEnt_h =
1318 myGraph_->k_numRowEntries_;
1321 lclTotalNumEntries =
1325 k_ptrs = packedRowOffsets;
1326 k_ptrs_const = k_ptrs;
1330 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1331 (static_cast<size_t> (k_ptrs.extent (0)) != lclNumRows + 1,
1333 "(unpacked branch) After packing k_ptrs, "
1334 "k_ptrs.extent(0) = " << k_ptrs.extent (0) <<
" != "
1335 "lclNumRows+1 = " << (lclNumRows+1) <<
".");
1336 const auto valToCheck = getEntryOnHost (k_ptrs, lclNumRows);
1337 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1338 (valToCheck != lclTotalNumEntries, std::logic_error,
1339 "(unpacked branch) After filling k_ptrs, "
1340 "k_ptrs(lclNumRows=" << lclNumRows <<
") = " << valToCheck
1341 <<
" != total number of entries on the calling process = "
1342 << lclTotalNumEntries <<
".");
1347 std::ostringstream os;
1348 os << *prefix <<
"Allocate packed local column indices: "
1349 << lclTotalNumEntries << endl;
1350 std::cerr << os.str ();
1352 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::lclInds", lclTotalNumEntries);
1354 std::ostringstream os;
1355 os << *prefix <<
"Allocate packed values: "
1356 << lclTotalNumEntries << endl;
1357 std::cerr << os.str ();
1359 k_vals = values_type (
"Tpetra::CrsMatrix::values", lclTotalNumEntries);
1371 using inds_packer_type = pack_functor<
1372 typename Graph::local_graph_device_type::entries_type::non_const_type,
1373 typename Graph::local_inds_dualv_type::t_dev::const_type,
1374 typename Graph::local_graph_device_type::row_map_type::non_const_type,
1375 typename Graph::local_graph_device_type::row_map_type>;
1376 inds_packer_type indsPacker (
1378 myGraph_->lclIndsUnpacked_wdv.getDeviceView(Access::ReadOnly),
1379 k_ptrs, curRowOffsets);
1381 using range_type = Kokkos::RangePolicy<exec_space, LocalOrdinal>;
1382 Kokkos::parallel_for
1383 (
"Tpetra::CrsMatrix pack column indices",
1384 range_type (0, lclNumRows), indsPacker);
1388 using vals_packer_type = pack_functor<
1389 typename values_type::non_const_type,
1390 typename values_type::const_type,
1391 typename row_map_type::non_const_type,
1392 typename row_map_type::const_type>;
1393 vals_packer_type valsPacker (
1395 this->valuesUnpacked_wdv.getDeviceView(Access::ReadOnly),
1396 k_ptrs, curRowOffsets);
1397 Kokkos::parallel_for (
"Tpetra::CrsMatrix pack values",
1398 range_type (0, lclNumRows), valsPacker);
1401 const char myPrefix[] =
"(\"Optimize Storage\""
1402 "=true branch) After packing, ";
1403 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1404 (k_ptrs.extent (0) == 0, std::logic_error, myPrefix
1405 <<
"k_ptrs.extent(0) = 0. This probably means that "
1406 "rowPtrsUnpacked_ was never allocated.");
1407 if (k_ptrs.extent (0) != 0) {
1408 const size_t numOffsets (k_ptrs.extent (0));
1409 const auto valToCheck =
1410 getEntryOnHost (k_ptrs, numOffsets - 1);
1411 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1412 (
size_t (valToCheck) != k_vals.extent (0),
1413 std::logic_error, myPrefix <<
1414 "k_ptrs(" << (numOffsets-1) <<
") = " << valToCheck <<
1415 " != k_vals.extent(0) = " << k_vals.extent (0) <<
".");
1416 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1417 (
size_t (valToCheck) != k_inds.extent (0),
1418 std::logic_error, myPrefix <<
1419 "k_ptrs(" << (numOffsets-1) <<
") = " << valToCheck <<
1420 " != k_inds.extent(0) = " << k_inds.extent (0) <<
".");
1424 myGraph_->setRowPtrsPacked(k_ptrs_const);
1425 myGraph_->lclIndsPacked_wdv =
1432 myGraph_->rowPtrsPacked_dev_ = myGraph_->rowPtrsUnpacked_dev_;
1433 myGraph_->rowPtrsPacked_host_ = myGraph_->rowPtrsUnpacked_host_;
1434 myGraph_->packedUnpackedRowPtrsMatch_ =
true;
1435 myGraph_->lclIndsPacked_wdv = myGraph_->lclIndsUnpacked_wdv;
1436 valuesPacked_wdv = valuesUnpacked_wdv;
1439 std::ostringstream os;
1440 os << *prefix <<
"Storage already packed: rowPtrsUnpacked_: "
1441 << myGraph_->getRowPtrsUnpackedHost().extent(0) <<
", lclIndsUnpacked_wdv: "
1442 << myGraph_->lclIndsUnpacked_wdv.extent(0) <<
", valuesUnpacked_wdv: "
1443 << valuesUnpacked_wdv.extent(0) << endl;
1444 std::cerr << os.str();
1448 const char myPrefix[] =
1449 "(\"Optimize Storage\"=false branch) ";
1450 auto rowPtrsUnpackedHost = myGraph_->getRowPtrsUnpackedHost();
1451 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1452 (myGraph_->rowPtrsUnpacked_dev_.extent (0) == 0, std::logic_error, myPrefix
1453 <<
"myGraph->rowPtrsUnpacked_dev_.extent(0) = 0. This probably means "
1454 "that rowPtrsUnpacked_ was never allocated.");
1455 if (myGraph_->rowPtrsUnpacked_dev_.extent (0) != 0) {
1456 const size_t numOffsets = rowPtrsUnpackedHost.extent (0);
1457 const auto valToCheck = rowPtrsUnpackedHost(numOffsets - 1);
1458 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1459 (
size_t (valToCheck) != valuesPacked_wdv.extent (0),
1460 std::logic_error, myPrefix <<
1461 "k_ptrs_const(" << (numOffsets-1) <<
") = " << valToCheck
1462 <<
" != valuesPacked_wdv.extent(0) = "
1463 << valuesPacked_wdv.extent (0) <<
".");
1464 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1465 (
size_t (valToCheck) != myGraph_->lclIndsPacked_wdv.extent (0),
1466 std::logic_error, myPrefix <<
1467 "k_ptrs_const(" << (numOffsets-1) <<
") = " << valToCheck
1468 <<
" != myGraph_->lclIndsPacked.extent(0) = "
1469 << myGraph_->lclIndsPacked_wdv.extent (0) <<
".");
1475 const char myPrefix[] =
"After packing, ";
1476 auto rowPtrsPackedHost = myGraph_->getRowPtrsPackedHost();
1477 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1478 (
size_t (rowPtrsPackedHost.extent (0)) !=
size_t (lclNumRows + 1),
1479 std::logic_error, myPrefix <<
"myGraph_->rowPtrsPacked_host_.extent(0) = "
1480 << rowPtrsPackedHost.extent (0) <<
" != lclNumRows+1 = " <<
1481 (lclNumRows+1) <<
".");
1482 if (rowPtrsPackedHost.extent (0) != 0) {
1483 const size_t numOffsets (rowPtrsPackedHost.extent (0));
1484 const size_t valToCheck = rowPtrsPackedHost(numOffsets-1);
1485 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1486 (valToCheck !=
size_t (valuesPacked_wdv.extent (0)),
1487 std::logic_error, myPrefix <<
"k_ptrs_const(" <<
1488 (numOffsets-1) <<
") = " << valToCheck
1489 <<
" != valuesPacked_wdv.extent(0) = "
1490 << valuesPacked_wdv.extent (0) <<
".");
1491 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1492 (valToCheck !=
size_t (myGraph_->lclIndsPacked_wdv.extent (0)),
1493 std::logic_error, myPrefix <<
"k_ptrs_const(" <<
1494 (numOffsets-1) <<
") = " << valToCheck
1495 <<
" != myGraph_->lclIndsPacked_wdvk_inds.extent(0) = "
1496 << myGraph_->lclIndsPacked_wdv.extent (0) <<
".");
1504 const bool defaultOptStorage =
1505 ! isStaticGraph () || staticGraph_->isStorageOptimized ();
1506 const bool requestOptimizedStorage =
1507 (! params.is_null () &&
1508 params->get (
"Optimize Storage", defaultOptStorage)) ||
1509 (params.is_null () && defaultOptStorage);
1514 if (requestOptimizedStorage) {
1519 std::ostringstream os;
1520 os << *prefix <<
"Optimizing storage: free k_numRowEntries_: "
1521 << myGraph_->k_numRowEntries_.extent(0) << endl;
1522 std::cerr << os.str();
1525 myGraph_->k_numRowEntries_ = row_entries_type ();
1530 myGraph_->rowPtrsUnpacked_dev_ = myGraph_->rowPtrsPacked_dev_;
1531 myGraph_->rowPtrsUnpacked_host_ = myGraph_->rowPtrsPacked_host_;
1532 myGraph_->packedUnpackedRowPtrsMatch_ =
true;
1533 myGraph_->lclIndsUnpacked_wdv = myGraph_->lclIndsPacked_wdv;
1534 valuesUnpacked_wdv = valuesPacked_wdv;
1536 myGraph_->storageStatus_ = Details::STORAGE_1D_PACKED;
1537 this->storageStatus_ = Details::STORAGE_1D_PACKED;
1541 std::ostringstream os;
1542 os << *prefix <<
"User requested NOT to optimize storage"
1544 std::cerr << os.str();
1549 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1554 using ::Tpetra::Details::ProfilingRegion;
1555 using Teuchos::ArrayRCP;
1556 using Teuchos::Array;
1557 using Teuchos::null;
1561 using row_map_type =
typename Graph::local_graph_device_type::row_map_type;
1562 using non_const_row_map_type =
typename row_map_type::non_const_type;
1563 using values_type =
typename local_matrix_device_type::values_type;
1564 ProfilingRegion regionFLM(
"Tpetra::CrsMatrix::fillLocalMatrix");
1565 const size_t lclNumRows = getLocalNumRows();
1568 std::unique_ptr<std::string> prefix;
1570 prefix = this->createPrefix(
"CrsMatrix",
"fillLocalMatrix");
1571 std::ostringstream os;
1572 os << *prefix <<
"lclNumRows: " << lclNumRows << endl;
1573 std::cerr << os.str ();
1585 size_t nodeNumEntries = staticGraph_->getLocalNumEntries ();
1586 size_t nodeNumAllocated = staticGraph_->getLocalAllocationSize ();
1587 row_map_type k_rowPtrs = staticGraph_->rowPtrsPacked_dev_;
1589 row_map_type k_ptrs;
1595 bool requestOptimizedStorage =
true;
1596 const bool default_OptimizeStorage =
1597 ! isStaticGraph() || staticGraph_->isStorageOptimized();
1598 if (! params.is_null() &&
1599 ! params->get(
"Optimize Storage", default_OptimizeStorage)) {
1600 requestOptimizedStorage =
false;
1607 if (! staticGraph_->isStorageOptimized () &&
1608 requestOptimizedStorage) {
1610 (
true, std::runtime_error,
"You requested optimized storage "
1611 "by setting the \"Optimize Storage\" flag to \"true\" in "
1612 "the ParameterList, or by virtue of default behavior. "
1613 "However, the associated CrsGraph was filled separately and "
1614 "requested not to optimize storage. Therefore, the "
1615 "CrsMatrix cannot optimize storage.");
1616 requestOptimizedStorage =
false;
1641 if (nodeNumEntries != nodeNumAllocated) {
1643 std::ostringstream os;
1644 os << *prefix <<
"Unpacked 1-D storage: numEnt="
1645 << nodeNumEntries <<
", allocSize=" << nodeNumAllocated
1647 std::cerr << os.str();
1652 std::ostringstream os;
1653 os << *prefix <<
"Allocate packed row offsets: "
1654 << (lclNumRows+1) << endl;
1655 std::cerr << os.str();
1657 non_const_row_map_type tmpk_ptrs (
"Tpetra::CrsGraph::ptr",
1662 size_t lclTotalNumEntries = 0;
1665 typename row_entries_type::const_type numRowEnt_h =
1666 staticGraph_->k_numRowEntries_;
1668 lclTotalNumEntries =
1675 std::ostringstream os;
1676 os << *prefix <<
"Allocate packed values: "
1677 << lclTotalNumEntries << endl;
1678 std::cerr << os.str ();
1680 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1684 typename values_type::non_const_type,
1685 typename values_type::const_type,
1686 typename row_map_type::non_const_type,
1687 typename row_map_type::const_type> valsPacker
1688 (k_vals, valuesUnpacked_wdv.getDeviceView(Access::ReadOnly),
1689 tmpk_ptrs, k_rowPtrs);
1692 using range_type = Kokkos::RangePolicy<exec_space, LocalOrdinal>;
1693 Kokkos::parallel_for (
"Tpetra::CrsMatrix pack values",
1694 range_type (0, lclNumRows), valsPacker);
1698 valuesPacked_wdv = valuesUnpacked_wdv;
1700 std::ostringstream os;
1701 os << *prefix <<
"Storage already packed: "
1702 <<
"valuesUnpacked_wdv: " << valuesUnpacked_wdv.extent(0) << endl;
1703 std::cerr << os.str();
1708 if (requestOptimizedStorage) {
1711 valuesUnpacked_wdv = valuesPacked_wdv;
1713 this->storageStatus_ = Details::STORAGE_1D_PACKED;
1717 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1722 const typename crs_graph_type::SLocalGlobalViews& newInds,
1723 const Teuchos::ArrayView<impl_scalar_type>& oldRowVals,
1724 const Teuchos::ArrayView<const impl_scalar_type>& newRowVals,
1725 const ELocalGlobal lg,
1726 const ELocalGlobal I)
1728 const size_t oldNumEnt = rowInfo.numEntries;
1729 const size_t numInserted = graph.insertIndices (rowInfo, newInds, lg, I);
1735 if (numInserted > 0) {
1736 const size_t startOffset = oldNumEnt;
1737 memcpy ((
void*) &oldRowVals[startOffset], &newRowVals[0],
1738 numInserted *
sizeof (impl_scalar_type));
1742 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1746 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1747 const Teuchos::ArrayView<const Scalar>& values,
1751 const char tfecfFuncName[] =
"insertLocalValues: ";
1753 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1754 (! this->isFillActive (), std::runtime_error,
1755 "Fill is not active. After calling fillComplete, you must call "
1756 "resumeFill before you may insert entries into the matrix again.");
1757 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1758 (this->isStaticGraph (), std::runtime_error,
1759 "Cannot insert indices with static graph; use replaceLocalValues() "
1763 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1764 (graph.
colMap_.is_null (), std::runtime_error,
1765 "Cannot insert local indices without a column map.");
1766 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1768 std::runtime_error,
"Graph indices are global; use "
1769 "insertGlobalValues().");
1770 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1771 (values.size () != indices.size (), std::runtime_error,
1772 "values.size() = " << values.size ()
1773 <<
" != indices.size() = " << indices.size () <<
".");
1774 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1775 ! graph.
rowMap_->isNodeLocalElement (lclRow), std::runtime_error,
1776 "Local row index " << lclRow <<
" does not belong to this process.");
1778 if (! graph.indicesAreAllocated ()) {
1782 this->allocateValues (LocalIndices, GraphNotYetAllocated, verbose);
1785 #ifdef HAVE_TPETRA_DEBUG
1786 const size_t numEntriesToAdd =
static_cast<size_t> (indices.size ());
1791 using Teuchos::toString;
1794 Teuchos::Array<LocalOrdinal> badColInds;
1795 bool allInColMap =
true;
1796 for (
size_t k = 0; k < numEntriesToAdd; ++k) {
1798 allInColMap =
false;
1799 badColInds.push_back (indices[k]);
1802 if (! allInColMap) {
1803 std::ostringstream os;
1804 os <<
"You attempted to insert entries in owned row " << lclRow
1805 <<
", at the following column indices: " << toString (indices)
1807 os <<
"Of those, the following indices are not in the column Map on "
1808 "this process: " << toString (badColInds) <<
"." << endl <<
"Since "
1809 "the matrix has a column Map already, it is invalid to insert "
1810 "entries at those locations.";
1811 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1812 (
true, std::invalid_argument, os.str ());
1815 #endif // HAVE_TPETRA_DEBUG
1819 auto valsView = this->getValuesViewHostNonConst(rowInfo);
1821 auto fun = [&](
size_t const k,
size_t const ,
size_t const offset) {
1822 valsView[offset] += values[k]; };
1823 std::function<void(size_t const, size_t const, size_t const)> cb(std::ref(fun));
1824 graph.insertLocalIndicesImpl(lclRow, indices, cb);
1825 }
else if (CM ==
INSERT) {
1826 auto fun = [&](
size_t const k,
size_t const ,
size_t const offset) {
1827 valsView[offset] = values[k]; };
1828 std::function<void(size_t const, size_t const, size_t const)> cb(std::ref(fun));
1829 graph.insertLocalIndicesImpl(lclRow, indices, cb);
1831 std::ostringstream os;
1832 os <<
"You attempted to use insertLocalValues with CombineMode " <<
combineModeToString(CM)
1833 <<
"but this has not been implemented." << endl;
1834 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1835 (
true, std::invalid_argument, os.str ());
1839 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1843 const LocalOrdinal numEnt,
1844 const Scalar vals[],
1845 const LocalOrdinal cols[],
1848 Teuchos::ArrayView<const LocalOrdinal> colsT (cols, numEnt);
1849 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
1850 this->insertLocalValues (localRow, colsT, valsT, CM);
1853 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1858 const GlobalOrdinal gblColInds[],
1860 const size_t numInputEnt)
1862 #ifdef HAVE_TPETRA_DEBUG
1863 const char tfecfFuncName[] =
"insertGlobalValuesImpl: ";
1865 const size_t curNumEnt = rowInfo.numEntries;
1866 #endif // HAVE_TPETRA_DEBUG
1868 if (! graph.indicesAreAllocated ()) {
1871 using ::Tpetra::Details::Behavior;
1872 const bool verbose = Behavior::verbose(
"CrsMatrix");
1873 this->allocateValues (GlobalIndices, GraphNotYetAllocated, verbose);
1878 rowInfo = graph.
getRowInfo (rowInfo.localRow);
1881 auto valsView = this->getValuesViewHostNonConst(rowInfo);
1882 auto fun = [&](
size_t const k,
size_t const ,
size_t const offset){
1883 valsView[offset] += vals[k];
1885 std::function<void(size_t const, size_t const, size_t const)> cb(std::ref(fun));
1886 #ifdef HAVE_TPETRA_DEBUG
1892 #ifdef HAVE_TPETRA_DEBUG
1893 size_t newNumEnt = curNumEnt + numInserted;
1894 const size_t chkNewNumEnt =
1896 if (chkNewNumEnt != newNumEnt) {
1897 std::ostringstream os;
1898 os << std::endl <<
"newNumEnt = " << newNumEnt
1899 <<
" != graph.getNumEntriesInLocalRow(" << rowInfo.localRow
1900 <<
") = " << chkNewNumEnt <<
"." << std::endl
1901 <<
"\torigNumEnt: " << origNumEnt << std::endl
1902 <<
"\tnumInputEnt: " << numInputEnt << std::endl
1903 <<
"\tgblColInds: [";
1904 for (
size_t k = 0; k < numInputEnt; ++k) {
1905 os << gblColInds[k];
1906 if (k +
size_t (1) < numInputEnt) {
1910 os <<
"]" << std::endl
1912 for (
size_t k = 0; k < numInputEnt; ++k) {
1914 if (k +
size_t (1) < numInputEnt) {
1918 os <<
"]" << std::endl;
1920 if (this->supportsRowViews ()) {
1921 values_host_view_type vals2;
1922 if (this->isGloballyIndexed ()) {
1923 global_inds_host_view_type gblColInds2;
1924 const GlobalOrdinal gblRow =
1925 graph.
rowMap_->getGlobalElement (rowInfo.localRow);
1927 Tpetra::Details::OrdinalTraits<GlobalOrdinal>::invalid ()) {
1928 os <<
"Local row index " << rowInfo.localRow <<
" is invalid!"
1932 bool getViewThrew =
false;
1934 this->getGlobalRowView (gblRow, gblColInds2, vals2);
1936 catch (std::exception& e) {
1937 getViewThrew =
true;
1938 os <<
"getGlobalRowView threw exception:" << std::endl
1939 << e.what () << std::endl;
1941 if (! getViewThrew) {
1942 os <<
"\tNew global column indices: ";
1943 for (
size_t jjj = 0; jjj < gblColInds2.extent(0); jjj++)
1944 os << gblColInds2[jjj] <<
" ";
1946 os <<
"\tNew values: ";
1947 for (
size_t jjj = 0; jjj < vals2.extent(0); jjj++)
1948 os << vals2[jjj] <<
" ";
1953 else if (this->isLocallyIndexed ()) {
1954 local_inds_host_view_type lclColInds2;
1955 this->getLocalRowView (rowInfo.localRow, lclColInds2, vals2);
1956 os <<
"\tNew local column indices: ";
1957 for (
size_t jjj = 0; jjj < lclColInds2.extent(0); jjj++)
1958 os << lclColInds2[jjj] <<
" ";
1960 os <<
"\tNew values: ";
1961 for (
size_t jjj = 0; jjj < vals2.extent(0); jjj++)
1962 os << vals2[jjj] <<
" ";
1967 os <<
"Please report this bug to the Tpetra developers.";
1968 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1969 (
true, std::logic_error, os.str ());
1971 #endif // HAVE_TPETRA_DEBUG
1974 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1978 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
1979 const Teuchos::ArrayView<const Scalar>& values)
1981 using Teuchos::toString;
1984 typedef LocalOrdinal LO;
1985 typedef GlobalOrdinal GO;
1986 typedef Tpetra::Details::OrdinalTraits<LO> OTLO;
1987 typedef typename Teuchos::ArrayView<const GO>::size_type size_type;
1988 const char tfecfFuncName[] =
"insertGlobalValues: ";
1990 #ifdef HAVE_TPETRA_DEBUG
1991 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
1992 (values.size () != indices.size (), std::runtime_error,
1993 "values.size() = " << values.size () <<
" != indices.size() = "
1994 << indices.size () <<
".");
1995 #endif // HAVE_TPETRA_DEBUG
1999 const map_type& rowMap = * (this->getCrsGraphRef ().rowMap_);
2002 if (lclRow == OTLO::invalid ()) {
2009 this->insertNonownedGlobalValues (gblRow, indices, values);
2012 if (this->isStaticGraph ()) {
2014 const int myRank = rowMap.getComm ()->getRank ();
2015 const int numProcs = rowMap.getComm ()->getSize ();
2016 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2017 (
true, std::runtime_error,
2018 "The matrix was constructed with a constant (\"static\") graph, "
2019 "yet the given global row index " << gblRow <<
" is in the row "
2020 "Map on the calling process (with rank " << myRank <<
", of " <<
2021 numProcs <<
" process(es)). In this case, you may not insert "
2022 "new entries into rows owned by the calling process.");
2026 const IST*
const inputVals =
2027 reinterpret_cast<const IST*
> (values.getRawPtr ());
2028 const GO*
const inputGblColInds = indices.getRawPtr ();
2029 const size_t numInputEnt = indices.size ();
2038 if (! graph.
colMap_.is_null ()) {
2044 #ifdef HAVE_TPETRA_DEBUG
2045 Teuchos::Array<GO> badColInds;
2046 #endif // HAVE_TPETRA_DEBUG
2047 const size_type numEntriesToInsert = indices.size ();
2048 bool allInColMap =
true;
2049 for (size_type k = 0; k < numEntriesToInsert; ++k) {
2051 allInColMap =
false;
2052 #ifdef HAVE_TPETRA_DEBUG
2053 badColInds.push_back (indices[k]);
2056 #endif // HAVE_TPETRA_DEBUG
2059 if (! allInColMap) {
2060 std::ostringstream os;
2061 os <<
"You attempted to insert entries in owned row " << gblRow
2062 <<
", at the following column indices: " << toString (indices)
2064 #ifdef HAVE_TPETRA_DEBUG
2065 os <<
"Of those, the following indices are not in the column Map "
2066 "on this process: " << toString (badColInds) <<
"." << endl
2067 <<
"Since the matrix has a column Map already, it is invalid "
2068 "to insert entries at those locations.";
2070 os <<
"At least one of those indices is not in the column Map "
2071 "on this process." << endl <<
"It is invalid to insert into "
2072 "columns not in the column Map on the process that owns the "
2074 #endif // HAVE_TPETRA_DEBUG
2075 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2076 (
true, std::invalid_argument, os.str ());
2080 this->insertGlobalValuesImpl (graph, rowInfo, inputGblColInds,
2081 inputVals, numInputEnt);
2086 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2090 const LocalOrdinal numEnt,
2091 const Scalar vals[],
2092 const GlobalOrdinal inds[])
2094 Teuchos::ArrayView<const GlobalOrdinal> indsT (inds, numEnt);
2095 Teuchos::ArrayView<const Scalar> valsT (vals, numEnt);
2096 this->insertGlobalValues (globalRow, indsT, valsT);
2100 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2104 const GlobalOrdinal gblRow,
2105 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2106 const Teuchos::ArrayView<const Scalar>& values,
2109 typedef impl_scalar_type IST;
2110 typedef LocalOrdinal LO;
2111 typedef GlobalOrdinal GO;
2112 typedef Tpetra::Details::OrdinalTraits<LO> OTLO;
2113 const char tfecfFuncName[] =
"insertGlobalValuesFiltered: ";
2116 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2117 (values.size () != indices.size (), std::runtime_error,
2118 "values.size() = " << values.size () <<
" != indices.size() = "
2119 << indices.size () <<
".");
2124 const map_type& rowMap = * (this->getCrsGraphRef ().rowMap_);
2125 const LO lclRow = rowMap.getLocalElement (gblRow);
2126 if (lclRow == OTLO::invalid ()) {
2133 this->insertNonownedGlobalValues (gblRow, indices, values);
2136 if (this->isStaticGraph ()) {
2138 const int myRank = rowMap.getComm ()->getRank ();
2139 const int numProcs = rowMap.getComm ()->getSize ();
2140 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2141 (
true, std::runtime_error,
2142 "The matrix was constructed with a constant (\"static\") graph, "
2143 "yet the given global row index " << gblRow <<
" is in the row "
2144 "Map on the calling process (with rank " << myRank <<
", of " <<
2145 numProcs <<
" process(es)). In this case, you may not insert "
2146 "new entries into rows owned by the calling process.");
2149 crs_graph_type& graph = * (this->myGraph_);
2150 const IST*
const inputVals =
2151 reinterpret_cast<const IST*
> (values.getRawPtr ());
2152 const GO*
const inputGblColInds = indices.getRawPtr ();
2153 const size_t numInputEnt = indices.size ();
2154 RowInfo rowInfo = graph.getRowInfo (lclRow);
2156 if (!graph.colMap_.is_null() && graph.isLocallyIndexed()) {
2163 const map_type& colMap = * (graph.colMap_);
2164 size_t curOffset = 0;
2165 while (curOffset < numInputEnt) {
2169 Teuchos::Array<LO> lclIndices;
2170 size_t endOffset = curOffset;
2171 for ( ; endOffset < numInputEnt; ++endOffset) {
2172 auto lclIndex = colMap.getLocalElement(inputGblColInds[endOffset]);
2173 if (lclIndex != OTLO::invalid())
2174 lclIndices.push_back(lclIndex);
2181 const LO numIndInSeq = (endOffset - curOffset);
2182 if (numIndInSeq != 0) {
2183 this->insertLocalValues(lclRow, lclIndices(), values(curOffset, numIndInSeq));
2189 const bool invariant = endOffset == numInputEnt ||
2190 colMap.getLocalElement (inputGblColInds[endOffset]) == OTLO::invalid ();
2191 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2192 (! invariant, std::logic_error, std::endl <<
"Invariant failed!");
2194 curOffset = endOffset + 1;
2197 else if (! graph.colMap_.is_null ()) {
2198 const map_type& colMap = * (graph.colMap_);
2199 size_t curOffset = 0;
2200 while (curOffset < numInputEnt) {
2204 size_t endOffset = curOffset;
2205 for ( ; endOffset < numInputEnt &&
2206 colMap.getLocalElement (inputGblColInds[endOffset]) != OTLO::invalid ();
2212 const LO numIndInSeq = (endOffset - curOffset);
2213 if (numIndInSeq != 0) {
2214 rowInfo = graph.getRowInfo(lclRow);
2215 this->insertGlobalValuesImpl (graph, rowInfo,
2216 inputGblColInds + curOffset,
2217 inputVals + curOffset,
2224 const bool invariant = endOffset == numInputEnt ||
2225 colMap.getLocalElement (inputGblColInds[endOffset]) == OTLO::invalid ();
2226 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
2227 (! invariant, std::logic_error, std::endl <<
"Invariant failed!");
2229 curOffset = endOffset + 1;
2233 this->insertGlobalValuesImpl (graph, rowInfo, inputGblColInds,
2234 inputVals, numInputEnt);
2239 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2241 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2242 insertGlobalValuesFilteredChecked(
2243 const GlobalOrdinal gblRow,
2244 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2245 const Teuchos::ArrayView<const Scalar>& values,
2246 const char*
const prefix,
2254 insertGlobalValuesFiltered(gblRow, indices, values, debug);
2256 catch(std::exception& e) {
2257 std::ostringstream os;
2259 const size_t maxNumToPrint =
2261 os << *prefix <<
": insertGlobalValuesFiltered threw an "
2262 "exception: " << e.what() << endl
2263 <<
"Global row index: " << gblRow << endl;
2271 os <<
": insertGlobalValuesFiltered threw an exception: "
2274 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, os.str());
2278 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2284 const LocalOrdinal inds[],
2286 const LocalOrdinal numElts)
2288 typedef LocalOrdinal LO;
2289 typedef GlobalOrdinal GO;
2290 const bool sorted = graph.
isSorted ();
2300 for (LO j = 0; j < numElts; ++j) {
2301 const LO lclColInd = inds[j];
2302 const size_t offset =
2303 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2304 lclColInd, hint, sorted);
2305 if (offset != rowInfo.numEntries) {
2306 rowVals[offset] = newVals[j];
2313 if (graph.
colMap_.is_null ()) {
2314 return Teuchos::OrdinalTraits<LO>::invalid ();
2322 for (LO j = 0; j < numElts; ++j) {
2324 if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
2325 const size_t offset =
2326 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2327 gblColInd, hint, sorted);
2328 if (offset != rowInfo.numEntries) {
2329 rowVals[offset] = newVals[j];
2348 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2352 const Teuchos::ArrayView<const LocalOrdinal>& lclCols,
2353 const Teuchos::ArrayView<const Scalar>& vals)
2355 typedef LocalOrdinal LO;
2357 const LO numInputEnt =
static_cast<LO
> (lclCols.size ());
2358 if (static_cast<LO> (vals.size ()) != numInputEnt) {
2359 return Teuchos::OrdinalTraits<LO>::invalid ();
2361 const LO*
const inputInds = lclCols.getRawPtr ();
2362 const Scalar*
const inputVals = vals.getRawPtr ();
2363 return this->replaceLocalValues (localRow, numInputEnt,
2364 inputVals, inputInds);
2367 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2373 const Kokkos::View<const local_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
2374 const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals)
2377 const LO numInputEnt = inputInds.extent(0);
2378 if (numInputEnt != static_cast<LO>(inputVals.extent(0))) {
2379 return Teuchos::OrdinalTraits<LO>::invalid();
2381 const Scalar*
const inVals =
2382 reinterpret_cast<const Scalar*
>(inputVals.data());
2383 return this->replaceLocalValues(localRow, numInputEnt,
2384 inVals, inputInds.data());
2387 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2391 const LocalOrdinal numEnt,
2392 const Scalar inputVals[],
2393 const LocalOrdinal inputCols[])
2396 typedef LocalOrdinal LO;
2398 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2400 return Teuchos::OrdinalTraits<LO>::invalid ();
2405 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2408 return static_cast<LO
> (0);
2410 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
2411 const IST*
const inVals =
reinterpret_cast<const IST*
> (inputVals);
2412 return this->replaceLocalValuesImpl (curRowVals.data (), graph, rowInfo,
2413 inputCols, inVals, numEnt);
2416 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2422 const GlobalOrdinal inds[],
2424 const LocalOrdinal numElts)
2426 Teuchos::ArrayView<const GlobalOrdinal> indsT(inds, numElts);
2428 [&](
size_t const k,
size_t const ,
size_t const offset) {
2429 rowVals[offset] = newVals[k];
2431 std::function<void(size_t const, size_t const, size_t const)> cb(std::ref(fun));
2435 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2439 const Teuchos::ArrayView<const GlobalOrdinal>& inputGblColInds,
2440 const Teuchos::ArrayView<const Scalar>& inputVals)
2442 typedef LocalOrdinal LO;
2444 const LO numInputEnt =
static_cast<LO
> (inputGblColInds.size ());
2445 if (static_cast<LO> (inputVals.size ()) != numInputEnt) {
2446 return Teuchos::OrdinalTraits<LO>::invalid ();
2448 return this->replaceGlobalValues (globalRow, numInputEnt,
2449 inputVals.getRawPtr (),
2450 inputGblColInds.getRawPtr ());
2453 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2457 const LocalOrdinal numEnt,
2458 const Scalar inputVals[],
2459 const GlobalOrdinal inputGblColInds[])
2462 typedef LocalOrdinal LO;
2464 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2466 return Teuchos::OrdinalTraits<LO>::invalid ();
2471 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2474 return static_cast<LO
> (0);
2477 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
2478 const IST*
const inVals =
reinterpret_cast<const IST*
> (inputVals);
2479 return this->replaceGlobalValuesImpl (curRowVals.data (), graph, rowInfo,
2480 inputGblColInds, inVals, numEnt);
2483 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2489 const Kokkos::View<const global_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
2490 const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals)
2499 const LO numInputEnt =
static_cast<LO
>(inputInds.extent(0));
2500 if (static_cast<LO>(inputVals.extent(0)) != numInputEnt) {
2501 return Teuchos::OrdinalTraits<LO>::invalid();
2503 const Scalar*
const inVals =
2504 reinterpret_cast<const Scalar*
>(inputVals.data());
2505 return this->replaceGlobalValues(globalRow, numInputEnt, inVals,
2509 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2515 const GlobalOrdinal inds[],
2517 const LocalOrdinal numElts,
2520 typedef LocalOrdinal LO;
2521 typedef GlobalOrdinal GO;
2523 const bool sorted = graph.
isSorted ();
2532 if (graph.
colMap_.is_null ()) {
2543 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
2545 for (LO j = 0; j < numElts; ++j) {
2547 if (lclColInd != LINV) {
2548 const size_t offset =
2549 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2550 lclColInd, hint, sorted);
2551 if (offset != rowInfo.numEntries) {
2553 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2556 rowVals[offset] += newVals[j];
2569 for (LO j = 0; j < numElts; ++j) {
2570 const GO gblColInd = inds[j];
2571 const size_t offset =
2572 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2573 gblColInd, hint, sorted);
2574 if (offset != rowInfo.numEntries) {
2576 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2579 rowVals[offset] += newVals[j];
2593 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2597 const Teuchos::ArrayView<const GlobalOrdinal>& inputGblColInds,
2598 const Teuchos::ArrayView<const Scalar>& inputVals,
2601 typedef LocalOrdinal LO;
2603 const LO numInputEnt =
static_cast<LO
> (inputGblColInds.size ());
2604 if (static_cast<LO> (inputVals.size ()) != numInputEnt) {
2605 return Teuchos::OrdinalTraits<LO>::invalid ();
2607 return this->sumIntoGlobalValues (gblRow, numInputEnt,
2608 inputVals.getRawPtr (),
2609 inputGblColInds.getRawPtr (),
2613 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2617 const LocalOrdinal numInputEnt,
2618 const Scalar inputVals[],
2619 const GlobalOrdinal inputGblColInds[],
2623 typedef LocalOrdinal LO;
2624 typedef GlobalOrdinal GO;
2626 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2628 return Teuchos::OrdinalTraits<LO>::invalid ();
2633 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2638 using Teuchos::ArrayView;
2639 ArrayView<const GO> inputGblColInds_av(
2640 numInputEnt == 0 ?
nullptr : inputGblColInds,
2642 ArrayView<const Scalar> inputVals_av(
2643 numInputEnt == 0 ?
nullptr :
2644 inputVals, numInputEnt);
2649 this->insertNonownedGlobalValues (gblRow, inputGblColInds_av,
2660 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
2661 const IST*
const inVals =
reinterpret_cast<const IST*
> (inputVals);
2662 return this->sumIntoGlobalValuesImpl (curRowVals.data (), graph, rowInfo,
2663 inputGblColInds, inVals,
2664 numInputEnt, atomic);
2668 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2672 const LocalOrdinal numInputEnt,
2673 const impl_scalar_type inputVals[],
2674 const LocalOrdinal inputCols[],
2675 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2678 using Tpetra::Details::OrdinalTraits;
2679 typedef LocalOrdinal LO;
2681 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2683 return Teuchos::OrdinalTraits<LO>::invalid ();
2685 const crs_graph_type& graph = * (this->staticGraph_);
2686 const RowInfo rowInfo = graph.getRowInfo (lclRow);
2688 if (rowInfo.localRow == OrdinalTraits<size_t>::invalid ()) {
2691 return static_cast<LO
> (0);
2693 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
2694 return this->transformLocalValues (curRowVals.data (), graph,
2695 rowInfo, inputCols, inputVals,
2696 numInputEnt, f, atomic);
2699 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2701 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2702 transformGlobalValues (
const GlobalOrdinal gblRow,
2703 const LocalOrdinal numInputEnt,
2704 const impl_scalar_type inputVals[],
2705 const GlobalOrdinal inputCols[],
2706 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2709 using Tpetra::Details::OrdinalTraits;
2710 typedef LocalOrdinal LO;
2712 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
2714 return OrdinalTraits<LO>::invalid ();
2716 const crs_graph_type& graph = * (this->staticGraph_);
2717 const RowInfo rowInfo = graph.getRowInfoFromGlobalRowIndex (gblRow);
2719 if (rowInfo.localRow == OrdinalTraits<size_t>::invalid ()) {
2722 return static_cast<LO
> (0);
2724 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
2725 return this->transformGlobalValues (curRowVals.data (), graph,
2726 rowInfo, inputCols, inputVals,
2727 numInputEnt, f, atomic);
2730 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2732 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2733 transformLocalValues (impl_scalar_type rowVals[],
2734 const crs_graph_type& graph,
2735 const RowInfo& rowInfo,
2736 const LocalOrdinal inds[],
2737 const impl_scalar_type newVals[],
2738 const LocalOrdinal numElts,
2739 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2742 typedef impl_scalar_type ST;
2743 typedef LocalOrdinal LO;
2744 typedef GlobalOrdinal GO;
2751 const bool sorted = graph.isSorted ();
2756 if (graph.isLocallyIndexed ()) {
2759 auto colInds = graph.getLocalIndsViewHost (rowInfo);
2761 for (LO j = 0; j < numElts; ++j) {
2762 const LO lclColInd = inds[j];
2763 const size_t offset =
2764 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2765 lclColInd, hint, sorted);
2766 if (offset != rowInfo.numEntries) {
2775 volatile ST*
const dest = &rowVals[offset];
2776 (void) atomic_binary_function_update (dest, newVals[j], f);
2780 rowVals[offset] = f (rowVals[offset], newVals[j]);
2787 else if (graph.isGloballyIndexed ()) {
2791 if (graph.colMap_.is_null ()) {
2798 const map_type& colMap = * (graph.colMap_);
2801 auto colInds = graph.getGlobalIndsViewHost (rowInfo);
2803 const GO GINV = Teuchos::OrdinalTraits<GO>::invalid ();
2804 for (LO j = 0; j < numElts; ++j) {
2805 const GO gblColInd = colMap.getGlobalElement (inds[j]);
2806 if (gblColInd != GINV) {
2807 const size_t offset =
2808 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2809 gblColInd, hint, sorted);
2810 if (offset != rowInfo.numEntries) {
2819 volatile ST*
const dest = &rowVals[offset];
2820 (void) atomic_binary_function_update (dest, newVals[j], f);
2824 rowVals[offset] = f (rowVals[offset], newVals[j]);
2839 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2841 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
2842 transformGlobalValues (impl_scalar_type rowVals[],
2843 const crs_graph_type& graph,
2844 const RowInfo& rowInfo,
2845 const GlobalOrdinal inds[],
2846 const impl_scalar_type newVals[],
2847 const LocalOrdinal numElts,
2848 std::function<impl_scalar_type (
const impl_scalar_type&,
const impl_scalar_type&) > f,
2851 typedef impl_scalar_type ST;
2852 typedef LocalOrdinal LO;
2853 typedef GlobalOrdinal GO;
2860 const bool sorted = graph.isSorted ();
2865 if (graph.isGloballyIndexed ()) {
2868 auto colInds = graph.getGlobalIndsViewHost (rowInfo);
2870 for (LO j = 0; j < numElts; ++j) {
2871 const GO gblColInd = inds[j];
2872 const size_t offset =
2873 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2874 gblColInd, hint, sorted);
2875 if (offset != rowInfo.numEntries) {
2884 volatile ST*
const dest = &rowVals[offset];
2885 (void) atomic_binary_function_update (dest, newVals[j], f);
2889 rowVals[offset] = f (rowVals[offset], newVals[j]);
2896 else if (graph.isLocallyIndexed ()) {
2900 if (graph.colMap_.is_null ()) {
2906 const map_type& colMap = * (graph.colMap_);
2909 auto colInds = graph.getLocalIndsViewHost (rowInfo);
2911 const LO LINV = Teuchos::OrdinalTraits<LO>::invalid ();
2912 for (LO j = 0; j < numElts; ++j) {
2913 const LO lclColInd = colMap.getLocalElement (inds[j]);
2914 if (lclColInd != LINV) {
2915 const size_t offset =
2916 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2917 lclColInd, hint, sorted);
2918 if (offset != rowInfo.numEntries) {
2927 volatile ST*
const dest = &rowVals[offset];
2928 (void) atomic_binary_function_update (dest, newVals[j], f);
2932 rowVals[offset] = f (rowVals[offset], newVals[j]);
2947 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2953 const LocalOrdinal inds[],
2955 const LocalOrdinal numElts,
2958 typedef LocalOrdinal LO;
2959 typedef GlobalOrdinal GO;
2961 const bool sorted = graph.
isSorted ();
2971 for (LO j = 0; j < numElts; ++j) {
2972 const LO lclColInd = inds[j];
2973 const size_t offset =
2974 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
2975 lclColInd, hint, sorted);
2976 if (offset != rowInfo.numEntries) {
2978 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
2981 rowVals[offset] += newVals[j];
2989 if (graph.
colMap_.is_null ()) {
2990 return Teuchos::OrdinalTraits<LO>::invalid ();
2998 for (LO j = 0; j < numElts; ++j) {
3000 if (gblColInd != Teuchos::OrdinalTraits<GO>::invalid ()) {
3001 const size_t offset =
3002 KokkosSparse::findRelOffset (colInds, rowInfo.numEntries,
3003 gblColInd, hint, sorted);
3004 if (offset != rowInfo.numEntries) {
3006 Kokkos::atomic_add (&rowVals[offset], newVals[j]);
3009 rowVals[offset] += newVals[j];
3029 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3033 const Teuchos::ArrayView<const LocalOrdinal>& indices,
3034 const Teuchos::ArrayView<const Scalar>& values,
3038 const LO numInputEnt =
static_cast<LO
>(indices.size());
3039 if (static_cast<LO>(values.size()) != numInputEnt) {
3040 return Teuchos::OrdinalTraits<LO>::invalid();
3042 const LO*
const inputInds = indices.getRawPtr();
3043 const scalar_type*
const inputVals = values.getRawPtr();
3044 return this->sumIntoLocalValues(localRow, numInputEnt,
3045 inputVals, inputInds, atomic);
3048 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3054 const Kokkos::View<const local_ordinal_type*, Kokkos::AnonymousSpace>& inputInds,
3055 const Kokkos::View<const impl_scalar_type*, Kokkos::AnonymousSpace>& inputVals,
3059 const LO numInputEnt =
static_cast<LO
>(inputInds.extent(0));
3060 if (static_cast<LO>(inputVals.extent(0)) != numInputEnt) {
3061 return Teuchos::OrdinalTraits<LO>::invalid();
3064 reinterpret_cast<const scalar_type*
>(inputVals.data());
3065 return this->sumIntoLocalValues(localRow, numInputEnt, inVals,
3066 inputInds.data(), atomic);
3069 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3073 const LocalOrdinal numEnt,
3074 const Scalar vals[],
3075 const LocalOrdinal cols[],
3079 typedef LocalOrdinal LO;
3081 if (! this->isFillActive () || this->staticGraph_.is_null ()) {
3083 return Teuchos::OrdinalTraits<LO>::invalid ();
3088 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
3091 return static_cast<LO
> (0);
3093 auto curRowVals = this->getValuesViewHostNonConst (rowInfo);
3094 const IST*
const inputVals =
reinterpret_cast<const IST*
> (vals);
3095 return this->sumIntoLocalValuesImpl (curRowVals.data (), graph, rowInfo,
3096 cols, inputVals, numEnt, atomic);
3099 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3101 values_dualv_type::t_host::const_type
3105 if (rowinfo.allocSize == 0 || valuesUnpacked_wdv.extent(0) == 0)
3106 return typename values_dualv_type::t_host::const_type ();
3108 return valuesUnpacked_wdv.getHostSubview(rowinfo.offset1D,
3113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3115 values_dualv_type::t_host
3119 if (rowinfo.allocSize == 0 || valuesUnpacked_wdv.extent(0) == 0)
3120 return typename values_dualv_type::t_host ();
3122 return valuesUnpacked_wdv.getHostSubview(rowinfo.offset1D,
3127 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3129 values_dualv_type::t_dev::const_type
3133 if (rowinfo.allocSize == 0 || valuesUnpacked_wdv.extent(0) == 0)
3134 return typename values_dualv_type::t_dev::const_type ();
3136 return valuesUnpacked_wdv.getDeviceSubview(rowinfo.offset1D,
3141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3143 values_dualv_type::t_dev
3147 if (rowinfo.allocSize == 0 || valuesUnpacked_wdv.extent(0) == 0)
3148 return typename values_dualv_type::t_dev ();
3150 return valuesUnpacked_wdv.getDeviceSubview(rowinfo.offset1D,
3156 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3160 nonconst_local_inds_host_view_type &indices,
3161 nonconst_values_host_view_type &values,
3162 size_t& numEntries)
const
3164 using Teuchos::ArrayView;
3165 using Teuchos::av_reinterpret_cast;
3166 const char tfecfFuncName[] =
"getLocalRowCopy: ";
3168 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3169 (! this->hasColMap (), std::runtime_error,
3170 "The matrix does not have a column Map yet. This means we don't have "
3171 "local indices for columns yet, so it doesn't make sense to call this "
3172 "method. If the matrix doesn't have a column Map yet, you should call "
3173 "fillComplete on it first.");
3175 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
3176 const size_t theNumEntries = rowinfo.numEntries;
3177 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3178 (static_cast<size_t> (indices.size ()) < theNumEntries ||
3179 static_cast<size_t> (values.size ()) < theNumEntries,
3180 std::runtime_error,
"Row with local index " << localRow <<
" has " <<
3181 theNumEntries <<
" entry/ies, but indices.size() = " <<
3182 indices.size () <<
" and values.size() = " << values.size () <<
".");
3183 numEntries = theNumEntries;
3185 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
3186 if (staticGraph_->isLocallyIndexed ()) {
3187 auto curLclInds = staticGraph_->getLocalIndsViewHost(rowinfo);
3188 auto curVals = getValuesViewHost(rowinfo);
3190 for (
size_t j = 0; j < theNumEntries; ++j) {
3191 values[j] = curVals[j];
3192 indices[j] = curLclInds(j);
3195 else if (staticGraph_->isGloballyIndexed ()) {
3197 const map_type& colMap = * (staticGraph_->colMap_);
3198 auto curGblInds = staticGraph_->getGlobalIndsViewHost(rowinfo);
3199 auto curVals = getValuesViewHost(rowinfo);
3201 for (
size_t j = 0; j < theNumEntries; ++j) {
3202 values[j] = curVals[j];
3210 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3214 nonconst_global_inds_host_view_type &indices,
3215 nonconst_values_host_view_type &values,
3216 size_t& numEntries)
const
3218 using Teuchos::ArrayView;
3219 using Teuchos::av_reinterpret_cast;
3220 const char tfecfFuncName[] =
"getGlobalRowCopy: ";
3223 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
3224 const size_t theNumEntries = rowinfo.numEntries;
3225 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3226 static_cast<size_t> (indices.size ()) < theNumEntries ||
3227 static_cast<size_t> (values.size ()) < theNumEntries,
3228 std::runtime_error,
"Row with global index " << globalRow <<
" has "
3229 << theNumEntries <<
" entry/ies, but indices.size() = " <<
3230 indices.size () <<
" and values.size() = " << values.size () <<
".");
3231 numEntries = theNumEntries;
3233 if (rowinfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid ()) {
3234 if (staticGraph_->isLocallyIndexed ()) {
3235 const map_type& colMap = * (staticGraph_->colMap_);
3236 auto curLclInds = staticGraph_->getLocalIndsViewHost(rowinfo);
3237 auto curVals = getValuesViewHost(rowinfo);
3239 for (
size_t j = 0; j < theNumEntries; ++j) {
3240 values[j] = curVals[j];
3244 else if (staticGraph_->isGloballyIndexed ()) {
3245 auto curGblInds = staticGraph_->getGlobalIndsViewHost(rowinfo);
3246 auto curVals = getValuesViewHost(rowinfo);
3248 for (
size_t j = 0; j < theNumEntries; ++j) {
3249 values[j] = curVals[j];
3250 indices[j] = curGblInds(j);
3257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3261 local_inds_host_view_type &indices,
3262 values_host_view_type &values)
const
3264 const char tfecfFuncName[] =
"getLocalRowView: ";
3266 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3267 isGloballyIndexed (), std::runtime_error,
"The matrix currently stores "
3268 "its indices as global indices, so you cannot get a view with local "
3269 "column indices. If the matrix has a column Map, you may call "
3270 "getLocalRowCopy() to get local column indices; otherwise, you may get "
3271 "a view with global column indices by calling getGlobalRowCopy().");
3273 const RowInfo rowInfo = staticGraph_->getRowInfo (localRow);
3274 if (rowInfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
3275 rowInfo.numEntries > 0) {
3276 indices = staticGraph_->lclIndsUnpacked_wdv.getHostSubview(
3280 values = valuesUnpacked_wdv.getHostSubview(rowInfo.offset1D,
3287 indices = local_inds_host_view_type();
3288 values = values_host_view_type();
3291 #ifdef HAVE_TPETRA_DEBUG
3292 const char suffix[] =
". This should never happen. Please report this "
3293 "bug to the Tpetra developers.";
3294 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3295 (static_cast<size_t> (indices.size ()) !=
3296 static_cast<size_t> (values.size ()), std::logic_error,
3297 "At the end of this method, for local row " << localRow <<
", "
3298 "indices.size() = " << indices.size () <<
" != values.size () = "
3299 << values.size () << suffix);
3300 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3301 (static_cast<size_t> (indices.size ()) !=
3302 static_cast<size_t> (rowInfo.numEntries), std::logic_error,
3303 "At the end of this method, for local row " << localRow <<
", "
3304 "indices.size() = " << indices.size () <<
" != rowInfo.numEntries = "
3305 << rowInfo.numEntries << suffix);
3306 const size_t expectedNumEntries = getNumEntriesInLocalRow (localRow);
3307 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3308 (rowInfo.numEntries != expectedNumEntries, std::logic_error,
"At the end "
3309 "of this method, for local row " << localRow <<
", rowInfo.numEntries = "
3310 << rowInfo.numEntries <<
" != getNumEntriesInLocalRow(localRow) = " <<
3311 expectedNumEntries << suffix);
3312 #endif // HAVE_TPETRA_DEBUG
3316 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3320 global_inds_host_view_type &indices,
3321 values_host_view_type &values)
const
3323 const char tfecfFuncName[] =
"getGlobalRowView: ";
3325 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3326 isLocallyIndexed (), std::runtime_error,
3327 "The matrix is locally indexed, so we cannot return a view of the row "
3328 "with global column indices. Use getGlobalRowCopy() instead.");
3333 staticGraph_->getRowInfoFromGlobalRowIndex (globalRow);
3334 if (rowInfo.localRow != Teuchos::OrdinalTraits<size_t>::invalid () &&
3335 rowInfo.numEntries > 0) {
3336 indices = staticGraph_->gblInds_wdv.getHostSubview(rowInfo.offset1D,
3339 values = valuesUnpacked_wdv.getHostSubview(rowInfo.offset1D,
3344 indices = global_inds_host_view_type();
3345 values = values_host_view_type();
3348 #ifdef HAVE_TPETRA_DEBUG
3349 const char suffix[] =
". This should never happen. Please report this "
3350 "bug to the Tpetra developers.";
3351 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3352 (static_cast<size_t> (indices.size ()) !=
3353 static_cast<size_t> (values.size ()), std::logic_error,
3354 "At the end of this method, for global row " << globalRow <<
", "
3355 "indices.size() = " << indices.size () <<
" != values.size () = "
3356 << values.size () << suffix);
3357 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3358 (static_cast<size_t> (indices.size ()) !=
3359 static_cast<size_t> (rowInfo.numEntries), std::logic_error,
3360 "At the end of this method, for global row " << globalRow <<
", "
3361 "indices.size() = " << indices.size () <<
" != rowInfo.numEntries = "
3362 << rowInfo.numEntries << suffix);
3363 const size_t expectedNumEntries = getNumEntriesInGlobalRow (globalRow);
3364 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3365 (rowInfo.numEntries != expectedNumEntries, std::logic_error,
"At the end "
3366 "of this method, for global row " << globalRow <<
", rowInfo.numEntries "
3367 "= " << rowInfo.numEntries <<
" != getNumEntriesInGlobalRow(globalRow) ="
3368 " " << expectedNumEntries << suffix);
3369 #endif // HAVE_TPETRA_DEBUG
3373 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3380 const size_t nlrs = staticGraph_->getLocalNumRows ();
3381 const size_t numEntries = staticGraph_->getLocalNumEntries ();
3382 if (! staticGraph_->indicesAreAllocated () ||
3383 nlrs == 0 || numEntries == 0) {
3388 auto vals = valuesPacked_wdv.getDeviceView(Access::ReadWrite);
3389 KokkosBlas::scal(vals, theAlpha, vals);
3394 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3405 const size_t numEntries = staticGraph_->getLocalNumEntries();
3406 if (! staticGraph_->indicesAreAllocated () || numEntries == 0) {
3414 Kokkos::fence(
"CrsMatrix::setAllToScalar");
3418 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3421 setAllValues (
const typename local_graph_device_type::row_map_type& rowPointers,
3422 const typename local_graph_device_type::entries_type::non_const_type& columnIndices,
3423 const typename local_matrix_device_type::values_type& values)
3426 ProfilingRegion region (
"Tpetra::CrsMatrix::setAllValues");
3427 const char tfecfFuncName[] =
"setAllValues: ";
3428 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3429 (columnIndices.size () != values.size (), std::invalid_argument,
3430 "columnIndices.size() = " << columnIndices.size () <<
" != values.size()"
3431 " = " << values.size () <<
".");
3432 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3433 (myGraph_.is_null (), std::runtime_error,
"myGraph_ must not be null.");
3436 myGraph_->setAllIndices (rowPointers, columnIndices);
3438 catch (std::exception &e) {
3439 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3440 (
true, std::runtime_error,
"myGraph_->setAllIndices() threw an "
3441 "exception: " << e.what ());
3448 auto lclGraph = myGraph_->getLocalGraphDevice ();
3449 const size_t numEnt = lclGraph.entries.extent (0);
3450 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3451 (lclGraph.row_map.extent (0) != rowPointers.extent (0) ||
3452 numEnt !=
static_cast<size_t> (columnIndices.extent (0)),
3453 std::logic_error,
"myGraph_->setAllIndices() did not correctly create "
3454 "local graph. Please report this bug to the Tpetra developers.");
3457 valuesUnpacked_wdv = valuesPacked_wdv;
3461 this->storageStatus_ = Details::STORAGE_1D_PACKED;
3463 checkInternalState ();
3466 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3472 ProfilingRegion region (
"Tpetra::CrsMatrix::setAllValues from KokkosSparse::CrsMatrix");
3474 auto graph = localDeviceMatrix.graph;
3477 auto rows = graph.row_map;
3478 auto columns = graph.entries;
3479 auto values = localDeviceMatrix.values;
3481 setAllValues(rows,columns,values);
3484 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3488 const Teuchos::ArrayRCP<LocalOrdinal>& ind,
3489 const Teuchos::ArrayRCP<Scalar>& val)
3491 using Kokkos::Compat::getKokkosViewDeepCopy;
3492 using Teuchos::ArrayRCP;
3493 using Teuchos::av_reinterpret_cast;
3496 typedef typename local_graph_device_type::row_map_type row_map_type;
3498 const char tfecfFuncName[] =
"setAllValues(ArrayRCP<size_t>, ArrayRCP<LO>, ArrayRCP<Scalar>): ";
3504 typename row_map_type::non_const_type ptrNative (
"ptr", ptr.size ());
3505 Kokkos::View<
const size_t*,
3506 typename row_map_type::array_layout,
3508 Kokkos::MemoryUnmanaged> ptrSizeT (ptr.getRawPtr (), ptr.size ());
3511 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3512 (ptrNative.extent (0) != ptrSizeT.extent (0),
3513 std::logic_error,
"ptrNative.extent(0) = " <<
3514 ptrNative.extent (0) <<
" != ptrSizeT.extent(0) = "
3515 << ptrSizeT.extent (0) <<
". Please report this bug to the "
3516 "Tpetra developers.");
3518 auto indIn = getKokkosViewDeepCopy<DT> (ind ());
3519 auto valIn = getKokkosViewDeepCopy<DT> (av_reinterpret_cast<IST> (val ()));
3520 this->setAllValues (ptrNative, indIn, valIn);
3523 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3528 const char tfecfFuncName[] =
"getLocalDiagOffsets: ";
3529 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3530 (staticGraph_.is_null (), std::runtime_error,
"The matrix has no graph.");
3537 const size_t lclNumRows = staticGraph_->getLocalNumRows ();
3538 if (static_cast<size_t> (offsets.size ()) < lclNumRows) {
3539 offsets.resize (lclNumRows);
3545 if (std::is_same<memory_space, Kokkos::HostSpace>::value) {
3550 Kokkos::MemoryUnmanaged> output_type;
3551 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3552 staticGraph_->getLocalDiagOffsets (offsetsOut);
3555 Kokkos::View<size_t*, device_type> offsetsTmp (
"diagOffsets", lclNumRows);
3556 staticGraph_->getLocalDiagOffsets (offsetsTmp);
3557 typedef Kokkos::View<
size_t*, Kokkos::HostSpace,
3558 Kokkos::MemoryUnmanaged> output_type;
3559 output_type offsetsOut (offsets.getRawPtr (), lclNumRows);
3565 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3570 using Teuchos::ArrayRCP;
3571 using Teuchos::ArrayView;
3572 using Teuchos::av_reinterpret_cast;
3573 const char tfecfFuncName[] =
"getLocalDiagCopy (1-arg): ";
3577 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3578 staticGraph_.is_null (), std::runtime_error,
3579 "This method requires that the matrix have a graph.");
3580 auto rowMapPtr = this->getRowMap ();
3581 if (rowMapPtr.is_null () || rowMapPtr->getComm ().is_null ()) {
3587 auto colMapPtr = this->getColMap ();
3588 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3589 (! this->hasColMap () || colMapPtr.is_null (), std::runtime_error,
3590 "This method requires that the matrix have a column Map.");
3591 const map_type& rowMap = * rowMapPtr;
3592 const map_type& colMap = * colMapPtr;
3593 const LO myNumRows =
static_cast<LO
> (this->getLocalNumRows ());
3595 #ifdef HAVE_TPETRA_DEBUG
3598 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3599 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3600 "The input Vector's Map must be compatible with the CrsMatrix's row "
3601 "Map. You may check this by using Map's isCompatible method: "
3602 "diag.getMap ()->isCompatible (A.getRowMap ());");
3603 #endif // HAVE_TPETRA_DEBUG
3605 if (this->isFillComplete ()) {
3608 const auto D_lcl_1d =
3609 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3611 const auto lclRowMap = rowMap.getLocalMap ();
3616 getLocalMatrixDevice ());
3624 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3629 Kokkos::MemoryUnmanaged>& offsets)
const
3631 typedef LocalOrdinal LO;
3633 #ifdef HAVE_TPETRA_DEBUG
3634 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3635 const map_type& rowMap = * (this->getRowMap ());
3638 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3639 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3640 "The input Vector's Map must be compatible with (in the sense of Map::"
3641 "isCompatible) the CrsMatrix's row Map.");
3642 #endif // HAVE_TPETRA_DEBUG
3652 const LO myNumRows =
static_cast<LO
> (this->getLocalNumRows ());
3655 Kokkos::subview (D_lcl, Kokkos::make_pair (LO (0), myNumRows), 0);
3657 KokkosSparse::getDiagCopy (D_lcl_1d, offsets,
3658 getLocalMatrixDevice ());
3661 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3665 const Teuchos::ArrayView<const size_t>& offsets)
const
3667 using LO = LocalOrdinal;
3668 using host_execution_space = Kokkos::DefaultHostExecutionSpace;
3671 #ifdef HAVE_TPETRA_DEBUG
3672 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
3673 const map_type& rowMap = * (this->getRowMap ());
3676 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3677 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
3678 "The input Vector's Map must be compatible with (in the sense of Map::"
3679 "isCompatible) the CrsMatrix's row Map.");
3680 #endif // HAVE_TPETRA_DEBUG
3692 auto lclVecHost1d = Kokkos::subview (lclVecHost, Kokkos::ALL (), 0);
3694 using host_offsets_view_type =
3695 Kokkos::View<
const size_t*, Kokkos::HostSpace,
3696 Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
3697 host_offsets_view_type h_offsets (offsets.getRawPtr (), offsets.size ());
3699 using range_type = Kokkos::RangePolicy<host_execution_space, LO>;
3700 const LO myNumRows =
static_cast<LO
> (this->getLocalNumRows ());
3701 const size_t INV = Tpetra::Details::OrdinalTraits<size_t>::invalid ();
3703 auto rowPtrsPackedHost = staticGraph_->getRowPtrsPackedHost();
3704 auto valuesPackedHost = valuesPacked_wdv.getHostView(Access::ReadOnly);
3705 Kokkos::parallel_for
3706 (
"Tpetra::CrsMatrix::getLocalDiagCopy",
3707 range_type (0, myNumRows),
3708 [&, INV, h_offsets] (
const LO lclRow) {
3709 lclVecHost1d(lclRow) = STS::zero ();
3710 if (h_offsets[lclRow] != INV) {
3711 auto curRowOffset = rowPtrsPackedHost (lclRow);
3712 lclVecHost1d(lclRow) =
3713 static_cast<IST
> (valuesPackedHost(curRowOffset+h_offsets[lclRow]));
3720 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3725 using ::Tpetra::Details::ProfilingRegion;
3726 using Teuchos::ArrayRCP;
3727 using Teuchos::ArrayView;
3728 using Teuchos::null;
3731 using Teuchos::rcpFromRef;
3733 const char tfecfFuncName[] =
"leftScale: ";
3735 ProfilingRegion region (
"Tpetra::CrsMatrix::leftScale");
3737 RCP<const vec_type> xp;
3738 if (this->getRangeMap ()->isSameAs (* (x.
getMap ()))) {
3741 auto exporter = this->getCrsGraphRef ().getExporter ();
3742 if (exporter.get () !=
nullptr) {
3743 RCP<vec_type> tempVec (
new vec_type (this->getRowMap ()));
3744 tempVec->doImport (x, *exporter,
REPLACE);
3748 xp = rcpFromRef (x);
3751 else if (this->getRowMap ()->isSameAs (* (x.
getMap ()))) {
3752 xp = rcpFromRef (x);
3755 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3756 (
true, std::invalid_argument,
"x's Map must be the same as "
3757 "either the row Map or the range Map of the CrsMatrix.");
3760 if (this->isFillComplete()) {
3761 auto x_lcl = xp->getLocalViewDevice (Access::ReadOnly);
3762 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
3765 x_lcl_1d,
false,
false);
3769 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3770 (
true, std::runtime_error,
"CrsMatrix::leftScale requires matrix to be"
3775 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3780 using ::Tpetra::Details::ProfilingRegion;
3781 using Teuchos::ArrayRCP;
3782 using Teuchos::ArrayView;
3783 using Teuchos::null;
3786 using Teuchos::rcpFromRef;
3788 const char tfecfFuncName[] =
"rightScale: ";
3790 ProfilingRegion region (
"Tpetra::CrsMatrix::rightScale");
3792 RCP<const vec_type> xp;
3793 if (this->getDomainMap ()->isSameAs (* (x.
getMap ()))) {
3796 auto importer = this->getCrsGraphRef ().getImporter ();
3797 if (importer.get () !=
nullptr) {
3798 RCP<vec_type> tempVec (
new vec_type (this->getColMap ()));
3799 tempVec->doImport (x, *importer,
REPLACE);
3803 xp = rcpFromRef (x);
3806 else if (this->getColMap ()->isSameAs (* (x.
getMap ()))) {
3807 xp = rcpFromRef (x);
3809 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3810 (
true, std::runtime_error,
"x's Map must be the same as "
3811 "either the domain Map or the column Map of the CrsMatrix.");
3814 if (this->isFillComplete()) {
3815 auto x_lcl = xp->getLocalViewDevice (Access::ReadOnly);
3816 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
3819 x_lcl_1d,
false,
false);
3823 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
3824 (
true, std::runtime_error,
"CrsMatrix::rightScale requires matrix to be"
3829 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3834 using Teuchos::ArrayView;
3835 using Teuchos::outArg;
3836 using Teuchos::REDUCE_SUM;
3837 using Teuchos::reduceAll;
3845 if (getLocalNumEntries() > 0) {
3846 if (isStorageOptimized ()) {
3849 const size_t numEntries = getLocalNumEntries ();
3850 auto values = valuesPacked_wdv.getHostView(Access::ReadOnly);
3851 for (
size_t k = 0; k < numEntries; ++k) {
3852 auto val = values[k];
3856 const mag_type val_abs = STS::abs (val);
3857 mySum += val_abs * val_abs;
3861 const LocalOrdinal numRows =
3862 static_cast<LocalOrdinal
> (this->getLocalNumRows ());
3863 for (LocalOrdinal r = 0; r < numRows; ++r) {
3864 const RowInfo rowInfo = myGraph_->getRowInfo (r);
3865 const size_t numEntries = rowInfo.numEntries;
3866 auto A_r = this->getValuesViewHost(rowInfo);
3867 for (
size_t k = 0; k < numEntries; ++k) {
3869 const mag_type val_abs = STS::abs (val);
3870 mySum += val_abs * val_abs;
3876 reduceAll<int, mag_type> (* (getComm ()), REDUCE_SUM,
3877 mySum, outArg (totalSum));
3878 return STM::sqrt (totalSum);
3881 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3886 const char tfecfFuncName[] =
"replaceColMap: ";
3890 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3891 myGraph_.is_null (), std::runtime_error,
3892 "This method does not work if the matrix has a const graph. The whole "
3893 "idea of a const graph is that you are not allowed to change it, but "
3894 "this method necessarily must modify the graph, since the graph owns "
3895 "the matrix's column Map.");
3896 myGraph_->replaceColMap (newColMap);
3899 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3903 const Teuchos::RCP<const map_type>& newColMap,
3904 const Teuchos::RCP<const import_type>& newImport,
3905 const bool sortEachRow)
3907 const char tfecfFuncName[] =
"reindexColumns: ";
3908 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3909 graph ==
nullptr && myGraph_.is_null (), std::invalid_argument,
3910 "The input graph is null, but the matrix does not own its graph.");
3912 crs_graph_type& theGraph = (graph ==
nullptr) ? *myGraph_ : *graph;
3913 const bool sortGraph =
false;
3917 if (sortEachRow && theGraph.isLocallyIndexed () && ! theGraph.isSorted ()) {
3918 const LocalOrdinal lclNumRows =
3919 static_cast<LocalOrdinal
> (theGraph.getLocalNumRows ());
3921 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
3923 const RowInfo rowInfo = theGraph.getRowInfo (row);
3924 auto lclColInds = theGraph.getLocalIndsViewHostNonConst (rowInfo);
3925 auto vals = this->getValuesViewHostNonConst (rowInfo);
3927 sort2 (lclColInds.data (),
3928 lclColInds.data () + rowInfo.numEntries,
3931 theGraph.indicesAreSorted_ =
true;
3935 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3940 const char tfecfFuncName[] =
"replaceDomainMap: ";
3941 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3942 myGraph_.is_null (), std::runtime_error,
3943 "This method does not work if the matrix has a const graph. The whole "
3944 "idea of a const graph is that you are not allowed to change it, but this"
3945 " method necessarily must modify the graph, since the graph owns the "
3946 "matrix's domain Map and Import objects.");
3947 myGraph_->replaceDomainMap (newDomainMap);
3950 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3954 Teuchos::RCP<const import_type>& newImporter)
3956 const char tfecfFuncName[] =
"replaceDomainMapAndImporter: ";
3957 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3958 myGraph_.is_null (), std::runtime_error,
3959 "This method does not work if the matrix has a const graph. The whole "
3960 "idea of a const graph is that you are not allowed to change it, but this"
3961 " method necessarily must modify the graph, since the graph owns the "
3962 "matrix's domain Map and Import objects.");
3963 myGraph_->replaceDomainMapAndImporter (newDomainMap, newImporter);
3966 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3971 const char tfecfFuncName[] =
"replaceRangeMap: ";
3972 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3973 myGraph_.is_null (), std::runtime_error,
3974 "This method does not work if the matrix has a const graph. The whole "
3975 "idea of a const graph is that you are not allowed to change it, but this"
3976 " method necessarily must modify the graph, since the graph owns the "
3977 "matrix's domain Map and Import objects.");
3978 myGraph_->replaceRangeMap (newRangeMap);
3981 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
3985 Teuchos::RCP<const export_type>& newExporter)
3987 const char tfecfFuncName[] =
"replaceRangeMapAndExporter: ";
3988 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3989 myGraph_.is_null (), std::runtime_error,
3990 "This method does not work if the matrix has a const graph. The whole "
3991 "idea of a const graph is that you are not allowed to change it, but this"
3992 " method necessarily must modify the graph, since the graph owns the "
3993 "matrix's domain Map and Import objects.");
3994 myGraph_->replaceRangeMapAndExporter (newRangeMap, newExporter);
3997 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4001 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
4002 const Teuchos::ArrayView<const Scalar>& values)
4004 using Teuchos::Array;
4005 typedef GlobalOrdinal GO;
4006 typedef typename Array<GO>::size_type size_type;
4008 const size_type numToInsert = indices.size ();
4011 std::pair<Array<GO>, Array<Scalar> >& curRow = nonlocals_[globalRow];
4012 Array<GO>& curRowInds = curRow.first;
4013 Array<Scalar>& curRowVals = curRow.second;
4014 const size_type newCapacity = curRowInds.size () + numToInsert;
4015 curRowInds.reserve (newCapacity);
4016 curRowVals.reserve (newCapacity);
4017 for (size_type k = 0; k < numToInsert; ++k) {
4018 curRowInds.push_back (indices[k]);
4019 curRowVals.push_back (values[k]);
4023 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4030 using Teuchos::Comm;
4031 using Teuchos::outArg;
4034 using Teuchos::REDUCE_MAX;
4035 using Teuchos::REDUCE_MIN;
4036 using Teuchos::reduceAll;
4040 typedef GlobalOrdinal GO;
4041 typedef typename Teuchos::Array<GO>::size_type size_type;
4042 const char tfecfFuncName[] =
"globalAssemble: ";
4043 ProfilingRegion regionGlobalAssemble (
"Tpetra::CrsMatrix::globalAssemble");
4045 const bool verbose = Behavior::verbose(
"CrsMatrix");
4046 std::unique_ptr<std::string> prefix;
4048 prefix = this->createPrefix(
"CrsMatrix",
"globalAssemble");
4049 std::ostringstream os;
4050 os << *prefix <<
"nonlocals_.size()=" << nonlocals_.size()
4052 std::cerr << os.str();
4054 RCP<const Comm<int> > comm = getComm ();
4056 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4057 (! isFillActive (), std::runtime_error,
"Fill must be active before "
4058 "you may call this method.");
4060 const size_t myNumNonlocalRows = nonlocals_.size ();
4067 const int iHaveNonlocalRows = (myNumNonlocalRows == 0) ? 0 : 1;
4068 int someoneHasNonlocalRows = 0;
4069 reduceAll<int, int> (*comm, REDUCE_MAX, iHaveNonlocalRows,
4070 outArg (someoneHasNonlocalRows));
4071 if (someoneHasNonlocalRows == 0) {
4085 RCP<const map_type> nonlocalRowMap;
4086 Teuchos::Array<size_t> numEntPerNonlocalRow (myNumNonlocalRows);
4088 Teuchos::Array<GO> myNonlocalGblRows (myNumNonlocalRows);
4089 size_type curPos = 0;
4090 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
4091 ++mapIter, ++curPos) {
4092 myNonlocalGblRows[curPos] = mapIter->first;
4095 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4096 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4103 sort2 (gblCols.begin (), gblCols.end (), vals.begin ());
4104 typename Teuchos::Array<GO>::iterator gblCols_newEnd;
4105 typename Teuchos::Array<Scalar>::iterator vals_newEnd;
4106 merge2 (gblCols_newEnd, vals_newEnd,
4107 gblCols.begin (), gblCols.end (),
4108 vals.begin (), vals.end ());
4109 gblCols.erase (gblCols_newEnd, gblCols.end ());
4110 vals.erase (vals_newEnd, vals.end ());
4111 numEntPerNonlocalRow[curPos] = gblCols.size ();
4122 GO myMinNonlocalGblRow = std::numeric_limits<GO>::max ();
4124 auto iter = std::min_element (myNonlocalGblRows.begin (),
4125 myNonlocalGblRows.end ());
4126 if (iter != myNonlocalGblRows.end ()) {
4127 myMinNonlocalGblRow = *iter;
4130 GO gblMinNonlocalGblRow = 0;
4131 reduceAll<int, GO> (*comm, REDUCE_MIN, myMinNonlocalGblRow,
4132 outArg (gblMinNonlocalGblRow));
4133 const GO indexBase = gblMinNonlocalGblRow;
4134 const global_size_t INV = Teuchos::OrdinalTraits<global_size_t>::invalid ();
4135 nonlocalRowMap = rcp (
new map_type (INV, myNonlocalGblRows (), indexBase, comm));
4144 std::ostringstream os;
4145 os << *prefix <<
"Create nonlocal matrix" << endl;
4146 std::cerr << os.str();
4148 RCP<crs_matrix_type> nonlocalMatrix =
4149 rcp (
new crs_matrix_type (nonlocalRowMap, numEntPerNonlocalRow ()));
4151 size_type curPos = 0;
4152 for (
auto mapIter = nonlocals_.begin (); mapIter != nonlocals_.end ();
4153 ++mapIter, ++curPos) {
4154 const GO gblRow = mapIter->first;
4156 Teuchos::Array<GO>& gblCols = (mapIter->second).first;
4157 Teuchos::Array<Scalar>& vals = (mapIter->second).second;
4159 nonlocalMatrix->insertGlobalValues (gblRow, gblCols (), vals ());
4171 auto origRowMap = this->getRowMap ();
4172 const bool origRowMapIsOneToOne = origRowMap->isOneToOne ();
4174 int isLocallyComplete = 1;
4176 if (origRowMapIsOneToOne) {
4178 std::ostringstream os;
4179 os << *prefix <<
"Original row Map is 1-to-1" << endl;
4180 std::cerr << os.str();
4182 export_type exportToOrig (nonlocalRowMap, origRowMap);
4184 isLocallyComplete = 0;
4187 std::ostringstream os;
4188 os << *prefix <<
"doExport from nonlocalMatrix" << endl;
4189 std::cerr << os.str();
4191 this->doExport (*nonlocalMatrix, exportToOrig,
Tpetra::ADD);
4196 std::ostringstream os;
4197 os << *prefix <<
"Original row Map is NOT 1-to-1" << endl;
4198 std::cerr << os.str();
4205 export_type exportToOneToOne (nonlocalRowMap, oneToOneRowMap);
4207 isLocallyComplete = 0;
4215 std::ostringstream os;
4216 os << *prefix <<
"Create & doExport into 1-to-1 matrix"
4218 std::cerr << os.str();
4220 crs_matrix_type oneToOneMatrix (oneToOneRowMap, 0);
4222 oneToOneMatrix.doExport(*nonlocalMatrix, exportToOneToOne,
4228 std::ostringstream os;
4229 os << *prefix <<
"Free nonlocalMatrix" << endl;
4230 std::cerr << os.str();
4232 nonlocalMatrix = Teuchos::null;
4236 std::ostringstream os;
4237 os << *prefix <<
"doImport from 1-to-1 matrix" << endl;
4238 std::cerr << os.str();
4240 import_type importToOrig (oneToOneRowMap, origRowMap);
4241 this->doImport (oneToOneMatrix, importToOrig,
Tpetra::ADD);
4249 std::ostringstream os;
4250 os << *prefix <<
"Free nonlocals_ (std::map)" << endl;
4251 std::cerr << os.str();
4253 decltype (nonlocals_) newNonlocals;
4254 std::swap (nonlocals_, newNonlocals);
4263 int isGloballyComplete = 0;
4264 reduceAll<int, int> (*comm, REDUCE_MIN, isLocallyComplete,
4265 outArg (isGloballyComplete));
4266 TEUCHOS_TEST_FOR_EXCEPTION
4267 (isGloballyComplete != 1, std::runtime_error,
"On at least one process, "
4268 "you called insertGlobalValues with a global row index which is not in "
4269 "the matrix's row Map on any process in its communicator.");
4272 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4277 if (! isStaticGraph ()) {
4278 myGraph_->resumeFill (params);
4280 #if KOKKOSKERNELS_VERSION >= 40299
4282 applyHelper.reset();
4284 fillComplete_ =
false;
4287 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4291 return getCrsGraphRef ().haveGlobalConstants ();
4294 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4299 const char tfecfFuncName[] =
"fillComplete(params): ";
4301 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4302 (this->getCrsGraph ().is_null (), std::logic_error,
4303 "getCrsGraph() returns null. This should not happen at this point. "
4304 "Please report this bug to the Tpetra developers.");
4314 Teuchos::RCP<const map_type> rangeMap = graph.
getRowMap ();
4315 Teuchos::RCP<const map_type> domainMap = rangeMap;
4316 this->fillComplete (domainMap, rangeMap, params);
4320 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4324 const Teuchos::RCP<const map_type>& rangeMap,
4325 const Teuchos::RCP<Teuchos::ParameterList>& params)
4329 using Teuchos::ArrayRCP;
4333 const char tfecfFuncName[] =
"fillComplete: ";
4334 ProfilingRegion regionFillComplete
4335 (
"Tpetra::CrsMatrix::fillComplete");
4336 const bool verbose = Behavior::verbose(
"CrsMatrix");
4337 std::unique_ptr<std::string> prefix;
4339 prefix = this->createPrefix(
"CrsMatrix",
"fillComplete(dom,ran,p)");
4340 std::ostringstream os;
4341 os << *prefix << endl;
4342 std::cerr << os.str ();
4345 "Tpetra::CrsMatrix::fillCompete",
4348 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4349 (! this->isFillActive () || this->isFillComplete (), std::runtime_error,
4350 "Matrix fill state must be active (isFillActive() "
4351 "must be true) before you may call fillComplete().");
4352 const int numProcs = this->getComm ()->getSize ();
4362 bool assertNoNonlocalInserts =
false;
4365 bool sortGhosts =
true;
4367 if (! params.is_null ()) {
4368 assertNoNonlocalInserts = params->get (
"No Nonlocal Changes",
4369 assertNoNonlocalInserts);
4370 if (params->isParameter (
"sort column map ghost gids")) {
4371 sortGhosts = params->get (
"sort column map ghost gids", sortGhosts);
4373 else if (params->isParameter (
"Sort column Map ghost GIDs")) {
4374 sortGhosts = params->get (
"Sort column Map ghost GIDs", sortGhosts);
4379 const bool needGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
4381 if (! this->myGraph_.is_null ()) {
4382 this->myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
4385 if (! this->getCrsGraphRef ().indicesAreAllocated ()) {
4386 if (this->hasColMap ()) {
4387 allocateValues(LocalIndices, GraphNotYetAllocated, verbose);
4390 allocateValues(GlobalIndices, GraphNotYetAllocated, verbose);
4395 if (needGlobalAssemble) {
4396 this->globalAssemble ();
4399 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4400 (numProcs == 1 && nonlocals_.size() > 0,
4401 std::runtime_error,
"Cannot have nonlocal entries on a serial run. "
4402 "An invalid entry (i.e., with row index not in the row Map) must have "
4403 "been submitted to the CrsMatrix.");
4406 if (this->isStaticGraph ()) {
4414 #ifdef HAVE_TPETRA_DEBUG
4432 const bool domainMapsMatch =
4433 this->staticGraph_->getDomainMap ()->isSameAs (*domainMap);
4434 const bool rangeMapsMatch =
4435 this->staticGraph_->getRangeMap ()->isSameAs (*rangeMap);
4437 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4438 (! domainMapsMatch, std::runtime_error,
4439 "The CrsMatrix's domain Map does not match the graph's domain Map. "
4440 "The graph cannot be changed because it was given to the CrsMatrix "
4441 "constructor as const. You can fix this by passing in the graph's "
4442 "domain Map and range Map to the matrix's fillComplete call.");
4444 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4445 (! rangeMapsMatch, std::runtime_error,
4446 "The CrsMatrix's range Map does not match the graph's range Map. "
4447 "The graph cannot be changed because it was given to the CrsMatrix "
4448 "constructor as const. You can fix this by passing in the graph's "
4449 "domain Map and range Map to the matrix's fillComplete call.");
4450 #endif // HAVE_TPETRA_DEBUG
4454 this->fillLocalMatrix (params);
4462 this->myGraph_->setDomainRangeMaps (domainMap, rangeMap);
4465 Teuchos::Array<int> remotePIDs (0);
4466 const bool mustBuildColMap = ! this->hasColMap ();
4467 if (mustBuildColMap) {
4468 this->myGraph_->makeColMap (remotePIDs);
4473 const std::pair<size_t, std::string> makeIndicesLocalResult =
4474 this->myGraph_->makeIndicesLocal(verbose);
4479 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4480 (makeIndicesLocalResult.first != 0, std::runtime_error,
4481 makeIndicesLocalResult.second);
4483 const bool sorted = this->myGraph_->isSorted ();
4484 const bool merged = this->myGraph_->isMerged ();
4485 this->sortAndMergeIndicesAndValues (sorted, merged);
4490 this->myGraph_->makeImportExport (remotePIDs, mustBuildColMap);
4494 this->fillLocalGraphAndMatrix (params);
4496 const bool callGraphComputeGlobalConstants = params.get () ==
nullptr ||
4497 params->get (
"compute global constants",
true);
4498 if (callGraphComputeGlobalConstants) {
4499 this->myGraph_->computeGlobalConstants ();
4502 this->myGraph_->computeLocalConstants ();
4504 this->myGraph_->fillComplete_ =
true;
4505 this->myGraph_->checkInternalState ();
4510 this->fillComplete_ =
true;
4513 "Tpetra::CrsMatrix::fillCompete",
"checkInternalState"
4515 this->checkInternalState ();
4519 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4523 const Teuchos::RCP<const map_type> & rangeMap,
4524 const Teuchos::RCP<const import_type>& importer,
4525 const Teuchos::RCP<const export_type>& exporter,
4526 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
4528 #ifdef HAVE_TPETRA_MMM_TIMINGS
4530 if(!params.is_null())
4531 label = params->get(
"Timer Label",label);
4532 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
4533 using Teuchos::TimeMonitor;
4535 Teuchos::TimeMonitor all(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-all")));
4538 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
4539 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! isFillActive() || isFillComplete(),
4540 std::runtime_error,
"Matrix fill state must be active (isFillActive() "
4541 "must be true) before calling fillComplete().");
4542 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4543 myGraph_.is_null (), std::logic_error,
"myGraph_ is null. This is not allowed.");
4546 #ifdef HAVE_TPETRA_MMM_TIMINGS
4547 Teuchos::TimeMonitor graph(*TimeMonitor::getNewTimer(prefix + std::string(
"eSFC-M-Graph")));
4550 myGraph_->expertStaticFillComplete (domainMap, rangeMap, importer, exporter,params);
4554 #ifdef HAVE_TPETRA_MMM_TIMINGS
4555 TimeMonitor fLGAM(*TimeMonitor::getNewTimer(prefix + std::string(
"eSFC-M-fLGAM")));
4558 fillLocalGraphAndMatrix (params);
4563 fillComplete_ =
true;
4566 #ifdef HAVE_TPETRA_DEBUG
4567 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(isFillActive(), std::logic_error,
4568 ": We're at the end of fillComplete(), but isFillActive() is true. "
4569 "Please report this bug to the Tpetra developers.");
4570 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(! isFillComplete(), std::logic_error,
4571 ": We're at the end of fillComplete(), but isFillActive() is true. "
4572 "Please report this bug to the Tpetra developers.");
4573 #endif // HAVE_TPETRA_DEBUG
4575 #ifdef HAVE_TPETRA_MMM_TIMINGS
4576 Teuchos::TimeMonitor cIS(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cIS")));
4579 checkInternalState();
4583 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4589 LocalOrdinal* beg = cols;
4590 LocalOrdinal* end = cols + rowLen;
4591 LocalOrdinal* newend = beg;
4593 LocalOrdinal* cur = beg + 1;
4597 while (cur != end) {
4598 if (*cur != *newend) {
4615 return newend - beg;
4618 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4623 using ::Tpetra::Details::ProfilingRegion;
4624 typedef LocalOrdinal LO;
4625 typedef typename Kokkos::View<LO*, device_type>::HostMirror::execution_space
4626 host_execution_space;
4627 typedef Kokkos::RangePolicy<host_execution_space, LO> range_type;
4628 const char tfecfFuncName[] =
"sortAndMergeIndicesAndValues: ";
4629 ProfilingRegion regionSAM (
"Tpetra::CrsMatrix::sortAndMergeIndicesAndValues");
4631 if (! sorted || ! merged) {
4632 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4633 (this->isStaticGraph (), std::runtime_error,
"Cannot sort or merge with "
4634 "\"static\" (const) graph, since the matrix does not own the graph.");
4635 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4636 (this->myGraph_.is_null (), std::logic_error,
"myGraph_ is null, but "
4637 "this matrix claims ! isStaticGraph(). "
4638 "Please report this bug to the Tpetra developers.");
4639 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
4640 (this->isStorageOptimized (), std::logic_error,
"It is invalid to call "
4641 "this method if the graph's storage has already been optimized. "
4642 "Please report this bug to the Tpetra developers.");
4645 const LO lclNumRows =
static_cast<LO
> (this->getLocalNumRows ());
4646 size_t totalNumDups = 0;
4651 auto vals_ = this->valuesUnpacked_wdv.getHostView(Access::ReadWrite);
4653 Kokkos::parallel_reduce (
"sortAndMergeIndicesAndValues", range_type (0, lclNumRows),
4654 [=] (
const LO lclRow,
size_t& numDups) {
4655 size_t rowBegin = rowBegins_(lclRow);
4656 size_t rowLen = rowLengths_(lclRow);
4657 LO* cols = cols_.data() + rowBegin;
4660 sort2 (cols, cols + rowLen, vals);
4663 size_t newRowLength = mergeRowIndicesAndValues (rowLen, cols, vals);
4664 rowLengths_(lclRow) = newRowLength;
4665 numDups += rowLen - newRowLength;
4678 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4689 using Teuchos::rcp_const_cast;
4690 using Teuchos::rcpFromRef;
4691 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4692 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
4698 if (alpha == ZERO) {
4701 }
else if (beta != ONE) {
4715 RCP<const import_type> importer = this->getGraph ()->getImporter ();
4716 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
4722 const bool Y_is_overwritten = (beta ==
ZERO);
4725 const bool Y_is_replicated =
4726 (! Y_in.
isDistributed () && this->getComm ()->getSize () != 1);
4734 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
4741 RCP<const MV> X_colMap;
4742 if (importer.is_null ()) {
4750 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in,
true);
4752 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
4757 X_colMap = rcpFromRef (X_in);
4761 ProfilingRegion regionImport (
"Tpetra::CrsMatrix::apply: Import");
4767 RCP<MV> X_colMapNonConst = getColumnMapMultiVector (X_in);
4770 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
4771 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
4778 RCP<MV> Y_rowMap = getRowMapMultiVector (Y_in);
4785 if (! exporter.is_null ()) {
4786 this->localApply (*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, ZERO);
4788 ProfilingRegion regionExport (
"Tpetra::CrsMatrix::apply: Export");
4794 if (Y_is_overwritten) {
4820 Y_rowMap = getRowMapMultiVector (Y_in,
true);
4827 this->localApply (*X_colMap, *Y_rowMap, Teuchos::NO_TRANS, alpha, beta);
4831 this->localApply (*X_colMap, Y_in, Teuchos::NO_TRANS, alpha, beta);
4839 if (Y_is_replicated) {
4840 ProfilingRegion regionReduce (
"Tpetra::CrsMatrix::apply: Reduce Y");
4845 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4850 const Teuchos::ETransp mode,
4855 using Teuchos::null;
4858 using Teuchos::rcp_const_cast;
4859 using Teuchos::rcpFromRef;
4860 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4863 if (alpha == ZERO) {
4876 else if (beta == ZERO) {
4894 RCP<const import_type> importer = this->getGraph ()->getImporter ();
4895 RCP<const export_type> exporter = this->getGraph ()->getExporter ();
4900 const bool Y_is_replicated = (! Y_in.
isDistributed () && this->getComm ()->getSize () != 1);
4901 const bool Y_is_overwritten = (beta ==
ZERO);
4902 if (Y_is_replicated && this->getComm ()->getRank () > 0) {
4908 X = rcp (
new MV (X_in, Teuchos::Copy));
4910 X = rcpFromRef (X_in);
4914 if (importer != Teuchos::null) {
4915 if (importMV_ != Teuchos::null && importMV_->getNumVectors() != numVectors) {
4918 if (importMV_ == null) {
4919 importMV_ = rcp (
new MV (this->getColMap (), numVectors));
4922 if (exporter != Teuchos::null) {
4923 if (exportMV_ != Teuchos::null && exportMV_->getNumVectors() != numVectors) {
4926 if (exportMV_ == null) {
4927 exportMV_ = rcp (
new MV (this->getRowMap (), numVectors));
4933 if (! exporter.is_null ()) {
4934 ProfilingRegion regionImport (
"Tpetra::CrsMatrix::apply (transpose): Import");
4935 exportMV_->doImport (X_in, *exporter,
INSERT);
4942 if (importer != Teuchos::null) {
4943 ProfilingRegion regionExport (
"Tpetra::CrsMatrix::apply (transpose): Export");
4950 importMV_->putScalar (ZERO);
4952 this->localApply (*X, *importMV_, mode, alpha, ZERO);
4954 if (Y_is_overwritten) {
4971 MV Y (Y_in, Teuchos::Copy);
4972 this->localApply (*X, Y, mode, alpha, beta);
4975 this->localApply (*X, Y_in, mode, alpha, beta);
4982 if (Y_is_replicated) {
4983 ProfilingRegion regionReduce (
"Tpetra::CrsMatrix::apply (transpose): Reduce Y");
4988 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
4993 const Teuchos::ETransp mode,
4994 const Scalar& alpha,
4995 const Scalar& beta)
const
4998 using Teuchos::NO_TRANS;
4999 ProfilingRegion regionLocalApply (
"Tpetra::CrsMatrix::localApply");
5006 const char tfecfFuncName[] =
"localApply: ";
5007 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5011 const bool transpose = (mode != Teuchos::NO_TRANS);
5012 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5014 getColMap ()->getLocalNumElements (), std::runtime_error,
5015 "NO_TRANS case: X has the wrong number of local rows. "
5017 "getColMap()->getLocalNumElements() = " <<
5018 getColMap ()->getLocalNumElements () <<
".");
5019 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5021 getRowMap ()->getLocalNumElements (), std::runtime_error,
5022 "NO_TRANS case: Y has the wrong number of local rows. "
5024 "getRowMap()->getLocalNumElements() = " <<
5025 getRowMap ()->getLocalNumElements () <<
".");
5026 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5028 getRowMap ()->getLocalNumElements (), std::runtime_error,
5029 "TRANS or CONJ_TRANS case: X has the wrong number of local "
5031 <<
" != getRowMap()->getLocalNumElements() = "
5032 << getRowMap ()->getLocalNumElements () <<
".");
5033 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5035 getColMap ()->getLocalNumElements (), std::runtime_error,
5036 "TRANS or CONJ_TRANS case: X has the wrong number of local "
5038 <<
" != getColMap()->getLocalNumElements() = "
5039 << getColMap ()->getLocalNumElements () <<
".");
5040 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5041 (! isFillComplete (), std::runtime_error,
"The matrix is not "
5042 "fill complete. You must call fillComplete() (possibly with "
5043 "domain and range Map arguments) without an intervening "
5044 "resumeFill() call before you may call this method.");
5045 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5047 std::runtime_error,
"X and Y must be constant stride.");
5052 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5053 (X_lcl.data () == Y_lcl.data () && X_lcl.data () !=
nullptr
5054 && X_lcl.extent(0) != 0,
5055 std::runtime_error,
"X and Y may not alias one another.");
5058 #if KOKKOSKERNELS_VERSION >= 40299
5059 auto A_lcl = getLocalMatrixDevice();
5061 if(!applyHelper.get()) {
5064 bool useMergePath =
false;
5065 #ifdef KOKKOSKERNELS_ENABLE_TPL_CUSPARSE
5071 if constexpr(std::is_same_v<execution_space, Kokkos::Cuda>) {
5072 LocalOrdinal nrows = getLocalNumRows();
5073 LocalOrdinal maxRowImbalance = 0;
5075 maxRowImbalance = getLocalMaxNumRowEntries() - (getLocalNumEntries() / nrows);
5078 useMergePath =
true;
5081 applyHelper = std::make_shared<ApplyHelper>(A_lcl.nnz(), A_lcl.graph.row_map,
5082 useMergePath ? KokkosSparse::SPMV_MERGE_PATH : KokkosSparse::SPMV_DEFAULT);
5086 const char* modeKK =
nullptr;
5089 case Teuchos::NO_TRANS:
5090 modeKK = KokkosSparse::NoTranspose;
break;
5091 case Teuchos::TRANS:
5092 modeKK = KokkosSparse::Transpose;
break;
5093 case Teuchos::CONJ_TRANS:
5094 modeKK = KokkosSparse::ConjugateTranspose;
break;
5096 throw std::invalid_argument(
"Tpetra::CrsMatrix::localApply: invalid mode");
5099 if(applyHelper->shouldUseIntRowptrs())
5101 auto A_lcl_int_rowptrs = applyHelper->getIntRowptrMatrix(A_lcl);
5103 &applyHelper->handle_int, modeKK,
5109 &applyHelper->handle, modeKK,
5113 LocalOrdinal nrows = getLocalNumRows();
5114 LocalOrdinal maxRowImbalance = 0;
5116 maxRowImbalance = getLocalMaxNumRowEntries() - (getLocalNumEntries() / nrows);
5118 auto matrix_lcl = getLocalMultiplyOperator();
5120 matrix_lcl->applyImbalancedRows (X_lcl, Y_lcl, mode, alpha, beta);
5122 matrix_lcl->apply (X_lcl, Y_lcl, mode, alpha, beta);
5126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5131 Teuchos::ETransp mode,
5136 const char fnName[] =
"Tpetra::CrsMatrix::apply";
5138 TEUCHOS_TEST_FOR_EXCEPTION
5139 (! isFillComplete (), std::runtime_error,
5140 fnName <<
": Cannot call apply() until fillComplete() "
5141 "has been called.");
5143 if (mode == Teuchos::NO_TRANS) {
5144 ProfilingRegion regionNonTranspose (fnName);
5145 this->applyNonTranspose (X, Y, alpha, beta);
5148 ProfilingRegion regionTranspose (
"Tpetra::CrsMatrix::apply (transpose)");
5149 this->applyTranspose (X, Y, mode, alpha, beta);
5154 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5156 Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node> >
5162 const char tfecfFuncName[] =
"convert: ";
5164 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5165 (! this->isFillComplete (), std::runtime_error,
"This matrix (the source "
5166 "of the conversion) is not fill complete. You must first call "
5167 "fillComplete() (possibly with the domain and range Map) without an "
5168 "intervening call to resumeFill(), before you may call this method.");
5170 RCP<output_matrix_type> newMatrix
5171 (
new output_matrix_type (this->getCrsGraph ()));
5175 copyConvert (newMatrix->getLocalMatrixDevice ().values,
5176 this->getLocalMatrixDevice ().values);
5180 newMatrix->fillComplete (this->getDomainMap (), this->getRangeMap ());
5186 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5193 const char tfecfFuncName[] =
"checkInternalState: ";
5194 const char err[] =
"Internal state is not consistent. "
5195 "Please report this bug to the Tpetra developers.";
5199 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5200 (staticGraph_.is_null (), std::logic_error, err);
5204 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5205 (! myGraph_.is_null () && myGraph_ != staticGraph_,
5206 std::logic_error, err);
5208 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5209 (isFillComplete () && ! staticGraph_->isFillComplete (),
5210 std::logic_error, err <<
" Specifically, the matrix is fill complete, "
5211 "but its graph is NOT fill complete.");
5214 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5215 (staticGraph_->indicesAreAllocated () &&
5216 staticGraph_->getLocalAllocationSize() > 0 &&
5217 staticGraph_->getLocalNumRows() > 0 &&
5218 valuesUnpacked_wdv.extent (0) == 0,
5219 std::logic_error, err);
5223 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5228 std::ostringstream os;
5230 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
5231 if (this->getObjectLabel () !=
"") {
5232 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
5234 if (isFillComplete ()) {
5235 os <<
"isFillComplete: true"
5236 <<
", global dimensions: [" << getGlobalNumRows () <<
", "
5237 << getGlobalNumCols () <<
"]"
5238 <<
", global number of entries: " << getGlobalNumEntries ()
5242 os <<
"isFillComplete: false"
5243 <<
", global dimensions: [" << getGlobalNumRows () <<
", "
5244 << getGlobalNumCols () <<
"]}";
5249 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5253 const Teuchos::EVerbosityLevel verbLevel)
const
5257 using Teuchos::ArrayView;
5258 using Teuchos::Comm;
5260 using Teuchos::TypeNameTraits;
5261 using Teuchos::VERB_DEFAULT;
5262 using Teuchos::VERB_NONE;
5263 using Teuchos::VERB_LOW;
5264 using Teuchos::VERB_MEDIUM;
5265 using Teuchos::VERB_HIGH;
5266 using Teuchos::VERB_EXTREME;
5268 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
5270 if (vl == VERB_NONE) {
5275 Teuchos::OSTab tab0 (out);
5277 RCP<const Comm<int> > comm = this->getComm();
5278 const int myRank = comm->getRank();
5279 const int numProcs = comm->getSize();
5281 for (
size_t dec=10; dec<getGlobalNumRows(); dec *= 10) {
5284 width = std::max<size_t> (width,
static_cast<size_t> (11)) + 2;
5294 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
5296 Teuchos::OSTab tab1 (out);
5299 if (this->getObjectLabel () !=
"") {
5300 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
5303 out <<
"Template parameters:" << endl;
5304 Teuchos::OSTab tab2 (out);
5305 out <<
"Scalar: " << TypeNameTraits<Scalar>::name () << endl
5306 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name () << endl
5307 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name () << endl
5308 <<
"Node: " << TypeNameTraits<Node>::name () << endl;
5310 if (isFillComplete()) {
5311 out <<
"isFillComplete: true" << endl
5312 <<
"Global dimensions: [" << getGlobalNumRows () <<
", "
5313 << getGlobalNumCols () <<
"]" << endl
5314 <<
"Global number of entries: " << getGlobalNumEntries () << endl
5315 << endl <<
"Global max number of entries in a row: "
5316 << getGlobalMaxNumRowEntries () << endl;
5319 out <<
"isFillComplete: false" << endl
5320 <<
"Global dimensions: [" << getGlobalNumRows () <<
", "
5321 << getGlobalNumCols () <<
"]" << endl;
5325 if (vl < VERB_MEDIUM) {
5331 out << endl <<
"Row Map:" << endl;
5333 if (getRowMap ().is_null ()) {
5335 out <<
"null" << endl;
5342 getRowMap ()->describe (out, vl);
5347 out <<
"Column Map: ";
5349 if (getColMap ().is_null ()) {
5351 out <<
"null" << endl;
5353 }
else if (getColMap () == getRowMap ()) {
5355 out <<
"same as row Map" << endl;
5361 getColMap ()->describe (out, vl);
5366 out <<
"Domain Map: ";
5368 if (getDomainMap ().is_null ()) {
5370 out <<
"null" << endl;
5372 }
else if (getDomainMap () == getRowMap ()) {
5374 out <<
"same as row Map" << endl;
5376 }
else if (getDomainMap () == getColMap ()) {
5378 out <<
"same as column Map" << endl;
5384 getDomainMap ()->describe (out, vl);
5389 out <<
"Range Map: ";
5391 if (getRangeMap ().is_null ()) {
5393 out <<
"null" << endl;
5395 }
else if (getRangeMap () == getDomainMap ()) {
5397 out <<
"same as domain Map" << endl;
5399 }
else if (getRangeMap () == getRowMap ()) {
5401 out <<
"same as row Map" << endl;
5407 getRangeMap ()->describe (out, vl);
5411 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5412 if (myRank == curRank) {
5413 out <<
"Process rank: " << curRank << endl;
5414 Teuchos::OSTab tab2 (out);
5415 if (! staticGraph_->indicesAreAllocated ()) {
5416 out <<
"Graph indices not allocated" << endl;
5419 out <<
"Number of allocated entries: "
5420 << staticGraph_->getLocalAllocationSize () << endl;
5422 out <<
"Number of entries: " << getLocalNumEntries () << endl
5423 <<
"Max number of entries per row: " << getLocalMaxNumRowEntries ()
5432 if (vl < VERB_HIGH) {
5437 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5438 if (myRank == curRank) {
5439 out << std::setw(width) <<
"Proc Rank"
5440 << std::setw(width) <<
"Global Row"
5441 << std::setw(width) <<
"Num Entries";
5442 if (vl == VERB_EXTREME) {
5443 out << std::setw(width) <<
"(Index,Value)";
5446 for (
size_t r = 0; r < getLocalNumRows (); ++r) {
5447 const size_t nE = getNumEntriesInLocalRow(r);
5448 GlobalOrdinal gid = getRowMap()->getGlobalElement(r);
5449 out << std::setw(width) << myRank
5450 << std::setw(width) << gid
5451 << std::setw(width) << nE;
5452 if (vl == VERB_EXTREME) {
5453 if (isGloballyIndexed()) {
5454 global_inds_host_view_type rowinds;
5455 values_host_view_type rowvals;
5456 getGlobalRowView (gid, rowinds, rowvals);
5457 for (
size_t j = 0; j < nE; ++j) {
5458 out <<
" (" << rowinds[j]
5459 <<
", " << rowvals[j]
5463 else if (isLocallyIndexed()) {
5464 local_inds_host_view_type rowinds;
5465 values_host_view_type rowvals;
5466 getLocalRowView (r, rowinds, rowvals);
5467 for (
size_t j=0; j < nE; ++j) {
5468 out <<
" (" << getColMap()->getGlobalElement(rowinds[j])
5469 <<
", " << rowvals[j]
5485 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5499 return (srcRowMat !=
nullptr);
5502 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5506 const typename crs_graph_type::padding_type& padding,
5512 using LO = local_ordinal_type;
5513 using row_ptrs_type =
5514 typename local_graph_device_type::row_map_type::non_const_type;
5515 using range_policy =
5516 Kokkos::RangePolicy<execution_space, Kokkos::IndexType<LO>>;
5517 const char tfecfFuncName[] =
"applyCrsPadding";
5518 const char suffix[] =
5519 ". Please report this bug to the Tpetra developers.";
5520 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::applyCrsPadding");
5522 std::unique_ptr<std::string> prefix;
5524 prefix = this->createPrefix(
"CrsMatrix", tfecfFuncName);
5525 std::ostringstream os;
5526 os << *prefix <<
"padding: ";
5529 std::cerr << os.str();
5531 const int myRank = ! verbose ? -1 : [&] () {
5532 auto map = this->getMap();
5533 if (map.is_null()) {
5536 auto comm = map->getComm();
5537 if (comm.is_null()) {
5540 return comm->getRank();
5544 if (! myGraph_->indicesAreAllocated()) {
5546 std::ostringstream os;
5547 os << *prefix <<
"Call allocateIndices" << endl;
5548 std::cerr << os.str();
5550 allocateValues(GlobalIndices, GraphNotYetAllocated, verbose);
5562 std::ostringstream os;
5563 os << *prefix <<
"Allocate row_ptrs_beg: "
5564 << myGraph_->getRowPtrsUnpackedHost().extent(0) << endl;
5565 std::cerr << os.str();
5567 using Kokkos::view_alloc;
5568 using Kokkos::WithoutInitializing;
5569 row_ptrs_type row_ptr_beg(view_alloc(
"row_ptr_beg", WithoutInitializing),
5570 myGraph_->rowPtrsUnpacked_dev_.extent(0));
5572 Kokkos::deep_copy(execution_space(),row_ptr_beg, myGraph_->rowPtrsUnpacked_dev_);
5574 const size_t N = row_ptr_beg.extent(0) == 0 ? size_t(0) :
5575 size_t(row_ptr_beg.extent(0) - 1);
5577 std::ostringstream os;
5578 os << *prefix <<
"Allocate row_ptrs_end: " << N << endl;
5579 std::cerr << os.str();
5581 row_ptrs_type row_ptr_end(
5582 view_alloc(
"row_ptr_end", WithoutInitializing), N);
5584 row_ptrs_type num_row_entries_d;
5586 const bool refill_num_row_entries =
5587 myGraph_->k_numRowEntries_.extent(0) != 0;
5589 if (refill_num_row_entries) {
5592 num_row_entries_d = create_mirror_view_and_copy(memory_space(),
5593 myGraph_->k_numRowEntries_);
5594 Kokkos::parallel_for
5595 (
"Fill end row pointers", range_policy(0, N),
5596 KOKKOS_LAMBDA (
const size_t i) {
5597 row_ptr_end(i) = row_ptr_beg(i) + num_row_entries_d(i);
5604 Kokkos::parallel_for
5605 (
"Fill end row pointers", range_policy(0, N),
5606 KOKKOS_LAMBDA (
const size_t i) {
5607 row_ptr_end(i) = row_ptr_beg(i+1);
5611 if (myGraph_->isGloballyIndexed()) {
5613 myGraph_->gblInds_wdv,
5614 valuesUnpacked_wdv, padding, myRank, verbose);
5615 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5616 const auto newColIndsLen = myGraph_->gblInds_wdv.extent(0);
5617 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5618 (newValuesLen != newColIndsLen, std::logic_error,
5619 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5620 <<
" != myGraph_->gblInds_wdv.extent(0)=" << newColIndsLen
5625 myGraph_->lclIndsUnpacked_wdv,
5626 valuesUnpacked_wdv, padding, myRank, verbose);
5627 const auto newValuesLen = valuesUnpacked_wdv.extent(0);
5628 const auto newColIndsLen = myGraph_->lclIndsUnpacked_wdv.extent(0);
5629 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5630 (newValuesLen != newColIndsLen, std::logic_error,
5631 ": After padding, valuesUnpacked_wdv.extent(0)=" << newValuesLen
5632 <<
" != myGraph_->lclIndsUnpacked_wdv.extent(0)=" << newColIndsLen
5636 if (refill_num_row_entries) {
5637 Kokkos::parallel_for
5638 (
"Fill num entries", range_policy(0, N),
5639 KOKKOS_LAMBDA (
const size_t i) {
5640 num_row_entries_d(i) = row_ptr_end(i) - row_ptr_beg(i);
5646 std::ostringstream os;
5647 os << *prefix <<
"Assign myGraph_->rowPtrsUnpacked_; "
5648 <<
"old size: " << myGraph_->rowPtrsUnpacked_host_.extent(0)
5649 <<
", new size: " << row_ptr_beg.extent(0) << endl;
5650 std::cerr << os.str();
5651 TEUCHOS_ASSERT( myGraph_->getRowPtrsUnpackedHost().extent(0) ==
5652 row_ptr_beg.extent(0) );
5654 myGraph_->setRowPtrsUnpacked(row_ptr_beg);
5657 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5659 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
5660 copyAndPermuteStaticGraph(
5661 const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& srcMat,
5662 const size_t numSameIDs,
5663 const LocalOrdinal permuteToLIDs[],
5664 const LocalOrdinal permuteFromLIDs[],
5665 const size_t numPermutes)
5667 using Details::ProfilingRegion;
5668 using Teuchos::Array;
5669 using Teuchos::ArrayView;
5671 using LO = LocalOrdinal;
5672 using GO = GlobalOrdinal;
5673 const char tfecfFuncName[] =
"copyAndPermuteStaticGraph";
5674 const char suffix[] =
5675 " Please report this bug to the Tpetra developers.";
5676 ProfilingRegion regionCAP
5677 (
"Tpetra::CrsMatrix::copyAndPermuteStaticGraph");
5681 std::unique_ptr<std::string> prefix;
5683 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5684 std::ostringstream os;
5685 os << *prefix <<
"Start" << endl;
5687 const char*
const prefix_raw =
5688 verbose ? prefix.get()->c_str() :
nullptr;
5690 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed ();
5695 const map_type& srcRowMap = * (srcMat.getRowMap ());
5696 nonconst_global_inds_host_view_type rowInds;
5697 nonconst_values_host_view_type rowVals;
5698 const LO numSameIDs_as_LID =
static_cast<LO
> (numSameIDs);
5699 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5703 const GO sourceGID = srcRowMap.getGlobalElement (sourceLID);
5704 const GO targetGID = sourceGID;
5706 ArrayView<const GO>rowIndsConstView;
5707 ArrayView<const Scalar> rowValsConstView;
5709 if (sourceIsLocallyIndexed) {
5710 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5711 if (rowLength > static_cast<size_t> (rowInds.size())) {
5712 Kokkos::resize(rowInds,rowLength);
5713 Kokkos::resize(rowVals,rowLength);
5717 nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(rowInds,std::make_pair((
size_t)0, rowLength));
5718 nonconst_values_host_view_type rowValsView = Kokkos::subview(rowVals,std::make_pair((
size_t)0, rowLength));
5723 size_t checkRowLength = 0;
5724 srcMat.getGlobalRowCopy (sourceGID, rowIndsView,
5725 rowValsView, checkRowLength);
5727 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5728 (rowLength != checkRowLength, std::logic_error,
"For "
5729 "global row index " << sourceGID <<
", the source "
5730 "matrix's getNumEntriesInGlobalRow returns a row length "
5731 "of " << rowLength <<
", but getGlobalRowCopy reports "
5732 "a row length of " << checkRowLength <<
"." << suffix);
5739 rowIndsConstView = Teuchos::ArrayView<const GO> (
5740 rowIndsView.data(), rowIndsView.extent(0),
5741 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5742 rowValsConstView = Teuchos::ArrayView<const Scalar> (
5743 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5744 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5749 global_inds_host_view_type rowIndsView;
5750 values_host_view_type rowValsView;
5751 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5756 rowIndsConstView = Teuchos::ArrayView<const GO> (
5757 rowIndsView.data(), rowIndsView.extent(0),
5758 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5759 rowValsConstView = Teuchos::ArrayView<const Scalar> (
5760 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5761 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5769 combineGlobalValues(targetGID, rowIndsConstView,
5771 prefix_raw, debug, verbose);
5775 std::ostringstream os;
5776 os << *prefix <<
"Do permutes" << endl;
5779 const map_type& tgtRowMap = * (this->getRowMap ());
5780 for (
size_t p = 0; p < numPermutes; ++p) {
5781 const GO sourceGID = srcRowMap.getGlobalElement (permuteFromLIDs[p]);
5782 const GO targetGID = tgtRowMap.getGlobalElement (permuteToLIDs[p]);
5784 ArrayView<const GO> rowIndsConstView;
5785 ArrayView<const Scalar> rowValsConstView;
5787 if (sourceIsLocallyIndexed) {
5788 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5789 if (rowLength > static_cast<size_t> (rowInds.size ())) {
5790 Kokkos::resize(rowInds,rowLength);
5791 Kokkos::resize(rowVals,rowLength);
5795 nonconst_global_inds_host_view_type rowIndsView = Kokkos::subview(rowInds,std::make_pair((
size_t)0, rowLength));
5796 nonconst_values_host_view_type rowValsView = Kokkos::subview(rowVals,std::make_pair((
size_t)0, rowLength));
5801 size_t checkRowLength = 0;
5802 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
5803 rowValsView, checkRowLength);
5805 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5806 (rowLength != checkRowLength, std::logic_error,
"For "
5807 "source matrix global row index " << sourceGID <<
", "
5808 "getNumEntriesInGlobalRow returns a row length of " <<
5809 rowLength <<
", but getGlobalRowCopy a row length of "
5810 << checkRowLength <<
"." << suffix);
5817 rowIndsConstView = Teuchos::ArrayView<const GO> (
5818 rowIndsView.data(), rowIndsView.extent(0),
5819 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5820 rowValsConstView = Teuchos::ArrayView<const Scalar> (
5821 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5822 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5827 global_inds_host_view_type rowIndsView;
5828 values_host_view_type rowValsView;
5829 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5834 rowIndsConstView = Teuchos::ArrayView<const GO> (
5835 rowIndsView.data(), rowIndsView.extent(0),
5836 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5837 rowValsConstView = Teuchos::ArrayView<const Scalar> (
5838 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5839 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5844 combineGlobalValues(targetGID, rowIndsConstView,
5846 prefix_raw, debug, verbose);
5850 std::ostringstream os;
5851 os << *prefix <<
"Done" << endl;
5855 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
5857 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
5858 copyAndPermuteNonStaticGraph(
5859 const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& srcMat,
5860 const size_t numSameIDs,
5861 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs_dv,
5862 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs_dv,
5863 const size_t numPermutes)
5865 using Details::ProfilingRegion;
5866 using Teuchos::Array;
5867 using Teuchos::ArrayView;
5869 using LO = LocalOrdinal;
5870 using GO = GlobalOrdinal;
5871 const char tfecfFuncName[] =
"copyAndPermuteNonStaticGraph";
5872 const char suffix[] =
5873 " Please report this bug to the Tpetra developers.";
5874 ProfilingRegion regionCAP
5875 (
"Tpetra::CrsMatrix::copyAndPermuteNonStaticGraph");
5879 std::unique_ptr<std::string> prefix;
5881 prefix = this->
createPrefix(
"CrsGraph", tfecfFuncName);
5882 std::ostringstream os;
5883 os << *prefix <<
"Start" << endl;
5885 const char*
const prefix_raw =
5886 verbose ? prefix.get()->c_str() :
nullptr;
5889 using row_graph_type = RowGraph<LO, GO, Node>;
5890 const row_graph_type& srcGraph = *(srcMat.getGraph());
5892 myGraph_->computeCrsPadding(srcGraph, numSameIDs,
5893 permuteToLIDs_dv, permuteFromLIDs_dv, verbose);
5894 applyCrsPadding(*padding, verbose);
5896 const bool sourceIsLocallyIndexed = srcMat.isLocallyIndexed ();
5901 const map_type& srcRowMap = * (srcMat.getRowMap ());
5902 const LO numSameIDs_as_LID =
static_cast<LO
> (numSameIDs);
5903 using gids_type = nonconst_global_inds_host_view_type;
5904 using vals_type = nonconst_values_host_view_type;
5907 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5911 const GO sourceGID = srcRowMap.getGlobalElement (sourceLID);
5912 const GO targetGID = sourceGID;
5914 ArrayView<const GO> rowIndsConstView;
5915 ArrayView<const Scalar> rowValsConstView;
5917 if (sourceIsLocallyIndexed) {
5919 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5920 if (rowLength > static_cast<size_t> (rowInds.extent(0))) {
5921 Kokkos::resize(rowInds,rowLength);
5922 Kokkos::resize(rowVals,rowLength);
5926 gids_type rowIndsView = Kokkos::subview(rowInds,std::make_pair((
size_t)0, rowLength));
5927 vals_type rowValsView = Kokkos::subview(rowVals,std::make_pair((
size_t)0, rowLength));
5932 size_t checkRowLength = 0;
5933 srcMat.getGlobalRowCopy (sourceGID, rowIndsView, rowValsView,
5936 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
5937 (rowLength != checkRowLength, std::logic_error,
": For "
5938 "global row index " << sourceGID <<
", the source "
5939 "matrix's getNumEntriesInGlobalRow returns a row length "
5940 "of " << rowLength <<
", but getGlobalRowCopy reports "
5941 "a row length of " << checkRowLength <<
"." << suffix);
5943 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
5944 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar *
>(rowValsView.data()), rowLength);
5947 global_inds_host_view_type rowIndsView;
5948 values_host_view_type rowValsView;
5949 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
5955 rowIndsConstView = Teuchos::ArrayView<const GO> (
5956 rowIndsView.data(), rowIndsView.extent(0),
5957 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5958 rowValsConstView = Teuchos::ArrayView<const Scalar> (
5959 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
5960 Teuchos::RCP_DISABLE_NODE_LOOKUP);
5966 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
5967 rowValsConstView, prefix_raw, debug, verbose);
5971 std::ostringstream os;
5972 os << *prefix <<
"Do permutes" << endl;
5974 const LO*
const permuteFromLIDs = permuteFromLIDs_dv.view_host().data();
5975 const LO*
const permuteToLIDs = permuteToLIDs_dv.view_host().data();
5977 const map_type& tgtRowMap = * (this->getRowMap ());
5978 for (
size_t p = 0; p < numPermutes; ++p) {
5979 const GO sourceGID = srcRowMap.getGlobalElement (permuteFromLIDs[p]);
5980 const GO targetGID = tgtRowMap.getGlobalElement (permuteToLIDs[p]);
5982 ArrayView<const GO> rowIndsConstView;
5983 ArrayView<const Scalar> rowValsConstView;
5985 if (sourceIsLocallyIndexed) {
5986 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5987 if (rowLength > static_cast<size_t> (rowInds.extent(0))) {
5988 Kokkos::resize(rowInds,rowLength);
5989 Kokkos::resize(rowVals,rowLength);
5993 gids_type rowIndsView = Kokkos::subview(rowInds,std::make_pair((
size_t)0, rowLength));
5994 vals_type rowValsView = Kokkos::subview(rowVals,std::make_pair((
size_t)0, rowLength));
5999 size_t checkRowLength = 0;
6000 srcMat.getGlobalRowCopy(sourceGID, rowIndsView,
6001 rowValsView, checkRowLength);
6003 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6004 (rowLength != checkRowLength, std::logic_error,
"For "
6005 "source matrix global row index " << sourceGID <<
", "
6006 "getNumEntriesInGlobalRow returns a row length of " <<
6007 rowLength <<
", but getGlobalRowCopy a row length of "
6008 << checkRowLength <<
"." << suffix);
6010 rowIndsConstView = Teuchos::ArrayView<const GO>(rowIndsView.data(), rowLength);
6011 rowValsConstView = Teuchos::ArrayView<const Scalar>(
reinterpret_cast<Scalar *
>(rowValsView.data()), rowLength);
6014 global_inds_host_view_type rowIndsView;
6015 values_host_view_type rowValsView;
6016 srcMat.getGlobalRowView(sourceGID, rowIndsView, rowValsView);
6022 rowIndsConstView = Teuchos::ArrayView<const GO> (
6023 rowIndsView.data(), rowIndsView.extent(0),
6024 Teuchos::RCP_DISABLE_NODE_LOOKUP);
6025 rowValsConstView = Teuchos::ArrayView<const Scalar> (
6026 reinterpret_cast<const Scalar*
>(rowValsView.data()), rowValsView.extent(0),
6027 Teuchos::RCP_DISABLE_NODE_LOOKUP);
6033 insertGlobalValuesFilteredChecked(targetGID, rowIndsConstView,
6034 rowValsConstView, prefix_raw, debug, verbose);
6038 std::ostringstream os;
6039 os << *prefix <<
"Done" << endl;
6043 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6048 const size_t numSameIDs,
6049 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteToLIDs,
6050 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& permuteFromLIDs,
6059 const char tfecfFuncName[] =
"copyAndPermute: ";
6060 ProfilingRegion regionCAP(
"Tpetra::CrsMatrix::copyAndPermute");
6062 const bool verbose = Behavior::verbose(
"CrsMatrix");
6063 std::unique_ptr<std::string> prefix;
6065 prefix = this->createPrefix(
"CrsMatrix",
"copyAndPermute");
6066 std::ostringstream os;
6067 os << *prefix << endl
6068 << *prefix <<
" numSameIDs: " << numSameIDs << endl
6069 << *prefix <<
" numPermute: " << permuteToLIDs.extent(0)
6078 <<
"isStaticGraph: " << (isStaticGraph() ?
"true" :
"false")
6080 std::cerr << os.str ();
6083 const auto numPermute = permuteToLIDs.extent (0);
6084 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6085 (numPermute != permuteFromLIDs.extent (0),
6086 std::invalid_argument,
"permuteToLIDs.extent(0) = "
6087 << numPermute <<
"!= permuteFromLIDs.extent(0) = "
6088 << permuteFromLIDs.extent (0) <<
".");
6093 const RMT& srcMat =
dynamic_cast<const RMT&
> (srcObj);
6094 if (isStaticGraph ()) {
6095 TEUCHOS_ASSERT( ! permuteToLIDs.need_sync_host () );
6096 auto permuteToLIDs_h = permuteToLIDs.view_host ();
6097 TEUCHOS_ASSERT( ! permuteFromLIDs.need_sync_host () );
6098 auto permuteFromLIDs_h = permuteFromLIDs.view_host ();
6100 copyAndPermuteStaticGraph(srcMat, numSameIDs,
6101 permuteToLIDs_h.data(),
6102 permuteFromLIDs_h.data(),
6106 copyAndPermuteNonStaticGraph(srcMat, numSameIDs, permuteToLIDs,
6107 permuteFromLIDs, numPermute);
6111 std::ostringstream os;
6112 os << *prefix <<
"Done" << endl;
6113 std::cerr << os.str();
6117 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6122 const Kokkos::DualView<const local_ordinal_type*, buffer_device_type>& exportLIDs,
6123 Kokkos::DualView<char*, buffer_device_type>& exports,
6124 Kokkos::DualView<size_t*, buffer_device_type> numPacketsPerLID,
6125 size_t& constantNumPackets)
6130 using Teuchos::outArg;
6131 using Teuchos::REDUCE_MAX;
6132 using Teuchos::reduceAll;
6134 typedef LocalOrdinal LO;
6135 typedef GlobalOrdinal GO;
6136 const char tfecfFuncName[] =
"packAndPrepare: ";
6137 ProfilingRegion regionPAP (
"Tpetra::CrsMatrix::packAndPrepare");
6139 const bool debug = Behavior::debug(
"CrsMatrix");
6140 const bool verbose = Behavior::verbose(
"CrsMatrix");
6143 Teuchos::RCP<const Teuchos::Comm<int> > pComm = this->getComm ();
6144 if (pComm.is_null ()) {
6147 const Teuchos::Comm<int>& comm = *pComm;
6148 const int myRank = comm.getSize ();
6150 std::unique_ptr<std::string> prefix;
6152 prefix = this->createPrefix(
"CrsMatrix",
"packAndPrepare");
6153 std::ostringstream os;
6154 os << *prefix <<
"Start" << endl
6164 std::cerr << os.str ();
6187 std::ostringstream msg;
6190 using crs_matrix_type = CrsMatrix<Scalar, LO, GO, Node>;
6191 const crs_matrix_type* srcCrsMat =
6192 dynamic_cast<const crs_matrix_type*
> (&source);
6193 if (srcCrsMat !=
nullptr) {
6195 std::ostringstream os;
6196 os << *prefix <<
"Source matrix same (CrsMatrix) type as target; "
6197 "calling packNew" << endl;
6198 std::cerr << os.str ();
6201 srcCrsMat->packNew (exportLIDs, exports, numPacketsPerLID,
6202 constantNumPackets);
6204 catch (std::exception& e) {
6206 msg <<
"Proc " << myRank <<
": " << e.what () << std::endl;
6210 using Kokkos::HostSpace;
6211 using Kokkos::subview;
6212 using exports_type = Kokkos::DualView<char*, buffer_device_type>;
6213 using range_type = Kokkos::pair<size_t, size_t>;
6216 std::ostringstream os;
6217 os << *prefix <<
"Source matrix NOT same (CrsMatrix) type as target"
6219 std::cerr << os.str ();
6222 const row_matrix_type* srcRowMat =
6223 dynamic_cast<const row_matrix_type*
> (&source);
6224 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6225 (srcRowMat ==
nullptr, std::invalid_argument,
6226 "The source object of the Import or Export operation is neither a "
6227 "CrsMatrix (with the same template parameters as the target object), "
6228 "nor a RowMatrix (with the same first four template parameters as the "
6239 TEUCHOS_ASSERT( ! exportLIDs.need_sync_host () );
6240 auto exportLIDs_h = exportLIDs.view_host ();
6241 Teuchos::ArrayView<const LO> exportLIDs_av (exportLIDs_h.data (),
6242 exportLIDs_h.size ());
6246 Teuchos::Array<char> exports_a;
6252 numPacketsPerLID.clear_sync_state ();
6253 numPacketsPerLID.modify_host ();
6254 auto numPacketsPerLID_h = numPacketsPerLID.view_host ();
6255 Teuchos::ArrayView<size_t> numPacketsPerLID_av (numPacketsPerLID_h.data (),
6256 numPacketsPerLID_h.size ());
6261 srcRowMat->pack (exportLIDs_av, exports_a, numPacketsPerLID_av,
6262 constantNumPackets);
6264 catch (std::exception& e) {
6266 msg <<
"Proc " << myRank <<
": " << e.what () << std::endl;
6270 const size_t newAllocSize =
static_cast<size_t> (exports_a.size ());
6271 if (static_cast<size_t> (exports.extent (0)) < newAllocSize) {
6272 const std::string oldLabel = exports.d_view.label ();
6273 const std::string newLabel = (oldLabel ==
"") ?
"exports" : oldLabel;
6274 exports = exports_type (newLabel, newAllocSize);
6279 exports.modify_host();
6281 auto exports_h = exports.view_host ();
6282 auto exports_h_sub = subview (exports_h, range_type (0, newAllocSize));
6286 typedef typename exports_type::t_host::execution_space HES;
6287 typedef Kokkos::Device<HES, HostSpace> host_device_type;
6288 Kokkos::View<const char*, host_device_type>
6289 exports_a_kv (exports_a.getRawPtr (), newAllocSize);
6296 reduceAll<int, int> (comm, REDUCE_MAX, lclBad, outArg (gblBad));
6299 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6300 (
true, std::logic_error,
"packNew() or pack() threw an exception on "
6301 "one or more participating processes.");
6305 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
6306 (lclBad != 0, std::logic_error,
"packNew threw an exception on one "
6307 "or more participating processes. Here is this process' error "
6308 "message: " << msg.str ());
6312 std::ostringstream os;
6313 os << *prefix <<
"packAndPrepare: Done!" << endl
6323 std::cerr << os.str ();
6327 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
6329 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
6330 packRow (
char exports[],
6331 const size_t offset,
6332 const size_t numEnt,