42 #ifndef TPETRA_MULTIVECTOR_DEF_HPP 43 #define TPETRA_MULTIVECTOR_DEF_HPP 54 #include "Tpetra_Vector.hpp" 55 #include "Tpetra_Details_MultiVectorDistObjectKernels.hpp" 56 #include "Tpetra_KokkosRefactor_Details_MultiVectorDistObjectKernels.hpp" 58 #include "KokkosCompat_View.hpp" 59 #include "Kokkos_MV_GEMM.hpp" 60 #include "Kokkos_Blas1_MV.hpp" 61 #include "Kokkos_Random.hpp" 63 #ifdef HAVE_TPETRA_INST_FLOAT128 66 template<
class Generator>
67 struct rand<Generator, __float128> {
68 static KOKKOS_INLINE_FUNCTION __float128 max ()
70 return static_cast<__float128
> (1.0);
72 static KOKKOS_INLINE_FUNCTION __float128
77 const __float128 scalingFactor =
78 static_cast<__float128
> (std::numeric_limits<double>::min ()) /
79 static_cast<__float128> (2.0);
80 const __float128 higherOrderTerm =
static_cast<__float128
> (gen.drand ());
81 const __float128 lowerOrderTerm =
82 static_cast<__float128
> (gen.drand ()) * scalingFactor;
83 return higherOrderTerm + lowerOrderTerm;
85 static KOKKOS_INLINE_FUNCTION __float128
86 draw (Generator& gen,
const __float128& range)
89 const __float128 scalingFactor =
90 static_cast<__float128
> (std::numeric_limits<double>::min ()) /
91 static_cast<__float128> (2.0);
92 const __float128 higherOrderTerm =
93 static_cast<__float128
> (gen.drand (range));
94 const __float128 lowerOrderTerm =
95 static_cast<__float128
> (gen.drand (range)) * scalingFactor;
96 return higherOrderTerm + lowerOrderTerm;
98 static KOKKOS_INLINE_FUNCTION __float128
99 draw (Generator& gen,
const __float128& start,
const __float128& end)
102 const __float128 scalingFactor =
103 static_cast<__float128
> (std::numeric_limits<double>::min ()) /
104 static_cast<__float128> (2.0);
105 const __float128 higherOrderTerm =
106 static_cast<__float128
> (gen.drand (start, end));
107 const __float128 lowerOrderTerm =
108 static_cast<__float128
> (gen.drand (start, end)) * scalingFactor;
109 return higherOrderTerm + lowerOrderTerm;
113 #endif // HAVE_TPETRA_INST_FLOAT128 129 template<
class ST,
class LO,
class GO,
class NT>
131 allocDualView (
const size_t lclNumRows,
const size_t numCols,
const bool zeroOut =
true)
134 const char* label =
"MV::DualView";
138 return dual_view_type (label, lclNumRows, numCols);
141 return dual_view_type (label, lclNumRows, numCols);
149 typename dual_view_type::t_dev d_view (Kokkos::ViewAllocateWithoutInitializing (label), lclNumRows, numCols);
150 #ifdef HAVE_TPETRA_DEBUG 158 const ST nan = Kokkos::Details::ArithTraits<ST>::nan ();
159 KokkosBlas::fill (d_view, nan);
160 #endif // HAVE_TPETRA_DEBUG 161 typename dual_view_type::t_host h_view = Kokkos::create_mirror_view (d_view);
166 dual_view_type dv (d_view, h_view);
167 dv.template modify<typename dual_view_type::t_dev::memory_space> ();
169 return dual_view_type (d_view, h_view);
178 template<
class T,
class ExecSpace>
179 struct MakeUnmanagedView {
188 typedef typename Kokkos::Impl::if_c<
189 Kokkos::Impl::VerifyExecutionCanAccessMemorySpace<ExecSpace, Kokkos::HostSpace>::value,
190 typename ExecSpace::device_type,
191 typename Kokkos::DualView<T*, ExecSpace>::host_mirror_space>::type host_exec_space;
192 typedef Kokkos::LayoutLeft array_layout;
193 typedef Kokkos::View<T*, array_layout, host_exec_space,
194 Kokkos::MemoryUnmanaged> view_type;
196 static view_type getView (
const Teuchos::ArrayView<T>& x_in)
198 const size_t numEnt =
static_cast<size_t> (x_in.size ());
202 return view_type (x_in.getRawPtr (), numEnt);
210 template<
class DualViewType>
212 takeSubview (
const DualViewType& X,
214 #ifdef KOKKOS_USING_EXPERIMENTAL_VIEW
215 const Kokkos::Experimental::Impl::ALL_t&,
219 const std::pair<size_t, size_t>& colRng)
221 if (X.dimension_0 () == 0 && X.dimension_1 () != 0) {
222 return DualViewType (
"MV::DualView", 0, colRng.second - colRng.first);
225 return subview (X, Kokkos::ALL (), colRng);
232 template<
class DualViewType>
234 takeSubview (
const DualViewType& X,
235 const std::pair<size_t, size_t>& rowRng,
236 const std::pair<size_t, size_t>& colRng)
238 if (X.dimension_0 () == 0 && X.dimension_1 () != 0) {
239 return DualViewType (
"MV::DualView", 0, colRng.second - colRng.first);
242 return subview (X, rowRng, colRng);
250 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
252 MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
253 vectorIndexOutOfRange (
const size_t VectorIndex)
const {
254 return (VectorIndex < 1 && VectorIndex != 0) || VectorIndex >= getNumVectors();
257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
258 MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic>::
263 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
266 const size_t numVecs,
267 const bool zeroOut) :
270 TEUCHOS_TEST_FOR_EXCEPTION(
271 numVecs < 1, std::invalid_argument,
"Tpetra::MultiVector::MultiVector" 272 "(map,numVecs,zeroOut): numVecs = " << numVecs <<
" < 1.");
274 view_ = allocDualView<Scalar, LocalOrdinal, GlobalOrdinal, Node> (lclNumRows, numVecs, zeroOut);
278 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
287 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
290 const Teuchos::DataAccess copyOrView) :
297 const char tfecfFuncName[] =
"MultiVector(const MultiVector&, " 298 "const Teuchos::DataAccess): ";
300 if (copyOrView == Teuchos::Copy) {
304 this->
view_ = cpy.view_;
308 else if (copyOrView == Teuchos::View) {
311 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
312 true, std::invalid_argument,
"Second argument 'copyOrView' has an " 313 "invalid value " << copyOrView <<
". Valid values include " 314 "Teuchos::Copy = " << Teuchos::Copy <<
" and Teuchos::View = " 315 << Teuchos::View <<
".");
319 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
322 const dual_view_type& view) :
327 const char tfecfFuncName[] =
"MultiVector(map,view): ";
333 const size_t LDA = (
origView_.dimension_1 () > 1) ? stride[1] :
336 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
337 LDA < lclNumRows, std::invalid_argument,
"The input Kokkos::DualView's " 338 "column stride LDA = " << LDA <<
" < getLocalLength() = " << lclNumRows
339 <<
". This may also mean that the input view's first dimension (number " 340 "of rows = " << view.dimension_0 () <<
") does not not match the number " 341 "of entries in the Map on the calling process.");
345 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
348 const typename dual_view_type::t_dev& d_view) :
351 using Teuchos::ArrayRCP;
353 const char tfecfFuncName[] =
"MultiVector(map,d_view): ";
358 d_view.stride (stride);
359 const size_t LDA = (d_view.dimension_1 () > 1) ? stride[1] :
360 d_view.dimension_0 ();
362 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
363 LDA < lclNumRows, std::invalid_argument,
"The input Kokkos::View's " 364 "column stride LDA = " << LDA <<
" < getLocalLength() = " << lclNumRows
365 <<
". This may also mean that the input view's first dimension (number " 366 "of rows = " << d_view.dimension_0 () <<
") does not not match the " 367 "number of entries in the Map on the calling process.");
376 view_.template modify<execution_space> ();
380 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
383 const dual_view_type& view,
384 const dual_view_type& origView) :
389 const char tfecfFuncName[] =
"MultiVector(map,view,origView): ";
395 const size_t LDA = (
origView_.dimension_1 () > 1) ? stride[1] :
398 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
399 LDA < lclNumRows, std::invalid_argument,
"The input Kokkos::DualView's " 400 "column stride LDA = " << LDA <<
" < getLocalLength() = " << lclNumRows
401 <<
". This may also mean that the input origView's first dimension (number " 402 "of rows = " << origView.dimension_0 () <<
") does not not match the number " 403 "of entries in the Map on the calling process.");
407 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
410 const dual_view_type& view,
411 const Teuchos::ArrayView<const size_t>& whichVectors) :
418 using Kokkos::subview;
419 const char tfecfFuncName[] =
"MultiVector(map,view,whichVectors): ";
421 const size_t lclNumRows = map.is_null () ? size_t (0) :
422 map->getNodeNumElements ();
429 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
430 view.dimension_1 () != 0 &&
static_cast<size_t> (view.dimension_0 ()) < lclNumRows,
431 std::invalid_argument,
"view.dimension_0() = " << view.dimension_0 ()
432 <<
" < map->getNodeNumElements() = " << lclNumRows <<
".");
433 if (whichVectors.size () != 0) {
434 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
435 view.dimension_1 () != 0 && view.dimension_1 () == 0,
436 std::invalid_argument,
"view.dimension_1() = 0, but whichVectors.size()" 437 " = " << whichVectors.size () <<
" > 0.");
438 size_t maxColInd = 0;
439 typedef Teuchos::ArrayView<const size_t>::size_type size_type;
440 for (size_type k = 0; k < whichVectors.size (); ++k) {
441 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
442 whichVectors[k] == Teuchos::OrdinalTraits<size_t>::invalid (),
443 std::invalid_argument,
"whichVectors[" << k <<
"] = " 444 "Teuchos::OrdinalTraits<size_t>::invalid().");
445 maxColInd = std::max (maxColInd, whichVectors[k]);
447 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
448 view.dimension_1 () != 0 &&
static_cast<size_t> (view.dimension_1 ()) <= maxColInd,
449 std::invalid_argument,
"view.dimension_1() = " << view.dimension_1 ()
450 <<
" <= max(whichVectors) = " << maxColInd <<
".");
457 const size_t LDA = (
origView_.dimension_1 () > 1) ? stride[1] :
459 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
460 LDA < lclNumRows, std::invalid_argument,
461 "LDA = " << LDA <<
" < this->getLocalLength() = " << lclNumRows <<
".");
463 if (whichVectors.size () == 1) {
472 const std::pair<size_t, size_t> colRng (whichVectors[0],
473 whichVectors[0] + 1);
481 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
484 const dual_view_type& view,
485 const dual_view_type& origView,
486 const Teuchos::ArrayView<const size_t>& whichVectors) :
493 using Kokkos::subview;
494 const char tfecfFuncName[] =
"MultiVector(map,view,origView,whichVectors): ";
503 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
504 view.dimension_1 () != 0 &&
static_cast<size_t> (view.dimension_0 ()) < lclNumRows,
505 std::invalid_argument,
"view.dimension_0() = " << view.dimension_0 ()
506 <<
" < map->getNodeNumElements() = " << lclNumRows <<
".");
507 if (whichVectors.size () != 0) {
508 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
509 view.dimension_1 () != 0 && view.dimension_1 () == 0,
510 std::invalid_argument,
"view.dimension_1() = 0, but whichVectors.size()" 511 " = " << whichVectors.size () <<
" > 0.");
512 size_t maxColInd = 0;
513 typedef Teuchos::ArrayView<const size_t>::size_type size_type;
514 for (size_type k = 0; k < whichVectors.size (); ++k) {
515 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
516 whichVectors[k] == Teuchos::OrdinalTraits<size_t>::invalid (),
517 std::invalid_argument,
"whichVectors[" << k <<
"] = " 518 "Teuchos::OrdinalTraits<size_t>::invalid().");
519 maxColInd = std::max (maxColInd, whichVectors[k]);
521 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
522 view.dimension_1 () != 0 &&
static_cast<size_t> (view.dimension_1 ()) <= maxColInd,
523 std::invalid_argument,
"view.dimension_1() = " << view.dimension_1 ()
524 <<
" <= max(whichVectors) = " << maxColInd <<
".");
530 const size_t LDA = (
origView_.dimension_1 () > 1) ? stride[1] :
532 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
533 LDA < lclNumRows, std::invalid_argument,
"Input DualView's column stride" 534 " = " << LDA <<
" < this->getLocalLength() = " << lclNumRows <<
".");
536 if (whichVectors.size () == 1) {
545 const std::pair<size_t, size_t> colRng (whichVectors[0],
546 whichVectors[0] + 1);
554 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
557 const Teuchos::ArrayView<const Scalar>& data,
559 const size_t numVecs) :
562 using Kokkos::subview;
563 using Teuchos::ArrayView;
564 using Teuchos::av_reinterpret_cast;
566 typedef LocalOrdinal LO;
567 typedef GlobalOrdinal GO;
568 typedef typename dual_view_type::host_mirror_space HMS;
569 typedef MakeUnmanagedView<const IST, device_type> view_getter_type;
570 typedef typename view_getter_type::view_type in_view_type;
571 typedef Kokkos::View<IST*, Kokkos::LayoutLeft, HMS> out_view_type;
572 const char tfecfFuncName[] =
"MultiVector(map,data,LDA,numVecs): ";
577 const size_t lclNumRows =
578 map.is_null () ? size_t (0) : map->getNodeNumElements ();
579 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
580 (LDA < lclNumRows, std::invalid_argument,
"LDA = " << LDA <<
" < " 581 "map->getNodeNumElements() = " << lclNumRows <<
".");
583 const size_t minNumEntries = LDA * (numVecs - 1) + lclNumRows;
584 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC
585 (static_cast<size_t> (data.size ()) < minNumEntries,
586 std::invalid_argument,
"Input Teuchos::ArrayView does not have enough " 587 "entries, given the input Map and number of vectors in the MultiVector." 588 " data.size() = " << data.size () <<
" < (LDA*(numVecs-1)) + " 589 "map->getNodeNumElements () = " << minNumEntries <<
".");
591 view_ = allocDualView<Scalar, LO, GO, Node> (lclNumRows, numVecs);
592 view_.template modify<HMS> ();
594 ArrayView<const IST> X_in_av = av_reinterpret_cast<
const IST> (data);
595 in_view_type X_in = view_getter_type::getView (X_in_av);
596 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
597 for (
size_t j = 0; j < numVecs; ++j) {
598 const std::pair<size_t, size_t> rng (j*LDA, j*LDA + lclNumRows);
599 in_view_type X_j_in = subview (X_in, rng);
600 out_view_type X_j_out = subview (view_.h_view, rowRng, j);
606 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
608 MultiVector (
const Teuchos::RCP<const map_type>& map,
609 const Teuchos::ArrayView<
const ArrayView<const Scalar> >& ArrayOfPtrs,
610 const size_t numVecs) :
613 using Kokkos::subview;
614 using Teuchos::ArrayView;
615 using Teuchos::av_reinterpret_cast;
617 typedef LocalOrdinal LO;
618 typedef GlobalOrdinal GO;
619 typedef typename dual_view_type::host_mirror_space HMS;
620 typedef MakeUnmanagedView<const IST, device_type> view_getter_type;
621 typedef typename view_getter_type::view_type in_view_type;
622 typedef Kokkos::View<IST*, Kokkos::LayoutLeft, HMS> out_view_type;
623 const char tfecfFuncName[] =
"MultiVector(map,ArrayOfPtrs,numVecs): ";
625 const size_t lclNumRows =
626 map.is_null () ? size_t (0) : map->getNodeNumElements ();
627 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
628 numVecs < 1 || numVecs != static_cast<size_t> (ArrayOfPtrs.size ()),
630 "ArrayOfPtrs.size() must be strictly positive and as large as ArrayOfPtrs.");
631 for (
size_t j = 0; j < numVecs; ++j) {
632 ArrayView<const Scalar> X_j_av = ArrayOfPtrs[j];
633 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
634 static_cast<size_t> (X_j_av.size ()) < lclNumRows,
635 std::invalid_argument,
"ArrayOfPtrs[" << j <<
"].size() = " 636 << X_j_av.size () <<
" < map->getNodeNumElements() = " << lclNumRows
640 view_ = allocDualView<Scalar, LO, GO, Node> (lclNumRows, numVecs);
641 view_.template modify<HMS> ();
643 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
644 for (
size_t j = 0; j < numVecs; ++j) {
645 ArrayView<const IST> X_j_av =
646 av_reinterpret_cast<
const IST> (ArrayOfPtrs[j]);
647 in_view_type X_j_in (X_j_av.getRawPtr (), lclNumRows);
648 out_view_type X_j_out = subview (view_.h_view, rowRng, j);
654 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
658 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
664 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
669 if (this->
getMap ().is_null ()) {
670 return static_cast<size_t> (0);
672 return this->
getMap ()->getNodeNumElements ();
676 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
681 if (this->
getMap ().is_null ()) {
682 return static_cast<size_t> (0);
684 return this->
getMap ()->getGlobalNumElements ();
688 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
698 const size_t LDA = (
origView_.dimension_1 () > 1) ? stride[1] :
origView_.dimension_0 ();
702 return static_cast<size_t> (0);
706 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
715 const this_type* src =
dynamic_cast<const this_type*
> (&sourceObj);
729 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
736 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
741 const Kokkos::View<const LocalOrdinal*, execution_space>& permuteToLIDs,
742 const Kokkos::View<const LocalOrdinal*, execution_space>& permuteFromLIDs)
744 using Teuchos::ArrayRCP;
745 using Teuchos::ArrayView;
747 using Kokkos::Compat::getKokkosViewDeepCopy;
748 using Kokkos::subview;
750 typename dual_view_type::array_layout,
754 const char tfecfFuncName[] =
"copyAndPermute";
756 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
757 permuteToLIDs.size() != permuteFromLIDs.size(), std::runtime_error,
758 ": permuteToLIDs and permuteFromLIDs must have the same size." 759 << std::endl <<
"permuteToLIDs.size() = " << permuteToLIDs.size ()
760 <<
" != permuteFromLIDs.size() = " << permuteFromLIDs.size () <<
".");
763 const MV& sourceMV =
dynamic_cast<const MV&
> (sourceObj);
787 if (numSameIDs > 0) {
788 const std::pair<size_t, size_t> rows (0, numSameIDs);
789 for (
size_t j = 0; j < numCols; ++j) {
791 const size_t srcCol =
792 sourceMV.isConstantStride () ? j : sourceMV.whichVectors_[j];
793 col_dual_view_type dst_j =
794 subview (
view_, rows, dstCol);
795 col_dual_view_type src_j =
796 subview (sourceMV.view_, rows, srcCol);
811 if (permuteFromLIDs.size() == 0 || permuteToLIDs.size() == 0)
815 KokkosRefactor::Details::permute_array_multi_column(
817 sourceMV.getKokkosView(),
823 KokkosRefactor::Details::permute_array_multi_column_variable_stride(
825 sourceMV.getKokkosView(),
829 getKokkosViewDeepCopy<execution_space> (sourceMV.whichVectors_ ()),
834 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
838 const Kokkos::View<const local_ordinal_type*, execution_space> &exportLIDs,
839 Kokkos::View<impl_scalar_type*, execution_space> &exports,
840 const Kokkos::View<size_t*, execution_space> &numExportPacketsPerLID,
841 size_t& constantNumPackets,
844 using Teuchos::Array;
845 using Teuchos::ArrayView;
846 using Kokkos::Compat::getKokkosViewDeepCopy;
851 if (exportLIDs.size () == 0) {
856 const MV& sourceMV =
dynamic_cast<const MV&
> (sourceObj);
860 (void) numExportPacketsPerLID;
882 constantNumPackets = numCols;
884 const size_t numExportLIDs = exportLIDs.size ();
885 const size_t newExportsSize = numCols * numExportLIDs;
886 if (exports.size () != newExportsSize) {
887 Kokkos::Compat::realloc (exports, newExportsSize);
898 if (sourceMV.isConstantStride ()) {
899 KokkosRefactor::Details::pack_array_single_column(
901 sourceMV.getKokkosView (),
906 KokkosRefactor::Details::pack_array_single_column(
908 sourceMV.getKokkosView (),
910 sourceMV.whichVectors_[0]);
914 if (sourceMV.isConstantStride ()) {
915 KokkosRefactor::Details::pack_array_multi_column(
917 sourceMV.getKokkosView (),
922 KokkosRefactor::Details::pack_array_multi_column_variable_stride(
924 sourceMV.getKokkosView (),
926 getKokkosViewDeepCopy<execution_space> (sourceMV.whichVectors_ ()),
933 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
936 unpackAndCombineNew (
const Kokkos::View<const local_ordinal_type*, execution_space> &importLIDs,
937 const Kokkos::View<const impl_scalar_type*, execution_space> &imports,
938 const Kokkos::View<size_t*, execution_space> &numPacketsPerLID,
939 size_t constantNumPackets,
943 using Teuchos::ArrayView;
944 using Kokkos::Compat::getKokkosViewDeepCopy;
945 const char tfecfFuncName[] =
"unpackAndCombine";
948 if (importLIDs.size () == 0) {
954 (void) numPacketsPerLID;
966 #ifdef HAVE_TPETRA_DEBUG 967 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
968 static_cast<size_t> (imports.size()) !=
getNumVectors()*importLIDs.size(),
970 ": 'imports' buffer size must be consistent with the amount of data to " 971 "be sent. " << std::endl <<
"imports.size() = " << imports.size()
972 <<
" != getNumVectors()*importLIDs.size() = " <<
getNumVectors() <<
"*" 973 << importLIDs.size() <<
" = " <<
getNumVectors() * importLIDs.size()
976 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
977 constantNumPackets == static_cast<size_t> (0), std::runtime_error,
978 ": constantNumPackets input argument must be nonzero.");
980 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
981 static_cast<size_t> (numVecs) != static_cast<size_t> (constantNumPackets),
982 std::runtime_error,
": constantNumPackets must equal numVecs.");
983 #endif // HAVE_TPETRA_DEBUG 985 if (numVecs > 0 && importLIDs.size () > 0) {
993 KokkosRefactor::Details::unpack_array_multi_column(
997 KokkosRefactor::Details::InsertOp(),
1001 KokkosRefactor::Details::unpack_array_multi_column_variable_stride(
1006 KokkosRefactor::Details::InsertOp(),
1010 else if (CM ==
ADD) {
1012 KokkosRefactor::Details::unpack_array_multi_column(
1016 KokkosRefactor::Details::AddOp(),
1020 KokkosRefactor::Details::unpack_array_multi_column_variable_stride(
1025 KokkosRefactor::Details::AddOp(),
1031 KokkosRefactor::Details::unpack_array_multi_column(
1035 KokkosRefactor::Details::AbsMaxOp(),
1039 KokkosRefactor::Details::unpack_array_multi_column_variable_stride(
1044 KokkosRefactor::Details::AbsMaxOp(),
1049 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1051 std::invalid_argument,
": Invalid CombineMode: " << CM <<
". Valid " 1052 "CombineMode values are ADD, REPLACE, INSERT, and ABSMAX.");
1057 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1063 return static_cast<size_t> (
view_.dimension_1 ());
1071 template<
class RV,
class XMV>
1073 lclDotImpl (
const RV& dotsOut,
1076 const size_t lclNumRows,
1077 const size_t numVecs,
1078 const Teuchos::ArrayView<const size_t>& whichVecsX,
1079 const Teuchos::ArrayView<const size_t>& whichVecsY,
1080 const bool constantStrideX,
1081 const bool constantStrideY)
1084 using Kokkos::subview;
1085 typedef typename RV::non_const_value_type
dot_type;
1086 #ifdef HAVE_TPETRA_DEBUG 1087 const char prefix[] =
"Tpetra::MultiVector::lclDotImpl: ";
1088 #endif // HAVE_TPETRA_DEBUG 1090 static_assert (Kokkos::Impl::is_view<RV>::value,
1091 "Tpetra::MultiVector::lclDotImpl: " 1092 "The first argument dotsOut is not a Kokkos::View.");
1093 static_assert (RV::rank == 1,
"Tpetra::MultiVector::lclDotImpl: " 1094 "The first argument dotsOut must have rank 1.");
1095 static_assert (Kokkos::Impl::is_view<XMV>::value,
1096 "Tpetra::MultiVector::lclDotImpl: The type of the 2nd and " 1097 "3rd arguments (X_lcl and Y_lcl) is not a Kokkos::View.");
1098 static_assert (XMV::rank == 2,
"Tpetra::MultiVector::lclDotImpl: " 1099 "X_lcl and Y_lcl must have rank 2.");
1104 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
1105 const std::pair<size_t, size_t> colRng (0, numVecs);
1106 RV theDots = subview (dotsOut, colRng);
1107 XMV X = subview (X_lcl, rowRng, Kokkos::ALL());
1108 XMV Y = subview (Y_lcl, rowRng, Kokkos::ALL());
1110 #ifdef HAVE_TPETRA_DEBUG 1111 if (lclNumRows != 0) {
1112 TEUCHOS_TEST_FOR_EXCEPTION
1113 (X.dimension_0 () != lclNumRows, std::logic_error, prefix <<
1114 "X.dimension_0() = " << X.dimension_0 () <<
" != lclNumRows " 1115 "= " << lclNumRows <<
". " 1116 "Please report this bug to the Tpetra developers.");
1117 TEUCHOS_TEST_FOR_EXCEPTION
1118 (Y.dimension_0 () != lclNumRows, std::logic_error, prefix <<
1119 "Y.dimension_0() = " << Y.dimension_0 () <<
" != lclNumRows " 1120 "= " << lclNumRows <<
". " 1121 "Please report this bug to the Tpetra developers.");
1125 TEUCHOS_TEST_FOR_EXCEPTION
1127 (X.dimension_0 () != lclNumRows || X.dimension_1 () != numVecs),
1128 std::logic_error, prefix <<
"X is " << X.dimension_0 () <<
" x " <<
1129 X.dimension_1 () <<
" (constant stride), which differs from the " 1130 "local dimensions " << lclNumRows <<
" x " << numVecs <<
". " 1131 "Please report this bug to the Tpetra developers.");
1132 TEUCHOS_TEST_FOR_EXCEPTION
1133 (! constantStrideX &&
1134 (X.dimension_0 () != lclNumRows || X.dimension_1 () < numVecs),
1135 std::logic_error, prefix <<
"X is " << X.dimension_0 () <<
" x " <<
1136 X.dimension_1 () <<
" (NOT constant stride), but the local " 1137 "dimensions are " << lclNumRows <<
" x " << numVecs <<
". " 1138 "Please report this bug to the Tpetra developers.");
1139 TEUCHOS_TEST_FOR_EXCEPTION
1141 (Y.dimension_0 () != lclNumRows || Y.dimension_1 () != numVecs),
1142 std::logic_error, prefix <<
"Y is " << Y.dimension_0 () <<
" x " <<
1143 Y.dimension_1 () <<
" (constant stride), which differs from the " 1144 "local dimensions " << lclNumRows <<
" x " << numVecs <<
". " 1145 "Please report this bug to the Tpetra developers.");
1146 TEUCHOS_TEST_FOR_EXCEPTION
1147 (! constantStrideY &&
1148 (Y.dimension_0 () != lclNumRows || Y.dimension_1 () < numVecs),
1149 std::logic_error, prefix <<
"Y is " << Y.dimension_0 () <<
" x " <<
1150 Y.dimension_1 () <<
" (NOT constant stride), but the local " 1151 "dimensions are " << lclNumRows <<
" x " << numVecs <<
". " 1152 "Please report this bug to the Tpetra developers.");
1154 #endif // HAVE_TPETRA_DEBUG 1156 if (lclNumRows == 0) {
1157 const dot_type zero = Kokkos::Details::ArithTraits<dot_type>::zero ();
1161 if (constantStrideX && constantStrideY) {
1162 if(X.dimension_1() == 1) {
1163 typename RV::non_const_value_type result =
1164 KokkosBlas::dot (Kokkos::subview(X,Kokkos::ALL(),0),
1165 Kokkos::subview(Y,Kokkos::ALL(),0));
1168 KokkosBlas::dot (theDots, X, Y);
1176 for (
size_t k = 0; k < numVecs; ++k) {
1177 const size_t X_col = constantStrideX ? k : whichVecsX[k];
1178 const size_t Y_col = constantStrideY ? k : whichVecsY[k];
1179 KokkosBlas::dot (subview (theDots, k), subview (X, ALL (), X_col),
1180 subview (Y, ALL (), Y_col));
1188 gblDotImpl (
const RV& dotsOut,
1189 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1190 const bool distributed)
1192 using Teuchos::REDUCE_MAX;
1193 using Teuchos::REDUCE_SUM;
1194 using Teuchos::reduceAll;
1195 typedef typename RV::non_const_value_type
dot_type;
1197 const size_t numVecs = dotsOut.dimension_0 ();
1212 if (distributed && ! comm.is_null ()) {
1218 typename RV::non_const_type lclDots (Kokkos::ViewAllocateWithoutInitializing (
"tmp"), numVecs);
1220 const dot_type*
const lclSum = lclDots.ptr_on_device ();
1221 dot_type*
const gblSum = dotsOut.ptr_on_device ();
1222 const int nv =
static_cast<int> (numVecs);
1223 reduceAll<int, dot_type> (*comm, REDUCE_SUM, nv, lclSum, gblSum);
1228 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1232 const Kokkos::View<dot_type*, device_type>& dots)
const 1234 using Kokkos::create_mirror_view;
1235 using Kokkos::subview;
1236 using Teuchos::Comm;
1237 using Teuchos::null;
1240 typedef Kokkos::View<dot_type*, device_type> RV;
1241 const char tfecfFuncName[] =
"Tpetra::MultiVector::dot: ";
1248 const size_t numDots =
static_cast<size_t> (dots.dimension_0 ());
1250 #ifdef HAVE_TPETRA_DEBUG 1252 const bool compat = this->
getMap ()->isCompatible (* (A.
getMap ()));
1253 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1254 ! compat, std::invalid_argument,
"Tpetra::MultiVector::dot: *this is " 1255 "not compatible with the input MultiVector A. We only test for this " 1256 "in a debug build.");
1258 #endif // HAVE_TPETRA_DEBUG 1267 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1269 "MultiVectors do not have the same local length. " 1270 "this->getLocalLength() = " << lclNumRows <<
" != " 1272 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1274 "MultiVectors must have the same number of columns (vectors). " 1275 "this->getNumVectors() = " << numVecs <<
" != " 1277 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1278 numDots != numVecs, std::runtime_error,
1279 "The output array 'dots' must have the same number of entries as the " 1280 "number of columns (vectors) in *this and A. dots.dimension_0() = " <<
1281 numDots <<
" != this->getNumVectors() = " << numVecs <<
".");
1283 const std::pair<size_t, size_t> colRng (0, numVecs);
1284 RV dotsOut = subview (dots, colRng);
1285 RCP<const Comm<int> > comm = this->
getMap ().is_null () ? null :
1286 this->
getMap ()->getComm ();
1290 const bool oneMemorySpace =
1291 Kokkos::Impl::is_same<
typename dual_view_type::t_dev::memory_space,
1292 typename dual_view_type::t_host::memory_space>::value;
1293 if (! oneMemorySpace && A.
view_.modified_host() > A.
view_.modified_device()) {
1296 typedef typename dual_view_type::t_host XMV;
1299 this->
view_.template sync<typename XMV::memory_space> ();
1300 lclDotImpl<RV, XMV> (dotsOut,
view_.h_view, A.
view_.h_view,
1301 lclNumRows, numVecs,
1304 typename RV::HostMirror dotsOutHost = create_mirror_view (dotsOut);
1306 gblDotImpl<typename RV::HostMirror> (dotsOutHost, comm,
1312 typedef typename dual_view_type::t_dev XMV;
1315 this->
view_.template sync<typename XMV::memory_space> ();
1316 lclDotImpl<RV, XMV> (dotsOut,
view_.d_view, A.
view_.d_view,
1317 lclNumRows, numVecs,
1325 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1329 const Teuchos::ArrayView<dot_type>& dots)
const 1331 typedef Kokkos::View<dot_type*, device_type> dev_dots_view_type;
1332 typedef MakeUnmanagedView<dot_type, device_type> view_getter_type;
1333 typedef typename view_getter_type::view_type host_dots_view_type;
1335 const size_t numDots =
static_cast<size_t> (dots.size ());
1336 host_dots_view_type dotsHostView (dots.getRawPtr (), numDots);
1337 dev_dots_view_type dotsDevView (
"MV::dot tmp", numDots);
1338 this->
dot (A, dotsDevView);
1343 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1346 norm2 (
const Teuchos::ArrayView<mag_type>& norms)
const 1348 typedef Kokkos::View<mag_type*, device_type> dev_norms_view_type;
1349 typedef MakeUnmanagedView<mag_type, device_type> view_getter_type;
1350 typedef typename view_getter_type::view_type host_norms_view_type;
1352 const size_t numNorms =
static_cast<size_t> (norms.size ());
1353 host_norms_view_type normsHostView (norms.getRawPtr (), numNorms);
1354 dev_norms_view_type normsDevView (
"MV::norm2 tmp", numNorms);
1355 this->
norm2 (normsDevView);
1360 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1363 norm2 (
const Kokkos::View<mag_type*, device_type>& norms)
const 1369 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1370 void TPETRA_DEPRECATED
1373 const Teuchos::ArrayView<mag_type> &norms)
const 1376 using Kokkos::subview;
1377 using Teuchos::Comm;
1378 using Teuchos::null;
1380 using Teuchos::reduceAll;
1381 using Teuchos::REDUCE_SUM;
1382 typedef Kokkos::Details::ArithTraits<impl_scalar_type> ATS;
1383 typedef Kokkos::Details::ArithTraits<mag_type> ATM;
1384 typedef Kokkos::View<mag_type*, device_type> norms_view_type;
1385 const char tfecfFuncName[] =
"normWeighted: ";
1388 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1389 static_cast<size_t> (norms.size ()) != numVecs, std::runtime_error,
1390 "norms.size() = " << norms.size () <<
" != this->getNumVectors() = " 1394 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1395 ! OneW && weights.
getNumVectors () != numVecs, std::runtime_error,
1396 "The input MultiVector of weights must contain either one column, " 1397 "or must have the same number of columns as *this. " 1399 <<
" and this->getNumVectors() = " << numVecs <<
".");
1401 #ifdef HAVE_TPETRA_DEBUG 1402 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1403 ! this->
getMap ()->isCompatible (*weights.
getMap ()), std::runtime_error,
1404 "MultiVectors do not have compatible Maps:" << std::endl
1405 <<
"this->getMap(): " << std::endl << *this->
getMap()
1406 <<
"weights.getMap(): " << std::endl << *weights.
getMap() << std::endl);
1409 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
1411 "MultiVectors do not have the same local length.");
1412 #endif // HAVE_TPETRA_DEBUG 1414 norms_view_type lclNrms (
"lclNrms", numVecs);
1416 view_.template sync<device_type> ();
1417 weights.
view_.template sync<device_type> ();
1419 typename dual_view_type::t_dev X_lcl = this->
view_.d_view;
1420 typename dual_view_type::t_dev W_lcl = weights.
view_.d_view;
1423 KokkosBlas::nrm2w_squared (lclNrms, X_lcl, W_lcl);
1426 for (
size_t j = 0; j < numVecs; ++j) {
1429 const size_t W_col = OneW ?
static_cast<size_t> (0) :
1431 KokkosBlas::nrm2w_squared (subview (lclNrms, j),
1432 subview (X_lcl, ALL (), X_col),
1433 subview (W_lcl, ALL (), W_col));
1439 RCP<const Comm<int> > comm = this->
getMap ().is_null () ?
1440 Teuchos::null : this->
getMap ()->getComm ();
1444 reduceAll<int, mag_type> (*comm, REDUCE_SUM,
static_cast<int> (numVecs),
1445 lclNrms.ptr_on_device (), norms.getRawPtr ());
1446 for (
size_t k = 0; k < numVecs; ++k) {
1447 norms[k] = ATM::sqrt (norms[k] * OneOverN);
1451 typename norms_view_type::HostMirror lclNrms_h =
1452 Kokkos::create_mirror_view (lclNrms);
1454 for (
size_t k = 0; k < numVecs; ++k) {
1455 norms[k] = ATM::sqrt (ATS::magnitude (lclNrms_h(k)) * OneOverN);
1461 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1464 norm1 (
const Teuchos::ArrayView<mag_type>& norms)
const 1466 typedef Kokkos::View<mag_type*, device_type> dev_norms_view_type;
1467 typedef MakeUnmanagedView<mag_type, device_type> view_getter_type;
1468 typedef typename view_getter_type::view_type host_norms_view_type;
1470 const size_t numNorms =
static_cast<size_t> (norms.size ());
1471 host_norms_view_type normsHostView (norms.getRawPtr (), numNorms);
1472 dev_norms_view_type normsDevView (
"MV::norm1 tmp", numNorms);
1473 this->
norm1 (normsDevView);
1478 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1481 norm1 (
const Kokkos::View<mag_type*, device_type>& norms)
const 1487 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1490 normInf (
const Teuchos::ArrayView<mag_type>& norms)
const 1492 typedef Kokkos::View<mag_type*, device_type> dev_norms_view_type;
1493 typedef MakeUnmanagedView<mag_type, device_type> view_getter_type;
1494 typedef typename view_getter_type::view_type host_norms_view_type;
1496 const size_t numNorms =
static_cast<size_t> (norms.size ());
1497 host_norms_view_type normsHostView (norms.getRawPtr (), numNorms);
1498 dev_norms_view_type normsDevView (
"MV::normInf tmp", numNorms);
1504 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1507 normInf (
const Kokkos::View<mag_type*, device_type>& norms)
const 1521 template<
class RV,
class XMV>
1523 lclNormImpl (
const RV& normsOut,
1525 const size_t lclNumRows,
1526 const size_t numVecs,
1527 const Teuchos::ArrayView<const size_t>& whichVecs,
1528 const bool constantStride,
1532 using Kokkos::subview;
1533 typedef typename RV::non_const_value_type
mag_type;
1535 static_assert (Kokkos::Impl::is_view<RV>::value,
1536 "Tpetra::MultiVector::lclNormImpl: " 1537 "The first argument RV is not a Kokkos::View.");
1538 static_assert (RV::rank == 1,
"Tpetra::MultiVector::lclNormImpl: " 1539 "The first argument normsOut must have rank 1.");
1540 static_assert (Kokkos::Impl::is_view<XMV>::value,
1541 "Tpetra::MultiVector::lclNormImpl: " 1542 "The second argument X_lcl is not a Kokkos::View.");
1543 static_assert (XMV::rank == 2,
"Tpetra::MultiVector::lclNormImpl: " 1544 "The second argument X_lcl must have rank 2.");
1549 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
1550 const std::pair<size_t, size_t> colRng (0, numVecs);
1551 RV theNorms = subview (normsOut, colRng);
1552 XMV X = subview (X_lcl, rowRng, Kokkos::ALL());
1557 TEUCHOS_TEST_FOR_EXCEPTION(
1558 lclNumRows != 0 && constantStride && ( \
1559 ( X.dimension_0 () != lclNumRows ) ||
1560 ( X.dimension_1 () != numVecs ) ),
1561 std::logic_error,
"Constant Stride X's dimensions are " << X.dimension_0 () <<
" x " 1562 << X.dimension_1 () <<
", which differ from the local dimensions " 1563 << lclNumRows <<
" x " << numVecs <<
". Please report this bug to " 1564 "the Tpetra developers.");
1566 TEUCHOS_TEST_FOR_EXCEPTION(
1567 lclNumRows != 0 && !constantStride && ( \
1568 ( X.dimension_0 () != lclNumRows ) ||
1569 ( X.dimension_1 () < numVecs ) ),
1570 std::logic_error,
"Strided X's dimensions are " << X.dimension_0 () <<
" x " 1571 << X.dimension_1 () <<
", which are incompatible with the local dimensions " 1572 << lclNumRows <<
" x " << numVecs <<
". Please report this bug to " 1573 "the Tpetra developers.");
1575 if (lclNumRows == 0) {
1576 const mag_type zeroMag = Kokkos::Details::ArithTraits<mag_type>::zero ();
1580 if (constantStride) {
1581 if (whichNorm == IMPL_NORM_INF) {
1582 KokkosBlas::nrmInf (theNorms, X);
1584 else if (whichNorm == IMPL_NORM_ONE) {
1585 KokkosBlas::nrm1 (theNorms, X);
1587 else if (whichNorm == IMPL_NORM_TWO) {
1588 KokkosBlas::nrm2_squared (theNorms, X);
1591 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
1600 for (
size_t k = 0; k < numVecs; ++k) {
1601 const size_t X_col = constantStride ? k : whichVecs[k];
1602 if (whichNorm == IMPL_NORM_INF) {
1603 KokkosBlas::nrmInf (subview (theNorms, k),
1604 subview (X, ALL (), X_col));
1606 else if (whichNorm == IMPL_NORM_ONE) {
1607 KokkosBlas::nrm1 (subview (theNorms, k),
1608 subview (X, ALL (), X_col));
1610 else if (whichNorm == IMPL_NORM_TWO) {
1611 KokkosBlas::nrm2_squared (subview (theNorms, k),
1612 subview (X, ALL (), X_col));
1615 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
1624 gblNormImpl (
const RV& normsOut,
1625 const Teuchos::RCP<
const Teuchos::Comm<int> >& comm,
1626 const bool distributed,
1629 using Teuchos::REDUCE_MAX;
1630 using Teuchos::REDUCE_SUM;
1631 using Teuchos::reduceAll;
1632 typedef typename RV::non_const_value_type
mag_type;
1634 const size_t numVecs = normsOut.dimension_0 ();
1649 if (distributed && ! comm.is_null ()) {
1655 RV lclNorms (
"MV::normImpl lcl", numVecs);
1657 const mag_type*
const lclSum = lclNorms.ptr_on_device ();
1658 mag_type*
const gblSum = normsOut.ptr_on_device ();
1659 const int nv =
static_cast<int> (numVecs);
1660 if (whichNorm == IMPL_NORM_INF) {
1661 reduceAll<int, mag_type> (*comm, REDUCE_MAX, nv, lclSum, gblSum);
1663 reduceAll<int, mag_type> (*comm, REDUCE_SUM, nv, lclSum, gblSum);
1667 if (whichNorm == IMPL_NORM_TWO) {
1673 const bool inHostMemory =
1674 Kokkos::Impl::is_same<
typename RV::memory_space,
1675 typename RV::host_mirror_space::memory_space>::value;
1677 for (
size_t j = 0; j < numVecs; ++j) {
1678 normsOut(j) = Kokkos::Details::ArithTraits<mag_type>::sqrt (normsOut(j));
1686 KokkosBlas::Impl::SquareRootFunctor<RV> f (normsOut);
1687 Kokkos::parallel_for (numVecs, f);
1694 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1697 normImpl (
const Kokkos::View<mag_type*, device_type>& norms,
1700 using Kokkos::create_mirror_view;
1701 using Kokkos::subview;
1702 using Teuchos::Comm;
1703 using Teuchos::null;
1706 typedef Kokkos::View<mag_type*, device_type> RV;
1713 const size_t numNorms =
static_cast<size_t> (norms.dimension_0 ());
1714 TEUCHOS_TEST_FOR_EXCEPTION(
1715 numNorms < numVecs, std::runtime_error,
"Tpetra::MultiVector::normImpl: " 1716 "'norms' must have at least as many entries as the number of vectors in " 1717 "*this. norms.dimension_0() = " << numNorms <<
" < this->getNumVectors()" 1718 " = " << numVecs <<
".");
1720 const std::pair<size_t, size_t> colRng (0, numVecs);
1721 RV normsOut = subview (norms, colRng);
1723 EWhichNormImpl lclNormType;
1724 if (whichNorm == NORM_ONE) {
1725 lclNormType = IMPL_NORM_ONE;
1726 }
else if (whichNorm == NORM_TWO) {
1727 lclNormType = IMPL_NORM_TWO;
1729 lclNormType = IMPL_NORM_INF;
1732 RCP<const Comm<int> > comm = this->
getMap ().is_null () ? null :
1733 this->
getMap ()->getComm ();
1737 const bool oneMemorySpace =
1738 Kokkos::Impl::is_same<
typename dual_view_type::t_dev::memory_space,
1739 typename dual_view_type::t_host::memory_space>::value;
1740 if (! oneMemorySpace &&
view_.modified_host() >
view_.modified_device()) {
1743 typedef typename dual_view_type::t_host XMV;
1744 lclNormImpl<RV, XMV> (normsOut,
view_.h_view, lclNumRows, numVecs,
1747 typename RV::HostMirror normsOutHost = create_mirror_view (normsOut);
1749 gblNormImpl<typename RV::HostMirror> (normsOutHost, comm,
1756 typedef typename dual_view_type::t_dev XMV;
1757 lclNormImpl<RV, XMV> (normsOut,
view_.d_view, lclNumRows, numVecs,
1760 gblNormImpl<RV> (normsOut, comm, this->
isDistributed (), lclNormType);
1764 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1767 meanValue (
const Teuchos::ArrayView<impl_scalar_type>& means)
const 1771 using Kokkos::subview;
1772 using Teuchos::Comm;
1774 using Teuchos::reduceAll;
1775 using Teuchos::REDUCE_SUM;
1776 typedef Kokkos::Details::ArithTraits<impl_scalar_type> ATS;
1780 const size_t numMeans =
static_cast<size_t> (means.size ());
1782 TEUCHOS_TEST_FOR_EXCEPTION(
1783 numMeans != numVecs, std::runtime_error,
1784 "Tpetra::MultiVector::meanValue: means.size() = " << numMeans
1785 <<
" != this->getNumVectors() = " << numVecs <<
".");
1787 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
1788 const std::pair<size_t, size_t> colRng (0, numVecs);
1793 typedef Kokkos::View<impl_scalar_type*, device_type> local_view_type;
1795 typename local_view_type::HostMirror::array_layout,
1797 Kokkos::MemoryTraits<Kokkos::Unmanaged> > host_local_view_type;
1798 host_local_view_type meansOut (means.getRawPtr (), numMeans);
1800 RCP<const Comm<int> > comm = this->
getMap ().is_null () ? Teuchos::null :
1801 this->
getMap ()->getComm ();
1805 const bool oneMemorySpace =
1806 Kokkos::Impl::is_same<
typename dual_view_type::t_dev::memory_space,
1807 typename dual_view_type::t_host::memory_space>::value;
1808 if (! oneMemorySpace &&
view_.modified_host() >
view_.modified_device()) {
1810 typename dual_view_type::t_host X_lcl =
1811 subview (this->
view_.h_view, rowRng, Kokkos::ALL());
1814 typename local_view_type::HostMirror lclSums (
"MV::meanValue tmp", numVecs);
1816 KokkosBlas::sum (lclSums, X_lcl);
1819 for (
size_t j = 0; j < numVecs; ++j) {
1821 KokkosBlas::sum (subview (lclSums, j), subview (X_lcl, ALL (), col));
1829 reduceAll (*comm, REDUCE_SUM, static_cast<int> (numVecs),
1830 lclSums.ptr_on_device (), meansOut.ptr_on_device ());
1838 typename dual_view_type::t_dev X_lcl =
1839 subview (this->
view_.d_view, rowRng, Kokkos::ALL());
1842 local_view_type lclSums (
"MV::meanValue tmp", numVecs);
1844 KokkosBlas::sum (lclSums, X_lcl);
1847 for (
size_t j = 0; j < numVecs; ++j) {
1849 KokkosBlas::sum (subview (lclSums, j), subview (X_lcl, ALL (), col));
1858 reduceAll (*comm, REDUCE_SUM, static_cast<int> (numVecs),
1859 lclSums.ptr_on_device (), meansOut.ptr_on_device ());
1869 const impl_scalar_type OneOverN =
1871 for (
size_t k = 0; k < numMeans; ++k) {
1872 meansOut(k) = meansOut(k) * OneOverN;
1877 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1883 using Kokkos::subview;
1885 typedef Kokkos::Details::ArithTraits<IST> ATS;
1886 typedef Kokkos::Random_XorShift64_Pool<typename device_type::execution_space> pool_type;
1887 typedef typename pool_type::generator_type generator_type;
1898 const uint64_t myRank =
1899 static_cast<uint64_t
> (this->
getMap ()->getComm ()->getRank ());
1900 uint64_t seed64 =
static_cast<uint64_t
> (std::rand ()) + myRank + 17311uLL;
1901 unsigned int seed =
static_cast<unsigned int> (seed64&0xffffffff);
1903 pool_type rand_pool (seed);
1904 const IST max = Kokkos::rand<generator_type, IST>::max ();
1905 const IST min = ATS::is_signed ? IST (-max) : ATS::zero ();
1908 Kokkos::fill_random (
view_.d_view, rand_pool, min, max);
1909 view_.template modify<device_type> ();
1913 view_.template sync<device_type> ();
1914 typedef Kokkos::View<IST*, device_type> view_type;
1915 for (
size_t k = 0; k < numVecs; ++k) {
1917 view_type X_k = subview (
view_.d_view, ALL (), col);
1918 Kokkos::fill_random (X_k, rand_pool, min, max);
1920 view_.template modify<device_type> ();
1925 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
1932 using Kokkos::subview;
1933 typedef typename dual_view_type::t_dev::device_type DMS;
1934 typedef typename dual_view_type::t_host::device_type HMS;
1939 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
1940 const std::pair<size_t, size_t> colRng (0, numVecs);
1944 if (
view_.modified_device() >=
view_.modified_host()) {
1949 typedef typename dual_view_type::t_dev mv_view_type;
1952 typename mv_view_type::array_layout, DMS> vec_view_type;
1954 this->
template modify<DMS> ();
1956 subview (this->
getDualView ().
template view<DMS> (),
1957 rowRng, Kokkos::ALL());
1960 subview (X, ALL (), static_cast<size_t> (0));
1967 for (
size_t k = 0; k < numVecs; ++k) {
1969 vec_view_type X_k = subview (X, ALL (), col);
1975 typedef typename dual_view_type::t_host mv_view_type;
1977 typename mv_view_type::array_layout, HMS> vec_view_type;
1979 this->
template modify<HMS> ();
1981 subview (this->
getDualView ().
template view<HMS> (),
1982 rowRng, Kokkos::ALL());
1985 subview (X, ALL (), static_cast<size_t> (0));
1992 for (
size_t k = 0; k < numVecs; ++k) {
1994 vec_view_type X_k = subview (X, ALL (), col);
2002 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2007 using Teuchos::ArrayRCP;
2008 using Teuchos::Comm;
2015 TEUCHOS_TEST_FOR_EXCEPTION(
2017 "Tpetra::MultiVector::replaceMap: This method does not currently work " 2018 "if the MultiVector is a column view of another MultiVector (that is, if " 2019 "isConstantStride() == false).");
2049 #ifdef HAVE_TEUCHOS_DEBUG 2063 #endif // HAVE_TEUCHOS_DEBUG 2065 if (this->
getMap ().is_null ()) {
2070 TEUCHOS_TEST_FOR_EXCEPTION(
2071 newMap.is_null (), std::invalid_argument,
2072 "Tpetra::MultiVector::replaceMap: both current and new Maps are null. " 2073 "This probably means that the input Map is incorrect.");
2077 const size_t newNumRows = newMap->getNodeNumElements ();
2078 const size_t origNumRows =
view_.dimension_0 ();
2081 if (origNumRows != newNumRows ||
view_.dimension_1 () != numCols) {
2082 view_ = allocDualView<Scalar, LocalOrdinal, GlobalOrdinal, Node> (newNumRows, numCols);
2085 else if (newMap.is_null ()) {
2088 const size_t newNumRows =
static_cast<size_t> (0);
2090 view_ = allocDualView<Scalar, LocalOrdinal, GlobalOrdinal, Node> (newNumRows, numCols);
2093 this->
map_ = newMap;
2096 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2102 using Kokkos::subview;
2105 if (theAlpha == Kokkos::Details::ArithTraits<impl_scalar_type>::one ()) {
2110 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
2111 const std::pair<size_t, size_t> colRng (0, numVecs);
2113 typedef typename dual_view_type::t_dev dev_view_type;
2114 typedef typename dual_view_type::t_host host_view_type;
2123 const bool oneMemorySpace =
2124 Kokkos::Impl::is_same<
typename dev_view_type::memory_space,
2125 typename host_view_type::memory_space>::value;
2126 if (! oneMemorySpace &&
view_.modified_host() >
view_.modified_device()) {
2127 auto Y_lcl = subview (this->
view_.h_view, rowRng, Kokkos::ALL());
2130 KokkosBlas::scal (Y_lcl, theAlpha, Y_lcl);
2133 for (
size_t k = 0; k < numVecs; ++k) {
2135 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2136 KokkosBlas::scal (Y_k, theAlpha, Y_k);
2141 auto Y_lcl = subview (this->
view_.d_view, rowRng, Kokkos::ALL());
2144 KokkosBlas::scal (Y_lcl, theAlpha, Y_lcl);
2147 for (
size_t k = 0; k < numVecs; ++k) {
2149 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2150 KokkosBlas::scal (Y_k, theAlpha, Y_k);
2157 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2160 scale (
const Teuchos::ArrayView<const Scalar>& alphas)
2163 const size_t numAlphas =
static_cast<size_t> (alphas.size ());
2164 TEUCHOS_TEST_FOR_EXCEPTION(
2165 numAlphas != numVecs, std::invalid_argument,
"Tpetra::MultiVector::" 2166 "scale: alphas.size() = " << numAlphas <<
" != this->getNumVectors() = " 2170 typedef Kokkos::DualView<impl_scalar_type*, device_type> k_alphas_type ;
2171 k_alphas_type k_alphas (
"alphas::tmp", numAlphas);
2172 k_alphas.template modify<typename k_alphas_type::host_mirror_space> ();
2173 for (
size_t i = 0; i < numAlphas; ++i) {
2176 k_alphas.template sync<typename k_alphas_type::memory_space> ();
2178 this->
scale (k_alphas.d_view);
2181 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2184 scale (
const Kokkos::View<const impl_scalar_type*, device_type>& alphas)
2187 using Kokkos::subview;
2191 TEUCHOS_TEST_FOR_EXCEPTION(
2192 static_cast<size_t> (alphas.dimension_0 ()) != numVecs,
2193 std::invalid_argument,
"Tpetra::MultiVector::scale(alphas): " 2194 "alphas.dimension_0() = " << alphas.dimension_0 ()
2195 <<
" != this->getNumVectors () = " << numVecs <<
".");
2196 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
2197 const std::pair<size_t, size_t> colRng (0, numVecs);
2199 typedef typename dual_view_type::t_dev dev_view_type;
2200 typedef typename dual_view_type::t_host host_view_type;
2208 const bool oneMemorySpace =
2209 Kokkos::Impl::is_same<
typename dev_view_type::memory_space,
2210 typename host_view_type::memory_space>::value;
2211 if (! oneMemorySpace &&
2212 this->
view_.modified_host() > this->
view_.modified_device()) {
2217 typename input_view_type::HostMirror alphas_h =
2218 Kokkos::create_mirror_view (alphas);
2221 auto Y_lcl = subview (this->
view_.h_view, rowRng, Kokkos::ALL());
2224 KokkosBlas::scal (Y_lcl, alphas_h, Y_lcl);
2227 for (
size_t k = 0; k < numVecs; ++k) {
2229 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2232 KokkosBlas::scal (Y_k, alphas_h(k), Y_k);
2237 auto Y_lcl = subview (this->
view_.d_view, rowRng, Kokkos::ALL());
2240 KokkosBlas::scal (Y_lcl, alphas, Y_lcl);
2243 for (
size_t k = 0; k < numVecs; ++k) {
2245 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2251 KokkosBlas::scal (Y_k, alphas(k), Y_k);
2257 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2264 using Kokkos::subview;
2265 const char tfecfFuncName[] =
"scale: ";
2270 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2272 "this->getLocalLength() = " << lclNumRows <<
" != A.getLocalLength() = " 2274 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2276 "this->getNumVectors() = " << numVecs <<
" != A.getNumVectors() = " 2280 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
2281 const std::pair<size_t, size_t> colRng (0, numVecs);
2283 typedef typename dual_view_type::t_dev dev_view_type;
2284 typedef typename dual_view_type::t_host host_view_type;
2288 const bool oneMemorySpace =
2289 Kokkos::Impl::is_same<
typename dev_view_type::memory_space,
2290 typename host_view_type::memory_space>::value;
2291 if (! oneMemorySpace && A.
view_.modified_host() > A.
view_.modified_device()) {
2295 this->
view_.template sync<typename host_view_type::memory_space> ();
2296 this->
view_.template modify<typename host_view_type::memory_space> ();
2297 auto Y_lcl = subview (this->
view_.h_view, rowRng, Kokkos::ALL());
2298 auto X_lcl = subview (A.
view_.h_view, rowRng, Kokkos::ALL());
2301 KokkosBlas::scal (Y_lcl, theAlpha, X_lcl);
2305 for (
size_t k = 0; k < numVecs; ++k) {
2308 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2309 auto X_k = subview (X_lcl, ALL (), X_col);
2311 KokkosBlas::scal (Y_k, theAlpha, X_k);
2318 this->
view_.template sync<typename dev_view_type::memory_space> ();
2319 this->
view_.template modify<typename dev_view_type::memory_space> ();
2320 auto Y_lcl = subview (this->
view_.d_view, rowRng, Kokkos::ALL());
2321 auto X_lcl = subview (A.
view_.d_view, rowRng, Kokkos::ALL());
2324 KokkosBlas::scal (Y_lcl, theAlpha, X_lcl);
2328 for (
size_t k = 0; k < numVecs; ++k) {
2331 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2332 auto X_k = subview (X_lcl, ALL (), X_col);
2334 KokkosBlas::scal (Y_k, theAlpha, X_k);
2341 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2346 const char tfecfFuncName[] =
"reciprocal";
2348 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2350 ": MultiVectors do not have the same local length. " 2353 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2355 ": MultiVectors do not have the same number of columns (vectors). " 2357 <<
" != A.getNumVectors() = " << A.
getNumVectors () <<
".");
2364 view_.template sync<device_type> ();
2365 view_.template modify<device_type> ();
2366 KokkosBlas::reciprocal (
view_.d_view, A.
view_.d_view);
2370 using Kokkos::subview;
2371 typedef Kokkos::View<impl_scalar_type*, device_type> view_type;
2373 view_.template sync<device_type> ();
2374 view_.template modify<device_type> ();
2378 A.
view_.template sync<device_type> ();
2379 A.
view_.template modify<device_type> ();
2381 for (
size_t k = 0; k < numVecs; ++k) {
2383 view_type vector_k = subview (
view_.d_view, ALL (), this_col);
2385 view_type vector_Ak = subview (A.
view_.d_view, ALL (), A_col);
2386 KokkosBlas::reciprocal(vector_k, vector_Ak);
2390 catch (std::runtime_error &e) {
2391 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
true, std::runtime_error,
2392 ": Caught exception from Kokkos: " << e.what () << std::endl);
2397 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2402 const char tfecfFuncName[] =
"abs";
2403 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2405 ": MultiVectors do not have the same local length. " 2408 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2410 ": MultiVectors do not have the same number of columns (vectors). " 2412 <<
" != A.getNumVectors() = " << A.
getNumVectors () <<
".");
2418 view_.template sync<device_type>();
2419 view_.template modify<device_type>();
2420 KokkosBlas::abs (
view_.d_view, A.
view_.d_view);
2424 using Kokkos::subview;
2425 typedef Kokkos::View<impl_scalar_type*, device_type> view_type;
2427 view_.template sync<device_type> ();
2428 view_.template modify<device_type> ();
2429 A.
view_.template sync<device_type> ();
2430 A.
view_.template modify<device_type> ();
2432 for (
size_t k=0; k < numVecs; ++k) {
2434 view_type vector_k = subview (
view_.d_view, ALL (), this_col);
2436 view_type vector_Ak = subview (A.
view_.d_view, ALL (), A_col);
2437 KokkosBlas::abs (vector_k, vector_Ak);
2443 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2451 using Kokkos::subview;
2452 const char tfecfFuncName[] =
"update: ";
2457 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2459 "this->getLocalLength() = " << lclNumRows <<
" != A.getLocalLength() = " 2461 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2463 "this->getNumVectors() = " << numVecs <<
" != A.getNumVectors() = " 2468 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
2469 const std::pair<size_t, size_t> colRng (0, numVecs);
2471 typedef typename dual_view_type::t_dev dev_view_type;
2472 typedef typename dual_view_type::t_host host_view_type;
2476 const bool oneMemorySpace =
2477 Kokkos::Impl::is_same<
typename dev_view_type::memory_space,
2478 typename host_view_type::memory_space>::value;
2479 if (! oneMemorySpace && A.
view_.modified_host() > A.
view_.modified_device()) {
2483 this->
view_.template sync<typename host_view_type::memory_space> ();
2484 this->
view_.template modify<typename host_view_type::memory_space> ();
2485 auto Y_lcl = subview (this->
view_.h_view, rowRng, Kokkos::ALL());
2486 auto X_lcl = subview (A.
view_.h_view, rowRng, Kokkos::ALL());
2489 KokkosBlas::axpby (theAlpha, X_lcl, theBeta, Y_lcl);
2493 for (
size_t k = 0; k < numVecs; ++k) {
2496 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2497 auto X_k = subview (X_lcl, ALL (), X_col);
2499 KokkosBlas::axpby (theAlpha, X_k, theBeta, Y_k);
2506 this->
view_.template sync<typename dev_view_type::memory_space> ();
2507 this->
view_.template modify<typename dev_view_type::memory_space> ();
2508 auto Y_lcl = subview (this->
view_.d_view, rowRng, Kokkos::ALL());
2509 auto X_lcl = subview (A.
view_.d_view, rowRng, Kokkos::ALL());
2512 KokkosBlas::axpby (theAlpha, X_lcl, theBeta, Y_lcl);
2516 for (
size_t k = 0; k < numVecs; ++k) {
2519 auto Y_k = subview (Y_lcl, ALL (), Y_col);
2520 auto X_k = subview (X_lcl, ALL (), X_col);
2522 KokkosBlas::axpby (theAlpha, X_k, theBeta, Y_k);
2528 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2535 const Scalar& gamma)
2538 using Kokkos::subview;
2539 const char tfecfFuncName[] =
"update(alpha,A,beta,B,gamma): ";
2542 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2544 "The input MultiVector A has " << A.
getLocalLength () <<
" local " 2545 "row(s), but this MultiVector has " << lclNumRows <<
" local " 2547 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2549 "The input MultiVector B has " << B.
getLocalLength () <<
" local " 2550 "row(s), but this MultiVector has " << lclNumRows <<
" local " 2553 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2555 "The input MultiVector A has " << A.
getNumVectors () <<
" column(s), " 2556 "but this MultiVector has " << numVecs <<
" column(s).");
2557 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2559 "The input MultiVector B has " << B.
getNumVectors () <<
" column(s), " 2560 "but this MultiVector has " << numVecs <<
" column(s).");
2570 this->
view_.template sync<typename dual_view_type::t_dev::memory_space> ();
2571 A.
view_.template sync<typename dual_view_type::t_dev::memory_space> ();
2572 B.
view_.template sync<typename dual_view_type::t_dev::memory_space> ();
2575 this->
template modify<typename dual_view_type::t_dev::memory_space> ();
2577 const std::pair<size_t, size_t> rowRng (0, lclNumRows);
2578 const std::pair<size_t, size_t> colRng (0, numVecs);
2583 auto C_lcl = subview (this->
view_.d_view, rowRng, Kokkos::ALL());
2584 auto A_lcl = subview (A.
view_.d_view, rowRng, Kokkos::ALL());
2585 auto B_lcl = subview (B.
view_.d_view, rowRng, Kokkos::ALL());
2588 KokkosBlas::update (theAlpha, A_lcl, theBeta, B_lcl, theGamma, C_lcl);
2593 for (
size_t k = 0; k < numVecs; ++k) {
2597 KokkosBlas::update (theAlpha, subview (A_lcl, rowRng, A_col),
2598 theBeta, subview (B_lcl, rowRng, B_col),
2599 theGamma, subview (C_lcl, rowRng, this_col));
2604 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2605 Teuchos::ArrayRCP<const Scalar>
2610 using Kokkos::subview;
2611 typedef typename dual_view_type::host_mirror_space host_type;
2612 typedef typename dual_view_type::t_host host_view_type;
2618 view_.template sync<host_type> ();
2621 host_view_type hostView =
view_.template view<host_type> ();
2623 host_view_type hostView_j;
2626 const std::pair<size_t, size_t> colRng (colStart, colStart+1);
2627 hostView_j = subview (hostView, ALL (), colRng);
2630 Teuchos::ArrayRCP<const impl_scalar_type> dataAsArcp =
2631 Kokkos::Compat::persistingView (hostView_j, 0,
getLocalLength ());
2633 #ifdef HAVE_TPETRA_DEBUG 2634 TEUCHOS_TEST_FOR_EXCEPTION(
2635 static_cast<size_t>(hostView_j.dimension_0 ()) < static_cast<size_t>(dataAsArcp.size ()), std::logic_error,
2636 "Tpetra::MultiVector::getData: hostView_j.dimension_0() = " 2637 << hostView_j.dimension_0 () <<
" < dataAsArcp.size() = " 2638 << dataAsArcp.size () <<
". " 2639 "Please report this bug to the Tpetra developers.");
2640 #endif // HAVE_TPETRA_DEBUG 2642 return Teuchos::arcp_reinterpret_cast<
const Scalar> (dataAsArcp);
2645 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2646 Teuchos::ArrayRCP<Scalar>
2651 using Kokkos::subview;
2652 typedef typename dual_view_type::host_mirror_space host_type;
2653 typedef typename dual_view_type::t_host host_view_type;
2659 view_.template sync<host_type> ();
2662 host_view_type hostView =
view_.template view<host_type> ();
2664 host_view_type hostView_j;
2666 hostView_j = subview (hostView, ALL (), Kokkos::pair<int,int>(j,j+1));
2674 view_.template modify<host_type> ();
2677 Teuchos::ArrayRCP<impl_scalar_type> dataAsArcp =
2678 Kokkos::Compat::persistingView (hostView_j, 0,
getLocalLength ());
2680 #ifdef HAVE_TPETRA_DEBUG 2681 TEUCHOS_TEST_FOR_EXCEPTION(
2682 static_cast<size_t>(hostView_j.dimension_0 ()) < static_cast<size_t>(dataAsArcp.size ()), std::logic_error,
2683 "Tpetra::MultiVector::getDataNonConst: hostView_j.dimension_0() = " 2684 << hostView_j.dimension_0 () <<
" < dataAsArcp.size() = " 2685 << dataAsArcp.size () <<
". " 2686 "Please report this bug to the Tpetra developers.");
2687 #endif // HAVE_TPETRA_DEBUG 2689 return Teuchos::arcp_reinterpret_cast<Scalar> (dataAsArcp);
2692 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2697 if (
this != &source) {
2698 base_type::operator= (source);
2720 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2721 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2723 subCopy (
const Teuchos::ArrayView<const size_t>& cols)
const 2730 bool contiguous =
true;
2731 const size_t numCopyVecs =
static_cast<size_t> (cols.size ());
2732 for (
size_t j = 1; j < numCopyVecs; ++j) {
2733 if (cols[j] != cols[j-1] + static_cast<size_t> (1)) {
2738 if (contiguous && numCopyVecs > 0) {
2739 return this->
subCopy (Teuchos::Range1D (cols[0], cols[numCopyVecs-1]));
2742 RCP<const MV> X_sub = this->
subView (cols);
2743 RCP<MV> Y (
new MV (this->
getMap (), numCopyVecs,
false));
2749 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2750 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2757 RCP<const MV> X_sub = this->
subView (colRng);
2758 RCP<MV> Y (
new MV (this->
getMap (), static_cast<size_t> (colRng.size ()),
false));
2763 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2770 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2777 template <
class Scalar,
class LO,
class GO,
class Node, const
bool classic>
2781 const size_t offset) :
2785 using Kokkos::subview;
2789 const char prefix[] =
"Tpetra::MultiVector constructor (offsetView): ";
2794 const int myRank = X.
getMap ()->getComm ()->getRank ();
2795 TEUCHOS_TEST_FOR_EXCEPTION(
2797 prefix <<
"Invalid input Map. The input Map owns " << newNumRows <<
2798 " entries on process " << myRank <<
". offset = " << offset <<
". " 2800 " rows on this process.");
2803 #ifdef HAVE_TPETRA_DEBUG 2804 const size_t strideBefore =
2810 #endif // HAVE_TPETRA_DEBUG 2812 const std::pair<size_t, size_t> rowRng (offset, offset + newNumRows);
2816 dual_view_type newView = subview (X.
origView_, rowRng, ALL ());
2822 if (newView.dimension_0 () == 0 &&
2823 newView.dimension_1 () != X.
view_.dimension_1 ()) {
2824 newView = allocDualView<Scalar, LO, GO, Node> (size_t (0),
2832 #ifdef HAVE_TPETRA_DEBUG 2835 static_cast<size_t> (0);
2841 const size_t strideRet = subViewMV.isConstantStride () ?
2842 subViewMV.getStride () :
2843 static_cast<size_t> (0);
2844 const size_t lclNumRowsRet = subViewMV.getLocalLength ();
2845 const size_t numColsRet = subViewMV.getNumVectors ();
2847 const char suffix[] =
". This should never happen. Please report this " 2848 "bug to the Tpetra developers.";
2850 TEUCHOS_TEST_FOR_EXCEPTION(
2852 std::logic_error, prefix <<
"Returned MultiVector has a number of rows " 2853 "different than the number of local indices in the input Map. " 2854 "lclNumRowsRet: " << lclNumRowsRet <<
", subMap.getNodeNumElements(): " 2856 TEUCHOS_TEST_FOR_EXCEPTION(
2857 strideBefore != strideAfter || lclNumRowsBefore != lclNumRowsAfter ||
2858 numColsBefore != numColsAfter || hostPtrBefore != hostPtrAfter,
2859 std::logic_error, prefix <<
"Original MultiVector changed dimensions, " 2860 "stride, or host pointer after taking offset view. strideBefore: " <<
2861 strideBefore <<
", strideAfter: " << strideAfter <<
", lclNumRowsBefore: " 2862 << lclNumRowsBefore <<
", lclNumRowsAfter: " << lclNumRowsAfter <<
2863 ", numColsBefore: " << numColsBefore <<
", numColsAfter: " <<
2864 numColsAfter <<
", hostPtrBefore: " << hostPtrBefore <<
", hostPtrAfter: " 2865 << hostPtrAfter << suffix);
2866 TEUCHOS_TEST_FOR_EXCEPTION(
2867 strideBefore != strideRet, std::logic_error, prefix <<
"Returned " 2868 "MultiVector has different stride than original MultiVector. " 2869 "strideBefore: " << strideBefore <<
", strideRet: " << strideRet <<
2870 ", numColsBefore: " << numColsBefore <<
", numColsRet: " << numColsRet
2872 TEUCHOS_TEST_FOR_EXCEPTION(
2873 numColsBefore != numColsRet, std::logic_error,
2874 prefix <<
"Returned MultiVector has a different number of columns than " 2875 "original MultiVector. numColsBefore: " << numColsBefore <<
", " 2876 "numColsRet: " << numColsRet << suffix);
2877 #endif // HAVE_TPETRA_DEBUG 2882 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2883 Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2886 const size_t offset)
const 2889 return Teuchos::rcp (
new MV (*
this, *subMap, offset));
2892 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2893 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2896 const size_t offset)
2899 return Teuchos::rcp (
new MV (*
this, *subMap, offset));
2902 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2903 Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2905 subView (
const Teuchos::ArrayView<const size_t>& cols)
const 2907 using Teuchos::Array;
2911 const size_t numViewCols =
static_cast<size_t> (cols.size ());
2912 TEUCHOS_TEST_FOR_EXCEPTION(
2913 numViewCols < 1, std::runtime_error,
"Tpetra::MultiVector::subView" 2914 "(const Teuchos::ArrayView<const size_t>&): The input array cols must " 2915 "contain at least one entry, but cols.size() = " << cols.size ()
2920 bool contiguous =
true;
2921 for (
size_t j = 1; j < numViewCols; ++j) {
2922 if (cols[j] != cols[j-1] + static_cast<size_t> (1)) {
2928 if (numViewCols == 0) {
2930 return rcp (
new MV (this->
getMap (), numViewCols));
2933 return this->
subView (Teuchos::Range1D (cols[0], cols[numViewCols-1]));
2941 Array<size_t> newcols (cols.size ());
2942 for (
size_t j = 0; j < numViewCols; ++j) {
2950 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
2951 Teuchos::RCP<const MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
2956 using Kokkos::subview;
2957 using Teuchos::Array;
2961 const char tfecfFuncName[] =
"subView(Range1D): ";
2968 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2969 static_cast<size_t> (colRng.size ()) > numVecs, std::runtime_error,
2970 "colRng.size() = " << colRng.size () <<
" > this->getNumVectors() = " 2972 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
2973 numVecs != 0 && colRng.size () != 0 &&
2974 (colRng.lbound () <
static_cast<Teuchos::Ordinal
> (0) ||
2975 static_cast<size_t> (colRng.ubound ()) >= numVecs),
2976 std::invalid_argument,
"Nonempty input range [" << colRng.lbound () <<
2977 "," << colRng.ubound () <<
"] exceeds the valid range of column indices " 2978 "[0, " << numVecs <<
"].");
2980 RCP<const MV> X_ret;
2990 const std::pair<size_t, size_t> rows (0, lclNumRows);
2991 if (colRng.size () == 0) {
2992 const std::pair<size_t, size_t> cols (0, 0);
2993 dual_view_type X_sub = takeSubview (this->
view_, ALL (), cols);
2999 const std::pair<size_t, size_t> cols (colRng.lbound (),
3000 colRng.ubound () + 1);
3001 dual_view_type X_sub = takeSubview (this->
view_, ALL (), cols);
3005 if (static_cast<size_t> (colRng.size ()) == static_cast<size_t> (1)) {
3008 const std::pair<size_t, size_t> col (
whichVectors_[0] + colRng.lbound (),
3010 dual_view_type X_sub = takeSubview (
view_, ALL (), col);
3014 Array<size_t> which (
whichVectors_.begin () + colRng.lbound (),
3021 #ifdef HAVE_TPETRA_DEBUG 3022 using Teuchos::Comm;
3023 using Teuchos::outArg;
3024 using Teuchos::REDUCE_MIN;
3025 using Teuchos::reduceAll;
3027 RCP<const Comm<int> > comm = this->
getMap ().is_null () ? Teuchos::null :
3028 this->
getMap ()->getComm ();
3029 if (! comm.is_null ()) {
3033 if (X_ret.is_null ()) {
3036 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
3037 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3038 lclSuccess != 1, std::logic_error,
"X_ret (the subview of this " 3039 "MultiVector; the return value of this method) is null on some MPI " 3040 "process in this MultiVector's communicator. This should never " 3041 "happen. Please report this bug to the Tpetra developers.");
3043 if (! X_ret.is_null () &&
3044 X_ret->getNumVectors () !=
static_cast<size_t> (colRng.size ())) {
3047 reduceAll<int, int> (*comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
3048 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3049 lclSuccess != 1, std::logic_error,
3050 "X_ret->getNumVectors() != colRng.size(), on at least one MPI process " 3051 "in this MultiVector's communicator. This should never happen. " 3052 "Please report this bug to the Tpetra developers.");
3054 #endif // HAVE_TPETRA_DEBUG 3060 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3061 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
3066 return Teuchos::rcp_const_cast<MV> (this->
subView (cols));
3070 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3071 Teuchos::RCP<MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
3076 return Teuchos::rcp_const_cast<MV> (this->
subView (colRng));
3080 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3081 Teuchos::RCP<const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
3086 using Kokkos::subview;
3090 #ifdef HAVE_TPETRA_DEBUG 3091 const char tfecfFuncName[] =
"getVector(NonConst): ";
3092 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3093 this->vectorIndexOutOfRange (j), std::runtime_error,
"Input index j (== " 3094 << j <<
") exceeds valid range [0, " << this->
getNumVectors ()
3096 #endif // HAVE_TPETRA_DEBUG 3098 static_cast<size_t> (j) :
3100 const std::pair<size_t, size_t> rng (jj, jj+1);
3101 return rcp (
new V (this->
getMap (),
3102 takeSubview (this->
view_, ALL (), rng),
3107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3108 Teuchos::RCP<Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node, classic> >
3113 return Teuchos::rcp_const_cast<V> (this->
getVector (j));
3117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3120 get1dCopy (
const Teuchos::ArrayView<Scalar>& A,
const size_t LDA)
const 3122 using Kokkos::subview;
3124 typedef MakeUnmanagedView<IST, device_type> view_getter_type;
3125 typedef typename view_getter_type::view_type input_col_type;
3127 typedef typename dual_view_type::t_host host_view_type;
3128 typedef typename dual_view_type::t_dev dev_view_type;
3129 typedef Kokkos::View<IST*,
3130 typename host_view_type::array_layout,
3131 typename host_view_type::memory_space> host_col_type;
3132 typedef Kokkos::View<IST*,
3133 typename dev_view_type::array_layout,
3134 typename dev_view_type::memory_space> dev_col_type;
3135 const char tfecfFuncName[] =
"get1dCopy: ";
3139 const std::pair<size_t, size_t> rowRange (0, numRows);
3140 const std::pair<size_t, size_t> colRange (0, numCols);
3142 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3143 LDA < numRows, std::runtime_error,
3144 "LDA = " << LDA <<
" < numRows = " << numRows <<
".");
3145 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3146 numRows > static_cast<size_t> (0) &&
3147 numCols > static_cast<size_t> (0) &&
3148 static_cast<size_t> (A.size ()) < LDA * (numCols - 1) + numRows,
3150 "A.size() = " << A.size () <<
", but its size must be at least " 3151 << (LDA * (numCols - 1) + numRows) <<
" to hold all the entries.");
3161 for (
size_t j = 0; j < numCols; ++j) {
3162 const size_t srcCol =
3164 const size_t dstCol = j;
3165 IST*
const dstColRaw =
3166 reinterpret_cast<IST*
> (A.getRawPtr () + LDA * dstCol);
3167 input_col_type dstColView (dstColRaw, numRows);
3171 if (
view_.modified_host() >
view_.modified_device()) {
3172 host_col_type srcColView =
3173 subview (
view_.h_view, rowRange, srcCol);
3174 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3175 dstColView.dimension_0 () != srcColView.dimension_0 (),
3176 std::logic_error,
": srcColView and dstColView have different " 3177 "dimensions. Please report this bug to the Tpetra developers.");
3181 dev_col_type srcColView =
3182 subview (
view_.d_view, rowRange, srcCol);
3183 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3184 dstColView.dimension_0 () != srcColView.dimension_0 (),
3185 std::logic_error,
": srcColView and dstColView have different " 3186 "dimensions. Please report this bug to the Tpetra developers.");
3193 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3196 get2dCopy (
const Teuchos::ArrayView<
const Teuchos::ArrayView<Scalar> >& ArrayOfPtrs)
const 3199 const char tfecfFuncName[] =
"get2dCopy: ";
3203 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3204 static_cast<size_t> (ArrayOfPtrs.size ()) != numCols,
3205 std::runtime_error,
"Input array of pointers must contain as many " 3206 "entries (arrays) as the MultiVector has columns. ArrayOfPtrs.size() = " 3207 << ArrayOfPtrs.size () <<
" != getNumVectors() = " << numCols <<
".");
3209 if (numRows != 0 && numCols != 0) {
3211 for (
size_t j = 0; j < numCols; ++j) {
3212 const size_t dstLen =
static_cast<size_t> (ArrayOfPtrs[j].size ());
3213 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3214 dstLen < numRows, std::invalid_argument,
"Array j = " << j <<
" of " 3215 "the input array of arrays is not long enough to fit all entries in " 3216 "that column of the MultiVector. ArrayOfPtrs[j].size() = " << dstLen
3217 <<
" < getLocalLength() = " << numRows <<
".");
3221 for (
size_t j = 0; j < numCols; ++j) {
3223 const size_t LDA =
static_cast<size_t> (ArrayOfPtrs[j].size ());
3224 X_j->get1dCopy (ArrayOfPtrs[j], LDA);
3229 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3230 Teuchos::ArrayRCP<const Scalar>
3235 return Teuchos::null;
3237 TEUCHOS_TEST_FOR_EXCEPTION(
3239 "get1dView: This MultiVector does not have constant stride, so it is " 3240 "not possible to view its data as a single array. You may check " 3241 "whether a MultiVector has constant stride by calling " 3242 "isConstantStride().");
3246 typedef typename dual_view_type::host_mirror_space host_type;
3247 view_.template sync<host_type> ();
3250 Teuchos::ArrayRCP<const impl_scalar_type> dataAsArcp =
3251 Kokkos::Compat::persistingView (
view_.template view<host_type> ());
3252 return Teuchos::arcp_reinterpret_cast<
const Scalar> (dataAsArcp);
3256 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3257 Teuchos::ArrayRCP<Scalar>
3262 return Teuchos::null;
3264 TEUCHOS_TEST_FOR_EXCEPTION(
3266 "get1dViewNonConst: This MultiVector does not have constant stride, so " 3267 "it is not possible to view its data as a single array. You may check " 3268 "whether a MultiVector has constant stride by calling " 3269 "isConstantStride().");
3273 typedef typename dual_view_type::host_mirror_space host_type;
3274 view_.template sync<host_type> ();
3277 Teuchos::ArrayRCP<impl_scalar_type> dataAsArcp =
3278 Kokkos::Compat::persistingView (
view_.template view<host_type> ());
3279 return Teuchos::arcp_reinterpret_cast<Scalar> (dataAsArcp);
3283 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3284 Teuchos::ArrayRCP<Teuchos::ArrayRCP<Scalar> >
3288 using Teuchos::ArrayRCP;
3290 typename dual_view_type::array_layout,
device_type> col_dual_view_type;
3293 ArrayRCP<ArrayRCP<Scalar> > views (numCols);
3294 for (
size_t j = 0; j < numCols; ++j) {
3296 col_dual_view_type X_col =
3297 Kokkos::subview (
view_, Kokkos::ALL (), col);
3298 ArrayRCP<impl_scalar_type> X_col_arcp =
3299 Kokkos::Compat::persistingView (X_col.d_view);
3300 views[j] = Teuchos::arcp_reinterpret_cast<Scalar> (X_col_arcp);
3305 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3306 Teuchos::ArrayRCP<Teuchos::ArrayRCP<const Scalar> >
3310 using Teuchos::ArrayRCP;
3312 typename dual_view_type::array_layout,
device_type> col_dual_view_type;
3315 ArrayRCP<ArrayRCP<const Scalar> > views (numCols);
3316 for (
size_t j = 0; j < numCols; ++j) {
3318 col_dual_view_type X_col =
3319 Kokkos::subview (
view_, Kokkos::ALL (), col);
3320 ArrayRCP<const impl_scalar_type> X_col_arcp =
3321 Kokkos::Compat::persistingView (X_col.d_view);
3322 views[j] = Teuchos::arcp_reinterpret_cast<
const Scalar> (X_col_arcp);
3327 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3331 Teuchos::ETransp transB,
3332 const Scalar& alpha,
3337 using Teuchos::CONJ_TRANS;
3338 using Teuchos::NO_TRANS;
3339 using Teuchos::TRANS;
3342 using Teuchos::rcpFromRef;
3343 typedef Kokkos::Details::ArithTraits<impl_scalar_type> ATS;
3344 typedef Teuchos::ScalarTraits<Scalar> STS;
3346 const char errPrefix[] =
"Tpetra::MultiVector::multiply: ";
3374 TEUCHOS_TEST_FOR_EXCEPTION(
3375 ATS::is_complex && (transA == TRANS || transB == TRANS),
3376 std::invalid_argument, errPrefix <<
"Transpose without conjugation " 3377 "(transA == TRANS || transB == TRANS) is not supported for complex Scalar " 3380 transA = (transA == NO_TRANS ? NO_TRANS : CONJ_TRANS);
3381 transB = (transB == NO_TRANS ? NO_TRANS : CONJ_TRANS);
3391 TEUCHOS_TEST_FOR_EXCEPTION(
3393 A_ncols != B_nrows, std::runtime_error, errPrefix <<
"Dimensions of " 3394 "*this, op(A), and op(B) must be consistent. Local part of *this is " 3396 <<
", A is " << A_nrows <<
" x " << A_ncols
3397 <<
", and B is " << B_nrows <<
" x " << B_ncols <<
".");
3403 const bool Case1 = C_is_local && A_is_local && B_is_local;
3405 const bool Case2 = C_is_local && ! A_is_local && ! B_is_local &&
3406 transA == CONJ_TRANS && transB == NO_TRANS;
3408 const bool Case3 = ! C_is_local && ! A_is_local && B_is_local &&
3412 TEUCHOS_TEST_FOR_EXCEPTION(
3413 ! Case1 && ! Case2 && ! Case3, std::runtime_error, errPrefix
3414 <<
"Multiplication of op(A) and op(B) into *this is not a " 3415 "supported use case.");
3417 if (beta != STS::zero () && Case2) {
3423 const int myRank = this->
getMap ()->getComm ()->getRank ();
3425 beta_local = STS::zero ();
3435 C_tmp = rcp (
new MV (*
this, Teuchos::Copy));
3437 C_tmp = rcp (
this,
false);
3440 RCP<const MV> A_tmp;
3442 A_tmp = rcp (
new MV (A, Teuchos::Copy));
3444 A_tmp = rcpFromRef (A);
3447 RCP<const MV> B_tmp;
3449 B_tmp = rcp (
new MV (B, Teuchos::Copy));
3451 B_tmp = rcpFromRef (B);
3454 TEUCHOS_TEST_FOR_EXCEPTION(
3455 ! C_tmp->isConstantStride () || ! B_tmp->isConstantStride () ||
3456 ! A_tmp->isConstantStride (), std::logic_error, errPrefix
3457 <<
"Failed to make temporary constant-stride copies of MultiVectors.");
3462 A_tmp->getDualView ().d_view, B_tmp->getDualView ().d_view,
3463 beta_local, C_tmp->getDualView ().d_view);
3469 A_tmp = Teuchos::null;
3470 B_tmp = Teuchos::null;
3478 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3487 using Kokkos::subview;
3488 const char tfecfFuncName[] =
"elementWiseMultiply: ";
3491 #ifdef HAVE_TPETRA_DEBUG 3492 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3495 "MultiVectors do not have the same local length.");
3496 #endif // HAVE_TPETRA_DEBUG 3497 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
3498 numVecs != B.
getNumVectors (), std::runtime_error,
"this->getNumVectors" 3499 "() = " << numVecs <<
" != B.getNumVectors() = " << B.
getNumVectors ()
3508 view_.template sync<device_type> ();
3509 view_.template modify<device_type> ();
3510 A.
view_.template sync<device_type> ();
3511 B.
view_.template sync<device_type> ();
3512 KokkosBlas::mult (scalarThis,
view_.d_view, scalarAB,
3513 subview (A.
view_.d_view, ALL (), 0),
3517 view_.template sync<device_type> ();
3518 view_.template modify<device_type> ();
3519 A.
view_.template sync<device_type> ();
3520 B.
view_.template sync<device_type> ();
3522 for (
size_t j = 0; j < numVecs; ++j) {
3526 KokkosBlas::mult (scalarThis,
3527 subview (
view_.d_view, ALL (), C_col),
3529 subview (A.
view_.d_view, ALL (), 0),
3530 subview (B.
view_.d_view, ALL (), B_col));
3535 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3541 using Kokkos::subview;
3542 using Teuchos::reduceAll;
3543 using Teuchos::REDUCE_SUM;
3544 typedef typename dual_view_type::t_dev device_view_type;
3545 typedef typename dual_view_type::host_mirror_space host_mirror_space;
3547 TEUCHOS_TEST_FOR_EXCEPTION(
3549 "Tpetra::MultiVector::reduce() should only be called with locally " 3550 "replicated or otherwise not distributed MultiVector objects.");
3551 const Teuchos::Comm<int>& comm = * (this->
getMap ()->getComm ());
3552 if (comm.getSize () == 1) {
3564 TEUCHOS_TEST_FOR_EXCEPTION(
3565 numLclRows > static_cast<size_t> (std::numeric_limits<int>::max ()),
3566 std::runtime_error,
"Tpetra::MultiVector::reduce: On Process " <<
3567 comm.getRank () <<
", the number of local rows " << numLclRows <<
3568 " does not fit in int.");
3578 device_view_type srcBuf;
3580 srcBuf =
view_.d_view;
3583 srcBuf = device_view_type (
"srcBuf", numLclRows, numCols);
3592 device_view_type tgtBuf (
"tgtBuf", numLclRows, numCols);
3594 const int reduceCount =
static_cast<int> (numLclRows * numCols);
3595 reduceAll<int, impl_scalar_type> (comm, REDUCE_SUM, reduceCount,
3596 srcBuf.ptr_on_device (),
3597 tgtBuf.ptr_on_device ());
3600 view_.template modify<execution_space> ();
3602 const std::pair<size_t, size_t> lclRowRange (0, numLclRows);
3603 device_view_type d_view =
3604 subview (
view_.d_view, lclRowRange, ALL ());
3610 for (
size_t j = 0; j < numCols; ++j) {
3611 device_view_type d_view_j =
3612 subview (d_view, ALL (), std::pair<int,int>(j,j+1));
3613 device_view_type tgtBuf_j =
3614 subview (tgtBuf, ALL (), std::pair<int,int>(j,j+1));
3626 view_.template sync<host_mirror_space> ();
3630 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3637 #ifdef HAVE_TPETRA_DEBUG 3638 const LocalOrdinal minLocalIndex = this->
getMap()->getMinLocalIndex();
3639 const LocalOrdinal maxLocalIndex = this->
getMap()->getMaxLocalIndex();
3640 TEUCHOS_TEST_FOR_EXCEPTION(
3641 MyRow < minLocalIndex || MyRow > maxLocalIndex,
3643 "Tpetra::MultiVector::replaceLocalValue: row index " << MyRow
3644 <<
" is invalid. The range of valid row indices on this process " 3645 << this->
getMap()->getComm()->getRank() <<
" is [" << minLocalIndex
3646 <<
", " << maxLocalIndex <<
"].");
3647 TEUCHOS_TEST_FOR_EXCEPTION(
3648 vectorIndexOutOfRange(VectorIndex),
3650 "Tpetra::MultiVector::replaceLocalValue: vector index " << VectorIndex
3651 <<
" of the multivector is invalid.");
3655 view_.h_view (MyRow, colInd) = ScalarValue;
3659 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3665 const bool atomic)
const 3667 #ifdef HAVE_TPETRA_DEBUG 3668 const LocalOrdinal minLocalIndex = this->
getMap()->getMinLocalIndex();
3669 const LocalOrdinal maxLocalIndex = this->
getMap()->getMaxLocalIndex();
3670 TEUCHOS_TEST_FOR_EXCEPTION(
3671 localRow < minLocalIndex || localRow > maxLocalIndex,
3673 "Tpetra::MultiVector::sumIntoLocalValue: row index " << localRow
3674 <<
" is invalid. The range of valid row indices on this process " 3675 << this->
getMap()->getComm()->getRank() <<
" is [" << minLocalIndex
3676 <<
", " << maxLocalIndex <<
"].");
3677 TEUCHOS_TEST_FOR_EXCEPTION(
3678 vectorIndexOutOfRange(col),
3680 "Tpetra::MultiVector::sumIntoLocalValue: vector index " << col
3681 <<
" of the multivector is invalid.");
3685 Kokkos::atomic_add (& (
view_.h_view(localRow, colInd)), value);
3688 view_.h_view (localRow, colInd) += value;
3693 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3702 const LocalOrdinal MyRow = this->
map_->getLocalElement (GlobalRow);
3703 #ifdef HAVE_TPETRA_DEBUG 3704 TEUCHOS_TEST_FOR_EXCEPTION(
3705 MyRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid (),
3707 "Tpetra::MultiVector::replaceGlobalValue: Global row index " << GlobalRow
3708 <<
"is not present on this process " 3709 << this->
getMap ()->getComm ()->getRank () <<
".");
3710 TEUCHOS_TEST_FOR_EXCEPTION(
3711 vectorIndexOutOfRange (VectorIndex), std::runtime_error,
3712 "Tpetra::MultiVector::replaceGlobalValue: Vector index " << VectorIndex
3713 <<
" of the multivector is invalid.");
3714 #endif // HAVE_TPETRA_DEBUG 3718 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3724 const bool atomic)
const 3728 const LocalOrdinal lclRow = this->
map_->getLocalElement (globalRow);
3729 #ifdef HAVE_TEUCHOS_DEBUG 3730 TEUCHOS_TEST_FOR_EXCEPTION(
3731 lclRow == Teuchos::OrdinalTraits<LocalOrdinal>::invalid (),
3733 "Tpetra::MultiVector::sumIntoGlobalValue: Global row index " << globalRow
3734 <<
"is not present on this process " 3735 << this->
getMap ()->getComm ()->getRank () <<
".");
3736 TEUCHOS_TEST_FOR_EXCEPTION(
3737 vectorIndexOutOfRange(col),
3739 "Tpetra::MultiVector::sumIntoGlobalValue: Vector index " << col
3740 <<
" of the multivector is invalid.");
3745 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3747 Teuchos::ArrayRCP<T>
3753 typename dual_view_type::array_layout,
3756 col_dual_view_type X_col =
3757 Kokkos::subview (
view_, Kokkos::ALL (), col);
3758 return Kokkos::Compat::persistingView (X_col.d_view);
3761 template <
class Scalar,
3763 class GlobalOrdinal,
3766 KokkosClassic::MultiVector<Scalar, Node>
3770 using Teuchos::ArrayRCP;
3771 typedef KokkosClassic::MultiVector<Scalar, Node> KMV;
3775 ArrayRCP<Scalar> data;
3777 data = Teuchos::null;
3781 Kokkos::Compat::persistingView (
view_.d_view) :
3783 data = Teuchos::arcp_reinterpret_cast<Scalar> (dataTmp);
3790 KMV kmv (this->
getMap ()->getNode ());
3797 template <
class Scalar,
3799 class GlobalOrdinal,
3808 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3814 std::ostringstream oss;
3815 oss << Teuchos::typeName (*
this) <<
" {" 3816 <<
"label: \"" << this->getObjectLabel () <<
"\"" 3820 if (isConstantStride ()) {
3821 oss <<
", columnStride: " <<
getStride ();
3827 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3831 const Teuchos::EVerbosityLevel verbLevel)
const 3833 using Teuchos::ArrayRCP;
3835 using Teuchos::VERB_DEFAULT;
3836 using Teuchos::VERB_NONE;
3837 using Teuchos::VERB_LOW;
3838 using Teuchos::VERB_MEDIUM;
3839 using Teuchos::VERB_HIGH;
3840 using Teuchos::VERB_EXTREME;
3845 const Teuchos::EVerbosityLevel vl =
3846 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
3848 RCP<const Teuchos::Comm<int> > comm = this->
getMap()->getComm();
3849 const int myImageID = comm->getRank();
3850 const int numImages = comm->getSize();
3852 if (vl != VERB_NONE) {
3854 Teuchos::OSTab tab0 (out);
3856 if (myImageID == 0) {
3857 out <<
"Tpetra::MultiVector:" << endl;
3858 Teuchos::OSTab tab1 (out);
3859 out <<
"Template parameters:" << endl;
3861 Teuchos::OSTab tab2 (out);
3862 out <<
"Scalar: " << Teuchos::TypeNameTraits<Scalar>::name () << endl
3863 <<
"LocalOrdinal: " << Teuchos::TypeNameTraits<LocalOrdinal>::name () << endl
3864 <<
"GlobalOrdinal: " << Teuchos::TypeNameTraits<GlobalOrdinal>::name () << endl
3865 <<
"Node: " << Teuchos::TypeNameTraits<Node>::name () << endl;
3867 out <<
"label: \"" << this->getObjectLabel () <<
"\"" << endl
3872 out <<
"columnStride: " <<
getStride () << endl;
3875 for (
int imageCtr = 0; imageCtr < numImages; ++imageCtr) {
3876 if (myImageID == imageCtr) {
3877 if (vl != VERB_LOW) {
3879 out <<
"Process " << myImageID <<
":" << endl;
3880 Teuchos::OSTab tab2 (out);
3885 if (vl != VERB_MEDIUM) {
3888 out <<
"columnStride: " <<
getStride() << endl;
3890 if (vl == VERB_EXTREME) {
3892 out <<
"values: " << endl;
3893 typename dual_view_type::t_host X = this->
getDualView ().h_view;
3899 if (j + 1 < getNumVectors ()) {
3903 if (i + 1 < getLocalLength ()) {
3920 #if TPETRA_USE_KOKKOS_DISTOBJECT 3921 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3929 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3937 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3945 #else // NOT TPETRA_USE_KOKKOS_DISTOBJECT 3947 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3953 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3959 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3965 #endif // TPETRA_USE_KOKKOS_DISTOBJECT 3967 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3975 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node, const
bool classic>
3980 using Kokkos::parallel_for;
3981 typedef LocalOrdinal LO;
3983 typedef typename dual_view_type::host_mirror_space::device_type HMDT;
3984 const bool debug =
false;
3986 TEUCHOS_TEST_FOR_EXCEPTION(
3989 "Tpetra::deep_copy: Global dimensions of the two Tpetra::MultiVector " 3990 "objects do not match. src has dimensions [" << src.
getGlobalLength ()
3991 <<
"," << src.
getNumVectors () <<
"], and *this has dimensions [" 3994 TEUCHOS_TEST_FOR_EXCEPTION(
3996 "Tpetra::deep_copy: The local row counts of the two Tpetra::MultiVector " 3997 "objects do not match. src has " << src.
getLocalLength () <<
" row(s) " 4000 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4001 std::cout <<
"*** MultiVector::assign: ";
4005 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4006 std::cout <<
"Both *this and src have constant stride" << std::endl;
4011 this->
template modify<DT> ();
4013 Details::localDeepCopyConstStride (this->
getDualView ().
template view<DT> (),
4015 this->
template sync<HMDT> ();
4018 this->
template modify<HMDT> ();
4020 Details::localDeepCopyConstStride (this->
getDualView ().
template view<HMDT> (),
4022 this->
template sync<DT> ();
4027 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4028 std::cout <<
"Only *this has constant stride";
4031 const LO numWhichVecs =
static_cast<LO
> (src.
whichVectors_.size ());
4032 const std::string whichVecsLabel (
"MV::deep_copy::whichVecs");
4038 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4039 std::cout <<
"; Copy from device version of src" << std::endl;
4045 typedef Kokkos::DualView<LO*, DT> whichvecs_type;
4046 whichvecs_type srcWhichVecs (whichVecsLabel, numWhichVecs);
4047 srcWhichVecs.template modify<HMDT> ();
4048 for (LO i = 0; i < numWhichVecs; ++i) {
4049 srcWhichVecs.h_view(i) =
static_cast<LO
> (src.
whichVectors_[i]);
4052 srcWhichVecs.template sync<DT> ();
4055 this->
template modify<DT> ();
4060 Details::localDeepCopy (this->
getDualView ().
template view<DT> (),
4062 true,
false, srcWhichVecs.d_view, srcWhichVecs.d_view);
4065 this->
template sync<HMDT> ();
4068 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4069 std::cout <<
"; Copy from host version of src" << std::endl;
4075 typedef Kokkos::View<LO*, HMDT> whichvecs_type;
4076 whichvecs_type srcWhichVecs (whichVecsLabel, numWhichVecs);
4077 for (LO i = 0; i < numWhichVecs; ++i) {
4083 Details::localDeepCopy (this->
getDualView ().
template view<HMDT> (),
4085 true,
false, srcWhichVecs, srcWhichVecs);
4087 this->
template sync<DT> ();
4092 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4093 std::cout <<
"Only src has constant stride" << std::endl;
4101 typedef Kokkos::DualView<LO*, DT> whichvecs_type;
4102 const std::string whichVecsLabel (
"MV::deep_copy::whichVecs");
4103 const LO numWhichVecs =
static_cast<LO
> (this->
whichVectors_.size ());
4104 whichvecs_type whichVecs (whichVecsLabel, numWhichVecs);
4105 whichVecs.template modify<HMDT> ();
4106 for (LO i = 0; i < numWhichVecs; ++i) {
4110 whichVecs.template sync<DT> ();
4113 Details::localDeepCopy (this->
getDualView ().
template view<DT> (),
4116 whichVecs.d_view, whichVecs.d_view);
4122 this->
template sync<HMDT> ();
4129 typedef Kokkos::View<LO*, HMDT> whichvecs_type;
4130 const LO numWhichVecs =
static_cast<LO
> (this->
whichVectors_.size ());
4131 whichvecs_type whichVecs (
"MV::deep_copy::whichVecs", numWhichVecs);
4132 for (LO i = 0; i < numWhichVecs; ++i) {
4137 Details::localDeepCopy (this->
getDualView ().
template view<HMDT> (),
4140 whichVecs, whichVecs);
4145 this->
template sync<DT> ();
4149 if (debug && this->
getMap ()->getComm ()->getRank () == 0) {
4150 std::cout <<
"Neither *this nor src has constant stride" << std::endl;
4159 const LO dstNumWhichVecs =
static_cast<LO
> (this->
whichVectors_.size ());
4160 Kokkos::DualView<LO*, DT> whichVecsDst (
"MV::deep_copy::whichVecsDst",
4162 whichVecsDst.template modify<HMDT> ();
4163 for (LO i = 0; i < dstNumWhichVecs; ++i) {
4164 whichVecsDst.h_view(i) =
static_cast<LO
> (this->
whichVectors_[i]);
4167 whichVecsDst.template sync<DT> ();
4173 const LO srcNumWhichVecs =
static_cast<LO
> (src.
whichVectors_.size ());
4174 Kokkos::DualView<LO*, DT> whichVecsSrc (
"MV::deep_copy::whichVecsSrc",
4176 whichVecsSrc.template modify<HMDT> ();
4177 for (LO i = 0; i < srcNumWhichVecs; ++i) {
4178 whichVecsSrc.h_view(i) =
static_cast<LO
> (src.
whichVectors_[i]);
4181 whichVecsSrc.template sync<DT> ();
4185 Details::localDeepCopy (this->
getDualView ().
template view<DT> (),
4188 whichVecsDst.d_view, whichVecsSrc.d_view);
4191 const LO dstNumWhichVecs =
static_cast<LO
> (this->
whichVectors_.size ());
4192 Kokkos::View<LO*, HMDT> whichVectorsDst (
"dstWhichVecs", dstNumWhichVecs);
4193 for (LO i = 0; i < dstNumWhichVecs; ++i) {
4198 const LO srcNumWhichVecs =
static_cast<LO
> (src.
whichVectors_.size ());
4199 Kokkos::View<LO*, HMDT> whichVectorsSrc (
"srcWhichVecs", srcNumWhichVecs);
4200 for (LO i = 0; i < srcNumWhichVecs; ++i) {
4206 Details::localDeepCopy (this->
getDualView ().
template view<HMDT> (),
4209 whichVectorsDst, whichVectorsSrc);
4216 this->
template sync<HMDT> ();
4223 template <
class Scalar,
class LO,
class GO,
class NT, const
bool classic>
4224 Teuchos::RCP<MultiVector<Scalar, LO, GO, NT, classic> >
4229 return Teuchos::rcp (
new MV (map, numVectors));
4232 template <
class ST,
class LO,
class GO,
class NT, const
bool classic>
4250 #define TPETRA_MULTIVECTOR_INSTANT(SCALAR,LO,GO,NODE) \ 4251 template class MultiVector< SCALAR , LO , GO , NODE >; \ 4252 template MultiVector< SCALAR , LO , GO , NODE > createCopy( const MultiVector< SCALAR , LO , GO , NODE >& src); \ 4253 template Teuchos::RCP<MultiVector< SCALAR , LO , GO , NODE > > createMultiVector (const Teuchos::RCP<const Map<LO, GO, NODE> >& map, size_t numVectors); 4255 #endif // TPETRA_MULTIVECTOR_DEF_HPP void scale(const Scalar &alpha)
Scale in place: this = alpha*this.
size_t getLocalLength() const
Local number of rows on the calling process.
EWhichNormImpl
Input argument for localNormImpl() (which see).
Namespace Tpetra contains the class and methods constituting the Tpetra library.
void get1dCopy(const Teuchos::ArrayView< Scalar > &A, const size_t LDA) const
Fill the given array with a copy of this multivector's local values.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > subViewNonConst(const Teuchos::Range1D &colRng)
Return a MultiVector with views of selected columns.
void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A, const Scalar &beta)
Update: this = beta*this + alpha*A.
friend void deep_copy(MultiVector< DS, DL, DG, DN, dstClassic > &dst, const MultiVector< SS, SL, SG, SN, srcClassic > &src)
Copy the contents of the MultiVector src into dst.
Teuchos::ArrayRCP< const Scalar > get1dView() const
Const persisting (1-D) view of this multivector's local values.
virtual void removeEmptyProcessesInPlace(const Teuchos::RCP< const map_type > &newMap)
Remove processes owning zero rows from the Map and their communicator.
MultiVector< ST, LO, GO, NT, classic > createCopy(const MultiVector< ST, LO, GO, NT, classic > &src)
Return a deep copy of the given MultiVector.
void normInf(const Kokkos::View< mag_type *, device_type > &norms) const
Compute the infinity-norm of each vector (column), storing the result in a device view...
void abs(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A)
Put element-wise absolute values of input Multi-vector in target: A = abs(this)
Teuchos::Array< size_t > whichVectors_
Indices of columns this multivector is viewing.
void replaceGlobalValue(GlobalOrdinal globalRow, size_t col, const impl_scalar_type &value) const
Replace value, using global (row) index.
One or more distributed dense vectors.
Teuchos::ArrayRCP< Scalar > getDataNonConst(size_t j)
View of the local values in a particular vector of this multivector.
Teuchos::ArrayRCP< T > getSubArrayRCP(Teuchos::ArrayRCP< T > arr, size_t j) const
Persisting view of j-th column in the given ArrayRCP.
void get2dCopy(const Teuchos::ArrayView< const Teuchos::ArrayView< Scalar > > &ArrayOfPtrs) const
Fill the given array with a copy of this multivector's local values.
Node::device_type device_type
The Kokkos device type.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > offsetViewNonConst(const Teuchos::RCP< const map_type > &subMap, const size_t offset)
Return a nonconst view of a subset of rows.
size_t getOrigNumLocalCols() const
"Original" number of columns in the (local) data.
void randomize()
Set all values in the multivector to pseudorandom numbers.
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.
bool isConstantStride() const
Whether this multivector has constant stride between columns.
virtual bool checkSizes(const SrcDistObject &sourceObj)
Whether data redistribution between sourceObj and this object is legal.
dual_view_type view_
The Kokkos::DualView containing the MultiVector's data.
Kokkos::Details::InnerProductSpaceTraits< impl_scalar_type >::dot_type dot_type
Type of an inner ("dot") product result.
void norm1(const Kokkos::View< mag_type *, device_type > &norms) const
Compute the one-norm of each vector (column), storing the result in a device view.
Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > getVector(const size_t j) const
Return a Vector which is a const view of column j.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > createMultiVector(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, const size_t numVectors)
Nonmember MultiVector "constructor": Create a MultiVector from a given Map.
void replaceMap(const Teuchos::RCP< const map_type > &map)
Replace the underlying Map in place.
void meanValue(const Teuchos::ArrayView< impl_scalar_type > &means) const
Compute mean (average) value of each column.
void reduce()
Sum values of a locally replicated multivector across all processes.
size_t global_size_t
Global size_t object.
Teuchos::ArrayRCP< Teuchos::ArrayRCP< const Scalar > > get2dView() const
Return const persisting pointers to values.
void norm2(const Kokkos::View< mag_type *, device_type > &norms) const
Compute the two-norm of each vector (column), storing the result in a device view.
Insert new values that don't currently exist.
void createViews() const
Hook for creating a const view.
dual_view_type getDualView() const
Get the Kokkos::DualView which implements local storage.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with the given verbosity level to a FancyOStream.
void GEMM(const char transA[], const char transB[], const CoefficientType &alpha, const ViewType1 &A, const ViewType2 &B, const CoefficientType &beta, const ViewType3 &C)
Small dense matrix-matrix multiply: C := alpha*A*B + beta*C
Teuchos::RCP< const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > offsetView(const Teuchos::RCP< const map_type > &subMap, const size_t offset) const
Return a const view of a subset of rows.
void sumIntoGlobalValue(const GlobalOrdinal globalRow, const size_t col, const impl_scalar_type &value, const bool atomic=useAtomicUpdatesByDefault) const
Add value to existing value, using global (row) index.
Node::execution_space execution_space
Type of the (new) Kokkos execution space.
void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, const Scalar &beta)
Matrix-matrix multiplication: this = beta*this + alpha*op(A)*op(B).
Sets up and executes a communication plan for a Tpetra DistObject.
CombineMode
Rule for combining data in an Import or Export.
Sum new values into existing values.
MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > & operator=(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &source)
Assign the contents of source to this multivector (deep copy).
Teuchos::RCP< const map_type > map_
The Map over which this object is distributed.
void createViewsNonConst(KokkosClassic::ReadWriteOption rwo)
Hook for creating a nonconst view.
bool isDistributed() const
Whether this is a globally distributed object.
Teuchos::ArrayRCP< const Scalar > getData(size_t j) const
Const view of the local values in a particular vector of this multivector.
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.
Teuchos::ArrayRCP< Scalar > get1dViewNonConst()
Nonconst persisting (1-D) view of this multivector's local values.
Replace existing values with new values.
size_t getStride() const
Stride between columns in the multivector.
void reciprocal(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A)
Put element-wise reciprocal values of input Multi-vector in target, this(i,j) = 1/A(i,j).
KokkosClassic::MultiVector< Scalar, Node > getLocalMV() const
A view of the underlying KokkosClassic::MultiVector object.
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
global_size_t getGlobalLength() const
Global number of rows in the multivector.
Kokkos::Details::ArithTraits< Scalar >::val_type impl_scalar_type
The type used internally in place of Scalar.
virtual ~MultiVector()
Destructor (virtual for memory safety of derived classes).
size_t getOrigNumLocalRows() const
"Original" number of rows in the (local) data.
void releaseViews() const
Hook for releasing views.
Kokkos::Details::ArithTraits< impl_scalar_type >::mag_type mag_type
Type of a norm result.
size_t getNodeNumElements() const
The number of elements belonging to the calling process.
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.
void TPETRA_DEPRECATED normWeighted(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &weights, const Teuchos::ArrayView< mag_type > &norms) const
Compute Weighted 2-norm (RMS Norm) of each column.
Class that provides GEMM for a particular Kokkos Device.
A distributed dense vector.
Stand-alone utility functions and macros.
virtual std::string description() const
A simple one-line description of this object.
void assign(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &src)
Copy the contents of src into *this (deep copy).
size_t getNumVectors() const
Number of columns in the multivector.
MultiVector()
Default constructor: makes a MultiVector with no rows or columns.
Teuchos::ArrayRCP< Teuchos::ArrayRCP< Scalar > > get2dViewNonConst()
Return non-const persisting pointers to values.
EWhichNorm
Input argument for normImpl() (which see).
virtual Teuchos::RCP< const map_type > getMap() const
The Map describing the parallel distribution of this object.
virtual size_t constantNumberOfPackets() const
Number of packets to send per LID.
void dot(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A, const Teuchos::ArrayView< dot_type > &dots) const
Compute the dot product of each corresponding pair of vectors (columns) in A and B.
Teuchos::RCP< Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > getVectorNonConst(const size_t j)
Return a Vector which is a nonconst view of column j.
Teuchos::RCP< const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > subView(const Teuchos::Range1D &colRng) const
Return a const MultiVector with const views of selected columns.
void elementWiseMultiply(Scalar scalarAB, const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &A, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > &B, Scalar scalarThis)
Multiply a Vector A elementwise by a MultiVector B.
Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node, classic > > subCopy(const Teuchos::Range1D &colRng) const
Return a MultiVector with copies of selected columns.
void normImpl(const Kokkos::View< mag_type *, device_type > &norms, const EWhichNorm whichNorm) const
Compute the norm of each vector (column), storing the result in a device View.
dual_view_type origView_
The "original view" of the MultiVector's data.
void sumIntoLocalValue(const LocalOrdinal localRow, const size_t col, const impl_scalar_type &value, const bool atomic=useAtomicUpdatesByDefault) const
Add value to existing value, using local (row) index.
void replaceLocalValue(LocalOrdinal localRow, size_t col, const impl_scalar_type &value) const
Replace value, using local (row) index.