42 #ifndef TPETRA_CRSMATRIX_DEF_HPP 43 #define TPETRA_CRSMATRIX_DEF_HPP 53 #include "Tpetra_RowMatrix.hpp" 57 #include "Teuchos_SerialDenseMatrix.hpp" 58 #include "Teuchos_as.hpp" 59 #include "Teuchos_ArrayRCP.hpp" 83 template<
class Scalar>
87 typedef Teuchos::ScalarTraits<Scalar> STS;
88 return std::max (STS::magnitude (x), STS::magnitude (y));
103 template <
class Ordinal,
class Scalar>
113 v (
Teuchos::ScalarTraits<Scalar>::zero ())
121 CrsIJV (Ordinal row, Ordinal col,
const Scalar &val) :
122 i (row), j (col), v (val)
130 bool operator< (const CrsIJV<Ordinal, Scalar>& rhs)
const {
138 return this->i < rhs.i;
161 template <
typename Ordinal,
typename Scalar>
162 class SerializationTraits<int, Tpetra::
Details::CrsIJV<Ordinal, Scalar> >
163 :
public DirectSerializationTraits<int, Tpetra::Details::CrsIJV<Ordinal, Scalar> >
169 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
171 CrsMatrix (
const Teuchos::RCP<const map_type>& rowMap,
172 size_t maxNumEntriesPerRow,
174 const RCP<Teuchos::ParameterList>& params) :
175 dist_object_type (rowMap),
177 Details::STORAGE_1D_UNPACKED :
178 Details::STORAGE_2D),
179 fillComplete_ (
false),
180 frobNorm_ (-STM::one ())
184 myGraph_ = rcp (
new crs_graph_type (rowMap, maxNumEntriesPerRow,
187 catch (std::exception& e) {
188 TEUCHOS_TEST_FOR_EXCEPTION(
189 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 190 "exception while allocating CrsGraph: " << e.what ());
192 staticGraph_ = myGraph_;
194 checkInternalState ();
197 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
200 const Teuchos::ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
202 const Teuchos::RCP<Teuchos::ParameterList>& params) :
207 fillComplete_ (false),
208 frobNorm_ (-STM::one ())
212 myGraph_ = rcp (
new Graph (rowMap, NumEntriesPerRowToAlloc, pftype, params));
214 catch (std::exception &e) {
215 TEUCHOS_TEST_FOR_EXCEPTION(
216 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 217 "exception while allocating CrsGraph: " << e.what ());
219 staticGraph_ = myGraph_;
224 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
227 const Teuchos::RCP<const map_type>& colMap,
228 size_t maxNumEntriesPerRow,
230 const Teuchos::RCP<Teuchos::ParameterList>& params) :
239 TEUCHOS_TEST_FOR_EXCEPTION(! staticGraph_.is_null(), std::logic_error,
240 "Tpetra::CrsMatrix ctor (row Map, col Map, maxNumEntriesPerRow, ...): " 241 "staticGraph_ is not null at the beginning of the constructor. " 242 "Please report this bug to the Tpetra developers.");
243 TEUCHOS_TEST_FOR_EXCEPTION(! myGraph_.is_null(), std::logic_error,
244 "Tpetra::CrsMatrix ctor (row Map, col Map, maxNumEntriesPerRow, ...): " 245 "myGraph_ is not null at the beginning of the constructor. " 246 "Please report this bug to the Tpetra developers.");
248 myGraph_ = rcp (
new Graph (rowMap, colMap, maxNumEntriesPerRow,
251 catch (std::exception &e) {
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 254 "exception while allocating CrsGraph: " << e.what ());
256 staticGraph_ = myGraph_;
261 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
264 const Teuchos::RCP<const map_type>& colMap,
265 const Teuchos::ArrayRCP<const size_t>& numEntPerRow,
267 const Teuchos::RCP<Teuchos::ParameterList>& params) :
277 myGraph_ = rcp (
new Graph (rowMap, colMap, numEntPerRow, pftype, params));
279 catch (std::exception &e) {
280 TEUCHOS_TEST_FOR_EXCEPTION(
281 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 282 "exception while allocating CrsGraph: " << e.what ());
284 staticGraph_ = myGraph_;
289 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
291 CrsMatrix (
const Teuchos::RCP<const crs_graph_type>& graph,
292 const Teuchos::RCP<Teuchos::ParameterList>& params) :
294 staticGraph_ (graph),
299 const char tfecfFuncName[] =
"CrsMatrix(graph[,params])";
300 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(staticGraph_.is_null (),
301 std::runtime_error,
": When calling the CrsMatrix constructor that " 302 "accepts a static graph, the pointer to the graph must not be null.");
304 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( ! staticGraph_->isFillComplete (),
305 std::runtime_error,
": The specified graph is not fill-complete. You " 306 "must invoke fillComplete() on the graph before using it to construct a " 307 "CrsMatrix. Note that calling resumeFill() makes the graph not fill-" 308 "complete, even if you had previously called fillComplete(). In that " 309 "case, you must call fillComplete() on the graph again.");
317 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
320 const Teuchos::RCP<const map_type>& colMap,
321 const typename local_matrix_type::row_map_type& rowPointers,
322 const typename local_graph_type::entries_type::non_const_type& columnIndices,
323 const typename local_matrix_type::values_type& values,
324 const Teuchos::RCP<Teuchos::ParameterList>& params) :
332 myGraph_ = rcp (
new Graph (rowMap, colMap, rowPointers,
333 columnIndices, params));
335 catch (std::exception &e) {
336 TEUCHOS_TEST_FOR_EXCEPTION(
337 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 338 "exception while allocating CrsGraph: " << e.what ());
340 staticGraph_ = myGraph_;
341 k_values1D_ = values;
346 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
348 CrsMatrix (
const Teuchos::RCP<const map_type>& rowMap,
349 const Teuchos::RCP<const map_type>& colMap,
350 const Teuchos::ArrayRCP<size_t> & rowPointers,
351 const Teuchos::ArrayRCP<LocalOrdinal> & columnIndices,
352 const Teuchos::ArrayRCP<Scalar> & values,
353 const Teuchos::RCP<Teuchos::ParameterList>& params) :
361 myGraph_ = rcp (
new Graph (rowMap, colMap, rowPointers,
362 columnIndices, params));
364 catch (std::exception &e) {
365 TEUCHOS_TEST_FOR_EXCEPTION(
366 true, std::runtime_error,
"Tpetra::CrsMatrix constructor: Caught " 367 "exception while allocating CrsGraph: " << e.what ());
369 staticGraph_ = myGraph_;
373 Teuchos::ArrayRCP<impl_scalar_type> vals =
375 k_values1D_ = Kokkos::Compat::getKokkosViewDeepCopy<device_type> (vals ());
380 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
383 const Teuchos::RCP<const map_type>& colMap,
385 const Teuchos::RCP<Teuchos::ParameterList>& params) :
392 using Teuchos::ArrayRCP;
396 const char tfecfFuncName[] =
"Tpetra::CrsMatrix(rowMap,colMap,lclMatrix,params): ";
399 myGraph_ = rcp (
new Graph (rowMap, colMap, lclMatrix.graph, params));
401 catch (std::exception &e) {
402 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
403 true, std::runtime_error,
"Caught exception while allocating " 404 "CrsGraph: " << e.what ());
406 staticGraph_ = myGraph_;
417 #ifdef HAVE_TPETRA_DEBUG 418 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillActive (), std::logic_error,
419 "We're at the end of fillComplete(), but isFillActive() is true. " 420 "Please report this bug to the Tpetra developers.");
421 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete (), std::logic_error,
422 "We're at the end of fillComplete(), but isFillComplete() is false. " 423 "Please report this bug to the Tpetra developers.");
424 #endif // HAVE_TPETRA_DEBUG 428 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
433 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
434 Teuchos::RCP<const Teuchos::Comm<int> >
440 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
447 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
454 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
461 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
468 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
475 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
482 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
489 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
496 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
503 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
510 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
517 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
524 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
531 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
538 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
545 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
552 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
556 return getCrsGraph ()->getNumEntriesInGlobalRow (globalRow);
559 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
563 return getCrsGraph ()->getNumEntriesInLocalRow (localRow);
566 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
570 return getCrsGraph ()->getGlobalMaxNumRowEntries ();
573 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
580 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
587 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
588 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
594 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
595 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
601 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
602 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
608 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
609 Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >
615 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
616 Teuchos::RCP<const RowGraph<LocalOrdinal, GlobalOrdinal, Node> >
619 if (staticGraph_ != Teuchos::null) {
625 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
626 Teuchos::RCP<const CrsGraph<LocalOrdinal, GlobalOrdinal, Node, classic> >
629 if (staticGraph_ != Teuchos::null) {
635 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
642 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
649 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
653 return myGraph_.is_null ();
656 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
663 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
670 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
675 #ifdef HAVE_TPETRA_DEBUG 679 if ((gas == GraphAlreadyAllocated) != staticGraph_->indicesAreAllocated()) {
680 const std::string err1 (
"allocateValues: The caller has asserted that " 682 const std::string err2 (
"already allocated, but the static graph says " 683 "that its indices are ");
684 const std::string err3 (
"already allocated. Please report this bug to " 685 "the Tpetra developers.");
686 TEUCHOS_TEST_FOR_EXCEPTION(gas == GraphAlreadyAllocated && ! staticGraph_->indicesAreAllocated(),
687 std::logic_error, err1 << err2 <<
"not " << err3);
688 TEUCHOS_TEST_FOR_EXCEPTION(gas != GraphAlreadyAllocated && staticGraph_->indicesAreAllocated(),
689 std::logic_error, err1 <<
"not " << err2 << err3);
697 TEUCHOS_TEST_FOR_EXCEPTION(
698 ! staticGraph_->indicesAreAllocated() && myGraph_.is_null(),
700 "allocateValues: The static graph says that its indices are not " 701 "allocated, but the graph is not owned by the matrix. Please report " 702 "this bug to the Tpetra developers.");
703 #endif // HAVE_TPETRA_DEBUG 705 if (gas == GraphNotYetAllocated) {
706 myGraph_->allocateIndices (lg);
717 const size_t lclNumRows = staticGraph_->getNodeNumRows ();
718 typename Graph::local_graph_type::row_map_type k_ptrs =
719 staticGraph_->k_rowPtrs_;
720 TEUCHOS_TEST_FOR_EXCEPTION(
721 k_ptrs.dimension_0 () != lclNumRows+1, std::logic_error,
722 "Tpetra::CrsMatrix::allocateValues: With StaticProfile, row offsets " 723 "array has length " << k_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 724 << (lclNumRows+1) <<
".");
729 const size_t lclTotalNumEntries = k_ptrs(lclNumRows);
732 typedef typename local_matrix_type::values_type values_type;
733 k_values1D_ = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
741 values2D_ = staticGraph_->template allocateValues2D<impl_scalar_type> ();
745 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
748 getAllValues (Teuchos::ArrayRCP<const size_t>& rowPointers,
749 Teuchos::ArrayRCP<const LocalOrdinal>& columnIndices,
750 Teuchos::ArrayRCP<const Scalar>& values)
const 752 const char tfecfFuncName[] =
"getAllValues: ";
753 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
754 columnIndices.size () != values.size (), std::runtime_error,
755 "Requires that columnIndices and values are the same size.");
757 RCP<const crs_graph_type> relevantGraph =
getCrsGraph ();
758 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
759 relevantGraph.is_null (), std::runtime_error,
760 "Requires that getCrsGraph() is not null.");
762 rowPointers = relevantGraph->getNodeRowPtrs ();
764 catch (std::exception &e) {
765 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
766 true, std::runtime_error,
767 "Caught exception while calling graph->getNodeRowPtrs(): " 771 columnIndices = relevantGraph->getNodePackedIndices ();
773 catch (std::exception &e) {
774 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
775 true, std::runtime_error,
776 "Caught exception while calling graph->getNodePackedIndices(): " 779 Teuchos::ArrayRCP<const impl_scalar_type> vals =
780 Kokkos::Compat::persistingView (k_values1D_);
781 values = Teuchos::arcp_reinterpret_cast<
const Scalar> (vals);
784 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
789 using Kokkos::create_mirror_view;
790 using Teuchos::arcp_const_cast;
791 using Teuchos::ArrayRCP;
795 typedef ArrayRCP<size_t>::size_type size_type;
796 typedef typename local_matrix_type::row_map_type row_map_type;
797 typedef typename Graph::t_numRowEntries_ row_entries_type;
798 typedef typename Graph::local_graph_type::entries_type::non_const_type lclinds_1d_type;
799 typedef typename local_matrix_type::values_type values_type;
803 TEUCHOS_TEST_FOR_EXCEPTION(
804 myGraph_.is_null (), std::logic_error,
"Tpetra::CrsMatrix::" 805 "fillLocalGraphAndMatrix (called from fillComplete or " 806 "expertStaticFillComplete): The nonconst graph (myGraph_) is null. This " 807 "means that the matrix has a const (a.k.a. \"static\") graph. This may " 808 "mean that fillComplete or expertStaticFillComplete has a bug, since it " 809 "should never call fillLocalGraphAndMatrix in that case. " 810 "Please report this bug to the Tpetra developers.");
820 typename row_map_type::non_const_type k_ptrs;
821 row_map_type k_ptrs_const;
822 lclinds_1d_type k_inds;
828 lclinds_1d_type k_lclInds1D_ = myGraph_->k_lclInds1D_;
834 row_entries_type k_numRowEnt = myGraph_->k_numRowEntries_;
835 typename row_entries_type::t_host h_numRowEnt = k_numRowEnt.h_view;
847 TEUCHOS_TEST_FOR_EXCEPTION(
848 static_cast<size_t> (k_numRowEnt.dimension_0 ()) != lclNumRows,
849 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix (called " 850 "from fillComplete or expertStaticFillComplete): For the " 851 "DynamicProfile branch, k_numRowEnt has the wrong length. " 852 "k_numRowEnt.dimension_0() = " << k_numRowEnt.dimension_0 ()
853 <<
" != getNodeNumRows() = " << lclNumRows <<
"");
861 size_t lclTotalNumEntries = 0;
863 typename row_map_type::non_const_type::HostMirror h_ptrs;
868 typename row_map_type::non_const_type packedRowOffsets (
"Tpetra::CrsGraph::ptr",
873 h_ptrs = create_mirror_view (packedRowOffsets);
875 for (size_type i = 0; i < static_cast<size_type> (lclNumRows); ++i) {
876 const size_t numEnt = h_numRowEnt(i);
877 lclTotalNumEntries += numEnt;
878 h_ptrs(i+1) = h_ptrs(i) + numEnt;
883 k_ptrs = packedRowOffsets;
884 k_ptrs_const = k_ptrs;
887 TEUCHOS_TEST_FOR_EXCEPTION(
888 static_cast<size_t> (k_ptrs.dimension_0 ()) != lclNumRows + 1,
889 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: In " 890 "DynamicProfile branch, after packing k_ptrs, k_ptrs.dimension_0()" 891 " = " << k_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 892 << (lclNumRows+1) <<
".");
893 TEUCHOS_TEST_FOR_EXCEPTION(
894 static_cast<size_t> (h_ptrs.dimension_0 ()) != lclNumRows + 1,
895 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: In " 896 "DynamicProfile branch, after packing h_ptrs, h_ptrs.dimension_0()" 897 " = " << h_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 898 << (lclNumRows+1) <<
".");
900 TEUCHOS_TEST_FOR_EXCEPTION(
901 k_ptrs(lclNumRows) != lclTotalNumEntries, std::logic_error,
902 "Tpetra::CrsMatrix::fillLocalGraphAndMatrix: In DynamicProfile branch, " 903 "after packing k_ptrs, k_ptrs(lclNumRows = " << lclNumRows <<
") = " <<
904 k_ptrs(lclNumRows) <<
" != total number of entries on the calling " 905 "process = " << lclTotalNumEntries <<
".");
908 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
909 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
912 typename lclinds_1d_type::HostMirror h_inds = create_mirror_view (k_inds);
913 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
916 ArrayRCP<Array<LocalOrdinal> > lclInds2D = myGraph_->lclInds2D_;
917 for (
size_t row = 0; row < lclNumRows; ++row) {
918 const size_t numEnt = h_numRowEnt(row);
919 std::copy (lclInds2D[row].begin(),
920 lclInds2D[row].begin() + numEnt,
921 h_inds.ptr_on_device() + h_ptrs(row));
922 std::copy (values2D_[row].begin(),
923 values2D_[row].begin() + numEnt,
924 h_vals.ptr_on_device() + h_ptrs(row));
931 if (k_ptrs.dimension_0 () != 0) {
932 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
933 TEUCHOS_TEST_FOR_EXCEPTION(
934 static_cast<size_t> (k_ptrs(numOffsets-1)) != k_vals.dimension_0 (),
935 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 936 "In DynamicProfile branch, after packing, k_ptrs(" << (numOffsets-1)
937 <<
") = " << k_ptrs(numOffsets-1) <<
" != k_vals.dimension_0() = " 938 << k_vals.dimension_0 () <<
".");
939 TEUCHOS_TEST_FOR_EXCEPTION(
940 static_cast<size_t> (k_ptrs(numOffsets-1)) != k_inds.dimension_0 (),
941 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 942 "In DynamicProfile branch, after packing, k_ptrs(" << (numOffsets-1)
943 <<
") = " << k_ptrs(numOffsets-1) <<
" != k_inds.dimension_0() = " 944 << k_inds.dimension_0 () <<
".");
954 typename Graph::local_graph_type::row_map_type curRowOffsets =
955 myGraph_->k_rowPtrs_;
956 TEUCHOS_TEST_FOR_EXCEPTION(
957 curRowOffsets.dimension_0 () == 0, std::logic_error,
958 "curRowOffsets has size zero, but shouldn't");
959 TEUCHOS_TEST_FOR_EXCEPTION(
960 curRowOffsets.dimension_0 () != lclNumRows + 1, std::logic_error,
961 "Tpetra::CrsMatrix::fillLocalGraphAndMatrix: curRowOffsets has size " 962 << curRowOffsets.dimension_0 () <<
" != lclNumRows + 1 = " 963 << (lclNumRows + 1) <<
".")
965 const size_t numOffsets = curRowOffsets.dimension_0 ();
967 TEUCHOS_TEST_FOR_EXCEPTION(
969 myGraph_->k_lclInds1D_.dimension_0 () != curRowOffsets(numOffsets - 1),
970 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 971 "numOffsets = " << numOffsets <<
" != 0 and " 972 "myGraph_->k_lclInds1D_.dimension_0() = " 973 << myGraph_->k_lclInds1D_.dimension_0 ()
974 <<
" != curRowOffsets(" << numOffsets <<
") = " 975 << curRowOffsets(numOffsets - 1) <<
".");
978 if (myGraph_->nodeNumEntries_ != myGraph_->nodeNumAllocated_) {
985 TEUCHOS_TEST_FOR_EXCEPTION(
986 static_cast<size_t> (k_numRowEnt.dimension_0 ()) != lclNumRows,
987 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix (called" 988 " from fillComplete or expertStaticFillComplete): In StaticProfile " 989 "unpacked branch, k_numRowEnt has the wrong length. " 990 "k_numRowEnt.dimension_0() = " << k_numRowEnt.dimension_0 ()
991 <<
" != getNodeNumRows() = " << lclNumRows <<
".");
993 if (curRowOffsets.dimension_0 () != 0) {
994 const size_t numOffsets =
995 static_cast<size_t> (curRowOffsets.dimension_0 ());
996 TEUCHOS_TEST_FOR_EXCEPTION(
997 curRowOffsets(numOffsets-1) != static_cast<size_t> (k_values1D_.dimension_0 ()),
998 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 999 "In StaticProfile branch, before allocating or packing, " 1000 "curRowOffsets(" << (numOffsets-1) <<
") = " 1001 << curRowOffsets(numOffsets - 1)
1002 <<
" != k_values1D_.dimension_0() = " 1003 << k_values1D_.dimension_0 () <<
".");
1004 TEUCHOS_TEST_FOR_EXCEPTION(
1005 static_cast<size_t> (curRowOffsets(numOffsets - 1)) !=
1006 myGraph_->k_lclInds1D_.dimension_0 (),
1007 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 1008 "In StaticProfile branch, before allocating or packing, " 1009 "curRowOffsets(" << (numOffsets-1) <<
") = " 1010 << curRowOffsets(numOffsets - 1)
1011 <<
" != myGraph_->k_lclInds1D_.dimension_0() = " 1012 << myGraph_->k_lclInds1D_.dimension_0 () <<
".");
1021 size_t lclTotalNumEntries = 0;
1023 typename row_map_type::non_const_type::HostMirror h_ptrs;
1028 typename row_map_type::non_const_type packedRowOffsets (
"Tpetra::CrsGraph::ptr",
1036 h_ptrs = create_mirror_view (packedRowOffsets);
1038 for (size_type i = 0; i < static_cast<size_type> (lclNumRows); ++i) {
1039 const size_t numEnt = h_numRowEnt(i);
1040 lclTotalNumEntries += numEnt;
1041 h_ptrs(i+1) = h_ptrs(i) + numEnt;
1046 k_ptrs = packedRowOffsets;
1047 k_ptrs_const = k_ptrs;
1050 TEUCHOS_TEST_FOR_EXCEPTION(
1051 static_cast<size_t> (k_ptrs.dimension_0 ()) != lclNumRows + 1,
1052 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: For " 1053 "the StaticProfile unpacked-but-pack branch, after packing k_ptrs, " 1054 "k_ptrs.dimension_0() = " << k_ptrs.dimension_0 () <<
" != " 1055 "lclNumRows+1 = " << (lclNumRows+1) <<
".");
1057 TEUCHOS_TEST_FOR_EXCEPTION(
1058 k_ptrs(lclNumRows) != lclTotalNumEntries, std::logic_error,
1059 "Tpetra::CrsMatrix::fillLocalGraphAndMatrix: In StaticProfile " 1060 "unpacked-but-pack branch, after filling k_ptrs, k_ptrs(lclNumRows=" 1061 << lclNumRows <<
") = " << k_ptrs(lclNumRows) <<
" != total number " 1062 "of entries on the calling process = " << lclTotalNumEntries <<
".");
1065 k_inds = lclinds_1d_type (
"Tpetra::CrsGraph::ind", lclTotalNumEntries);
1066 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1078 typedef pack_functor<
typename Graph::local_graph_type::entries_type::non_const_type,
1079 typename Graph::local_graph_type::row_map_type>
1081 inds_packer_type indsPacker (k_inds, myGraph_->k_lclInds1D_,
1082 k_ptrs, curRowOffsets);
1083 Kokkos::parallel_for (lclNumRows, indsPacker);
1087 typedef pack_functor<values_type, row_map_type> vals_packer_type;
1088 vals_packer_type valsPacker (k_vals, this->k_values1D_,
1089 k_ptrs, curRowOffsets);
1090 Kokkos::parallel_for (lclNumRows, valsPacker);
1092 TEUCHOS_TEST_FOR_EXCEPTION(
1093 k_ptrs.dimension_0 () == 0, std::logic_error,
"Tpetra::CrsMatrix::" 1094 "fillLocalGraphAndMatrix: In StaticProfile \"Optimize Storage\" = " 1095 "true branch, after packing, k_ptrs.dimension_0() = 0. This " 1096 "probably means that k_rowPtrs_ was never allocated.");
1097 if (k_ptrs.dimension_0 () != 0) {
1098 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
1099 TEUCHOS_TEST_FOR_EXCEPTION(
1100 static_cast<size_t> (k_ptrs(numOffsets - 1)) != k_vals.dimension_0 (),
1101 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 1102 "In StaticProfile \"Optimize Storage\"=true branch, after packing, " 1103 "k_ptrs(" << (numOffsets-1) <<
") = " << k_ptrs(numOffsets-1) <<
1104 " != k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1105 TEUCHOS_TEST_FOR_EXCEPTION(
1106 static_cast<size_t> (k_ptrs(numOffsets - 1)) != k_inds.dimension_0 (),
1107 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 1108 "In StaticProfile \"Optimize Storage\"=true branch, after packing, " 1109 "k_ptrs(" << (numOffsets-1) <<
") = " << k_ptrs(numOffsets-1) <<
1110 " != k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1114 k_ptrs_const = myGraph_->k_rowPtrs_;
1115 k_inds = myGraph_->k_lclInds1D_;
1116 k_vals = this->k_values1D_;
1118 TEUCHOS_TEST_FOR_EXCEPTION(
1119 k_ptrs_const.dimension_0 () == 0, std::logic_error,
"Tpetra::CrsMatrix::" 1120 "fillLocalGraphAndMatrix: In StaticProfile \"Optimize Storage\" = " 1121 "false branch, k_ptrs_const.dimension_0() = 0. This probably means that " 1122 "k_rowPtrs_ was never allocated.");
1123 if (k_ptrs_const.dimension_0 () != 0) {
1124 const size_t numOffsets =
static_cast<size_t> (k_ptrs_const.dimension_0 ());
1125 TEUCHOS_TEST_FOR_EXCEPTION(
1126 static_cast<size_t> (k_ptrs_const(numOffsets - 1)) != k_vals.dimension_0 (),
1127 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 1128 "In StaticProfile \"Optimize Storage\" = false branch, " 1129 "k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const(numOffsets - 1)
1130 <<
" != k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1131 TEUCHOS_TEST_FOR_EXCEPTION(
1132 static_cast<size_t> (k_ptrs_const(numOffsets - 1)) != k_inds.dimension_0 (),
1133 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: " 1134 "In StaticProfile \"Optimize Storage\" = false branch, " 1135 "k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const(numOffsets - 1)
1136 <<
" != k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1142 TEUCHOS_TEST_FOR_EXCEPTION(
1143 static_cast<size_t> (k_ptrs_const.dimension_0 ()) != lclNumRows + 1,
1144 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: After " 1145 "packing, k_ptrs_const.dimension_0() = " << k_ptrs_const.dimension_0 ()
1146 <<
" != lclNumRows+1 = " << (lclNumRows+1) <<
".");
1147 if (k_ptrs_const.dimension_0 () != 0) {
1148 const size_t numOffsets =
static_cast<size_t> (k_ptrs_const.dimension_0 ());
1149 TEUCHOS_TEST_FOR_EXCEPTION(
1150 static_cast<size_t> (k_ptrs_const(numOffsets - 1)) != k_vals.dimension_0 (),
1151 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: After " 1152 "packing, k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const(numOffsets-1)
1153 <<
" != k_vals.dimension_0() = " << k_vals.dimension_0 () <<
".");
1154 TEUCHOS_TEST_FOR_EXCEPTION(
1155 static_cast<size_t> (k_ptrs_const(numOffsets - 1)) != k_inds.dimension_0 (),
1156 std::logic_error,
"Tpetra::CrsMatrix::fillLocalGraphAndMatrix: After " 1157 "packing, k_ptrs_const(" << (numOffsets-1) <<
") = " << k_ptrs_const(numOffsets-1)
1158 <<
" != k_inds.dimension_0() = " << k_inds.dimension_0 () <<
".");
1165 const bool defaultOptStorage =
1167 const bool requestOptimizedStorage =
1168 (! params.is_null () && params->get (
"Optimize Storage", defaultOptStorage)) ||
1169 (params.is_null () && defaultOptStorage);
1176 if (requestOptimizedStorage) {
1182 myGraph_->lclInds2D_ = null;
1183 myGraph_->k_numRowEntries_ = row_entries_type ();
1186 this->values2D_ = null;
1189 myGraph_->k_rowPtrs_ = k_ptrs_const;
1190 myGraph_->k_lclInds1D_ = k_inds;
1191 this->k_values1D_ = k_vals;
1195 myGraph_->nodeNumAllocated_ = myGraph_->nodeNumEntries_;
1199 myGraph_->storageStatus_ = Details::STORAGE_1D_PACKED;
1211 myGraph_->lclGraph_ =
1217 myGraph_->lclGraph_);
1220 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1225 using Kokkos::create_mirror_view;
1226 using Teuchos::ArrayRCP;
1227 using Teuchos::Array;
1228 using Teuchos::null;
1231 typedef LocalOrdinal LO;
1232 typedef typename Graph::t_numRowEntries_ row_entries_type;
1233 typedef typename Graph::local_graph_type::row_map_type row_map_type;
1234 typedef typename row_map_type::non_const_type non_const_row_map_type;
1235 typedef typename local_matrix_type::values_type values_type;
1239 RCP<node_type> node = rowMap.
getNode ();
1250 ArrayRCP<Array<LO> > lclInds2D = staticGraph_->lclInds2D_;
1251 size_t nodeNumEntries = staticGraph_->nodeNumEntries_;
1252 size_t nodeNumAllocated = staticGraph_->nodeNumAllocated_;
1253 row_map_type k_rowPtrs_ = staticGraph_->lclGraph_.row_map;
1255 row_map_type k_ptrs;
1261 bool requestOptimizedStorage =
true;
1262 const bool default_OptimizeStorage =
1264 if (! params.is_null () && ! params->get (
"Optimize Storage", default_OptimizeStorage)) {
1265 requestOptimizedStorage =
false;
1272 if (! staticGraph_->isStorageOptimized () && requestOptimizedStorage) {
1274 "::fillLocalMatrix(): You requested optimized storage by setting the" 1275 "\"Optimize Storage\" flag to \"true\" in the parameter list, or by virtue" 1276 "of default behavior. However, the associated CrsGraph was filled separately" 1277 "and requested not to optimize storage. Therefore, the CrsMatrix cannot" 1278 "optimize storage.");
1279 requestOptimizedStorage =
false;
1286 row_entries_type k_numRowEnt = staticGraph_->k_numRowEntries_;
1287 typename row_entries_type::t_host h_numRowEnt = k_numRowEnt.h_view;
1312 size_t lclTotalNumEntries = 0;
1314 typename non_const_row_map_type::HostMirror h_ptrs;
1316 non_const_row_map_type packedRowOffsets (
"Tpetra::CrsGraph::ptr",
1321 h_ptrs = create_mirror_view (packedRowOffsets);
1323 for (
size_t i = 0; i < lclNumRows; ++i) {
1324 const size_t numEnt = h_numRowEnt(i);
1325 lclTotalNumEntries += numEnt;
1326 h_ptrs(i+1) = h_ptrs(i) + numEnt;
1329 k_ptrs = packedRowOffsets;
1332 TEUCHOS_TEST_FOR_EXCEPTION(
1333 static_cast<size_t> (k_ptrs.dimension_0 ()) != lclNumRows + 1,
1334 std::logic_error,
"Tpetra::CrsMatrix::fillLocalMatrix: In " 1335 "DynamicProfile branch, after packing k_ptrs, k_ptrs.dimension_0()" 1336 " = " << k_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 1337 << (lclNumRows+1) <<
".");
1338 TEUCHOS_TEST_FOR_EXCEPTION(
1339 static_cast<size_t> (h_ptrs.dimension_0 ()) != lclNumRows + 1,
1340 std::logic_error,
"Tpetra::CrsMatrix::fillLocalMatrix: In " 1341 "DynamicProfile branch, after packing h_ptrs, h_ptrs.dimension_0()" 1342 " = " << h_ptrs.dimension_0 () <<
" != (lclNumRows+1) = " 1343 << (lclNumRows+1) <<
".");
1345 TEUCHOS_TEST_FOR_EXCEPTION(
1346 k_ptrs(lclNumRows) != lclTotalNumEntries, std::logic_error,
1347 "Tpetra::CrsMatrix::fillLocalMatrix: In DynamicProfile branch, " 1348 "after packing k_ptrs, k_ptrs(lclNumRows = " << lclNumRows <<
") = " <<
1349 k_ptrs(lclNumRows) <<
" != total number of entries on the calling " 1350 "process = " << lclTotalNumEntries <<
".");
1353 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1355 typename values_type::HostMirror h_vals = create_mirror_view (k_vals);
1357 for (
size_t lclRow = 0; lclRow < lclNumRows; ++lclRow) {
1358 const size_t numEnt = h_numRowEnt(lclRow);
1359 std::copy (values2D_[lclRow].begin(),
1360 values2D_[lclRow].begin() + numEnt,
1361 h_vals.ptr_on_device() + h_ptrs(lclRow));
1367 if (k_ptrs.dimension_0 () != 0) {
1368 const size_t numOffsets =
static_cast<size_t> (k_ptrs.dimension_0 ());
1369 TEUCHOS_TEST_FOR_EXCEPTION(
1370 static_cast<size_t> (k_ptrs(numOffsets-1)) != k_vals.dimension_0 (),
1371 std::logic_error,
"Tpetra::CrsMatrix::fillLocalMatrix: " 1372 "In DynamicProfile branch, after packing, k_ptrs(" << (numOffsets-1)
1373 <<
") = " << k_ptrs(numOffsets-1) <<
" != k_vals.dimension_0() = " 1374 << k_vals.dimension_0 () <<
".");
1395 if (nodeNumEntries != nodeNumAllocated) {
1398 non_const_row_map_type tmpk_ptrs (
"Tpetra::CrsGraph::ptr",
1403 size_t lclTotalNumEntries = 0;
1409 typename non_const_row_map_type::HostMirror h_ptrs =
1410 create_mirror_view (tmpk_ptrs);
1412 for (
size_t i = 0; i < lclNumRows; ++i) {
1413 const size_t numEnt = h_numRowEnt(i);
1414 lclTotalNumEntries += numEnt;
1415 h_ptrs(i+1) = h_ptrs(i) + numEnt;
1422 k_vals = values_type (
"Tpetra::CrsMatrix::val", lclTotalNumEntries);
1425 typedef pack_functor<values_type, row_map_type> packer_type;
1426 packer_type valsPacker (k_vals, k_values1D_, tmpk_ptrs, k_rowPtrs_);
1427 Kokkos::parallel_for (lclNumRows, valsPacker);
1430 k_vals = k_values1D_;
1435 if (requestOptimizedStorage) {
1439 k_values1D_ = k_vals;
1450 staticGraph_->getLocalGraph ());
1453 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1457 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1458 const Teuchos::ArrayView<const Scalar>& values)
1460 using Teuchos::Array;
1461 using Teuchos::ArrayView;
1462 using Teuchos::av_reinterpret_cast;
1463 using Teuchos::toString;
1465 const char tfecfFuncName[] =
"insertLocalValues";
1467 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillActive (), std::runtime_error,
1468 ": Fill is not active. After calling fillComplete, you must call " 1469 "resumeFill before you may insert entries into the matrix again.");
1470 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isStaticGraph (), std::runtime_error,
1471 " cannot insert indices with static graph; use replaceLocalValues() instead.");
1472 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(myGraph_->isGloballyIndexed(),
1473 std::runtime_error,
": graph indices are global; use insertGlobalValues().");
1474 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
hasColMap (), std::runtime_error,
1475 " cannot insert local indices without a column map.");
1476 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(values.size() != indices.size(),
1477 std::runtime_error,
": values.size() must equal indices.size().");
1478 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1479 !
getRowMap()->isNodeLocalElement(localRow), std::runtime_error,
1480 ": Local row index " << localRow <<
" does not belong to this process.");
1482 if (! myGraph_->indicesAreAllocated ()) {
1486 catch (std::exception& e) {
1487 TEUCHOS_TEST_FOR_EXCEPTION(
1488 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1489 "allocateValues(LocalIndices,GraphNotYetAllocated) threw an " 1490 "exception: " << e.what ());
1494 const size_t numEntriesToAdd =
static_cast<size_t> (indices.size ());
1495 #ifdef HAVE_TPETRA_DEBUG 1502 Array<LocalOrdinal> badColInds;
1503 bool allInColMap =
true;
1504 for (
size_t k = 0; k < numEntriesToAdd; ++k) {
1506 allInColMap =
false;
1507 badColInds.push_back (indices[k]);
1510 if (! allInColMap) {
1511 std::ostringstream os;
1512 os <<
"Tpetra::CrsMatrix::insertLocalValues: You attempted to insert " 1513 "entries in owned row " << localRow <<
", at the following column " 1514 "indices: " << toString (indices) <<
"." << endl;
1515 os <<
"Of those, the following indices are not in the column Map on " 1516 "this process: " << toString (badColInds) <<
"." << endl <<
"Since " 1517 "the matrix has a column Map already, it is invalid to insert " 1518 "entries at those locations.";
1519 TEUCHOS_TEST_FOR_EXCEPTION(! allInColMap, std::invalid_argument, os.str ());
1522 #endif // HAVE_TPETRA_DEBUG 1524 #ifdef HAVE_TPETRA_DEBUG 1527 rowInfo = myGraph_->getRowInfo (localRow);
1528 }
catch (std::exception& e) {
1529 TEUCHOS_TEST_FOR_EXCEPTION(
1530 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1531 "myGraph_->getRowInfo threw an exception: " << e.what ());
1534 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1535 #endif // HAVE_TPETRA_DEBUG 1537 const size_t curNumEntries = rowInfo.numEntries;
1538 const size_t newNumEntries = curNumEntries + numEntriesToAdd;
1539 if (newNumEntries > rowInfo.allocSize) {
1540 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1542 ": new indices exceed statically allocated graph structure.");
1546 rowInfo = myGraph_->template updateLocalAllocAndValues<impl_scalar_type> (rowInfo,
1548 values2D_[localRow]);
1549 }
catch (std::exception& e) {
1550 TEUCHOS_TEST_FOR_EXCEPTION(
1551 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1552 "myGraph_->updateGlobalAllocAndValues threw an exception: " 1556 typename Graph::SLocalGlobalViews indsView;
1557 indsView.linds = indices;
1559 #ifdef HAVE_TPETRA_DEBUG 1560 ArrayView<impl_scalar_type> valsView;
1563 }
catch (std::exception& e) {
1564 TEUCHOS_TEST_FOR_EXCEPTION(
1565 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1566 "getViewNonConst threw an exception: " << e.what ());
1569 ArrayView<impl_scalar_type> valsView = this->
getViewNonConst (rowInfo);
1570 #endif // HAVE_TPETRA_DEBUG 1572 ArrayView<const impl_scalar_type> valsIn =
1575 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, indsView,
1579 }
catch (std::exception& e) {
1580 TEUCHOS_TEST_FOR_EXCEPTION(
1581 true, std::runtime_error,
"Tpetra::CrsMatrix::insertLocalValues: " 1582 "myGraph_->insertIndicesAndValues threw an exception: " 1586 #ifdef HAVE_TPETRA_DEBUG 1587 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
1588 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1589 chkNewNumEntries != newNumEntries, std::logic_error,
1590 ": The row should have " << newNumEntries <<
" entries after insert, but " 1591 "instead has " << chkNewNumEntries <<
". Please report this bug to the " 1592 "Tpetra developers.");
1593 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isLocallyIndexed(), std::logic_error,
1594 ": At end of insertLocalValues(), this CrsMatrix is not locally indexed. " 1595 "Please report this bug to the Tpetra developers.");
1596 #endif // HAVE_TPETRA_DEBUG 1599 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1603 const Teuchos::ArrayView<const LocalOrdinal>& indices,
1604 const Teuchos::ArrayView<const Scalar>& values)
1606 using Teuchos::Array;
1607 using Teuchos::ArrayView;
1608 using Teuchos::av_reinterpret_cast;
1609 const char tfecfFuncName[] =
"insertLocalValues: ";
1611 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillActive (), std::runtime_error,
1612 "Requires that fill is active.");
1613 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isStaticGraph (), std::runtime_error,
1614 "Cannot insert indices with static graph; use replaceLocalValues() instead.");
1615 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(myGraph_->isGloballyIndexed(),
1616 std::runtime_error,
"Graph indices are global; use insertGlobalValues().");
1617 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1618 !
hasColMap (), std::runtime_error,
"The matrix has no column Map yet, " 1619 "so you cannot insert local indices. If you created the matrix without " 1620 "a column Map (or without a fill-complete graph), you must call " 1621 "fillComplete to create the column Map, before you may work with local " 1623 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1624 values.size () != indices.size (), std::runtime_error,
"values.size() = " 1625 << values.size () <<
" != indices.size() = " << indices.size ()<<
".");
1626 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1627 !
getRowMap()->isNodeLocalElement (localRow), std::runtime_error,
1628 "Local row index " << localRow <<
" does not belong to this process.");
1629 if (! myGraph_->indicesAreAllocated ()) {
1634 Array<LocalOrdinal> f_inds (indices);
1635 ArrayView<const impl_scalar_type> valsIn =
1637 Array<impl_scalar_type> f_vals (valsIn);
1638 const size_t numFilteredEntries =
1639 myGraph_->template filterLocalIndicesAndValues<impl_scalar_type> (f_inds (),
1641 if (numFilteredEntries > 0) {
1642 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1643 const size_t curNumEntries = rowInfo.numEntries;
1644 const size_t newNumEntries = curNumEntries + numFilteredEntries;
1645 if (newNumEntries > rowInfo.allocSize) {
1646 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1648 ": new indices exceed statically allocated graph structure. " 1649 "newNumEntries (" << newNumEntries <<
" > rowInfo.allocSize (" 1650 << rowInfo.allocSize <<
").");
1653 myGraph_->template updateLocalAllocAndValues<impl_scalar_type> (rowInfo,
1655 values2D_[localRow]);
1657 typename Graph::SLocalGlobalViews inds_view;
1658 inds_view.linds = f_inds (0, numFilteredEntries);
1659 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
1661 f_vals, LocalIndices,
1663 #ifdef HAVE_TPETRA_DEBUG 1664 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
1665 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
1666 std::logic_error,
": Internal logic error. Please contact Tpetra team.");
1667 #endif // HAVE_TPETRA_DEBUG 1669 #ifdef HAVE_TPETRA_DEBUG 1670 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isLocallyIndexed(), std::logic_error,
1671 ": At end of insertLocalValues(), this CrsMatrix is not locally indexed. " 1672 "Please report this bug to the Tpetra developers.");
1673 #endif // HAVE_TPETRA_DEBUG 1677 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1681 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
1682 const Teuchos::ArrayView<const Scalar>& values)
1684 using Teuchos::Array;
1685 using Teuchos::ArrayView;
1686 using Teuchos::av_reinterpret_cast;
1687 using Teuchos::toString;
1689 typedef LocalOrdinal LO;
1690 typedef GlobalOrdinal GO;
1691 typedef typename ArrayView<const GO>::size_type size_type;
1692 const char tfecfFuncName[] =
"insertGlobalValues: ";
1694 #ifdef HAVE_TPETRA_DEBUG 1695 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1696 values.size () != indices.size (), std::runtime_error,
1697 "values.size() = " << values.size() <<
" != indices.size() = " 1698 << indices.size() <<
".");
1699 #endif // HAVE_TPETRA_DEBUG 1701 const LO localRow =
getRowMap ()->getLocalElement (globalRow);
1703 if (localRow == OTL::invalid ()) {
1704 insertNonownedGlobalValues (globalRow, indices, values);
1709 std::ostringstream err;
1710 const int myRank =
getRowMap ()->getComm ()->getRank ();
1711 const int numProcs =
getRowMap ()->getComm ()->getSize ();
1713 err <<
"The matrix was constructed with a constant (\"static\") graph, " 1714 "yet the given global row index " << globalRow <<
" is in the row " 1715 "Map on the calling process (with rank " << myRank <<
", of " <<
1716 numProcs <<
" process(es)). In this case, you may not insert new " 1717 "entries into rows owned by the calling process.";
1719 if (!
getRowMap ()->isNodeGlobalElement (globalRow)) {
1720 err <<
" Furthermore, GID->LID conversion with the row Map claims that " 1721 "the global row index is owned on the calling process, yet " 1722 "getRowMap()->isNodeGlobalElement(globalRow) returns false. That's" 1723 " weird! This might indicate a Map bug. Please report this to the" 1724 " Tpetra developers.";
1726 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1730 if (! myGraph_->indicesAreAllocated ()) {
1734 catch (std::exception& e) {
1735 TEUCHOS_TEST_FOR_EXCEPTION(
1736 true, std::runtime_error,
"Tpetra::CrsMatrix::insertGlobalValues: " 1737 "allocateValues(GlobalIndices,GraphNotYetAllocated) threw an " 1738 "exception: " << e.what ());
1742 const size_type numEntriesToInsert = indices.size ();
1755 #ifdef HAVE_TPETRA_DEBUG 1756 Array<GO> badColInds;
1757 #endif // HAVE_TPETRA_DEBUG 1758 bool allInColMap =
true;
1759 for (size_type k = 0; k < numEntriesToInsert; ++k) {
1761 allInColMap =
false;
1762 #ifdef HAVE_TPETRA_DEBUG 1763 badColInds.push_back (indices[k]);
1766 #endif // HAVE_TPETRA_DEBUG 1769 if (! allInColMap) {
1770 std::ostringstream os;
1771 os <<
"You attempted to insert entries in owned row " << globalRow
1772 <<
", at the following column indices: " << toString (indices)
1774 #ifdef HAVE_TPETRA_DEBUG 1775 os <<
"Of those, the following indices are not in the column Map on " 1776 "this process: " << toString (badColInds) <<
"." << endl <<
"Since " 1777 "the matrix has a column Map already, it is invalid to insert " 1778 "entries at those locations.";
1780 os <<
"At least one of those indices is not in the column Map on this " 1781 "process." << endl <<
"It is invalid to insert into columns not in " 1782 "the column Map on the process that owns the row.";
1783 #endif // HAVE_TPETRA_DEBUG 1784 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1785 ! allInColMap, std::invalid_argument, os.str ());
1789 typename Graph::SLocalGlobalViews inds_view;
1790 ArrayView<const impl_scalar_type> vals_view;
1792 inds_view.ginds = indices;
1795 #ifdef HAVE_TPETRA_DEBUG 1798 rowInfo = myGraph_->getRowInfo (localRow);
1799 }
catch (std::exception& e) {
1800 TEUCHOS_TEST_FOR_EXCEPTION(
1801 true, std::runtime_error,
"myGraph_->getRowInfo(localRow=" << localRow
1802 <<
") threw an exception: " << e.what ());
1805 RowInfo rowInfo = myGraph_->getRowInfo (localRow);
1806 #endif // HAVE_TPETRA_DEBUG 1808 const size_t curNumEntries = rowInfo.numEntries;
1809 const size_t newNumEntries =
1810 curNumEntries +
static_cast<size_t> (numEntriesToInsert);
1811 if (newNumEntries > rowInfo.allocSize) {
1812 TEUCHOS_TEST_FOR_EXCEPTION(
1814 std::runtime_error,
"Tpetra::CrsMatrix::insertGlobalValues: new " 1815 "indices exceed statically allocated graph structure. curNumEntries" 1816 " (" << curNumEntries <<
") + numEntriesToInsert (" <<
1817 numEntriesToInsert <<
") > allocSize (" << rowInfo.allocSize <<
").");
1822 myGraph_->template updateGlobalAllocAndValues<impl_scalar_type> (rowInfo,
1824 values2D_[localRow]);
1825 }
catch (std::exception& e) {
1826 TEUCHOS_TEST_FOR_EXCEPTION(
1827 true, std::runtime_error,
"myGraph_->updateGlobalAllocAndValues" 1828 "(...) threw an exception: " << e.what ());
1836 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
1839 GlobalIndices, GlobalIndices);
1845 myGraph_->template insertIndicesAndValues<impl_scalar_type> (rowInfo, inds_view,
1848 GlobalIndices, LocalIndices);
1851 catch (std::exception& e) {
1852 TEUCHOS_TEST_FOR_EXCEPTION(
1853 true, std::runtime_error,
"myGraph_->insertIndicesAndValues(...) " 1854 "threw an exception: " << e.what ());
1857 #ifdef HAVE_TPETRA_DEBUG 1858 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (localRow);
1859 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
1860 std::logic_error,
": There should be a total of " << newNumEntries
1861 <<
" entries in the row, but the graph now reports " << chkNewNumEntries
1862 <<
" entries. Please report this bug to the Tpetra developers.");
1863 #endif // HAVE_TPETRA_DEBUG 1868 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1872 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
1873 const Teuchos::ArrayView<const Scalar>& values)
1875 using Teuchos::Array;
1876 using Teuchos::ArrayView;
1877 using Teuchos::av_reinterpret_cast;
1878 typedef LocalOrdinal LO;
1879 typedef GlobalOrdinal GO;
1881 const char tfecfFuncName[] =
"insertGlobalValuesFiltered: ";
1891 #ifdef HAVE_TPETRA_DEBUG 1892 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1893 values.size () != indices.size (), std::runtime_error,
1894 "values.size() = " << values.size() <<
" != indices.size() = " 1895 << indices.size() <<
".");
1896 #endif // HAVE_TPETRA_DEBUG 1898 ArrayView<const ST> valsIn = av_reinterpret_cast<
const ST> (values);
1899 const LO lrow =
getRowMap ()->getLocalElement (globalRow);
1901 if (lrow != Teuchos::OrdinalTraits<LO>::invalid ()) {
1904 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1906 "The matrix was constructed with a static graph. In that case, " 1907 "it is forbidden to insert new entries into rows owned by the " 1908 "calling process.");
1909 if (! myGraph_->indicesAreAllocated ()) {
1912 typename Graph::SLocalGlobalViews inds_view;
1913 ArrayView<const ST> vals_view;
1918 Array<GO> filtered_indices;
1919 Array<ST> filtered_values;
1923 filtered_indices.assign (indices.begin (), indices.end ());
1924 filtered_values.assign (valsIn.begin (), valsIn.end ());
1925 const size_t numFilteredEntries =
1926 myGraph_->template filterGlobalIndicesAndValues<ST> (filtered_indices (),
1927 filtered_values ());
1928 inds_view.ginds = filtered_indices (0, numFilteredEntries);
1929 vals_view = filtered_values (0, numFilteredEntries);
1932 inds_view.ginds = indices;
1935 const size_t numFilteredEntries = vals_view.size ();
1937 if (numFilteredEntries > 0) {
1938 RowInfo rowInfo = myGraph_->getRowInfo (lrow);
1939 const size_t curNumEntries = rowInfo.numEntries;
1940 const size_t newNumEntries = curNumEntries + numFilteredEntries;
1941 if (newNumEntries > rowInfo.allocSize) {
1942 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1944 "New indices exceed statically allocated graph structure.");
1947 rowInfo = myGraph_->template updateGlobalAllocAndValues<ST> (rowInfo,
1955 myGraph_->template insertIndicesAndValues<ST> (rowInfo, inds_view,
1958 GlobalIndices, GlobalIndices);
1964 myGraph_->template insertIndicesAndValues<ST> (rowInfo, inds_view,
1967 GlobalIndices, LocalIndices);
1969 #ifdef HAVE_TPETRA_DEBUG 1971 const size_t chkNewNumEntries = myGraph_->getNumEntriesInLocalRow (lrow);
1972 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(chkNewNumEntries != newNumEntries,
1973 std::logic_error,
": There should be a total of " << newNumEntries
1974 <<
" entries in the row, but the graph now reports " << chkNewNumEntries
1975 <<
" entries. Please report this bug to the Tpetra developers.");
1977 #endif // HAVE_TPETRA_DEBUG 1981 insertNonownedGlobalValues (globalRow, indices, values);
1986 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1990 const Teuchos::ArrayView<const LocalOrdinal> &indices,
1991 const Teuchos::ArrayView<const Scalar>& values)
const 1993 using Kokkos::MemoryUnmanaged;
1996 typedef LocalOrdinal LO;
1998 typedef typename View<LO*, DD>::HostMirror::device_type HD;
1999 typedef View<const ST*, HD, MemoryUnmanaged> ISVT;
2000 typedef View<const LO*, HD, MemoryUnmanaged> LIVT;
2004 return Teuchos::OrdinalTraits<LO>::invalid ();
2007 const RowInfo rowInfo = staticGraph_->getRowInfo (localRow);
2008 if (rowInfo.localRow == Teuchos::OrdinalTraits<size_t>::invalid ()) {
2011 return static_cast<LO
> (0);
2014 auto curVals = this->getRowViewNonConst (rowInfo);
2015 typedef typename std::remove_const<typename std::remove_reference<decltype (curVals)>::type>::type OSVT;
2016 const ST* valsRaw =
reinterpret_cast<const ST*
> (values.getRawPtr ());
2017 ISVT valsIn (valsRaw, values.size ());
2018 LIVT indsIn (indices.getRawPtr (), indices.size ());
2019 return staticGraph_->template replaceLocalValues<OSVT, LIVT, ISVT> (rowInfo,
2026 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2030 const Kokkos::View<
const GlobalOrdinal*,
device_type,
2031 Kokkos::MemoryUnmanaged>& inputInds,
2033 Kokkos::MemoryUnmanaged>& inputVals)
const 2045 return this->
template transformGlobalValues<BF, DD> (globalRow, inputInds,
2046 inputVals, BF (),
false);
2050 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2054 const Teuchos::ArrayView<const GlobalOrdinal>& inputInds,
2055 const Teuchos::ArrayView<const Scalar>& inputVals)
const 2057 using Kokkos::MemoryUnmanaged;
2060 typedef GlobalOrdinal GO;
2062 typedef typename View<GO*, DD>::HostMirror::device_type HD;
2068 const ST*
const rawInputVals =
2069 reinterpret_cast<const ST*
> (inputVals.getRawPtr ());
2071 View<const ST*, HD, MemoryUnmanaged> inputValsK (rawInputVals,
2073 View<const GO*, HD, MemoryUnmanaged> inputIndsK (inputInds.getRawPtr (),
2078 return this->
template transformGlobalValues<BF, HD> (globalRow, inputIndsK,
2079 inputValsK, BF (),
false);
2083 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2087 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
2088 const Teuchos::ArrayView<const Scalar>& values,
2091 using Kokkos::MemoryUnmanaged;
2094 typedef LocalOrdinal LO;
2095 typedef GlobalOrdinal GO;
2097 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2101 return Teuchos::OrdinalTraits<LO>::invalid ();
2108 const LO lrow = this->staticGraph_.is_null () ?
2109 myGraph_->rowMap_->getLocalElement (globalRow) :
2110 staticGraph_->rowMap_->getLocalElement (globalRow);
2113 if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) {
2118 this->insertNonownedGlobalValues (globalRow, indices, values);
2125 return static_cast<LO
> (indices.size ());
2128 if (staticGraph_.is_null ()) {
2129 return Teuchos::OrdinalTraits<LO>::invalid ();
2131 const RowInfo rowInfo = this->staticGraph_->getRowInfo (lrow);
2133 auto curVals = this->getRowViewNonConst (rowInfo);
2134 const ST* valsRaw =
reinterpret_cast<const ST*
> (values.getRawPtr ());
2135 View<const ST*, HD, MemoryUnmanaged> valsIn (valsRaw, values.size ());
2136 View<const GO*, HD, MemoryUnmanaged> indsIn (indices.getRawPtr (),
2138 return staticGraph_->template sumIntoGlobalValues<ST, HD, DD> (rowInfo,
2146 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2150 const Teuchos::ArrayView<const LocalOrdinal>& indices,
2151 const Teuchos::ArrayView<const Scalar>& values,
2152 const bool atomic)
const 2154 using Kokkos::MemoryUnmanaged;
2157 typedef LocalOrdinal LO;
2159 typedef typename View<LO*, DD>::HostMirror::device_type HD;
2161 typedef View<const ST*, HD, MemoryUnmanaged> IVT;
2162 typedef View<const LO*, HD, MemoryUnmanaged> IIT;
2164 const ST* valsRaw =
reinterpret_cast<const ST*
> (values.getRawPtr ());
2165 IVT valsIn (valsRaw, values.size ());
2166 IIT indsIn (indices.getRawPtr (), indices.size ());
2167 return this->
template sumIntoLocalValues<IIT, IVT> (localRow, indsIn,
2171 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2172 Teuchos::ArrayView<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type>
2176 using Kokkos::MemoryUnmanaged;
2178 using Teuchos::ArrayView;
2180 typedef std::pair<size_t, size_t> range_type;
2182 if (k_values1D_.dimension_0 () != 0 && rowinfo.allocSize > 0) {
2183 #ifdef HAVE_TPETRA_DEBUG 2184 TEUCHOS_TEST_FOR_EXCEPTION(
2185 rowinfo.offset1D + rowinfo.allocSize > k_values1D_.dimension_0 (),
2186 std::range_error,
"Tpetra::CrsMatrix::getView: Invalid access " 2187 "to 1-D storage of values." << std::endl <<
"rowinfo.offset1D (" <<
2188 rowinfo.offset1D <<
") + rowinfo.allocSize (" << rowinfo.allocSize <<
2189 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2190 #endif // HAVE_TPETRA_DEBUG 2191 range_type range (rowinfo.offset1D, rowinfo.offset1D + rowinfo.allocSize);
2192 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
2199 subview_type sv = Kokkos::subview (subview_type (k_values1D_), range);
2200 const ST*
const sv_raw = (rowinfo.allocSize == 0) ? NULL : sv.ptr_on_device ();
2201 return ArrayView<const ST> (sv_raw, rowinfo.allocSize);
2203 else if (values2D_ != null) {
2204 return values2D_[rowinfo.localRow] ();
2207 return ArrayView<impl_scalar_type> ();
2211 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2212 Kokkos::View<const typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type*,
2214 Kokkos::MemoryUnmanaged>
2218 using Kokkos::MemoryUnmanaged;
2221 typedef View<const ST*, execution_space, MemoryUnmanaged> subview_type;
2222 typedef std::pair<size_t, size_t> range_type;
2224 if (k_values1D_.dimension_0 () != 0 && rowInfo.allocSize > 0) {
2225 #ifdef HAVE_TPETRA_DEBUG 2226 TEUCHOS_TEST_FOR_EXCEPTION(
2227 rowInfo.offset1D + rowInfo.allocSize > k_values1D_.dimension_0 (),
2228 std::range_error,
"Tpetra::CrsMatrix::getRowView: Invalid access " 2229 "to 1-D storage of values." << std::endl <<
"rowInfo.offset1D (" <<
2230 rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize <<
2231 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2232 #endif // HAVE_TPETRA_DEBUG 2233 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
2240 return Kokkos::subview (subview_type (k_values1D_), range);
2242 else if (values2D_ != null) {
2243 Teuchos::ArrayView<const ST> rowView = values2D_[rowInfo.localRow] ();
2244 return subview_type (rowView.getRawPtr (), rowView.size ());
2247 return subview_type ();
2251 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2252 Kokkos::View<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type*,
2253 typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::execution_space,
2254 Kokkos::MemoryUnmanaged>
2258 using Kokkos::MemoryUnmanaged;
2261 typedef View<ST*, execution_space, MemoryUnmanaged> subview_type;
2262 typedef std::pair<size_t, size_t> range_type;
2264 if (k_values1D_.dimension_0 () != 0 && rowInfo.allocSize > 0) {
2265 #ifdef HAVE_TPETRA_DEBUG 2266 TEUCHOS_TEST_FOR_EXCEPTION(
2267 rowInfo.offset1D + rowInfo.allocSize > k_values1D_.dimension_0 (),
2268 std::range_error,
"Tpetra::CrsMatrix::getRowViewNonConst: Invalid access " 2269 "to 1-D storage of values." << std::endl <<
"rowInfo.offset1D (" <<
2270 rowInfo.offset1D <<
") + rowInfo.allocSize (" << rowInfo.allocSize <<
2271 ") > k_values1D_.dimension_0() (" << k_values1D_.dimension_0 () <<
").");
2272 #endif // HAVE_TPETRA_DEBUG 2273 range_type range (rowInfo.offset1D, rowInfo.offset1D + rowInfo.allocSize);
2280 return Kokkos::subview (subview_type (k_values1D_), range);
2282 else if (values2D_ != null) {
2283 Teuchos::ArrayView<ST> rowView = values2D_[rowInfo.localRow] ();
2284 return subview_type (rowView.getRawPtr (), rowView.size ());
2287 return subview_type ();
2291 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2292 Teuchos::ArrayView<typename CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::impl_scalar_type>
2299 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2303 const Teuchos::ArrayView<LocalOrdinal>& indices,
2304 const Teuchos::ArrayView<Scalar>& values,
2305 size_t& numEntries)
const 2307 using Teuchos::ArrayView;
2308 using Teuchos::av_reinterpret_cast;
2309 typedef LocalOrdinal LO;
2310 typedef GlobalOrdinal GO;
2312 TEUCHOS_TEST_FOR_EXCEPTION(
2314 "Tpetra::CrsMatrix::getLocalRowCopy: The matrix is globally indexed and " 2315 "does not have a column Map yet. That means we don't have local indices " 2316 "for columns yet, so it doesn't make sense to call this method. If the " 2317 "matrix doesn't have a column Map yet, you should call fillComplete on " 2319 TEUCHOS_TEST_FOR_EXCEPTION(
2320 ! staticGraph_->hasRowInfo (), std::runtime_error,
2321 "Tpetra::CrsMatrix::getLocalRowCopy: The graph's row information was " 2322 "deleted at fillComplete().");
2324 if (! this->
getRowMap ()->isNodeLocalElement (localRow)) {
2330 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
2331 const size_t theNumEntries = rowinfo.numEntries;
2333 TEUCHOS_TEST_FOR_EXCEPTION(
2334 static_cast<size_t> (indices.size ()) < theNumEntries ||
2335 static_cast<size_t> (values.size ()) < theNumEntries,
2337 "Tpetra::CrsMatrix::getLocalRowCopy: The given row " << localRow
2338 <<
" has " << theNumEntries <<
" entries. One or both of the given " 2339 "ArrayViews are not long enough to store that many entries. indices " 2340 "can store " << indices.size() <<
" entries and values can store " 2341 << values.size() <<
" entries.");
2343 numEntries = theNumEntries;
2345 if (staticGraph_->isLocallyIndexed ()) {
2346 ArrayView<const LO> indrowview = staticGraph_->getLocalView (rowinfo);
2347 ArrayView<const Scalar> valrowview =
2348 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2349 std::copy (indrowview.begin (), indrowview.begin () + numEntries, indices.begin ());
2350 std::copy (valrowview.begin (), valrowview.begin () + numEntries, values.begin ());
2352 else if (staticGraph_->isGloballyIndexed ()) {
2353 ArrayView<const GO> indrowview = staticGraph_->getGlobalView (rowinfo);
2354 ArrayView<const Scalar> valrowview =
2355 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2356 std::copy (valrowview.begin (), valrowview.begin () + numEntries, values.begin ());
2359 for (
size_t j=0; j < numEntries; ++j) {
2368 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2372 const Teuchos::ArrayView<GlobalOrdinal>& indices,
2373 const Teuchos::ArrayView<Scalar>& values,
2374 size_t& numEntries)
const 2376 using Teuchos::ArrayView;
2377 using Teuchos::av_reinterpret_cast;
2378 typedef LocalOrdinal LO;
2379 typedef GlobalOrdinal GO;
2381 const char tfecfFuncName[] =
"getGlobalRowCopy: ";
2382 const LocalOrdinal lrow =
getRowMap ()->getLocalElement (globalRow);
2383 if (lrow == OTL::invalid ()) {
2389 const RowInfo rowinfo = staticGraph_->getRowInfo (lrow);
2390 const size_t theNumEntries = rowinfo.numEntries;
2392 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2393 static_cast<size_t> (indices.size ()) < theNumEntries ||
2394 static_cast<size_t> (values.size ()) < theNumEntries,
2396 "The given row " << globalRow <<
", corresponding to local row " << lrow
2397 <<
", has " << theNumEntries <<
" entries. One or both of the given " 2398 "ArrayView input arguments are not long enough to store that many " 2399 "entries. indices.size() = " << indices.size() <<
", values.size() = " 2400 << values.size () <<
", but the number of entries in the row is " 2401 << theNumEntries <<
".");
2404 numEntries = theNumEntries;
2406 if (staticGraph_->isGloballyIndexed ()) {
2407 ArrayView<const GO> indrowview = staticGraph_->getGlobalView (rowinfo);
2408 ArrayView<const Scalar> valrowview =
2409 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2410 std::copy (indrowview.begin (), indrowview.begin () + numEntries, indices.begin ());
2411 std::copy (valrowview.begin (), valrowview.begin () + numEntries, values.begin ());
2413 else if (staticGraph_->isLocallyIndexed ()) {
2414 ArrayView<const LO> indrowview = staticGraph_->getLocalView(rowinfo);
2415 ArrayView<const Scalar> valrowview =
2416 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2417 std::copy (valrowview.begin (), valrowview.begin () + numEntries, values.begin ());
2418 for (
size_t j = 0; j < numEntries; ++j) {
2419 indices[j] =
getColMap ()->getGlobalElement (indrowview[j]);
2423 #ifdef HAVE_TPETRA_DEBUG 2425 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2426 staticGraph_->indicesAreAllocated (), std::logic_error,
2427 "Internal logic error. Please contact Tpetra team.");
2428 #endif // HAVE_TPETRA_DEBUG 2433 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2437 Teuchos::ArrayView<const LocalOrdinal>& indices,
2438 Teuchos::ArrayView<const Scalar>& values)
const 2440 using Teuchos::ArrayView;
2441 using Teuchos::av_reinterpret_cast;
2442 typedef LocalOrdinal LO;
2444 const char tfecfFuncName[] =
"getLocalRowView: ";
2445 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2447 "its indices as global indices, so you cannot get a view with local " 2448 "column indices. If the matrix has a column Map, you may call " 2449 "getLocalRowCopy() to get local column indices; otherwise, you may get " 2450 "a view with global column indices by calling getGlobalRowCopy().");
2451 indices = Teuchos::null;
2452 values = Teuchos::null;
2453 #ifdef HAVE_TPETRA_DEBUG 2454 size_t numEntries = 0;
2455 #endif // HAVE_TPETRA_DEBUG 2456 if (
getRowMap ()->isNodeLocalElement (localRow)) {
2457 const RowInfo rowinfo = staticGraph_->getRowInfo (localRow);
2458 #ifdef HAVE_TPETRA_DEBUG 2459 numEntries = rowinfo.numEntries;
2460 #endif // HAVE_TPETRA_DEBUG 2461 if (rowinfo.numEntries > 0) {
2462 ArrayView<const LO> indTmp = staticGraph_->getLocalView (rowinfo);
2463 ArrayView<const Scalar> valTmp =
2464 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2465 indices = indTmp (0, rowinfo.numEntries);
2466 values = valTmp (0, rowinfo.numEntries);
2470 #ifdef HAVE_TPETRA_DEBUG 2471 const char suffix[] =
". This should never happen. Please report this " 2472 "bug to the Tpetra developers.";
2473 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2474 static_cast<size_t>(indices.size ()) != static_cast<size_t>(values.size ()), std::logic_error,
2475 "At the end of this method, for local row " << localRow <<
", " 2476 "indices.size() = " << indices.size () <<
" != values.size () = " 2477 << values.size () << suffix);
2478 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2479 static_cast<size_t>(indices.size ()) != static_cast<size_t>(numEntries), std::logic_error,
2480 "At the end of this method, for local row " << localRow <<
", " 2481 "indices.size() = " << indices.size () <<
" != numEntries = " 2482 << numEntries << suffix);
2484 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2485 numEntries != expectedNumEntries, std::logic_error,
2486 "At the end of this method, for local row " << localRow <<
", numEntries" 2487 " = " << numEntries <<
" != getNumEntriesInLocalRow(localRow)" 2488 " = "<< expectedNumEntries << suffix);
2489 #endif // HAVE_TPETRA_DEBUG 2492 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2496 Teuchos::ArrayView<const GlobalOrdinal>& indices,
2497 Teuchos::ArrayView<const Scalar>& values)
const 2499 using Teuchos::ArrayView;
2500 using Teuchos::av_reinterpret_cast;
2501 typedef GlobalOrdinal GO;
2502 const char tfecfFuncName[] =
"getGlobalRowView: ";
2504 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2506 "The matrix is locally indexed, so we cannot return a view of the row " 2507 "with global column indices. Use getGlobalRowCopy() instead.");
2508 indices = Teuchos::null;
2509 values = Teuchos::null;
2510 const LocalOrdinal lrow =
getRowMap ()->getLocalElement (globalRow);
2511 if (lrow != Teuchos::OrdinalTraits<LocalOrdinal>::invalid ()) {
2514 const RowInfo rowinfo = staticGraph_->getRowInfo(lrow);
2515 if (rowinfo.numEntries > 0) {
2516 ArrayView<const GO> indTmp = staticGraph_->getGlobalView (rowinfo);
2517 ArrayView<const Scalar> valTmp =
2518 av_reinterpret_cast<
const Scalar> (this->
getView (rowinfo));
2519 indices = indTmp (0, rowinfo.numEntries);
2520 values = valTmp (0, rowinfo.numEntries);
2523 #ifdef HAVE_TPETRA_DEBUG 2524 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2526 indices.size () != values.size (),
2528 "Violated stated post-conditions. Please contact Tpetra team.");
2529 #endif // HAVE_TPETRA_DEBUG 2532 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2537 typedef LocalOrdinal LO;
2538 typedef Kokkos::SparseRowView<local_matrix_type> row_view_type;
2539 typedef typename Teuchos::Array<Scalar>::size_type size_type;
2540 const char tfecfFuncName[] =
"scale: ";
2543 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2545 "Fill must be active before you may call this method. " 2546 "Please call resumeFill() to make fill active.");
2548 const size_t nlrs = staticGraph_->getNodeNumRows ();
2549 const size_t numAlloc = staticGraph_->getNodeAllocationSize ();
2550 const size_t numEntries = staticGraph_->getNodeNumEntries ();
2551 if (! staticGraph_->indicesAreAllocated () || nlrs == 0 ||
2552 numAlloc == 0 || numEntries == 0) {
2558 for (LO lclRow = 0; lclRow < lclNumRows; ++lclRow) {
2559 row_view_type row_i =
lclMatrix_.template row<typename row_view_type::size_type> (lclRow);
2560 for (LO k = 0; k < row_i.length; ++k) {
2562 row_i.value (k) *= theAlpha;
2567 for (
size_t row = 0; row < nlrs; ++row) {
2569 Teuchos::ArrayView<impl_scalar_type> rowVals = values2D_[row] ();
2570 for (size_type k = 0; k < numEnt; ++k) {
2571 rowVals[k] *= theAlpha;
2578 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2583 const char tfecfFuncName[] =
"setAllToScalar: ";
2585 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2587 "Fill must be active before you may call this method. " 2588 "Please call resumeFill() to make fill active.");
2594 const size_t nlrs = staticGraph_->getNodeNumRows(),
2595 numAlloc = staticGraph_->getNodeAllocationSize(),
2596 numEntries = staticGraph_->getNodeNumEntries();
2597 if (! staticGraph_->indicesAreAllocated () || numAlloc == 0 || numEntries == 0) {
2601 const ProfileType profType = staticGraph_->getProfileType ();
2609 for (
size_t row = 0; row < nlrs; ++row) {
2610 std::fill (values2D_[row].begin (), values2D_[row].end (), theAlpha);
2616 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2620 const typename local_graph_type::entries_type::non_const_type& columnIndices,
2621 const typename local_matrix_type::values_type& values)
2623 const char tfecfFuncName[] =
"setAllValues";
2624 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2625 columnIndices.size () != values.size (), std::runtime_error,
2626 ": columnIndices and values must have the same size. columnIndices.size() = " 2627 << columnIndices.size () <<
" != values.size() = " << values.size () <<
".");
2628 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2629 myGraph_.is_null (), std::runtime_error,
": myGraph_ must not be null.");
2632 myGraph_->setAllIndices (rowPointers, columnIndices);
2634 catch (std::exception &e) {
2635 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2636 true, std::runtime_error,
": Caught exception while calling myGraph_->" 2637 "setAllIndices(): " << e.what ());
2639 k_values1D_ = values;
2643 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2647 const Teuchos::ArrayRCP<LocalOrdinal>& columnIndices,
2648 const Teuchos::ArrayRCP<Scalar>& values)
2650 const char tfecfFuncName[] =
"setAllValues: ";
2651 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2652 columnIndices.size () != values.size (), std::runtime_error,
2653 "columnIndices.size() = " << columnIndices.size () <<
" != " 2654 "values.size() = " << values.size () <<
".");
2655 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2656 myGraph_.is_null (), std::runtime_error,
"myGraph_ must not be null.");
2659 myGraph_->setAllIndices (rowPointers, columnIndices);
2661 catch (std::exception &e) {
2662 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2663 true, std::runtime_error,
"Caught exception while calling myGraph_->" 2664 "setAllIndices(): " << e.what ());
2666 Teuchos::ArrayRCP<impl_scalar_type> vals =
2668 k_values1D_ = Kokkos::Compat::getKokkosViewDeepCopy<device_type> (vals ());
2672 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2677 using Teuchos::ArrayRCP;
2678 using Teuchos::ArrayView;
2679 typedef LocalOrdinal LO;
2680 const char tfecfFuncName[] =
"getLocalDiagOffsets";
2682 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2684 ": This method requires that the matrix have a column Map.");
2685 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2686 staticGraph_.is_null (), std::runtime_error,
2687 ": This method requires that the matrix have a graph.");
2693 if (static_cast<LO> (offsets.size ()) != myNumRows) {
2694 offsets.resize (myNumRows);
2697 #ifdef HAVE_TPETRA_DEBUG 2698 bool allRowMapDiagEntriesInColMap =
true;
2699 bool allDiagEntriesFound =
true;
2700 #endif // HAVE_TPETRA_DEBUG 2704 for (LO r = 0; r < myNumRows; ++r) {
2708 #ifdef HAVE_TPETRA_DEBUG 2709 if (rlid == Teuchos::OrdinalTraits<LO>::invalid ()) {
2710 allRowMapDiagEntriesInColMap =
false;
2712 #endif // HAVE_TPETRA_DEBUG 2714 if (rlid != Teuchos::OrdinalTraits<LO>::invalid ()) {
2715 const RowInfo rowinfo = staticGraph_->getRowInfo (r);
2716 if (rowinfo.numEntries > 0) {
2717 offsets[r] = staticGraph_->findLocalIndex (rowinfo, rlid);
2720 offsets[r] = Teuchos::OrdinalTraits<size_t>::invalid ();
2721 #ifdef HAVE_TPETRA_DEBUG 2722 allDiagEntriesFound =
false;
2723 #endif // HAVE_TPETRA_DEBUG 2728 #ifdef HAVE_TPETRA_DEBUG 2729 using Teuchos::reduceAll;
2732 const bool localSuccess =
2733 allRowMapDiagEntriesInColMap && allDiagEntriesFound;
2734 int localResults[3];
2735 localResults[0] = allRowMapDiagEntriesInColMap ? 1 : 0;
2736 localResults[1] = allDiagEntriesFound ? 1 : 0;
2740 ! localSuccess ?
getComm ()->getRank () :
getComm ()->getSize ();
2741 int globalResults[3];
2742 globalResults[0] = 0;
2743 globalResults[1] = 0;
2744 globalResults[2] = 0;
2745 reduceAll<int, int> (* (
getComm ()), Teuchos::REDUCE_MIN,
2746 3, localResults, globalResults);
2747 if (globalResults[0] == 0 || globalResults[1] == 0) {
2748 std::ostringstream os;
2750 globalResults[0] == 0 && globalResults[1] == 0;
2751 os <<
": At least one process (including Process " << globalResults[2]
2752 <<
") had the following issue" << (both ?
"s" :
"") <<
":" << endl;
2753 if (globalResults[0] == 0) {
2754 os <<
" - The column Map does not contain at least one diagonal entry " 2755 "of the matrix." << endl;
2757 if (globalResults[1] == 0) {
2758 os <<
" - There is a row on that / those process(es) that does not " 2759 "contain a diagonal entry." << endl;
2761 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::runtime_error, os.str());
2763 #endif // HAVE_TPETRA_DEBUG 2766 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2771 using Teuchos::ArrayRCP;
2772 using Teuchos::ArrayView;
2773 using Teuchos::av_reinterpret_cast;
2774 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
2776 typedef typename vec_type::dual_view_type dual_view_type;
2777 typedef typename dual_view_type::host_mirror_space::execution_space host_execution_space;
2779 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2781 "This method requires that the matrix have a column Map.");
2782 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2783 staticGraph_.is_null (), std::runtime_error,
2784 "This method requires that the matrix have a graph.");
2788 #ifdef HAVE_TPETRA_DEBUG 2791 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2792 ! dvec.
getMap ()->isCompatible (rowMap), std::runtime_error,
2793 ": The input Vector's Map must be compatible with the CrsMatrix's row " 2794 "Map. You may check this by using Map's isCompatible method: " 2795 "dvec.getMap ()->isCompatible (A.getRowMap ());");
2796 #endif // HAVE_TPETRA_DEBUG 2802 lclVec.template modify<host_execution_space> ();
2803 typedef typename dual_view_type::t_host host_view_type;
2804 host_view_type lclVecHost = lclVec.h_view;
2810 typename host_view_type::array_layout,
2811 typename host_view_type::device_type,
2812 typename host_view_type::memory_traits>
2814 host_view_1d_type lclVecHost1d =
2815 Kokkos::subview (lclVecHost, Kokkos::ALL (), 0);
2818 const LocalOrdinal myNumRows =
2820 for (LocalOrdinal r = 0; r < myNumRows; ++r) {
2821 lclVecHost1d(r) = STS::zero ();
2822 const GlobalOrdinal rgid = rowMap.getGlobalElement (r);
2825 if (rlid != Teuchos::OrdinalTraits<LocalOrdinal>::invalid ()) {
2826 const RowInfo rowinfo = staticGraph_->getRowInfo (r);
2827 if (rowinfo.numEntries > 0) {
2828 const size_t j = staticGraph_->findLocalIndex (rowinfo, rlid);
2829 if (j != Teuchos::OrdinalTraits<size_t>::invalid ()) {
2833 ArrayView<const impl_scalar_type> view = this->
getView (rowinfo);
2834 lclVecHost1d(r) = view[j];
2839 lclVec.template sync<execution_space> ();
2842 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2846 const Teuchos::ArrayView<const size_t>& offsets)
const 2848 using Teuchos::ArrayRCP;
2849 using Teuchos::ArrayView;
2851 typedef typename vec_type::dual_view_type dual_view_type;
2852 typedef typename dual_view_type::host_mirror_space::execution_space host_execution_space;
2854 #ifdef HAVE_TPETRA_DEBUG 2855 const char tfecfFuncName[] =
"getLocalDiagCopy: ";
2859 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2860 ! diag.
getMap ()->isCompatible (rowMap), std::runtime_error,
2861 "The input Vector's Map must be compatible with (in the sense of Map::" 2862 "isCompatible) the CrsMatrix's row Map.");
2863 #endif // HAVE_TPETRA_DEBUG 2869 lclVec.template modify<host_execution_space> ();
2870 typedef typename dual_view_type::t_host host_view_type;
2871 host_view_type lclVecHost = lclVec.h_view;
2877 typename host_view_type::array_layout,
2878 typename host_view_type::device_type,
2879 typename host_view_type::memory_traits>
2881 host_view_1d_type lclVecHost1d =
2882 Kokkos::subview (lclVecHost, Kokkos::ALL (), 0);
2884 Kokkos::View<const size_t*, Kokkos::HostSpace, Kokkos::MemoryTraits<Kokkos::Unmanaged> >
2885 h_offsets(&offsets[0],offsets.size());
2887 const LocalOrdinal myNumRows =
2889 typedef Kokkos::RangePolicy<host_execution_space, LocalOrdinal> policy_type;
2890 Kokkos::parallel_for (policy_type (0, myNumRows), [&] (
const LocalOrdinal& lclRow) {
2891 lclVecHost1d(lclRow) = STS::zero ();
2892 if (h_offsets[lclRow] != Teuchos::OrdinalTraits<size_t>::invalid ()) {
2899 auto curRow =
lclMatrix_.template rowConst<size_t>(lclRow);
2900 lclVecHost1d(lclRow) =
static_cast<impl_scalar_type
> (curRow.value(h_offsets[lclRow]));
2903 lclVec.template sync<execution_space> ();
2907 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2912 using Teuchos::ArrayRCP;
2913 using Teuchos::ArrayView;
2914 using Teuchos::null;
2917 using Teuchos::rcpFromRef;
2919 const char tfecfFuncName[] =
"leftScale";
2923 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2925 ": matrix must be fill complete.");
2926 RCP<const vec_type> xp;
2934 RCP<vec_type> tempVec = rcp (
new vec_type (
getRowMap ()));
2939 xp = rcpFromRef (x);
2943 xp = rcpFromRef (x);
2946 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::invalid_argument,
": The " 2947 "input scaling vector x's Map must be the same as either the row Map or " 2948 "the range Map of the CrsMatrix.");
2950 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
2951 ArrayView<impl_scalar_type> rowValues = null;
2953 const LocalOrdinal lclNumRows =
2955 for (LocalOrdinal i = 0; i < lclNumRows; ++i) {
2956 const RowInfo rowinfo = staticGraph_->getRowInfo (i);
2959 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
2960 rowValues[j] *= scaleValue;
2965 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2970 using Teuchos::ArrayRCP;
2971 using Teuchos::ArrayView;
2972 using Teuchos::null;
2975 using Teuchos::rcpFromRef;
2977 const char tfecfFuncName[] =
"rightScale: ";
2981 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2982 !
isFillComplete (), std::runtime_error,
"Matrix must be fill complete.");
2983 RCP<const vec_type> xp;
2989 RCP<vec_type> tempVec = rcp (
new vec_type (
getColMap ()));
2994 xp = rcpFromRef (x);
2998 xp = rcpFromRef (x);
3000 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3001 true, std::runtime_error,
"The vector x must have the same Map as " 3002 "either the row Map or the range Map.");
3005 ArrayRCP<const Scalar> vectorVals = xp->getData (0);
3006 ArrayView<impl_scalar_type> rowValues = null;
3008 const LocalOrdinal lclNumRows =
3010 for (LocalOrdinal i = 0; i < lclNumRows; ++i) {
3011 const RowInfo rowinfo = staticGraph_->getRowInfo (i);
3013 ArrayView<const LocalOrdinal> colInds;
3015 for (
size_t j = 0; j < rowinfo.numEntries; ++j) {
3021 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3026 using Teuchos::outArg;
3027 using Teuchos::REDUCE_SUM;
3028 using Teuchos::reduceAll;
3029 typedef typename Teuchos::ArrayRCP<const impl_scalar_type>::size_type size_type;
3037 if (frobNorm == -STM::one ()) {
3043 const size_type numEntries =
3045 for (size_type k = 0; k < numEntries; ++k) {
3051 const mag_type val_abs = STS::abs (val);
3052 mySum += val_abs * val_abs;
3056 const LocalOrdinal numRows =
3058 for (LocalOrdinal r = 0; r < numRows; ++r) {
3059 const RowInfo rowInfo = myGraph_->getRowInfo (r);
3060 const size_type numEntries =
3061 static_cast<size_type
> (rowInfo.numEntries);
3062 ArrayView<const impl_scalar_type> A_r =
3063 this->
getView (rowInfo).view (0, numEntries);
3064 for (size_type k = 0; k < numEntries; ++k) {
3066 const mag_type val_abs = STS::abs (val);
3067 mySum += val_abs * val_abs;
3073 reduceAll<int, mag_type> (* (
getComm ()), REDUCE_SUM,
3074 mySum, outArg (totalSum));
3075 frobNorm = STM::sqrt (totalSum);
3086 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3091 const char tfecfFuncName[] =
"replaceColMap: ";
3095 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3096 myGraph_.is_null (), std::runtime_error,
3097 "This method does not work if the matrix has a const graph. The whole " 3098 "idea of a const graph is that you are not allowed to change it, but " 3099 "this method necessarily must modify the graph, since the graph owns " 3100 "the matrix's column Map.");
3101 myGraph_->replaceColMap (newColMap);
3104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3108 const Teuchos::RCP<const map_type>& newColMap,
3109 const Teuchos::RCP<const import_type>& newImport,
3110 const bool sortEachRow)
3112 const char tfecfFuncName[] =
"reindexColumns: ";
3113 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3114 graph == NULL && myGraph_.is_null (), std::invalid_argument,
3115 "The input graph is NULL, but the matrix does not own its graph.");
3118 const bool sortGraph =
false;
3120 if (sortEachRow && theGraph.isLocallyIndexed () && ! theGraph.isSorted ()) {
3124 const LocalOrdinal lclNumRows =
3125 static_cast<LocalOrdinal
> (theGraph.getNodeNumRows ());
3126 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
3127 const RowInfo rowInfo = theGraph.getRowInfo (row);
3128 Teuchos::ArrayView<impl_scalar_type> rv = this->
getViewNonConst (rowInfo);
3129 theGraph.template sortRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
3131 theGraph.indicesAreSorted_ =
true;
3135 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3139 Teuchos::RCP<const import_type>& newImporter)
3141 const char tfecfFuncName[] =
"replaceDomainMapAndImporter: ";
3142 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3143 myGraph_.is_null (), std::runtime_error,
3144 "This method does not work if the matrix has a const graph. The whole " 3145 "idea of a const graph is that you are not allowed to change it, but this" 3146 " method necessarily must modify the graph, since the graph owns the " 3147 "matrix's domain Map and Import objects.");
3148 myGraph_->replaceDomainMapAndImporter (newDomainMap, newImporter);
3151 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3155 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
3156 const Teuchos::ArrayView<const Scalar>& values)
3158 using Teuchos::Array;
3159 typedef GlobalOrdinal GO;
3160 typedef typename Array<GO>::size_type size_type;
3162 const size_type numToInsert = indices.size ();
3165 std::pair<Array<GO>, Array<Scalar> >& curRow =
nonlocals_[globalRow];
3166 Array<GO>& curRowInds = curRow.first;
3167 Array<Scalar>& curRowVals = curRow.second;
3168 const size_type newCapacity = curRowInds.size () + numToInsert;
3169 curRowInds.reserve (newCapacity);
3170 curRowVals.reserve (newCapacity);
3171 for (size_type k = 0; k < numToInsert; ++k) {
3172 curRowInds.push_back (indices[k]);
3173 curRowVals.push_back (values[k]);
3177 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3182 using Teuchos::arcp;
3183 using Teuchos::Array;
3184 using Teuchos::ArrayRCP;
3185 using Teuchos::ArrayView;
3186 using Teuchos::CommRequest;
3187 using Teuchos::gatherAll;
3188 using Teuchos::isend;
3189 using Teuchos::ireceive;
3190 using Teuchos::null;
3191 using Teuchos::outArg;
3193 using Teuchos::rcpFromRef;
3194 using Teuchos::REDUCE_MAX;
3195 using Teuchos::reduceAll;
3196 using Teuchos::SerialDenseMatrix;
3197 using Teuchos::tuple;
3198 using Teuchos::waitAll;
3199 using std::make_pair;
3201 typedef GlobalOrdinal GO;
3202 typedef typename Array<GO>::size_type size_type;
3205 typedef std::map<GO, pair<Array<GO>, Array<Scalar> > > nonlocals_map_type;
3206 typedef typename nonlocals_map_type::const_iterator nonlocals_iter_type;
3208 const char tfecfFuncName[] =
"globalAssemble";
3209 const Teuchos::Comm<int>& comm = * (
getComm ());
3210 const int numImages = comm.getSize ();
3211 const int myImageID = comm.getRank ();
3213 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3214 !
isFillActive (), std::runtime_error,
": requires that fill is active.");
3220 size_t MyNonlocals =
static_cast<size_t> (
nonlocals_.size ());
3221 size_t MaxGlobalNonlocals = 0;
3222 reduceAll<int, size_t> (comm, REDUCE_MAX, MyNonlocals,
3223 outArg (MaxGlobalNonlocals));
3224 if (MaxGlobalNonlocals == 0) {
3242 Array<pair<int,GlobalOrdinal> > IdsAndRows;
3243 std::map<GlobalOrdinal,int> NLR2Id;
3244 SerialDenseMatrix<int,char> globalNeighbors;
3245 Array<int> sendIDs, recvIDs;
3249 std::set<GlobalOrdinal> setOfRows;
3250 for (nonlocals_iter_type iter =
nonlocals_.begin ();
3252 setOfRows.insert (iter->first);
3255 Array<GlobalOrdinal> NLRs (setOfRows.size ());
3256 std::copy (setOfRows.begin (), setOfRows.end (), NLRs.begin ());
3259 Array<int> NLRIds (NLRs.size ());
3262 getRowMap ()->getRemoteIndexList (NLRs (), NLRIds ());
3265 reduceAll<int, int> (comm, REDUCE_MAX, lclerr, outArg (gblerr));
3266 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3267 gblerr, std::runtime_error,
": non-local entries correspond to " 3274 IdsAndRows.reserve (NLRs.size ());
3275 Array<char> localNeighbors (numImages, 0);
3276 typename Array<GO>::const_iterator nlr;
3277 typename Array<int>::const_iterator id;
3278 for (nlr = NLRs.begin (),
id = NLRIds.begin ();
3279 nlr != NLRs.end (); ++nlr, ++id) {
3281 localNeighbors[*id] = 1;
3282 IdsAndRows.push_back (make_pair (*
id, *nlr));
3284 for (
int j = 0; j < numImages; ++j) {
3285 if (localNeighbors[j]) {
3286 sendIDs.push_back (j);
3290 std::sort (IdsAndRows.begin (), IdsAndRows.end ());
3304 globalNeighbors.shapeUninitialized (numImages, numImages);
3305 gatherAll (comm, numImages, localNeighbors.getRawPtr (),
3306 numImages*numImages, globalNeighbors.values ());
3318 for (
int j=0; j<numImages; ++j) {
3319 if (globalNeighbors (myImageID, j)) {
3320 recvIDs.push_back (j);
3323 const size_t numRecvs = recvIDs.size ();
3328 Array<Details::CrsIJV<GlobalOrdinal, Scalar> > IJVSendBuffer;
3329 Array<size_t> sendSizes (sendIDs.size(), 0);
3330 size_t numSends = 0;
3331 for (
typename Array<pair<int, GlobalOrdinal> >::const_iterator IdAndRow = IdsAndRows.begin();
3332 IdAndRow != IdsAndRows.end(); ++IdAndRow)
3334 const int id = IdAndRow->first;
3335 const GO row = IdAndRow->second;
3338 if (sendIDs[numSends] !=
id) {
3340 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3341 sendIDs[numSends] !=
id, std::logic_error,
3342 ": internal logic error. Contact Tpetra team.");
3346 pair<Array<GO>, Array<Scalar> >& nonlocalsRow =
nonlocals_[row];
3347 ArrayView<const GO> nonlocalsRow_colInds = nonlocalsRow.first ();
3348 ArrayView<const Scalar> nonlocalsRow_values = nonlocalsRow.second ();
3349 const size_type numNonlocalsRow = nonlocalsRow_colInds.size ();
3351 for (size_type k = 0; k < numNonlocalsRow; ++k) {
3352 const Scalar val = nonlocalsRow_values[k];
3353 const GO col = nonlocalsRow_colInds[k];
3355 sendSizes[numSends]++;
3358 if (IdsAndRows.size () > 0) {
3361 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3362 static_cast<size_type> (numSends) != sendIDs.size(),
3363 std::logic_error,
": internal logic error. Contact Tpetra team.");
3373 Array<RCP<CommRequest<int> > > sendRequests;
3374 for (
size_t s = 0; s < numSends ; ++s) {
3376 sendRequests.push_back (isend<int, size_t> (comm, rcpFromRef (sendSizes[s]), sendIDs[s]));
3379 Array<RCP<CommRequest<int> > > recvRequests;
3380 Array<size_t> recvSizes (numRecvs);
3381 for (
size_t r = 0; r < numRecvs; ++r) {
3384 recvRequests.push_back (ireceive<int, size_t> (comm, rcpFromRef (recvSizes[r]), recvIDs[r]));
3387 if (! sendRequests.empty ()) {
3388 waitAll (comm, sendRequests ());
3390 if (! recvRequests.empty ()) {
3391 waitAll (comm, recvRequests ());
3394 sendRequests.clear ();
3395 recvRequests.clear ();
3401 Array<ArrayView<Details::CrsIJV<GO, Scalar> > > sendBuffers (numSends, null);
3404 for (
size_t s=0; s<numSends; ++s) {
3405 sendBuffers[s] = IJVSendBuffer (cur, sendSizes[s]);
3406 cur += sendSizes[s];
3410 for (
size_t s = 0; s < numSends; ++s) {
3413 ArrayRCP<Details::CrsIJV<GO, Scalar> > tmparcp =
3414 arcp (sendBuffers[s].getRawPtr (), 0, sendBuffers[s].size (),
false);
3419 size_t totalRecvSize = std::accumulate (recvSizes.begin (), recvSizes.end (), 0);
3420 Array<Details::CrsIJV<GO, Scalar> > IJVRecvBuffer (totalRecvSize);
3422 Array<ArrayView<Details::CrsIJV<GO, Scalar> > > recvBuffers (numRecvs, null);
3425 for (
size_t r = 0; r < numRecvs; ++r) {
3426 recvBuffers[r] = IJVRecvBuffer (cur, recvSizes[r]);
3427 cur += recvSizes[r];
3431 for (
size_t r = 0; r < numRecvs ; ++r) {
3434 ArrayRCP<Details::CrsIJV<GO, Scalar> > tmparcp =
3435 arcp (recvBuffers[r].getRawPtr (), 0, recvBuffers[r].size (),
false);
3436 recvRequests.push_back (ireceive (comm, tmparcp, recvIDs[r]));
3439 if (! sendRequests.empty ()) {
3440 waitAll (comm, sendRequests ());
3442 if (! recvRequests.empty ()) {
3443 waitAll (comm, recvRequests ());
3446 sendRequests.clear ();
3447 recvRequests.clear ();
3457 typedef typename Array<Details::CrsIJV<GO, Scalar> >::const_iterator ijv_iter_type;
3459 for (ijv_iter_type ijv = IJVRecvBuffer.begin ();
3460 ijv != IJVRecvBuffer.end (); ++ijv) {
3465 for (ijv_iter_type ijv = IJVRecvBuffer.begin ();
3466 ijv != IJVRecvBuffer.end (); ++ijv) {
3470 catch (std::runtime_error &e) {
3471 std::ostringstream outmsg;
3472 outmsg << e.what() << std::endl
3473 <<
"caught in globalAssemble() in " << __FILE__ <<
":" << __LINE__
3475 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, outmsg.str());
3483 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3489 myGraph_->resumeFill (params);
3495 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3509 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3524 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3529 TEUCHOS_TEST_FOR_EXCEPTION(
3530 getCrsGraph ().is_null (), std::logic_error,
"Tpetra::CrsMatrix::" 3531 "fillComplete(params): getCrsGraph() returns null. " 3532 "This should not happen at this point. " 3533 "Please report this bug to the Tpetra developers.");
3543 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3546 fillComplete (
const Teuchos::RCP<const map_type>& domainMap,
3547 const Teuchos::RCP<const map_type>& rangeMap,
3548 const Teuchos::RCP<Teuchos::ParameterList>& params)
3550 using Teuchos::ArrayRCP;
3553 const char tfecfFuncName[] =
"fillComplete";
3555 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3557 std::runtime_error,
": Matrix fill state must be active (isFillActive() " 3558 "must be true) before you may call fillComplete().");
3559 const int numProcs =
getComm ()->getSize ();
3567 bool assertNoNonlocalInserts =
false;
3570 bool sortGhosts =
true;
3572 if (! params.is_null ()) {
3573 assertNoNonlocalInserts = params->get (
"No Nonlocal Changes",
3574 assertNoNonlocalInserts);
3575 if (params->isParameter (
"sort column map ghost gids")) {
3576 sortGhosts = params->get (
"sort column map ghost gids", sortGhosts);
3578 else if (params->isParameter (
"Sort column Map ghost GIDs")) {
3579 sortGhosts = params->get (
"Sort column Map ghost GIDs", sortGhosts);
3584 const bool needGlobalAssemble = ! assertNoNonlocalInserts && numProcs > 1;
3586 if (! myGraph_.is_null ()) {
3587 myGraph_->sortGhostsAssociatedWithEachProcessor_ = sortGhosts;
3601 if (needGlobalAssemble) {
3605 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3607 std::runtime_error,
": cannot have nonlocal entries on a serial run. " 3608 "An invalid entry (i.e., with row index not in the row Map) must have " 3609 "been submitted to the CrsMatrix.");
3630 const bool domainMapsMatch = staticGraph_->getDomainMap ()->isSameAs (*domainMap);
3631 const bool rangeMapsMatch = staticGraph_->getRangeMap ()->isSameAs (*rangeMap);
3633 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3634 ! domainMapsMatch, std::runtime_error,
3635 ": The CrsMatrix's domain Map does not match the graph's domain Map. " 3636 "The graph cannot be changed because it was given to the CrsMatrix " 3637 "constructor as const. You can fix this by passing in the graph's " 3638 "domain Map and range Map to the matrix's fillComplete call.");
3640 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3641 ! rangeMapsMatch, std::runtime_error,
3642 ": The CrsMatrix's range Map does not match the graph's range Map. " 3643 "The graph cannot be changed because it was given to the CrsMatrix " 3644 "constructor as const. You can fix this by passing in the graph's " 3645 "domain Map and range Map to the matrix's fillComplete call.");
3652 myGraph_->setDomainRangeMaps (domainMap, rangeMap);
3655 if (! myGraph_->hasColMap ()) {
3656 myGraph_->makeColMap ();
3661 myGraph_->makeIndicesLocal ();
3663 if (! myGraph_->isSorted ()) {
3666 if (! myGraph_->isMerged ()) {
3670 myGraph_->makeImportExport ();
3671 myGraph_->computeGlobalConstants ();
3672 myGraph_->fillComplete_ =
true;
3673 myGraph_->checkInternalState ();
3677 if (myGraph_.is_null ()) {
3698 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3702 const Teuchos::RCP<const map_type> & rangeMap,
3703 const Teuchos::RCP<const import_type>& importer,
3704 const Teuchos::RCP<const export_type>& exporter,
3705 const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
3707 #ifdef HAVE_TPETRA_MMM_TIMINGS 3709 if(!params.is_null())
3710 label = params->get(
"Timer Label",label);
3711 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
3712 using Teuchos::TimeMonitor;
3713 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-Graph"))));
3716 const char tfecfFuncName[] =
"expertStaticFillComplete: ";
3718 std::runtime_error,
"Matrix fill state must be active (isFillActive() " 3719 "must be true) before calling fillComplete().");
3720 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3721 myGraph_.is_null (), std::logic_error,
"myGraph_ is null. This is not allowed.");
3725 myGraph_->expertStaticFillComplete (domainMap, rangeMap, importer, exporter,params);
3727 #ifdef HAVE_TPETRA_MMM_TIMINGS 3728 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cGC"))));
3733 #ifdef HAVE_TPETRA_MMM_TIMINGS 3734 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-fLGAM"))));
3746 #ifdef HAVE_TPETRA_DEBUG 3747 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
isFillActive(), std::logic_error,
3748 ": We're at the end of fillComplete(), but isFillActive() is true. " 3749 "Please report this bug to the Tpetra developers.");
3750 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!
isFillComplete(), std::logic_error,
3751 ": We're at the end of fillComplete(), but isFillActive() is true. " 3752 "Please report this bug to the Tpetra developers.");
3753 #endif // HAVE_TPETRA_DEBUG 3755 #ifdef HAVE_TPETRA_MMM_TIMINGS 3756 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"ESFC-M-cIS"))));
3762 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3767 TEUCHOS_TEST_FOR_EXCEPTION(
3768 isStaticGraph (), std::runtime_error,
"Tpetra::CrsMatrix::sortEntries: " 3769 "Cannot sort with static graph.");
3770 if (! myGraph_->isSorted ()) {
3771 const LocalOrdinal lclNumRows =
3773 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
3774 const RowInfo rowInfo = myGraph_->getRowInfo (row);
3775 Teuchos::ArrayView<impl_scalar_type> rv = this->
getViewNonConst (rowInfo);
3776 myGraph_->template sortRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
3779 myGraph_->indicesAreSorted_ =
true;
3783 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3788 TEUCHOS_TEST_FOR_EXCEPTION(
3790 "mergeRedundantEntries: Cannot merge with static graph.");
3791 if (! myGraph_->isMerged ()) {
3792 const LocalOrdinal lclNumRows =
3794 for (LocalOrdinal row = 0; row < lclNumRows; ++row) {
3795 const RowInfo rowInfo = myGraph_->getRowInfo (row);
3796 Teuchos::ArrayView<impl_scalar_type> rv = this->
getViewNonConst (rowInfo);
3797 myGraph_->template mergeRowIndicesAndValues<impl_scalar_type> (rowInfo, rv);
3799 myGraph_->noRedundancies_ =
true;
3803 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3811 using Teuchos::null;
3814 using Teuchos::rcp_const_cast;
3815 using Teuchos::rcpFromRef;
3816 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
3817 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
3823 if (alpha == ZERO) {
3826 }
else if (beta != ONE) {
3840 RCP<const import_type> importer = this->
getGraph ()->getImporter ();
3841 RCP<const export_type> exporter = this->
getGraph ()->getExporter ();
3847 const bool Y_is_overwritten = (beta ==
ZERO);
3858 if (Y_is_replicated && this->
getComm ()->getRank () > 0) {
3865 RCP<const MV> X_colMap;
3866 if (importer.is_null ()) {
3876 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
3881 X_colMap = rcpFromRef (X_in);
3892 X_colMapNonConst->doImport (X_in, *importer,
INSERT);
3893 X_colMap = rcp_const_cast<
const MV> (X_colMapNonConst);
3905 if (! exporter.is_null ()) {
3906 this->
template localMultiply<Scalar, Scalar> (*X_colMap, *Y_rowMap,
3913 if (Y_is_overwritten) {
3945 this->
template localMultiply<Scalar, Scalar> (*X_colMap,
3952 this->
template localMultiply<Scalar, Scalar> (*X_colMap, Y_in,
3962 if (Y_is_replicated) {
3967 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3972 const Teuchos::ETransp mode,
3976 using Teuchos::null;
3979 using Teuchos::rcp_const_cast;
3980 using Teuchos::rcpFromRef;
3981 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
3984 if (alpha == ZERO) {
4007 RCP<const import_type> importer = this->
getGraph ()->getImporter ();
4008 RCP<const export_type> exporter = this->
getGraph ()->getExporter ();
4014 const bool Y_is_overwritten = (beta ==
ZERO);
4015 if (Y_is_replicated && this->
getComm ()->getRank () > 0) {
4021 X = rcp (
new MV (X_in, Teuchos::Copy));
4023 X = rcpFromRef (X_in);
4027 if (importer != null) {
4035 if (exporter != null) {
4046 if (! exporter.is_null ()) {
4054 if (importer != null) {
4062 this->
template localMultiply<Scalar, Scalar> (*X, *
importMV_, mode,
4064 if (Y_is_overwritten) {
4081 MV Y (Y_in, Teuchos::Copy);
4082 this->
template localMultiply<Scalar, Scalar> (*X, Y, mode, alpha, beta);
4085 this->
template localMultiply<Scalar, Scalar> (*X, Y_in, mode, alpha, beta);
4092 if (Y_is_replicated) {
4097 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4102 Teuchos::ETransp mode,
4106 TEUCHOS_TEST_FOR_EXCEPTION(
4108 "Tpetra::CrsMatrix::apply(): Cannot call apply() until fillComplete() " 4109 "has been called.");
4111 if (mode == Teuchos::NO_TRANS) {
4119 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4132 const Scalar& dampingFactor,
4134 const int numSweeps)
const 4139 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4145 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
4146 const Scalar& dampingFactor,
4148 const int numSweeps)
const 4150 using Teuchos::null;
4153 using Teuchos::rcp_const_cast;
4154 using Teuchos::rcpFromRef;
4157 TEUCHOS_TEST_FOR_EXCEPTION(
4159 "Tpetra::CrsMatrix::gaussSeidel: cannot call this method until " 4160 "fillComplete() has been called.");
4161 TEUCHOS_TEST_FOR_EXCEPTION(
4163 std::invalid_argument,
4164 "Tpetra::CrsMatrix::gaussSeidel: The number of sweeps must be , " 4165 "nonnegative but you provided numSweeps = " << numSweeps <<
" < 0.");
4169 KokkosClassic::ESweepDirection localDirection;
4170 if (direction == Forward) {
4171 localDirection = KokkosClassic::Forward;
4173 else if (direction == Backward) {
4174 localDirection = KokkosClassic::Backward;
4176 else if (direction == Symmetric) {
4178 localDirection = KokkosClassic::Forward;
4181 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
4182 "Tpetra::CrsMatrix::gaussSeidel: The 'direction' enum does not have " 4183 "any of its valid values: Forward, Backward, or Symmetric.");
4186 if (numSweeps == 0) {
4193 RCP<const import_type> importer = this->
getGraph()->getImporter();
4194 RCP<const export_type> exporter = this->
getGraph()->getExporter();
4195 TEUCHOS_TEST_FOR_EXCEPTION(
4196 ! exporter.is_null (), std::runtime_error,
4197 "Tpetra's gaussSeidel implementation requires that the row, domain, " 4198 "and range Maps be the same. This cannot be the case, because the " 4199 "matrix has a nontrivial Export object.");
4202 RCP<const map_type> rangeMap = this->
getRangeMap ();
4203 RCP<const map_type> rowMap = this->
getGraph ()->getRowMap ();
4204 RCP<const map_type> colMap = this->
getGraph ()->getColMap ();
4206 #ifdef HAVE_TEUCHOS_DEBUG 4211 TEUCHOS_TEST_FOR_EXCEPTION(
4212 ! X.
getMap ()->isSameAs (*domainMap),
4214 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4215 "multivector X be in the domain Map of the matrix.");
4216 TEUCHOS_TEST_FOR_EXCEPTION(
4217 ! B.
getMap ()->isSameAs (*rangeMap),
4219 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4220 "B be in the range Map of the matrix.");
4221 TEUCHOS_TEST_FOR_EXCEPTION(
4222 ! D.
getMap ()->isSameAs (*rowMap),
4224 "Tpetra::CrsMatrix::gaussSeidel requires that the input " 4225 "D be in the row Map of the matrix.");
4226 TEUCHOS_TEST_FOR_EXCEPTION(
4227 ! rowMap->isSameAs (*rangeMap),
4229 "Tpetra::CrsMatrix::gaussSeidel requires that the row Map and the " 4230 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
4231 TEUCHOS_TEST_FOR_EXCEPTION(
4232 ! domainMap->isSameAs (*rangeMap),
4234 "Tpetra::CrsMatrix::gaussSeidel requires that the domain Map and " 4235 "the range Map of the matrix be the same.");
4241 #endif // HAVE_TEUCHOS_DEBUG 4251 B_in = rcpFromRef (B);
4260 B_in = rcp_const_cast<
const MV> (B_in_nonconst);
4265 "gaussSeidel: The current implementation of the Gauss-Seidel kernel " 4266 "requires that X and B both have constant stride. Since B does not " 4267 "have constant stride, we had to make a copy. This is a limitation of " 4268 "the current implementation and not your fault, but we still report it " 4269 "as an efficiency warning for your information.");
4278 RCP<MV> X_domainMap;
4280 bool copiedInput =
false;
4282 if (importer.is_null ()) {
4284 X_domainMap = rcpFromRef (X);
4285 X_colMap = X_domainMap;
4286 copiedInput =
false;
4293 X_domainMap = X_colMap;
4298 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the " 4299 "Gauss-Seidel kernel requires that X and B both have constant " 4300 "stride. Since X does not have constant stride, we had to make a " 4301 "copy. This is a limitation of the current implementation and not " 4302 "your fault, but we still report it as an efficiency warning for " 4303 "your information.");
4308 X_domainMap = rcpFromRef (X);
4312 X_colMap = X_domainMap->offsetViewNonConst (colMap, 0);
4322 X_colMap->doImport (X, *importer,
INSERT);
4323 copiedInput =
false;
4331 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
4332 X_colMap->doImport (X, *importer,
INSERT);
4336 "Tpetra::CrsMatrix::gaussSeidel: The current implementation of the " 4337 "Gauss-Seidel kernel requires that X and B both have constant stride. " 4338 "Since X does not have constant stride, we had to make a copy. " 4339 "This is a limitation of the current implementation and not your fault, " 4340 "but we still report it as an efficiency warning for your information.");
4344 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
4345 if (! importer.is_null () && sweep > 0) {
4347 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4351 if (direction != Symmetric) {
4352 if (rowIndices.is_null ()) {
4353 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4358 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4365 const bool doImportBetweenDirections =
false;
4366 if (rowIndices.is_null ()) {
4367 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4369 KokkosClassic::Forward);
4374 if (doImportBetweenDirections) {
4376 if (! importer.is_null ()) {
4377 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4380 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4382 KokkosClassic::Backward);
4385 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4388 KokkosClassic::Forward);
4389 if (doImportBetweenDirections) {
4391 if (! importer.is_null ()) {
4392 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4395 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4398 KokkosClassic::Backward);
4408 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4414 const Scalar& dampingFactor,
4416 const int numSweeps,
4417 const bool zeroInitialGuess)
const 4420 numSweeps, zeroInitialGuess);
4423 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4429 const Teuchos::ArrayView<LocalOrdinal>& rowIndices,
4430 const Scalar& dampingFactor,
4432 const int numSweeps,
4433 const bool zeroInitialGuess)
const 4435 using Teuchos::null;
4438 using Teuchos::rcpFromRef;
4439 using Teuchos::rcp_const_cast;
4441 const char prefix[] =
"Tpetra::CrsMatrix::(reordered)gaussSeidelCopy: ";
4442 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
4444 TEUCHOS_TEST_FOR_EXCEPTION(
4446 prefix <<
"The matrix is not fill complete.");
4447 TEUCHOS_TEST_FOR_EXCEPTION(
4448 numSweeps < 0, std::invalid_argument,
4449 prefix <<
"The number of sweeps must be nonnegative, " 4450 "but you provided numSweeps = " << numSweeps <<
" < 0.");
4454 KokkosClassic::ESweepDirection localDirection;
4455 if (direction == Forward) {
4456 localDirection = KokkosClassic::Forward;
4458 else if (direction == Backward) {
4459 localDirection = KokkosClassic::Backward;
4461 else if (direction == Symmetric) {
4463 localDirection = KokkosClassic::Forward;
4466 TEUCHOS_TEST_FOR_EXCEPTION(
4467 true, std::invalid_argument,
4468 prefix <<
"The 'direction' enum does not have any of its valid " 4469 "values: Forward, Backward, or Symmetric.");
4472 if (numSweeps == 0) {
4476 RCP<const import_type> importer = this->
getGraph ()->getImporter ();
4477 RCP<const export_type> exporter = this->
getGraph ()->getExporter ();
4478 TEUCHOS_TEST_FOR_EXCEPTION(
4479 ! exporter.is_null (), std::runtime_error,
4480 "This method's implementation currently requires that the matrix's row, " 4481 "domain, and range Maps be the same. This cannot be the case, because " 4482 "the matrix has a nontrivial Export object.");
4485 RCP<const map_type> rangeMap = this->
getRangeMap ();
4486 RCP<const map_type> rowMap = this->
getGraph ()->getRowMap ();
4487 RCP<const map_type> colMap = this->
getGraph ()->getColMap ();
4489 #ifdef HAVE_TEUCHOS_DEBUG 4494 TEUCHOS_TEST_FOR_EXCEPTION(
4495 ! X.
getMap ()->isSameAs (*domainMap), std::runtime_error,
4496 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4497 "multivector X be in the domain Map of the matrix.");
4498 TEUCHOS_TEST_FOR_EXCEPTION(
4499 ! B.
getMap ()->isSameAs (*rangeMap), std::runtime_error,
4500 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4501 "B be in the range Map of the matrix.");
4502 TEUCHOS_TEST_FOR_EXCEPTION(
4503 ! D.
getMap ()->isSameAs (*rowMap), std::runtime_error,
4504 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the input " 4505 "D be in the row Map of the matrix.");
4506 TEUCHOS_TEST_FOR_EXCEPTION(
4507 ! rowMap->isSameAs (*rangeMap), std::runtime_error,
4508 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the row Map and the " 4509 "range Map be the same (in the sense of Tpetra::Map::isSameAs).");
4510 TEUCHOS_TEST_FOR_EXCEPTION(
4511 ! domainMap->isSameAs (*rangeMap), std::runtime_error,
4512 "Tpetra::CrsMatrix::gaussSeidelCopy requires that the domain Map and " 4513 "the range Map of the matrix be the same.");
4519 #endif // HAVE_TEUCHOS_DEBUG 4530 RCP<MV> X_domainMap;
4531 bool copyBackOutput =
false;
4532 if (importer.is_null ()) {
4534 X_colMap = rcpFromRef (X);
4535 X_domainMap = rcpFromRef (X);
4541 if (zeroInitialGuess) {
4542 X_colMap->putScalar (ZERO);
4554 X_domainMap = X_colMap;
4555 if (! zeroInitialGuess) {
4558 }
catch (std::exception& e) {
4559 std::ostringstream os;
4560 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: " 4561 "deep_copy(*X_domainMap, X) threw an exception: " 4562 << e.what () <<
".";
4563 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
4566 copyBackOutput =
true;
4570 "gaussSeidelCopy: The current implementation of the Gauss-Seidel " 4571 "kernel requires that X and B both have constant stride. Since X " 4572 "does not have constant stride, we had to make a copy. This is a " 4573 "limitation of the current implementation and not your fault, but we " 4574 "still report it as an efficiency warning for your information.");
4579 X_domainMap = X_colMap->offsetViewNonConst (domainMap, 0);
4581 #ifdef HAVE_TPETRA_DEBUG 4585 if (X_colMap->getLocalLength () != 0 && X_domainMap->getLocalLength ()) {
4586 TEUCHOS_TEST_FOR_EXCEPTION(
4587 X_colMap_view.h_view.ptr_on_device () != X_domainMap_view.h_view.ptr_on_device (),
4588 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: " 4589 "Pointer to start of column Map view of X is not equal to pointer to " 4590 "start of (domain Map view of) X. This may mean that " 4591 "Tpetra::MultiVector::offsetViewNonConst is broken. " 4592 "Please report this bug to the Tpetra developers.");
4595 TEUCHOS_TEST_FOR_EXCEPTION(
4596 X_colMap_view.dimension_0 () < X_domainMap_view.dimension_0 () ||
4597 X_colMap->getLocalLength () < X_domainMap->getLocalLength (),
4598 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: " 4599 "X_colMap has fewer local rows than X_domainMap. " 4600 "X_colMap_view.dimension_0() = " << X_colMap_view.dimension_0 ()
4601 <<
", X_domainMap_view.dimension_0() = " 4602 << X_domainMap_view.dimension_0 ()
4603 <<
", X_colMap->getLocalLength() = " << X_colMap->getLocalLength ()
4604 <<
", and X_domainMap->getLocalLength() = " 4605 << X_domainMap->getLocalLength ()
4606 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 4607 "is broken. Please report this bug to the Tpetra developers.");
4609 TEUCHOS_TEST_FOR_EXCEPTION(
4610 X_colMap->getNumVectors () != X_domainMap->getNumVectors (),
4611 std::logic_error,
"Tpetra::CrsMatrix::gaussSeidelCopy: " 4612 "X_colMap has a different number of columns than X_domainMap. " 4613 "X_colMap->getNumVectors() = " << X_colMap->getNumVectors ()
4614 <<
" != X_domainMap->getNumVectors() = " 4615 << X_domainMap->getNumVectors ()
4616 <<
". This means that Tpetra::MultiVector::offsetViewNonConst " 4617 "is broken. Please report this bug to the Tpetra developers.");
4618 #endif // HAVE_TPETRA_DEBUG 4620 if (zeroInitialGuess) {
4622 X_colMap->putScalar (ZERO);
4632 X_colMap->doImport (X, *importer,
INSERT);
4634 copyBackOutput =
true;
4642 B_in = rcpFromRef (B);
4651 }
catch (std::exception& e) {
4652 std::ostringstream os;
4653 os <<
"Tpetra::CrsMatrix::reorderedGaussSeidelCopy: " 4654 "deep_copy(*B_in_nonconst, B) threw an exception: " 4655 << e.what () <<
".";
4656 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error, e.what ());
4658 B_in = rcp_const_cast<
const MV> (B_in_nonconst);
4663 "gaussSeidelCopy: The current implementation requires that B have " 4664 "constant stride. Since B does not have constant stride, we had to " 4665 "copy it into a separate constant-stride multivector. This is a " 4666 "limitation of the current implementation and not your fault, but we " 4667 "still report it as an efficiency warning for your information.");
4670 for (
int sweep = 0; sweep < numSweeps; ++sweep) {
4671 if (! importer.is_null () && sweep > 0) {
4674 X_colMap->doImport (*X_domainMap, *importer,
INSERT);
4678 if (direction != Symmetric) {
4679 if (rowIndices.is_null ()) {
4680 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4685 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4692 if (rowIndices.is_null ()) {
4693 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4695 KokkosClassic::Forward);
4701 this->
template localGaussSeidel<ST, ST> (*B_in, *X_colMap, D,
4703 KokkosClassic::Backward);
4706 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4709 KokkosClassic::Forward);
4710 this->
template reorderedLocalGaussSeidel<ST, ST> (*B_in, *X_colMap,
4713 KokkosClassic::Backward);
4719 if (copyBackOutput) {
4722 }
catch (std::exception& e) {
4723 TEUCHOS_TEST_FOR_EXCEPTION(
4724 true, std::runtime_error, prefix <<
"deep_copy(X, *X_domainMap) " 4725 "threw an exception: " << e.what ());
4730 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4732 Teuchos::RCP<CrsMatrix<T, LocalOrdinal, GlobalOrdinal, Node, classic> >
4736 using Teuchos::ArrayRCP;
4740 typedef typename out_mat_type::local_matrix_type out_lcl_mat_type;
4741 typedef typename out_lcl_mat_type::values_type out_vals_type;
4742 typedef ArrayRCP<size_t>::size_type size_type;
4743 const char tfecfFuncName[] =
"convert";
4745 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4746 !
isFillComplete (), std::runtime_error,
"This matrix (the source of " 4747 "the conversion) is not fill complete. You must first call " 4748 "fillComplete() (possibly with the domain and range Map) without an " 4749 "intervening call to resumeFill(), before you may call this method.");
4756 RCP<out_mat_type> newmat;
4760 newmat = rcp (
new out_mat_type (this->
getCrsGraph ()));
4764 const size_type numVals =
4765 static_cast<size_type
> (this->
lclMatrix_.values.dimension_0 ());
4771 out_vals_type newVals1D (
"Tpetra::CrsMatrix::val", numVals);
4772 for (size_type k = 0; k < numVals; ++k) {
4773 newVals1D(k) =
static_cast<T
> (this->k_values1D_(k));
4775 newmat->lclMatrix_ =
4776 out_lcl_mat_type (
"Tpetra::CrsMatrix::lclMatrix_",
4779 newmat->k_values1D_ = newVals1D;
4799 ArrayRCP<const size_t> ptr;
4800 ArrayRCP<const LocalOrdinal> ind;
4801 ArrayRCP<const Scalar> oldVal;
4802 this->getAllValues (ptr, ind, oldVal);
4804 RCP<const map_type> rowMap = this->
getRowMap ();
4805 RCP<const map_type> colMap = this->
getColMap ();
4809 const size_type numLocalRows =
4810 static_cast<size_type
> (rowMap->getNodeNumElements ());
4811 ArrayRCP<size_t> numEntriesPerRow (numLocalRows);
4812 for (size_type localRow = 0; localRow < numLocalRows; ++localRow) {
4813 numEntriesPerRow[localRow] =
4817 newmat = rcp (
new out_mat_type (rowMap, colMap, numEntriesPerRow,
4821 const size_type numVals = this->
lclMatrix_.values.dimension_0 ();
4822 ArrayRCP<T> newVals1D (numVals);
4824 for (size_type k = 0; k < numVals; ++k) {
4825 newVals1D[k] =
static_cast<T
> (this->k_values1D_(k));
4832 ArrayRCP<size_t> newPtr (ptr.size ());
4833 std::copy (ptr.begin (), ptr.end (), newPtr.begin ());
4834 ArrayRCP<LocalOrdinal> newInd (ind.size ());
4835 std::copy (ind.begin (), ind.end (), newInd.begin ());
4836 newmat->setAllValues (newPtr, newInd, newVals1D);
4842 RCP<const map_type> rangeMap = this->
getRangeMap ();
4843 RCP<const import_type> importer = this->
getCrsGraph ()->getImporter ();
4844 RCP<const export_type> exporter = this->
getCrsGraph ()->getExporter ();
4845 newmat->expertStaticFillComplete (domainMap, rangeMap, importer, exporter);
4852 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4857 #ifdef HAVE_TPETRA_DEBUG 4858 const char tfecfFuncName[] =
"checkInternalState: ";
4859 const char err[] =
"Internal state is not consistent. " 4860 "Please report this bug to the Tpetra developers.";
4864 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4865 staticGraph_.is_null (),
4866 std::logic_error, err);
4870 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4871 ! myGraph_.is_null () && myGraph_ != staticGraph_,
4872 std::logic_error, err);
4874 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4876 std::logic_error, err <<
" Specifically, the matrix is fill complete, " 4877 "but its graph is NOT fill complete.");
4879 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4881 std::logic_error, err);
4883 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4885 std::logic_error, err);
4887 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4889 std::logic_error, err);
4892 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4893 staticGraph_->indicesAreAllocated () &&
4894 staticGraph_->getNodeAllocationSize() > 0 &&
4895 staticGraph_->getNodeNumRows() > 0
4896 && values2D_.is_null () &&
4897 k_values1D_.dimension_0 () == 0,
4898 std::logic_error, err);
4900 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
4901 k_values1D_.dimension_0 () > 0 && values2D_ != null,
4902 std::logic_error, err <<
" Specifically, k_values1D_ is allocated (has " 4903 "size " << k_values1D_.dimension_0 () <<
" > 0) and values2D_ is also " 4904 "allocated. CrsMatrix is not suppose to have both a 1-D and a 2-D " 4905 "allocation at the same time.");
4909 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4914 std::ostringstream os;
4916 os <<
"Tpetra::CrsMatrix (Kokkos refactor): {";
4917 if (this->getObjectLabel () !=
"") {
4918 os <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
4921 os <<
"isFillComplete: true" 4928 os <<
"isFillComplete: false" 4935 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
4939 const Teuchos::EVerbosityLevel verbLevel)
const 4943 using Teuchos::Comm;
4945 using Teuchos::TypeNameTraits;
4946 using Teuchos::VERB_DEFAULT;
4947 using Teuchos::VERB_NONE;
4948 using Teuchos::VERB_LOW;
4949 using Teuchos::VERB_MEDIUM;
4950 using Teuchos::VERB_HIGH;
4951 using Teuchos::VERB_EXTREME;
4953 const Teuchos::EVerbosityLevel vl = (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
4955 if (vl == VERB_NONE) {
4959 Teuchos::OSTab tab0 (out);
4961 RCP<const Comm<int> > comm = this->
getComm();
4962 const int myRank = comm->getRank();
4963 const int numProcs = comm->getSize();
4968 width = std::max<size_t> (width,
static_cast<size_t> (11)) + 2;
4978 out <<
"Tpetra::CrsMatrix (Kokkos refactor):" << endl;
4980 Teuchos::OSTab tab1 (out);
4983 if (this->getObjectLabel () !=
"") {
4984 out <<
"Label: \"" << this->getObjectLabel () <<
"\", ";
4987 out <<
"Template parameters:" << endl;
4988 Teuchos::OSTab tab2 (out);
4989 out <<
"Scalar: " << TypeNameTraits<Scalar>::name () << endl
4990 <<
"LocalOrdinal: " << TypeNameTraits<LocalOrdinal>::name () << endl
4991 <<
"GlobalOrdinal: " << TypeNameTraits<GlobalOrdinal>::name () << endl
4992 <<
"Node: " << TypeNameTraits<Node>::name () << endl;
4995 out <<
"isFillComplete: true" << endl
5000 << endl <<
"Global max number of entries in a row: " 5004 out <<
"isFillComplete: false" << endl
5010 if (vl < VERB_MEDIUM) {
5016 out << endl <<
"Row Map:" << endl;
5020 out <<
"null" << endl;
5032 out <<
"Column Map: ";
5036 out <<
"null" << endl;
5040 out <<
"same as row Map" << endl;
5051 out <<
"Domain Map: ";
5055 out <<
"null" << endl;
5059 out <<
"same as row Map" << endl;
5063 out <<
"same as column Map" << endl;
5074 out <<
"Range Map: ";
5078 out <<
"null" << endl;
5082 out <<
"same as domain Map" << endl;
5086 out <<
"same as row Map" << endl;
5096 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5097 if (myRank == curRank) {
5098 out <<
"Process rank: " << curRank << endl;
5099 Teuchos::OSTab tab2 (out);
5100 if (! staticGraph_->indicesAreAllocated ()) {
5101 out <<
"Graph indices not allocated" << endl;
5104 out <<
"Number of allocated entries: " 5105 << staticGraph_->getNodeAllocationSize () << endl;
5120 if (vl < VERB_HIGH) {
5125 for (
int curRank = 0; curRank < numProcs; ++curRank) {
5126 if (myRank == curRank) {
5127 out << std::setw(width) <<
"Proc Rank" 5128 << std::setw(width) <<
"Global Row" 5129 << std::setw(width) <<
"Num Entries";
5130 if (vl == VERB_EXTREME) {
5131 out << std::setw(width) <<
"(Index,Value)";
5136 GlobalOrdinal gid =
getRowMap()->getGlobalElement(r);
5137 out << std::setw(width) << myRank
5138 << std::setw(width) << gid
5139 << std::setw(width) << nE;
5140 if (vl == VERB_EXTREME) {
5142 ArrayView<const GlobalOrdinal> rowinds;
5143 ArrayView<const Scalar> rowvals;
5145 for (
size_t j = 0; j < nE; ++j) {
5146 out <<
" (" << rowinds[j]
5147 <<
", " << rowvals[j]
5152 ArrayView<const LocalOrdinal> rowinds;
5153 ArrayView<const Scalar> rowvals;
5155 for (
size_t j=0; j < nE; ++j) {
5156 out <<
" (" <<
getColMap()->getGlobalElement(rowinds[j])
5157 <<
", " << rowvals[j]
5173 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5186 const row_matrix_type* srcRowMat =
5187 dynamic_cast<const row_matrix_type*
> (&source);
5188 return (srcRowMat != NULL);
5191 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5196 const Teuchos::ArrayView<const LocalOrdinal>& permuteToLIDs,
5197 const Teuchos::ArrayView<const LocalOrdinal>& permuteFromLIDs)
5199 using Teuchos::Array;
5200 using Teuchos::ArrayView;
5201 typedef LocalOrdinal LO;
5202 typedef GlobalOrdinal GO;
5205 const char tfecfFuncName[] =
"copyAndPermute: ";
5207 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5208 permuteToLIDs.size() != permuteFromLIDs.size(),
5209 std::invalid_argument,
"permuteToLIDs.size() = " << permuteToLIDs.size()
5210 <<
"!= permuteFromLIDs.size() = " << permuteFromLIDs.size() <<
".");
5215 const row_matrix_type& srcMat =
dynamic_cast<const row_matrix_type&
> (source);
5222 const map_type& srcRowMap = * (srcMat.getRowMap ());
5224 Array<Scalar> rowVals;
5225 const LO numSameIDs_as_LID =
static_cast<LO
> (numSameIDs);
5226 for (LO sourceLID = 0; sourceLID < numSameIDs_as_LID; ++sourceLID) {
5231 const GO targetGID = sourceGID;
5234 ArrayView<const GO> rowIndsConstView;
5235 ArrayView<const Scalar> rowValsConstView;
5237 if (sourceIsLocallyIndexed) {
5238 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5239 if (rowLength > static_cast<size_t> (rowInds.size())) {
5240 rowInds.resize (rowLength);
5241 rowVals.resize (rowLength);
5245 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
5246 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
5251 size_t checkRowLength = 0;
5252 srcMat.getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
5254 #ifdef HAVE_TPETRA_DEBUG 5255 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
5256 std::logic_error,
"For global row index " << sourceGID <<
", the source" 5257 " matrix's getNumEntriesInGlobalRow() method returns a row length of " 5258 << rowLength <<
", but the getGlobalRowCopy() method reports that " 5259 "the row length is " << checkRowLength <<
". Please report this bug " 5260 "to the Tpetra developers.");
5261 #endif // HAVE_TPETRA_DEBUG 5263 rowIndsConstView = rowIndsView.view (0, rowLength);
5264 rowValsConstView = rowValsView.view (0, rowLength);
5267 srcMat.getGlobalRowView (sourceGID, rowIndsConstView, rowValsConstView);
5274 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
REPLACE);
5280 combineGlobalValues (targetGID, rowIndsConstView, rowValsConstView,
INSERT);
5288 const size_t numPermuteToLIDs =
static_cast<size_t> (permuteToLIDs.size ());
5289 for (
size_t p = 0; p < numPermuteToLIDs; ++p) {
5294 ArrayView<const GO> rowIndsConstView;
5295 ArrayView<const Scalar> rowValsConstView;
5297 if (sourceIsLocallyIndexed) {
5298 const size_t rowLength = srcMat.getNumEntriesInGlobalRow (sourceGID);
5299 if (rowLength > static_cast<size_t> (rowInds.size ())) {
5300 rowInds.resize (rowLength);
5301 rowVals.resize (rowLength);
5305 ArrayView<GO> rowIndsView = rowInds.view (0, rowLength);
5306 ArrayView<Scalar> rowValsView = rowVals.view (0, rowLength);
5311 size_t checkRowLength = 0;
5312 srcMat.getGlobalRowCopy (sourceGID, rowIndsView, rowValsView, checkRowLength);
5314 #ifdef HAVE_TPETRA_DEBUG 5315 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(rowLength != checkRowLength,
5316 std::logic_error,
"For the source matrix's global row index " 5317 << sourceGID <<
", the source matrix's getNumEntriesInGlobalRow() " 5318 "method returns a row length of " << rowLength <<
", but the " 5319 "getGlobalRowCopy() method reports that the row length is " 5320 << checkRowLength <<
". Please report this bug to the Tpetra " 5322 #endif // HAVE_TPETRA_DEBUG 5324 rowIndsConstView = rowIndsView.view (0, rowLength);
5325 rowValsConstView = rowValsView.view (0, rowLength);
5328 srcMat.getGlobalRowView (sourceGID, rowIndsConstView, rowValsConstView);
5333 this->combineGlobalValues (targetGID, rowIndsConstView,
5337 this->combineGlobalValues (targetGID, rowIndsConstView,
5338 rowValsConstView,
INSERT);
5343 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5347 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
5348 Teuchos::Array<char>& exports,
5349 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5350 size_t& constantNumPackets,
5353 using Teuchos::Array;
5354 using Teuchos::ArrayView;
5355 using Teuchos::av_reinterpret_cast;
5356 typedef LocalOrdinal LO;
5357 typedef GlobalOrdinal GO;
5358 const char tfecfFuncName[] =
"packAndPrepare: ";
5381 const row_matrix_type* srcRowMat =
5382 dynamic_cast<const row_matrix_type*
> (&source);
5383 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5384 srcRowMat == NULL, std::invalid_argument,
5385 "The source object of the Import or Export operation is neither a " 5386 "CrsMatrix (with the same template parameters as the target object), " 5387 "nor a RowMatrix (with the same first four template parameters as the " 5389 #ifdef HAVE_TPETRA_DEBUG 5391 using Teuchos::reduceAll;
5392 std::ostringstream msg;
5395 srcRowMat->pack (exportLIDs, exports, numPacketsPerLID,
5396 constantNumPackets, distor);
5397 }
catch (std::exception& e) {
5402 const Teuchos::Comm<int>& comm = * (this->
getComm ());
5403 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
5404 lclBad, Teuchos::outArg (gblBad));
5406 const int myRank = comm.getRank ();
5407 const int numProcs = comm.getSize ();
5408 for (
int r = 0; r < numProcs; ++r) {
5409 if (r == myRank && lclBad != 0) {
5410 std::ostringstream os;
5411 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
5412 std::cerr << os.str ();
5418 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5419 true, std::logic_error,
"pack() threw an exception on one or " 5420 "more participating processes.");
5424 srcRowMat->pack (exportLIDs, exports, numPacketsPerLID,
5425 constantNumPackets, distor);
5426 #endif // HAVE_TPETRA_DEBUG 5429 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5432 packRow (
char*
const numEntOut,
5435 const size_t numEnt,
5436 const LocalOrdinal lclRow)
const 5438 using Teuchos::ArrayView;
5439 typedef LocalOrdinal LO;
5440 typedef GlobalOrdinal GO;
5442 const LO numEntLO =
static_cast<LO
> (numEnt);
5443 memcpy (numEntOut, &numEntLO,
sizeof (LO));
5448 ArrayView<const LO> indIn;
5449 ArrayView<const Scalar> valIn;
5454 for (
size_t k = 0; k < numEnt; ++k) {
5456 memcpy (indOut + k *
sizeof (GO), &gblIndIn,
sizeof (GO));
5458 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
5466 ArrayView<const GO> indIn;
5467 ArrayView<const Scalar> valIn;
5471 memcpy (indOut, indIn.getRawPtr (), numEnt *
sizeof (GO));
5472 memcpy (valOut, valIn.getRawPtr (), numEnt *
sizeof (Scalar));
5483 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5487 GlobalOrdinal*
const indInTmp,
5488 const size_t tmpSize,
5489 const char*
const valIn,
5490 const char*
const indIn,
5491 const size_t numEnt,
5492 const LocalOrdinal lclRow,
5495 if (tmpSize < numEnt || (numEnt != 0 && (valInTmp == NULL || indInTmp == NULL))) {
5498 memcpy (valInTmp, valIn, numEnt *
sizeof (Scalar));
5499 memcpy (indInTmp, indIn, numEnt *
sizeof (GlobalOrdinal));
5500 const GlobalOrdinal gblRow = this->
getRowMap ()->getGlobalElement (lclRow);
5501 Teuchos::ArrayView<Scalar> val ((numEnt == 0) ? NULL : valInTmp, numEnt);
5502 Teuchos::ArrayView<GlobalOrdinal> ind ((numEnt == 0) ? NULL : indInTmp, numEnt);
5503 this->combineGlobalValues (gblRow, ind, val, combineMode);
5508 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5512 size_t& totalNumEntries,
5513 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs)
const 5515 typedef LocalOrdinal LO;
5516 typedef GlobalOrdinal GO;
5517 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
5519 const size_type numExportLIDs = exportLIDs.size ();
5522 totalNumEntries = 0;
5523 for (size_type i = 0; i < numExportLIDs; ++i) {
5524 const LO lclRow = exportLIDs[i];
5528 if (curNumEntries == Teuchos::OrdinalTraits<size_t>::invalid ()) {
5531 totalNumEntries += curNumEntries;
5542 const size_t allocSize =
5543 static_cast<size_t> (numExportLIDs) *
sizeof (LO) +
5544 totalNumEntries * (
sizeof (Scalar) +
sizeof (GO));
5545 if (static_cast<size_t> (exports.size ()) < allocSize) {
5546 exports.resize (allocSize);
5550 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5553 pack (
const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
5554 Teuchos::Array<char>& exports,
5555 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5556 size_t& constantNumPackets,
5559 using Teuchos::Array;
5560 using Teuchos::ArrayView;
5561 using Teuchos::av_reinterpret_cast;
5563 typedef LocalOrdinal LO;
5564 typedef GlobalOrdinal GO;
5565 typedef typename ArrayView<const LO>::size_type size_type;
5566 const char tfecfFuncName[] =
"pack: ";
5568 const size_type numExportLIDs = exportLIDs.size ();
5569 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5570 numExportLIDs != numPacketsPerLID.size (), std::invalid_argument,
5571 "exportLIDs.size() = " << numExportLIDs <<
" != numPacketsPerLID.size()" 5572 " = " << numPacketsPerLID.size () <<
".");
5577 constantNumPackets = 0;
5582 size_t totalNumEntries = 0;
5583 allocatePackSpace (exports, totalNumEntries, exportLIDs);
5584 const size_t bufSize =
static_cast<size_t> (exports.size ());
5596 size_type firstBadIndex = 0;
5597 size_t firstBadOffset = 0;
5598 size_t firstBadNumBytes = 0;
5599 bool outOfBounds =
false;
5600 bool packErr =
false;
5602 char*
const exportsRawPtr = exports.getRawPtr ();
5604 for (size_type i = 0; i < numExportLIDs; ++i) {
5605 const LO lclRow = exportLIDs[i];
5610 numPacketsPerLID[i] = 0;
5613 char*
const numEntBeg = exportsRawPtr + offset;
5614 char*
const numEntEnd = numEntBeg +
sizeof (LO);
5615 char*
const valBeg = numEntEnd;
5616 char*
const valEnd = valBeg + numEnt *
sizeof (Scalar);
5617 char*
const indBeg = valEnd;
5618 const size_t numBytes =
sizeof (LO) +
5619 numEnt * (
sizeof (Scalar) +
sizeof (GO));
5620 if (offset > bufSize || offset + numBytes > bufSize) {
5622 firstBadOffset = offset;
5623 firstBadNumBytes = numBytes;
5627 packErr = ! packRow (numEntBeg, valBeg, indBeg, numEnt, lclRow);
5630 firstBadOffset = offset;
5631 firstBadNumBytes = numBytes;
5637 numPacketsPerLID[i] = numBytes;
5642 TEUCHOS_TEST_FOR_EXCEPTION(
5643 outOfBounds, std::logic_error,
"First invalid offset into 'exports' " 5644 "pack buffer at index i = " << firstBadIndex <<
". exportLIDs[i]: " 5645 << exportLIDs[firstBadIndex] <<
", bufSize: " << bufSize <<
", offset: " 5646 << firstBadOffset <<
", numBytes: " << firstBadNumBytes <<
".");
5647 TEUCHOS_TEST_FOR_EXCEPTION(
5648 packErr, std::logic_error,
"First error in packRow() at index i = " 5649 << firstBadIndex <<
". exportLIDs[i]: " << exportLIDs[firstBadIndex]
5650 <<
", bufSize: " << bufSize <<
", offset: " << firstBadOffset
5651 <<
", numBytes: " << firstBadNumBytes <<
".");
5654 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5658 const Teuchos::ArrayView<const GlobalOrdinal>& columnIndices,
5659 const Teuchos::ArrayView<const Scalar>& values,
5662 const char tfecfFuncName[] =
"combineGlobalValues: ";
5668 if (combineMode ==
ADD) {
5671 else if (combineMode ==
REPLACE) {
5674 else if (combineMode ==
ABSMAX) {
5677 this->
template transformGlobalValues<AbsMax<Scalar> > (globalRowIndex,
5681 else if (combineMode ==
INSERT) {
5682 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5684 "INSERT combine mode is not allowed if the matrix has a static graph " 5685 "(i.e., was constructed with the CrsMatrix constructor that takes a " 5686 "const CrsGraph pointer).");
5689 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5690 true, std::logic_error,
"Invalid combine mode; should never get " 5691 "here! Please report this bug to the Tpetra developers.");
5695 if (combineMode ==
ADD || combineMode ==
INSERT) {
5702 insertGlobalValuesFiltered (globalRowIndex, columnIndices, values);
5713 else if (combineMode ==
ABSMAX) {
5714 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5716 "ABSMAX combine mode when the matrix has a dynamic graph is not yet " 5719 else if (combineMode ==
REPLACE) {
5720 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5722 "REPLACE combine mode when the matrix has a dynamic graph is not yet " 5726 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5727 true, std::logic_error,
"Should never get here! Please report this " 5728 "bug to the Tpetra developers.");
5734 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5738 const Teuchos::ArrayView<const char>& imports,
5739 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5740 size_t constantNumPackets,
5744 #ifdef HAVE_TPETRA_DEBUG 5745 const char tfecfFuncName[] =
"unpackAndCombine: ";
5747 const char* validModeNames[4] = {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT"};
5748 const int numValidModes = 4;
5750 if (std::find (validModes, validModes+numValidModes, combineMode) ==
5751 validModes+numValidModes) {
5752 std::ostringstream os;
5753 os <<
"Invalid combine mode. Valid modes are {";
5754 for (
int k = 0; k < numValidModes; ++k) {
5755 os << validModeNames[k];
5756 if (k < numValidModes - 1) {
5761 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5762 true, std::invalid_argument, os.str ());
5766 using Teuchos::reduceAll;
5767 std::ostringstream msg;
5770 this->unpackAndCombineImpl (importLIDs, imports, numPacketsPerLID,
5771 constantNumPackets, distor, combineMode);
5772 }
catch (std::exception& e) {
5777 const Teuchos::Comm<int>& comm = * (this->
getComm ());
5778 reduceAll<int, int> (comm, Teuchos::REDUCE_MAX,
5779 lclBad, Teuchos::outArg (gblBad));
5781 const int myRank = comm.getRank ();
5782 const int numProcs = comm.getSize ();
5783 for (
int r = 0; r < numProcs; ++r) {
5784 if (r == myRank && lclBad != 0) {
5785 std::ostringstream os;
5786 os <<
"Proc " << myRank <<
": " << msg.str () << std::endl;
5787 std::cerr << os.str ();
5793 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5794 true, std::logic_error,
"unpackAndCombineImpl() threw an " 5795 "exception on one or more participating processes.");
5799 this->unpackAndCombineImpl (importLIDs, imports, numPacketsPerLID,
5800 constantNumPackets, distor, combineMode);
5801 #endif // HAVE_TPETRA_DEBUG 5804 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5808 const Teuchos::ArrayView<const char>& imports,
5809 const Teuchos::ArrayView<size_t>& numPacketsPerLID,
5810 size_t constantNumPackets,
5814 typedef LocalOrdinal LO;
5815 typedef GlobalOrdinal GO;
5816 typedef typename Teuchos::ArrayView<const LO>::size_type size_type;
5817 const char tfecfFuncName[] =
"unpackAndCombine: ";
5819 #ifdef HAVE_TPETRA_DEBUG 5821 const char* validModeNames[4] = {
"ADD",
"REPLACE",
"ABSMAX",
"INSERT"};
5822 const int numValidModes = 4;
5824 if (std::find (validModes, validModes+numValidModes, combineMode) ==
5825 validModes+numValidModes) {
5826 std::ostringstream os;
5827 os <<
"Invalid combine mode. Valid modes are {";
5828 for (
int k = 0; k < numValidModes; ++k) {
5829 os << validModeNames[k];
5830 if (k < numValidModes - 1) {
5835 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5836 true, std::invalid_argument, os.str ());
5838 #endif // HAVE_TPETRA_DEBUG 5840 const size_type numImportLIDs = importLIDs.size ();
5841 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5842 numImportLIDs != numPacketsPerLID.size (), std::invalid_argument,
5843 "importLIDs.size() = " << numImportLIDs <<
" != numPacketsPerLID.size()" 5844 <<
" = " << numPacketsPerLID.size () <<
".");
5850 size_type firstBadIndex = 0;
5851 size_t firstBadOffset = 0;
5852 size_t firstBadExpectedNumBytes = 0;
5853 size_t firstBadNumBytes = 0;
5854 LO firstBadNumEnt = 0;
5863 bool outOfBounds =
false;
5864 bool wrongNumBytes =
false;
5865 bool unpackErr =
false;
5867 const size_t bufSize =
static_cast<size_t> (imports.size ());
5868 const char*
const importsRawPtr = imports.getRawPtr ();
5877 Array<Scalar> valInTmp;
5879 for (size_type i = 0; i < numImportLIDs; ++i) {
5880 const LO lclRow = importLIDs[i];
5881 const size_t numBytes = numPacketsPerLID[i];
5884 const char*
const numEntBeg = importsRawPtr + offset;
5885 const char*
const numEntEnd = numEntBeg +
sizeof (LO);
5890 memcpy (&numEnt, numEntBeg,
sizeof (LO));
5892 const char*
const valBeg = numEntEnd;
5893 const char*
const valEnd =
5894 valBeg +
static_cast<size_t> (numEnt) *
sizeof (Scalar);
5895 const char*
const indBeg = valEnd;
5896 const size_t expectedNumBytes =
sizeof (LO) +
5897 static_cast<size_t> (numEnt) * (
sizeof (Scalar) +
sizeof (GO));
5899 if (expectedNumBytes > numBytes) {
5901 firstBadOffset = offset;
5902 firstBadExpectedNumBytes = expectedNumBytes;
5903 firstBadNumBytes = numBytes;
5904 firstBadNumEnt = numEnt;
5905 wrongNumBytes =
true;
5908 if (offset > bufSize || offset + numBytes > bufSize) {
5910 firstBadOffset = offset;
5911 firstBadExpectedNumBytes = expectedNumBytes;
5912 firstBadNumBytes = numBytes;
5913 firstBadNumEnt = numEnt;
5917 size_t tmpNumEnt =
static_cast<size_t> (valInTmp.size ());
5918 if (tmpNumEnt < static_cast<size_t> (numEnt) ||
5919 static_cast<size_t> (indInTmp.size ()) < static_cast<size_t> (numEnt)) {
5921 tmpNumEnt = std::max (static_cast<size_t> (numEnt), tmpNumEnt * 2);
5922 valInTmp.resize (tmpNumEnt);
5923 indInTmp.resize (tmpNumEnt);
5926 ! unpackRow (valInTmp.getRawPtr (), indInTmp.getRawPtr (), tmpNumEnt,
5927 valBeg, indBeg, numEnt, lclRow, combineMode);
5930 firstBadOffset = offset;
5931 firstBadExpectedNumBytes = expectedNumBytes;
5932 firstBadNumBytes = numBytes;
5933 firstBadNumEnt = numEnt;
5940 if (wrongNumBytes || outOfBounds || unpackErr) {
5941 std::ostringstream os;
5942 os <<
" importLIDs[i]: " << importLIDs[firstBadIndex]
5943 <<
", bufSize: " << bufSize
5944 <<
", offset: " << firstBadOffset
5945 <<
", numBytes: " << firstBadNumBytes
5946 <<
", expectedNumBytes: " << firstBadExpectedNumBytes
5947 <<
", numEnt: " << firstBadNumEnt;
5948 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5949 wrongNumBytes, std::logic_error,
"At index i = " << firstBadIndex
5950 <<
", expectedNumBytes > numBytes." << os.str ());
5951 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5952 outOfBounds, std::logic_error,
"First invalid offset into 'imports' " 5953 "unpack buffer at index i = " << firstBadIndex <<
"." << os.str ());
5954 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
5955 unpackErr, std::logic_error,
"First error in unpackRow() at index i = " 5956 << firstBadIndex <<
"." << os.str ());
5960 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
5961 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
5964 const bool force)
const 5966 using Teuchos::null;
5970 TEUCHOS_TEST_FOR_EXCEPTION(
5971 ! this->
hasColMap (), std::runtime_error,
"Tpetra::CrsMatrix::getColumn" 5972 "MapMultiVector: You may only call this method if the matrix has a " 5973 "column Map. If the matrix does not yet have a column Map, you should " 5974 "first call fillComplete (with domain and range Map if necessary).");
5978 TEUCHOS_TEST_FOR_EXCEPTION(
5980 "CrsMatrix::getColumnMapMultiVector: You may only call this method if " 5981 "this matrix's graph is fill complete.");
5984 RCP<const import_type> importer = this->
getGraph ()->getImporter ();
5985 RCP<const map_type> colMap = this->
getColMap ();
5998 if (! importer.is_null () || force) {
6000 X_colMap = rcp (
new MV (colMap, numVecs));
6017 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6018 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
6021 const bool force)
const 6023 using Teuchos::null;
6029 TEUCHOS_TEST_FOR_EXCEPTION(
6031 "CrsMatrix::getRowMapMultiVector: You may only call this method if this " 6032 "matrix's graph is fill complete.");
6035 RCP<const export_type> exporter = this->
getGraph ()->getExporter ();
6039 RCP<const map_type> rowMap = this->
getRowMap ();
6051 if (! exporter.is_null () || force) {
6053 Y_rowMap = rcp (
new MV (rowMap, numVecs));
6063 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6068 TEUCHOS_TEST_FOR_EXCEPTION(
6069 myGraph_.is_null (), std::logic_error,
"Tpetra::CrsMatrix::" 6070 "removeEmptyProcessesInPlace: This method does not work when the matrix " 6071 "was created with a constant graph (that is, when it was created using " 6072 "the version of its constructor that takes an RCP<const CrsGraph>). " 6073 "This is because the matrix is not allowed to modify the graph in that " 6074 "case, but removing empty processes requires modifying the graph.");
6075 myGraph_->removeEmptyProcessesInPlace (newMap);
6083 staticGraph_ = Teuchos::rcp_const_cast<
const Graph> (myGraph_);
6086 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6087 Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
6092 const Teuchos::RCP<const map_type>& domainMap,
6093 const Teuchos::RCP<const map_type>& rangeMap,
6094 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6096 using Teuchos::Array;
6097 using Teuchos::ArrayRCP;
6098 using Teuchos::ParameterList;
6101 using Teuchos::rcp_implicit_cast;
6102 using Teuchos::sublist;
6103 typedef LocalOrdinal LO;
6104 typedef GlobalOrdinal GO;
6108 const crs_matrix_type& B = *
this;
6109 const Scalar
ZERO = Teuchos::ScalarTraits<Scalar>::zero ();
6110 const Scalar ONE = Teuchos::ScalarTraits<Scalar>::one ();
6117 RCP<const map_type> A_rangeMap = A.
getRangeMap ();
6118 RCP<const map_type> B_domainMap = B.getDomainMap ();
6119 RCP<const map_type> B_rangeMap = B.getRangeMap ();
6121 RCP<const map_type> theDomainMap = domainMap;
6122 RCP<const map_type> theRangeMap = rangeMap;
6124 if (domainMap.is_null ()) {
6125 if (B_domainMap.is_null ()) {
6126 TEUCHOS_TEST_FOR_EXCEPTION(
6127 A_domainMap.is_null (), std::invalid_argument,
6128 "Tpetra::CrsMatrix::add: If neither A nor B have a domain Map, " 6129 "then you must supply a nonnull domain Map to this method.");
6130 theDomainMap = A_domainMap;
6132 theDomainMap = B_domainMap;
6135 if (rangeMap.is_null ()) {
6136 if (B_rangeMap.is_null ()) {
6137 TEUCHOS_TEST_FOR_EXCEPTION(
6138 A_rangeMap.is_null (), std::invalid_argument,
6139 "Tpetra::CrsMatrix::add: If neither A nor B have a range Map, " 6140 "then you must supply a nonnull range Map to this method.");
6141 theRangeMap = A_rangeMap;
6143 theRangeMap = B_rangeMap;
6147 #ifdef HAVE_TPETRA_DEBUG 6151 if (! A_domainMap.is_null () && ! A_rangeMap.is_null ()) {
6152 if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
6153 TEUCHOS_TEST_FOR_EXCEPTION(
6154 ! B_domainMap->isSameAs (*A_domainMap), std::invalid_argument,
6155 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a domain Map " 6156 "which is the same as (isSameAs) this RowMatrix's domain Map.");
6157 TEUCHOS_TEST_FOR_EXCEPTION(
6158 ! B_rangeMap->isSameAs (*A_rangeMap), std::invalid_argument,
6159 "Tpetra::CrsMatrix::add: The input RowMatrix A must have a range Map " 6160 "which is the same as (isSameAs) this RowMatrix's range Map.");
6161 TEUCHOS_TEST_FOR_EXCEPTION(
6162 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
6163 std::invalid_argument,
6164 "Tpetra::CrsMatrix::add: The input domain Map must be the same as " 6165 "(isSameAs) this RowMatrix's domain Map.");
6166 TEUCHOS_TEST_FOR_EXCEPTION(
6167 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
6168 std::invalid_argument,
6169 "Tpetra::CrsMatrix::add: The input range Map must be the same as " 6170 "(isSameAs) this RowMatrix's range Map.");
6173 else if (! B_domainMap.is_null () && ! B_rangeMap.is_null ()) {
6174 TEUCHOS_TEST_FOR_EXCEPTION(
6175 ! domainMap.is_null () && ! domainMap->isSameAs (*B_domainMap),
6176 std::invalid_argument,
6177 "Tpetra::CrsMatrix::add: The input domain Map must be the same as " 6178 "(isSameAs) this RowMatrix's domain Map.");
6179 TEUCHOS_TEST_FOR_EXCEPTION(
6180 ! rangeMap.is_null () && ! rangeMap->isSameAs (*B_rangeMap),
6181 std::invalid_argument,
6182 "Tpetra::CrsMatrix::add: The input range Map must be the same as " 6183 "(isSameAs) this RowMatrix's range Map.");
6186 TEUCHOS_TEST_FOR_EXCEPTION(
6187 domainMap.is_null () || rangeMap.is_null (), std::invalid_argument,
6188 "Tpetra::CrsMatrix::add: If neither A nor B have a domain and range " 6189 "Map, then you must supply a nonnull domain and range Map to this " 6192 #endif // HAVE_TPETRA_DEBUG 6197 bool callFillComplete =
true;
6198 RCP<ParameterList> constructorSublist;
6199 RCP<ParameterList> fillCompleteSublist;
6200 if (! params.is_null ()) {
6201 callFillComplete = params->get (
"Call fillComplete", callFillComplete);
6202 constructorSublist = sublist (params,
"Constructor parameters");
6203 fillCompleteSublist = sublist (params,
"fillComplete parameters");
6206 RCP<const map_type> A_rowMap = A.
getRowMap ();
6207 RCP<const map_type> B_rowMap = B.getRowMap ();
6208 RCP<const map_type> C_rowMap = B_rowMap;
6209 RCP<crs_matrix_type> C;
6216 if (A_rowMap->isSameAs (*B_rowMap)) {
6217 const LO localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
6218 ArrayRCP<size_t> C_maxNumEntriesPerRow (localNumRows, 0);
6221 if (alpha != ZERO) {
6222 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
6224 C_maxNumEntriesPerRow[localRow] += A_numEntries;
6229 for (LO localRow = 0; localRow < localNumRows; ++localRow) {
6230 const size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
6231 C_maxNumEntriesPerRow[localRow] += B_numEntries;
6235 if (constructorSublist.is_null ()) {
6236 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
6239 C = rcp (
new crs_matrix_type (C_rowMap, C_maxNumEntriesPerRow,
6250 if (constructorSublist.is_null ()) {
6254 constructorSublist));
6258 #ifdef HAVE_TPETRA_DEBUG 6259 TEUCHOS_TEST_FOR_EXCEPTION(C.is_null (), std::logic_error,
6260 "Tpetra::RowMatrix::add: C should not be null at this point. " 6261 "Please report this bug to the Tpetra developers.");
6262 #endif // HAVE_TPETRA_DEBUG 6269 if (alpha != ZERO) {
6270 const LO A_localNumRows =
static_cast<LO
> (A_rowMap->getNodeNumElements ());
6271 for (LO localRow = 0; localRow < A_localNumRows; ++localRow) {
6273 const GO globalRow = A_rowMap->getGlobalElement (localRow);
6274 if (A_numEntries > static_cast<size_t> (ind.size ())) {
6275 ind.resize (A_numEntries);
6276 val.resize (A_numEntries);
6278 ArrayView<GO> indView = ind (0, A_numEntries);
6279 ArrayView<Scalar> valView = val (0, A_numEntries);
6283 for (
size_t k = 0; k < A_numEntries; ++k) {
6284 valView[k] *= alpha;
6287 C->insertGlobalValues (globalRow, indView, valView);
6292 const LO B_localNumRows =
static_cast<LO
> (B_rowMap->getNodeNumElements ());
6293 for (LO localRow = 0; localRow < B_localNumRows; ++localRow) {
6294 size_t B_numEntries = B.getNumEntriesInLocalRow (localRow);
6295 const GO globalRow = B_rowMap->getGlobalElement (localRow);
6296 if (B_numEntries > static_cast<size_t> (ind.size ())) {
6297 ind.resize (B_numEntries);
6298 val.resize (B_numEntries);
6300 ArrayView<GO> indView = ind (0, B_numEntries);
6301 ArrayView<Scalar> valView = val (0, B_numEntries);
6302 B.getGlobalRowCopy (globalRow, indView, valView, B_numEntries);
6305 for (
size_t k = 0; k < B_numEntries; ++k) {
6309 C->insertGlobalValues (globalRow, indView, valView);
6313 if (callFillComplete) {
6314 if (fillCompleteSublist.is_null ()) {
6315 C->fillComplete (theDomainMap, theRangeMap);
6317 C->fillComplete (theDomainMap, theRangeMap, fillCompleteSublist);
6320 return rcp_implicit_cast<row_matrix_type> (C);
6323 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6327 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer,
6328 const Teuchos::RCP<const map_type>& domainMap,
6329 const Teuchos::RCP<const map_type>& rangeMap,
6330 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6332 using Teuchos::ArrayView;
6333 using Teuchos::Comm;
6334 using Teuchos::ParameterList;
6336 typedef LocalOrdinal LO;
6337 typedef GlobalOrdinal GO;
6342 #ifdef HAVE_TPETRA_MMM_TIMINGS 6344 if(!params.is_null())
6345 label = params->get(
"Timer Label",label);
6346 std::string prefix = std::string(
"Tpetra ")+ label + std::string(
": ");
6347 using Teuchos::TimeMonitor;
6348 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-1"))));
6358 TEUCHOS_TEST_FOR_EXCEPTION(
6359 xferAsImport == NULL && xferAsExport == NULL, std::invalid_argument,
6360 "Tpetra::CrsMatrix::transferAndFillComplete: The 'rowTransfer' input " 6361 "argument must be either an Import or an Export, and its template " 6362 "parameters must match the corresponding template parameters of the " 6367 const bool communication_needed = rowTransfer.getSourceMap ()->isDistributed ();
6373 bool reverseMode =
false;
6374 bool restrictComm =
false;
6375 RCP<ParameterList> matrixparams;
6376 if (! params.is_null ()) {
6377 reverseMode = params->get (
"Reverse Mode", reverseMode);
6378 restrictComm = params->get (
"Restrict Communicator", restrictComm);
6379 matrixparams = sublist (params,
"CrsMatrix");
6384 RCP<const map_type> MyRowMap = reverseMode ?
6385 rowTransfer.getSourceMap () : rowTransfer.getTargetMap ();
6386 RCP<const map_type> MyColMap;
6387 RCP<const map_type> MyDomainMap = ! domainMap.is_null () ?
6389 RCP<const map_type> MyRangeMap = ! rangeMap.is_null () ?
6391 RCP<const map_type> BaseRowMap = MyRowMap;
6392 RCP<const map_type> BaseDomainMap = MyDomainMap;
6400 if (! destMat.is_null ()) {
6411 const bool NewFlag = ! destMat->getGraph ()->isLocallyIndexed () &&
6412 ! destMat->getGraph ()->isGloballyIndexed ();
6413 TEUCHOS_TEST_FOR_EXCEPTION(
6414 ! NewFlag, std::invalid_argument,
"Tpetra::CrsMatrix::" 6415 "transferAndFillComplete: The input argument 'destMat' is only allowed " 6416 "to be nonnull, if its graph is empty (neither locally nor globally " 6425 TEUCHOS_TEST_FOR_EXCEPTION(
6426 ! destMat->getRowMap ()->isSameAs (*MyRowMap), std::invalid_argument,
6427 "Tpetra::CrsMatrix::transferAndFillComplete: The (row) Map of the " 6428 "input argument 'destMat' is not the same as the (row) Map specified " 6429 "by the input argument 'rowTransfer'.");
6430 TEUCHOS_TEST_FOR_EXCEPTION(
6431 ! destMat->checkSizes (*
this), std::invalid_argument,
6432 "Tpetra::CrsMatrix::transferAndFillComplete: You provided a nonnull " 6433 "destination matrix, but checkSizes() indicates that it is not a legal " 6434 "legal target for redistribution from the source matrix (*this). This " 6435 "may mean that they do not have the same dimensions.");
6449 TEUCHOS_TEST_FOR_EXCEPTION(
6450 ! (reverseMode ||
getRowMap ()->isSameAs (*rowTransfer.getSourceMap ())),
6451 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: " 6452 "rowTransfer->getSourceMap() must match this->getRowMap() in forward mode.");
6453 TEUCHOS_TEST_FOR_EXCEPTION(
6454 ! (! reverseMode ||
getRowMap ()->isSameAs (*rowTransfer.getTargetMap ())),
6455 std::invalid_argument,
"Tpetra::CrsMatrix::transferAndFillComplete: " 6456 "rowTransfer->getTargetMap() must match this->getRowMap() in reverse mode.");
6469 const size_t NumSameIDs = rowTransfer.getNumSameIDs();
6470 ArrayView<const LO> ExportLIDs = reverseMode ?
6471 rowTransfer.getRemoteLIDs () : rowTransfer.getExportLIDs ();
6472 ArrayView<const LO> RemoteLIDs = reverseMode ?
6473 rowTransfer.getExportLIDs () : rowTransfer.getRemoteLIDs ();
6474 ArrayView<const LO> PermuteToLIDs = reverseMode ?
6475 rowTransfer.getPermuteFromLIDs () : rowTransfer.getPermuteToLIDs ();
6476 ArrayView<const LO> PermuteFromLIDs = reverseMode ?
6477 rowTransfer.getPermuteToLIDs () : rowTransfer.getPermuteFromLIDs ();
6478 Distributor& Distor = rowTransfer.getDistributor ();
6481 Teuchos::Array<int> SourcePids;
6482 Teuchos::Array<int> TargetPids;
6483 int MyPID =
getComm ()->getRank ();
6486 RCP<const map_type> ReducedRowMap, ReducedColMap,
6487 ReducedDomainMap, ReducedRangeMap;
6488 RCP<const Comm<int> > ReducedComm;
6492 if (destMat.is_null ()) {
6493 destMat = rcp (
new this_type (MyRowMap, 0,
StaticProfile, matrixparams));
6500 ReducedRowMap = MyRowMap->removeEmptyProcesses ();
6501 ReducedComm = ReducedRowMap.is_null () ?
6503 ReducedRowMap->getComm ();
6504 destMat->removeEmptyProcessesInPlace (ReducedRowMap);
6506 ReducedDomainMap = MyRowMap.getRawPtr () == MyDomainMap.getRawPtr () ?
6508 MyDomainMap->replaceCommWithSubset (ReducedComm);
6509 ReducedRangeMap = MyRowMap.getRawPtr () == MyRangeMap.getRawPtr () ?
6511 MyRangeMap->replaceCommWithSubset (ReducedComm);
6514 MyRowMap = ReducedRowMap;
6515 MyDomainMap = ReducedDomainMap;
6516 MyRangeMap = ReducedRangeMap;
6519 if (! ReducedComm.is_null ()) {
6520 MyPID = ReducedComm->getRank ();
6527 ReducedComm = MyRowMap->getComm ();
6533 #ifdef HAVE_TPETRA_MMM_TIMINGS 6534 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ImportSetup"))));
6537 RCP<const import_type> MyImporter =
getGraph ()->getImporter ();
6539 if (! restrictComm && ! MyImporter.is_null () &&
6547 Import_Util::getPids (*MyImporter, SourcePids,
false);
6549 else if (MyImporter.is_null () && BaseDomainMap->isSameAs (*
getDomainMap ())) {
6551 SourcePids.resize (
getColMap ()->getNodeNumElements ());
6552 SourcePids.assign (
getColMap ()->getNodeNumElements (), MyPID);
6554 else if (BaseDomainMap->isSameAs (*BaseRowMap) &&
6557 IntVectorType TargetRow_pids (domainMap);
6558 IntVectorType SourceRow_pids (
getRowMap ());
6559 IntVectorType SourceCol_pids (
getColMap ());
6561 TargetRow_pids.putScalar (MyPID);
6562 if (! reverseMode && xferAsImport != NULL) {
6563 SourceRow_pids.doExport (TargetRow_pids, *xferAsImport,
INSERT);
6565 else if (reverseMode && xferAsExport != NULL) {
6566 SourceRow_pids.doExport (TargetRow_pids, *xferAsExport,
INSERT);
6568 else if (! reverseMode && xferAsExport != NULL) {
6569 SourceRow_pids.doImport (TargetRow_pids, *xferAsExport,
INSERT);
6571 else if (reverseMode && xferAsImport != NULL) {
6572 SourceRow_pids.doImport (TargetRow_pids, *xferAsImport,
INSERT);
6575 TEUCHOS_TEST_FOR_EXCEPTION(
6576 true, std::logic_error,
"Tpetra::CrsMatrix::" 6577 "transferAndFillComplete: Should never get here! " 6578 "Please report this bug to a Tpetra developer.");
6580 SourceCol_pids.doImport (SourceRow_pids, *MyImporter,
INSERT);
6581 SourcePids.resize (
getColMap ()->getNodeNumElements ());
6582 SourceCol_pids.get1dCopy (SourcePids ());
6585 TEUCHOS_TEST_FOR_EXCEPTION(
6586 true, std::invalid_argument,
"Tpetra::CrsMatrix::" 6587 "transferAndFillComplete: This method only allows either domainMap == " 6588 "getDomainMap (), or (domainMap == rowTransfer.getTargetMap () and " 6589 "getDomainMap () == getRowMap ()).");
6591 #ifdef HAVE_TPETRA_MMM_TIMINGS 6592 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Pack-2"))));
6615 size_t constantNumPackets = destMat->constantNumberOfPackets ();
6616 if (constantNumPackets == 0) {
6617 destMat->numExportPacketsPerLID_old_.resize (ExportLIDs.size ());
6618 destMat->numImportPacketsPerLID_old_.resize (RemoteLIDs.size ());
6625 const size_t rbufLen = RemoteLIDs.size() * constantNumPackets;
6626 if (static_cast<size_t> (destMat->imports_old_.size ()) != rbufLen) {
6627 destMat->imports_old_.resize (rbufLen);
6644 #ifdef HAVE_TPETRA_DEBUG 6646 using Teuchos::outArg;
6647 using Teuchos::REDUCE_MAX;
6648 using Teuchos::reduceAll;
6651 RCP<const Teuchos::Comm<int> > comm = this->
getComm ();
6652 const int myRank = comm->getRank ();
6653 const int numProcs = comm->getSize ();
6655 std::ostringstream os;
6658 Import_Util::packAndPrepareWithOwningPIDs (*
this, ExportLIDs,
6659 destMat->exports_old_,
6660 destMat->numExportPacketsPerLID_old_ (),
6661 constantNumPackets, Distor,
6664 catch (std::exception& e) {
6665 os <<
"Proc " << myRank <<
": " << e.what ();
6669 if (! comm.is_null ()) {
6670 reduceAll<int, int> (*comm, REDUCE_MAX, lclErr, outArg (gblErr));
6674 cerr <<
"packAndPrepareWithOwningPIDs threw an exception: " << endl;
6676 std::ostringstream err;
6677 for (
int r = 0; r < numProcs; ++r) {
6678 if (r == myRank && lclErr != 0) {
6679 cerr << os.str () << endl;
6686 TEUCHOS_TEST_FOR_EXCEPTION(
6687 true, std::logic_error,
"packAndPrepareWithOwningPIDs threw an " 6693 Import_Util::packAndPrepareWithOwningPIDs (*
this, ExportLIDs,
6694 destMat->exports_old_,
6695 destMat->numExportPacketsPerLID_old_ (),
6696 constantNumPackets, Distor,
6698 #endif // HAVE_TPETRA_DEBUG 6712 #ifdef HAVE_TPETRA_MMM_TIMINGS 6713 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Transfer"))));
6716 if (communication_needed) {
6718 if (constantNumPackets == 0) {
6719 Distor.doReversePostsAndWaits (destMat->numExportPacketsPerLID_old_ ().getConst (), 1,
6720 destMat->numImportPacketsPerLID_old_ ());
6721 size_t totalImportPackets = 0;
6722 for (
Array_size_type i = 0; i < destMat->numImportPacketsPerLID_old_.size (); ++i) {
6723 totalImportPackets += destMat->numImportPacketsPerLID_old_[i];
6725 destMat->imports_old_.resize (totalImportPackets);
6726 Distor.doReversePostsAndWaits (destMat->exports_old_ ().getConst (),
6727 destMat->numExportPacketsPerLID_old_ (),
6728 destMat->imports_old_ (),
6729 destMat->numImportPacketsPerLID_old_ ());
6732 Distor.doReversePostsAndWaits (destMat->exports_old_ ().getConst (),
6734 destMat->imports_old_ ());
6738 if (constantNumPackets == 0) {
6739 Distor.doPostsAndWaits (destMat->numExportPacketsPerLID_old_ ().getConst (), 1,
6740 destMat->numImportPacketsPerLID_old_ ());
6741 size_t totalImportPackets = 0;
6742 for (
Array_size_type i = 0; i < destMat->numImportPacketsPerLID_old_.size (); ++i) {
6743 totalImportPackets += destMat->numImportPacketsPerLID_old_[i];
6745 destMat->imports_old_.resize (totalImportPackets);
6746 Distor.doPostsAndWaits (destMat->exports_old_ ().getConst (),
6747 destMat->numExportPacketsPerLID_old_ (),
6748 destMat->imports_old_ (),
6749 destMat->numImportPacketsPerLID_old_ ());
6752 Distor.doPostsAndWaits (destMat->exports_old_ ().getConst (),
6754 destMat->imports_old_ ());
6772 #ifdef HAVE_TPETRA_MMM_TIMINGS 6773 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-1"))));
6776 Import_Util::unpackAndCombineWithOwningPIDsCount (*
this, RemoteLIDs,
6777 destMat->imports_old_ (),
6778 destMat->numImportPacketsPerLID_old_ (),
6779 constantNumPackets, Distor,
INSERT,
6780 NumSameIDs, PermuteToLIDs,
6782 size_t N = BaseRowMap->getNodeNumElements ();
6785 ArrayRCP<size_t> CSR_rowptr(N+1);
6786 ArrayRCP<GO> CSR_colind_GID;
6787 ArrayRCP<LO> CSR_colind_LID;
6788 ArrayRCP<Scalar> CSR_vals;
6789 CSR_colind_GID.resize (mynnz);
6790 CSR_vals.resize (mynnz);
6794 if (
typeid (LO) ==
typeid (GO)) {
6795 CSR_colind_LID = Teuchos::arcp_reinterpret_cast<LO> (CSR_colind_GID);
6798 CSR_colind_LID.resize (mynnz);
6816 Import_Util::unpackAndCombineIntoCrsArrays (*
this, RemoteLIDs, destMat->imports_old_ (),
6817 destMat->numImportPacketsPerLID_old_ (),
6818 constantNumPackets, Distor,
INSERT, NumSameIDs,
6819 PermuteToLIDs, PermuteFromLIDs, N, mynnz, MyPID,
6820 CSR_rowptr (), CSR_colind_GID (), CSR_vals (),
6821 SourcePids (), TargetPids);
6826 #ifdef HAVE_TPETRA_MMM_TIMINGS 6827 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-2"))));
6832 Teuchos::Array<int> RemotePids;
6833 Import_Util::lowCommunicationMakeColMapAndReindex (CSR_rowptr (),
6837 TargetPids, RemotePids,
6844 ReducedColMap = (MyRowMap.getRawPtr () == MyColMap.getRawPtr ()) ?
6846 MyColMap->replaceCommWithSubset (ReducedComm);
6847 MyColMap = ReducedColMap;
6851 destMat->replaceColMap (MyColMap);
6858 if (ReducedComm.is_null ()) {
6865 #ifdef HAVE_TPETRA_MMM_TIMINGS 6866 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC Unpack-3"))));
6868 Import_Util::sortCrsEntries (CSR_rowptr (),
6871 if ((! reverseMode && xferAsImport != NULL) ||
6872 (reverseMode && xferAsExport != NULL)) {
6873 Import_Util::sortCrsEntries (CSR_rowptr (),
6877 else if ((! reverseMode && xferAsExport != NULL) ||
6878 (reverseMode && xferAsImport != NULL)) {
6879 Import_Util::sortAndMergeCrsEntries (CSR_rowptr (),
6882 if (CSR_rowptr[N] != mynnz) {
6883 CSR_colind_LID.resize (CSR_rowptr[N]);
6884 CSR_vals.resize (CSR_rowptr[N]);
6888 TEUCHOS_TEST_FOR_EXCEPTION(
6889 true, std::logic_error,
"Tpetra::CrsMatrix::" 6890 "transferAndFillComplete: Should never get here! " 6891 "Please report this bug to a Tpetra developer.");
6902 destMat->setAllValues (CSR_rowptr, CSR_colind_LID, CSR_vals);
6908 Teuchos::ParameterList esfc_params;
6909 #ifdef HAVE_TPETRA_MMM_TIMINGS 6910 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC CreateImporter"))));
6912 RCP<import_type> MyImport = rcp (
new import_type (MyDomainMap, MyColMap, RemotePids));
6913 #ifdef HAVE_TPETRA_MMM_TIMINGS 6914 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(prefix + std::string(
"TAFC ESFC"))));
6916 esfc_params.set(
"Timer Label",prefix + std::string(
"TAFC"));
6919 destMat->expertStaticFillComplete (MyDomainMap, MyRangeMap, MyImport,Teuchos::null,rcp(&esfc_params,
false));
6922 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6927 const Teuchos::RCP<const map_type>& domainMap,
6928 const Teuchos::RCP<const map_type>& rangeMap,
6929 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6931 transferAndFillComplete (destMatrix, importer, domainMap, rangeMap, params);
6935 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
6940 const Teuchos::RCP<const map_type>& domainMap,
6941 const Teuchos::RCP<const map_type>& rangeMap,
6942 const Teuchos::RCP<Teuchos::ParameterList>& params)
const 6944 transferAndFillComplete (destMatrix, exporter, domainMap, rangeMap, params);
6955 #define TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR,LO,GO,NODE) \ 6957 template class CrsMatrix< SCALAR , LO , GO , NODE >; \ 6958 template RCP< CrsMatrix< SCALAR , LO , GO , NODE > > \ 6959 CrsMatrix< SCALAR , LO , GO , NODE >::convert< SCALAR > () const; 6961 #define TPETRA_CRSMATRIX_CONVERT_INSTANT(SO,SI,LO,GO,NODE) \ 6963 template RCP< CrsMatrix< SO , LO , GO , NODE > > \ 6964 CrsMatrix< SI , LO , GO , NODE >::convert< SO > () const; 6966 #define TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 6968 RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 6969 importAndFillCompleteCrsMatrix (const RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 6970 const Import<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6971 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6972 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& importer, \ 6973 const RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6974 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6975 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 6976 const RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6977 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6978 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 6979 const RCP<Teuchos::ParameterList>& params); 6981 #define TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 6983 RCP<CrsMatrix<SCALAR, LO, GO, NODE> > \ 6984 exportAndFillCompleteCrsMatrix (const RCP<const CrsMatrix<SCALAR, LO, GO, NODE> >& sourceMatrix, \ 6985 const Export<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6986 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6987 CrsMatrix<SCALAR, LO, GO, NODE>::node_type>& exporter, \ 6988 const RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6989 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6990 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& domainMap, \ 6991 const RCP<const Map<CrsMatrix<SCALAR, LO, GO, NODE>::local_ordinal_type, \ 6992 CrsMatrix<SCALAR, LO, GO, NODE>::global_ordinal_type, \ 6993 CrsMatrix<SCALAR, LO, GO, NODE>::node_type> >& rangeMap, \ 6994 const RCP<Teuchos::ParameterList>& params); 6996 #define TPETRA_CRSMATRIX_INSTANT(SCALAR, LO, GO ,NODE) \ 6997 TPETRA_CRSMATRIX_MATRIX_INSTANT(SCALAR, LO, GO, NODE) \ 6998 TPETRA_CRSMATRIX_IMPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) \ 6999 TPETRA_CRSMATRIX_EXPORT_AND_FILL_COMPLETE_INSTANT(SCALAR, LO, GO, NODE) 7001 #endif // TPETRA_CRSMATRIX_DEF_HPP void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
Kokkos::CrsMatrix< impl_scalar_type, LocalOrdinal, execution_space, void, typename local_graph_type::size_type > local_matrix_type
The specialization of Kokkos::CrsMatrix that represents the part of the sparse matrix on each MPI pro...
std::string description() const
A one-line description of this object.
Communication plan for data redistribution from a uniquely-owned to a (possibly) multiply-owned distr...
void reindexColumns(crs_graph_type *const graph, const Teuchos::RCP< const map_type > &newColMap, const Teuchos::RCP< const import_type > &newImport=Teuchos::null, const bool sortEachRow=true)
Reindex the column indices in place, and replace the column Map. Optionally, replace the Import objec...
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Teuchos::RCP< node_type > getNode() const
The Kokkos Node instance.
Functor for the the ABSMAX CombineMode of Import and Export operations.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
Teuchos::ArrayView< const impl_scalar_type > getView(RowInfo rowinfo) const
Constant view of all entries (including extra space) in the given row.
void getLocalRowCopy(LocalOrdinal localRow, const Teuchos::ArrayView< LocalOrdinal > &colInds, const Teuchos::ArrayView< Scalar > &vals, size_t &numEntries) const
Fill given arrays with a deep copy of the locally owned entries of the matrix in a given row...
virtual bool isLocallyIndexed() const =0
Whether matrix indices are locally indexed.
void checkInternalState() const
Check that this object's state is sane; throw if it's not.
void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to alpha.
Teuchos::RCP< const crs_graph_type > getCrsGraph() const
This matrix's graph, as a CrsGraph.
bool isNodeGlobalElement(GlobalOrdinal globalIndex) const
Whether the given global index is owned by this Map on the calling process.
Definition of the Tpetra::CrsGraph class.
void sortEntries()
Sort the entries of each row by their column indices.
Teuchos::RCP< const RowGraph< LocalOrdinal, GlobalOrdinal, Node > > getGraph() const
This matrix's graph, as a RowGraph.
virtual void copyAndPermute(const SrcDistObject &source, size_t numSameIDs, const Teuchos::ArrayView< const LocalOrdinal > &permuteToLIDs, const Teuchos::ArrayView< const LocalOrdinal > &permuteFromLIDs)
Perform copies and permutations that are local to this process.
void gaussSeidel(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps) const
"Hybrid" Jacobi + (Gauss-Seidel or SOR) on .
Teuchos::RCP< const map_type > getRowMap() const
The Map that describes the row distribution in this matrix.
LookupStatus
Return status of Map remote index lookup (getRemoteIndexList()).
void rightScale(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &x)
Details::EStorageStatus storageStatus_
Status of the matrix's storage, when not in a fill-complete state.
void getGlobalRowView(GlobalOrdinal GlobalRow, Teuchos::ArrayView< const GlobalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const
Get a constant, nonpersisting view of a row of this matrix, using global row and column indices...
Teuchos::RCP< const map_type > getColMap() const
The Map that describes the column distribution in this matrix.
One or more distributed dense vectors.
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0
The current number of entries on the calling process in the specified local row.
void mergeRedundantEntries()
Merge entries in each row with the same column indices.
bool isLocallyIndexed() const
Whether the matrix is locally indexed on the calling process.
virtual bool supportsRowViews() const
Return true if getLocalRowView() and getGlobalRowView() are valid for this object.
void deep_copy(MultiVector< DS, DL, DG, DN, dstClassic > &dst, const MultiVector< SS, SL, SG, SN, srcClassic > &src)
Copy the contents of the MultiVector src into dst.
bool isUpperTriangular() const
Indicates whether the matrix is upper triangular.
void replaceColMap(const Teuchos::RCP< const map_type > &newColMap)
Replace the matrix's column Map with the given Map.
virtual bool checkSizes(const SrcDistObject &source)
Compare the source and target (this) objects for compatibility.
void setAllValues(const typename local_matrix_type::row_map_type &rowPointers, const typename local_graph_type::entries_type::non_const_type &columnIndices, const typename local_matrix_type::values_type &values)
Sets the 1D pointer arrays of the graph.
size_t getNodeNumRows() const
The number of matrix rows owned by the calling process.
Node node_type
This class' fourth template parameter; the Kokkos device type.
bool fillComplete_
Whether the matrix is fill complete.
void deep_copy(const LittleBlock< ST2, LO > &dst, const LittleBlock< ST1, LO > &src, typename std::enable_if< std::is_convertible< ST1, ST2 >::value &&!std::is_const< ST2 >::value, int >::type *=NULL)
Copy the LittleBlock src into the LittleBlock dst.
Node::device_type device_type
The Kokkos device type.
bool isConstantStride() const
Whether this multivector has constant stride between columns.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRangeMap() const =0
The Map associated with the range of this operator, which must be compatible with Y...
std::enable_if< Kokkos::is_view< LocalIndicesViewType >::value &&Kokkos::is_view< ImplScalarViewType >::value &&std::is_same< typename LocalIndicesViewType::non_const_value_type, local_ordinal_type >::value &&std::is_same< typename ImplScalarViewType::non_const_value_type, impl_scalar_type >::value, LocalOrdinal >::type replaceLocalValues(const LocalOrdinal localRow, const typename UnmanagedView< LocalIndicesViewType >::type &inputInds, const typename UnmanagedView< ImplScalarViewType >::type &inputVals) const
Replace one or more entries' values, using local row and column indices.
GlobalOrdinal getIndexBase() const
The index base for global indices for this matrix.
void getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView< GlobalOrdinal > &Indices, const Teuchos::ArrayView< Scalar > &Values, size_t &NumEntries) const
Fill given arrays with a deep copy of the locally owned entries of the matrix in a given row...
Allocation information for a locally owned row in a CrsGraph or CrsMatrix.
Teuchos_Ordinal Array_size_type
Size type for Teuchos Array objects.
Teuchos::RCP< CrsMatrix< T, LocalOrdinal, GlobalOrdinal, Node, classic > > convert() const
Return another CrsMatrix with the same entries, but converted to a different Scalar type T...
CrsIJV()
Default constructor.
bool isNodeLocalElement(LocalOrdinal localIndex) const
Whether the given local index is valid for this Map on the calling process.
virtual void removeEmptyProcessesInPlace(const Teuchos::RCP< const map_type > &newMap)
Remove processes owning zero rows from the Maps and their communicator.
void insertLocalValues(const LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert one or more entries into the matrix, using local indices.
void getLocalRowView(LocalOrdinal LocalRow, Teuchos::ArrayView< const LocalOrdinal > &indices, Teuchos::ArrayView< const Scalar > &values) const
Get a constant, nonpersisting view of a row of this matrix, using local row and column indices...
device_type::execution_space execution_space
The Kokkos execution space.
mag_type frobNorm_
Cached Frobenius norm of the (global) matrix.
Implementation details of Tpetra.
Teuchos::ArrayView< impl_scalar_type > getViewNonConst(const RowInfo &rowinfo) const
Nonconst view of all entries (including extra space) in the given row.
void insertGlobalValues(const GlobalOrdinal globalRow, const Teuchos::ArrayView< const GlobalOrdinal > &cols, const Teuchos::ArrayView< const Scalar > &vals)
Insert one or more entries into the matrix, using global indices.
void reduce()
Sum values of a locally replicated multivector across all processes.
void fillLocalMatrix(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Fill data into the local matrix.
size_t global_size_t
Global size_t object.
LocalOrdinal sumIntoGlobalValues(const GlobalOrdinal globalRow, const Teuchos::ArrayView< const GlobalOrdinal > &cols, const Teuchos::ArrayView< const Scalar > &vals, const bool atomic=useAtomicUpdatesByDefault)
Sum into one or more sparse matrix entries, using global indices.
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, const Teuchos::ArrayView< GlobalOrdinal > &Indices, const Teuchos::ArrayView< Scalar > &Values, size_t &NumEntries) const =0
Get a copy of the given global row's entries.
Kokkos::StaticCrsGraph< LocalOrdinal, Kokkos::LayoutLeft, execution_space > local_graph_type
The type of the part of the sparse graph on each MPI process.
Traits class for "invalid" (flag) values of integer types that Tpetra uses as local ordinals or globa...
Insert new values that don't currently exist.
void getLocalDiagCopy(Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &diag) const
Get a copy of the diagonal entries of the matrix.
std::enable_if< Kokkos::is_view< LocalIndicesViewType >::value &&Kokkos::is_view< ImplScalarViewType >::value &&std::is_same< typename LocalIndicesViewType::non_const_value_type, local_ordinal_type >::value &&std::is_same< typename ImplScalarViewType::non_const_value_type, impl_scalar_type >::value, LocalOrdinal >::type sumIntoLocalValues(const LocalOrdinal localRow, const typename UnmanagedView< LocalIndicesViewType >::type &inputInds, const typename UnmanagedView< ImplScalarViewType >::type &inputVals, const bool atomic=useAtomicUpdatesByDefault) const
Sum into one or more sparse matrix entries, using local row and column indices.
dual_view_type getDualView() const
Get the Kokkos::DualView which implements local storage.
Teuchos::RCP< MV > getRowMapMultiVector(const MV &Y_rangeMap, const bool force=false) const
Create a (or fetch a cached) row Map MultiVector.
bool isLowerTriangular() const
Indicates whether the matrix is lower triangular.
ESweepDirection
Sweep direction for Gauss-Seidel or Successive Over-Relaxation (SOR).
void allocateValues(ELocalGlobal lg, GraphAllocationStatus gas)
Allocate values (and optionally indices) using the Node.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this matrix.
Communication plan for data redistribution from a (possibly) multiply-owned to a uniquely-owned distr...
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
bool isFillComplete() const
Whether the matrix is fill complete.
void scale(const Scalar &alpha)
Scale the matrix's values: this := alpha*this.
void resumeFill(const RCP< ParameterList > ¶ms=null)
Resume operations that may change the values or structure of the matrix.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getDomainMap() const =0
The Map associated with the domain of this operator, which must be compatible with X...
#define TPETRA_ABUSE_WARNING(throw_exception_test, Exception, msg)
Handle an abuse warning, according to HAVE_TPETRA_THROW_ABUSE_WARNINGS and HAVE_TPETRA_PRINT_ABUSE_WA...
void unpackAndCombine(const Teuchos::ArrayView< const LocalOrdinal > &importLIDs, const Teuchos::ArrayView< const char > &imports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t constantNumPackets, Distributor &distor, CombineMode combineMode)
Unpack the imported column indices and values, and combine into matrix.
Sets up and executes a communication plan for a Tpetra DistObject.
void reorderedGaussSeidel(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Teuchos::ArrayView< LocalOrdinal > &rowIndices, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps) const
Reordered "Hybrid" Jacobi + (Gauss-Seidel or SOR) on .
Struct representing a sparse matrix entry as an i,j,v triplet.
GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const
The global index corresponding to the given local index.
CombineMode
Rule for combining data in an Import or Export.
bool isStorageOptimized() const
Returns true if storage has been optimized.
Sum new values into existing values.
Teuchos::RCP< Node > getNode() const
Get this Map's Node object.
void exportAndFillComplete(Teuchos::RCP< CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > &destMatrix, const export_type &exporter, const Teuchos::RCP< const map_type > &domainMap=Teuchos::null, const Teuchos::RCP< const map_type > &rangeMap=Teuchos::null, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null) const
Export from this to the given destination matrix, and make the result fill complete.
Teuchos::RCP< const map_type > map_
The Map over which this object is distributed.
size_t getNodeNumDiags() const
Returns the number of local diagonal entries, based on global row/column index comparisons.
Utility functions for packing and unpacking sparse matrix entries.
bool isDistributed() const
Whether this is a globally distributed object.
ProfileType getProfileType() const
Returns true if the matrix was allocated with static data structures.
virtual ~CrsMatrix()
Destructor.
Replace old value with maximum of magnitudes of old and new values.
Abstract base class for objects that can be the source of an Import or Export operation.
global_size_t getGlobalNumRows() const
Number of global elements in the row map of this matrix.
global_size_t getGlobalNumDiags() const
Returns the number of global diagonal entries, based on global row/column index comparisons.
void replaceDomainMapAndImporter(const Teuchos::RCP< const map_type > &newDomainMap, Teuchos::RCP< const import_type > &newImporter)
Replace the current domain Map and Import with the given objects.
Teuchos::RCP< MV > getColumnMapMultiVector(const MV &X_domainMap, const bool force=false) const
Create a (or fetch a cached) column Map MultiVector.
Replace existing values with new values.
#define TPETRA_EFFICIENCY_WARNING(throw_exception_test, Exception, msg)
Print or throw an efficency warning.
Binary function that returns its second argument.
void reindexColumns(const Teuchos::RCP< const map_type > &newColMap, const Teuchos::RCP< const import_type > &newImport=Teuchos::null, const bool sortIndicesInEachRow=true)
Reindex the column indices in place, and replace the column Map. Optionally, replace the Import objec...
bool hasTransposeApply() const
Whether apply() allows applying the transpose or conjugate transpose.
Scalar v
Value of matrix entry.
void computeGlobalConstants()
Compute matrix properties that require collectives.
bool isFillActive() const
Whether the matrix is not fill complete.
bool isStaticGraph() const
Indicates that the graph is static, so that new entries cannot be added to this matrix.
Replace old values with zero.
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
global_size_t getGlobalNumCols() const
The number of global columns in the matrix.
bool hasColMap() const
Indicates whether the matrix has a well-defined column map.
mag_type getFrobeniusNorm() const
Compute and return the Frobenius norm of the matrix.
Kokkos::Details::ArithTraits< impl_scalar_type >::mag_type mag_type
Type of a norm result.
LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const
The local index corresponding to the given global index.
void leftScale(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &x)
A distributed graph accessed by rows (adjacency lists) and stored sparsely.
Teuchos::RCP< MV > importMV_
Column Map MultiVector used in apply() and gaussSeidel().
Ordinal i
(Global) row index
Kokkos::DualView< impl_scalar_type **, Kokkos::LayoutLeft, typename execution_space::execution_space > dual_view_type
Kokkos::DualView specialization used by this class.
Describes a parallel distribution of objects over processes.
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
global_size_t getGlobalNumEntries() const
The global number of entries in this matrix.
Internal functions and macros designed for use with Tpetra::Import and Tpetra::Export objects...
A read-only, row-oriented interface to a sparse matrix.
Scalar operator()(const Scalar &x, const Scalar &y)
Return the maximum of the magnitudes (absolute values) of x and y.
std::map< GlobalOrdinal, std::pair< Teuchos::Array< GlobalOrdinal >, Teuchos::Array< Scalar > > > nonlocals_
Nonlocal data added using insertGlobalValues().
A distributed dense vector.
Stand-alone utility functions and macros.
virtual void pack(const Teuchos::ArrayView< const LocalOrdinal > &exportLIDs, Teuchos::Array< char > &exports, const Teuchos::ArrayView< size_t > &numPacketsPerLID, size_t &constantNumPackets, Distributor &distor) const
Pack this object's data for an Import or Export.
void expertStaticFillComplete(const RCP< const map_type > &domainMap, const RCP< const map_type > &rangeMap, const RCP< const import_type > &importer=Teuchos::null, const RCP< const export_type > &exporter=Teuchos::null, const RCP< ParameterList > ¶ms=Teuchos::null)
Perform a fillComplete on a matrix that already has data.
Ordinal j
(Global) column index
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
size_t getNumVectors() const
Number of columns in the multivector.
void doExport(const SrcDistObject &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)
Export data into this object using an Export object ("forward mode").
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
void globalAssemble()
Communicate nonlocal contributions to other processes.
Kokkos::Details::ArithTraits< Scalar >::val_type impl_scalar_type
The type used internally in place of Scalar.
LocalOrdinal replaceGlobalValues(const GlobalOrdinal globalRow, const Kokkos::View< const GlobalOrdinal *, device_type, Kokkos::MemoryUnmanaged > &cols, const Kokkos::View< const impl_scalar_type *, device_type, Kokkos::MemoryUnmanaged > &vals) const
Replace one or more entries' values, using global indices.
void gaussSeidelCopy(MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps, const bool zeroInitialGuess) const
Version of gaussSeidel(), with fewer requirements on X.
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
void applyNonTranspose(const MV &X_in, MV &Y_in, Scalar alpha, Scalar beta) const
Special case of apply() for mode == Teuchos::NO_TRANS.
Teuchos::RCP< MV > exportMV_
Row Map MultiVector used in apply().
void getLocalDiagOffsets(Teuchos::ArrayRCP< size_t > &offsets) const
Get offsets of the diagonal entries in the matrix.
void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Compute a sparse matrix-MultiVector multiply.
size_t getNodeNumCols() const
The number of columns connected to the locally owned rows of this matrix.
void reorderedGaussSeidelCopy(MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &X, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &D, const Teuchos::ArrayView< LocalOrdinal > &rowIndices, const Scalar &dampingFactor, const ESweepDirection direction, const int numSweeps, const bool zeroInitialGuess) const
Version of reorderedGaussSeidel(), with fewer requirements on X.
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries on this node in the specified global row.
size_t getNodeNumEntries() const
The local number of entries in this matrix.
virtual Teuchos::RCP< RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > add(const Scalar &alpha, const RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap, const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms) const
Implementation of RowMatrix::add: return alpha*A + beta*this.
void clearGlobalConstants()
Clear matrix properties that require collectives.
CrsIJV(Ordinal row, Ordinal col, const Scalar &val)
Standard constructor.
bool isGloballyIndexed() const
Whether the matrix is globally indexed on the calling process.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
The Map that describes the distribution of rows over processes.
void applyTranspose(const MV &X_in, MV &Y_in, const Teuchos::ETransp mode, Scalar alpha, Scalar beta) const
Special case of apply() for mode != Teuchos::NO_TRANS.
void fillComplete(const RCP< const map_type > &domainMap, const RCP< const map_type > &rangeMap, const RCP< ParameterList > ¶ms=null)
Signal that data entry is complete, specifying domain and range maps.
local_matrix_type lclMatrix_
The local sparse matrix.
void importAndFillComplete(Teuchos::RCP< CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > &destMatrix, const import_type &importer, const Teuchos::RCP< const map_type > &domainMap, const Teuchos::RCP< const map_type > &rangeMap, const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null) const
Import from this to the given destination matrix, and make the result fill complete.
void fillLocalGraphAndMatrix(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Fill data into the local graph and matrix.
Teuchos::RCP< const map_type > getRangeMap() const
The range Map of this matrix.
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.