19 #ifndef AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
20 #define AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
22 #include <type_traits>
25 #include "Teuchos_CompilerCodeTweakMacros.hpp"
30 using Tpetra::MultiVector;
32 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
37 Node> >::MultiVecAdapter(
const Teuchos::RCP<multivec_t>& m )
41 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
51 Node> >::clone()
const
53 using MV = MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
54 Teuchos::RCP<MV> Y (
new MV (mv_->getMap(), mv_->getNumVectors(),
false));
55 Y->setCopyOrView (Teuchos::View);
59 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
60 typename MultiVecAdapter<
64 Node> >::multivec_t::impl_scalar_type *
69 Node> >::getMVPointer_impl()
const
71 TEUCHOS_TEST_FOR_EXCEPTION( this->getGlobalNumVectors() != 1,
72 std::invalid_argument,
73 "Amesos2_TpetraMultiVectorAdapter: getMVPointer_impl should only be called for case with a single vector and single MPI process" );
75 auto contig_local_view_2d = mv_->getLocalViewHost(Tpetra::Access::ReadWrite);
76 auto contig_local_view_1d = Kokkos::subview (contig_local_view_2d, Kokkos::ALL (), 0);
77 return contig_local_view_1d.data();
88 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
94 Node> >::get1dCopy(
const Teuchos::ArrayView<scalar_t>& av,
97 const Tpetra::Map<LocalOrdinal,
99 Node> > distribution_map,
104 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
105 const size_t num_vecs = getGlobalNumVectors ();
107 TEUCHOS_TEST_FOR_EXCEPTION(
108 distribution_map.getRawPtr () == NULL, std::invalid_argument,
109 "Amesos2::MultiVecAdapter::get1dCopy: distribution_map argument is null.");
110 TEUCHOS_TEST_FOR_EXCEPTION(
111 mv_.is_null (), std::logic_error,
112 "Amesos2::MultiVecAdapter::get1dCopy: mv_ is null.");
114 TEUCHOS_TEST_FOR_EXCEPTION(
115 this->getMap ().is_null (), std::logic_error,
116 "Amesos2::MultiVecAdapter::get1dCopy: this->getMap() returns null.");
118 #ifdef HAVE_AMESOS2_DEBUG
119 const size_t requested_vector_length = distribution_map->getLocalNumElements ();
120 TEUCHOS_TEST_FOR_EXCEPTION(
121 lda < requested_vector_length, std::invalid_argument,
122 "Amesos2::MultiVecAdapter::get1dCopy: On process " <<
123 distribution_map->getComm ()->getRank () <<
" of the distribution Map's "
124 "communicator, the given stride lda = " << lda <<
" is not large enough "
125 "for the local vector length " << requested_vector_length <<
".");
126 TEUCHOS_TEST_FOR_EXCEPTION(
127 as<size_t> (av.size ()) < as<size_t> ((num_vecs - 1) * lda + requested_vector_length),
128 std::invalid_argument,
"Amesos2::MultiVector::get1dCopy: MultiVector "
129 "storage not large enough given leading dimension and number of vectors." );
130 #endif // HAVE_AMESOS2_DEBUG
133 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
134 mv_->get1dCopy (av, lda);
141 RCP<const map_type> distMap;
142 if (exporter_.is_null () ||
143 ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
144 ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
150 distMap = rcp(
new map_type(*distribution_map));
152 exporter_ = rcp (
new export_type (this->getMap (), distMap));
155 distMap = exporter_->getTargetMap ();
158 multivec_t redist_mv (distMap, num_vecs);
161 redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
163 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
166 redist_mv.get1dCopy (av, lda);
172 auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
173 if ( redist_mv.isConstantStride() ) {
174 for (
size_t j = 0; j < num_vecs; ++j) {
175 auto av_j = av(lda*j, lda);
176 for (
size_t i = 0; i < lda; ++i ) {
177 av_j[i] = contig_local_view_2d(i,j);
186 const size_t lclNumRows = redist_mv.getLocalLength();
187 for (
size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
188 auto av_j = av(lda*j, lclNumRows);
189 auto X_lcl_j_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
190 auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
192 using val_type =
typename std::remove_const<typename decltype( X_lcl_j_1d )::value_type>::type;
193 Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
194 Kokkos::deep_copy (umavj, X_lcl_j_1d);
201 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
202 template <
typename KV>
205 MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>
206 >::get1dCopy_kokkos_view(
209 [[maybe_unused]]
size_t lda,
210 Teuchos::Ptr<
const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > distribution_map,
215 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
216 const size_t num_vecs = getGlobalNumVectors ();
218 TEUCHOS_TEST_FOR_EXCEPTION(
219 distribution_map.getRawPtr () == NULL, std::invalid_argument,
220 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: distribution_map argument is null.");
221 TEUCHOS_TEST_FOR_EXCEPTION(
222 mv_.is_null (), std::logic_error,
223 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: mv_ is null.");
225 TEUCHOS_TEST_FOR_EXCEPTION(
226 this->getMap ().is_null (), std::logic_error,
227 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: this->getMap() returns null.");
229 #ifdef HAVE_AMESOS2_DEBUG
230 const size_t requested_vector_length = distribution_map->getLocalNumElements ();
231 TEUCHOS_TEST_FOR_EXCEPTION(
232 lda < requested_vector_length, std::invalid_argument,
233 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: On process " <<
234 distribution_map->getComm ()->getRank () <<
" of the distribution Map's "
235 "communicator, the given stride lda = " << lda <<
" is not large enough "
236 "for the local vector length " << requested_vector_length <<
".");
240 #endif // HAVE_AMESOS2_DEBUG
243 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
244 if(mv_->isConstantStride()) {
247 deep_copy_only(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
251 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"Resolve handling for non-constant stride.");
252 TEUCHOS_UNREACHABLE_RETURN(
false);
260 RCP<const map_type> distMap;
261 if (exporter_.is_null () ||
262 ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
263 ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
269 distMap = rcp(
new map_type(*distribution_map));
271 exporter_ = rcp (
new export_type (this->getMap (), distMap));
274 distMap = exporter_->getTargetMap ();
277 multivec_t redist_mv (distMap, num_vecs);
280 redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
282 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
286 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
290 if(redist_mv.isConstantStride()) {
292 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
296 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"Kokkos adapter non-constant stride not imlemented.");
297 TEUCHOS_UNREACHABLE_RETURN(
false);
303 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
304 Teuchos::ArrayRCP<Scalar>
309 Node> >::get1dViewNonConst (
bool local)
315 TEUCHOS_TEST_FOR_EXCEPTION(
316 true, std::logic_error,
"Amesos2::MultiVecAdapter::get1dViewNonConst: "
370 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node>
376 Node> >::put1dData(
const Teuchos::ArrayView<const scalar_t>& new_data,
379 const Tpetra::Map<LocalOrdinal,
385 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
387 TEUCHOS_TEST_FOR_EXCEPTION(
388 source_map.getRawPtr () == NULL, std::invalid_argument,
389 "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
390 TEUCHOS_TEST_FOR_EXCEPTION(
391 mv_.is_null (), std::logic_error,
392 "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
394 TEUCHOS_TEST_FOR_EXCEPTION(
395 this->getMap ().is_null (), std::logic_error,
396 "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
398 const size_t num_vecs = getGlobalNumVectors ();
401 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
403 auto mv_view_to_modify_2d = mv_->getLocalViewHost(Tpetra::Access::OverwriteAll);
404 for (
size_t i = 0; i < lda; ++i ) {
405 mv_view_to_modify_2d(i,0) = new_data[i];
413 RCP<const map_type> srcMap;
414 if (importer_.is_null () ||
415 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
416 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
422 srcMap = rcp(
new map_type(*source_map));
423 importer_ = rcp (
new import_type (srcMap, this->getMap ()));
426 srcMap = importer_->getSourceMap ();
429 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
432 const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
433 mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
436 multivec_t redist_mv (srcMap, num_vecs);
437 if ( redist_mv.isConstantStride() ) {
438 auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
439 for (
size_t j = 0; j < num_vecs; ++j) {
440 auto av_j = new_data(lda*j, lda);
441 for (
size_t i = 0; i < lda; ++i ) {
442 contig_local_view_2d(i,j) = av_j[i];
451 const size_t lclNumRows = redist_mv.getLocalLength();
452 for (
size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
453 auto av_j = new_data(lda*j, lclNumRows);
454 auto X_lcl_j_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
455 auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
457 using val_type =
typename std::remove_const<typename decltype( X_lcl_j_1d )::value_type>::type;
458 Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
459 Kokkos::deep_copy (umavj, X_lcl_j_1d);
466 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
472 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node>
473 template <
typename KV>
479 Node> >::put1dData_kokkos_view(KV& kokkos_new_data,
482 const Tpetra::Map<LocalOrdinal,
488 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
490 TEUCHOS_TEST_FOR_EXCEPTION(
491 source_map.getRawPtr () == NULL, std::invalid_argument,
492 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: source_map argument is null.");
493 TEUCHOS_TEST_FOR_EXCEPTION(
494 mv_.is_null (), std::logic_error,
495 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: the internal MultiVector mv_ is null.");
497 TEUCHOS_TEST_FOR_EXCEPTION(
498 this->getMap ().is_null (), std::logic_error,
499 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: this->getMap() returns null.");
501 const size_t num_vecs = getGlobalNumVectors ();
504 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
509 auto mv_view_to_modify_2d = mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
511 deep_copy_only(mv_view_to_modify_2d, kokkos_new_data);
518 RCP<const map_type> srcMap;
519 if (importer_.is_null () ||
520 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
521 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
527 srcMap = rcp(
new map_type(*source_map));
528 importer_ = rcp (
new import_type (srcMap, this->getMap ()));
531 srcMap = importer_->getSourceMap ();
534 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
539 typedef typename multivec_t::dual_view_type::t_host::value_type tpetra_mv_view_type;
540 Kokkos::View<tpetra_mv_view_type**,
typename KV::array_layout,
541 Kokkos::HostSpace> convert_kokkos_new_data;
542 deep_copy_or_assign_view(convert_kokkos_new_data, kokkos_new_data);
543 #ifdef HAVE_TEUCHOS_COMPLEX
545 auto pData =
reinterpret_cast<Scalar*
>(convert_kokkos_new_data.data());
547 auto pData = convert_kokkos_new_data.data();
550 const multivec_t source_mv (srcMap, Teuchos::ArrayView<const scalar_t>(
551 pData, kokkos_new_data.size()), lda, num_vecs);
552 mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
555 multivec_t redist_mv (srcMap, num_vecs);
560 auto host_kokkos_new_data = Kokkos::create_mirror_view(kokkos_new_data);
561 Kokkos::deep_copy(host_kokkos_new_data, kokkos_new_data);
562 if ( redist_mv.isConstantStride() ) {
563 auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
564 for (
size_t j = 0; j < num_vecs; ++j) {
565 auto av_j = Kokkos::subview(host_kokkos_new_data, Kokkos::ALL, j);
566 for (
size_t i = 0; i < lda; ++i ) {
567 contig_local_view_2d(i,j) = av_j(i);
572 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::runtime_error,
"Kokkos adapter "
573 "CONTIGUOUS_AND_ROOTED not implemented for put1dData_kokkos_view "
574 "with non constant stride.");
580 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
586 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
587 template<
typename KV,
typename host_ordinal_type_array>
593 Node> >::gather (KV& kokkos_new_view,
594 host_ordinal_type_array &perm_g2l,
595 host_ordinal_type_array &recvCountRows,
596 host_ordinal_type_array &recvDisplRows,
599 auto nCols = this->mv_->getNumVectors();
600 auto nRows = this->mv_->getGlobalLength();
601 auto comm = this->mv_->getMap()->getComm();
602 auto myRank = comm->getRank();
604 Kokkos::resize(kokkos_new_view, nRows, nCols);
605 if (perm_g2l.extent(0) == nRows) {
606 Kokkos::resize(this->buf_, nRows, 1);
608 Kokkos::resize(this->buf_, 0, 1);
612 auto nRows_l = this->mv_->getLocalLength();
613 auto lclMV_d = this->mv_->getLocalViewDevice(Tpetra::Access::ReadOnly);
614 auto lclMV = Kokkos::create_mirror_view(lclMV_d);
615 Kokkos::deep_copy(lclMV, lclMV_d);
616 for (
size_t j=0; j<nCols; j++) {
618 const Scalar * sendbuf =
reinterpret_cast<const Scalar*
> (nRows_l > 0 ? &lclMV(0,j) : lclMV.data());
619 Scalar * recvbuf =
reinterpret_cast<Scalar*
> (this->buf_.extent(0) > 0 || myRank != 0 ? this->buf_.data() : &kokkos_new_view(0,j));
620 Teuchos::gatherv<int, Scalar> (sendbuf, nRows_l,
621 recvbuf, recvCountRows.data(), recvDisplRows.data(),
623 if (myRank == 0 && this->buf_.extent(0) > 0) {
624 for (global_size_t i=0; i<nRows; i++) kokkos_new_view(perm_g2l(i),j) = this->buf_(i,0);
631 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
632 template<
typename KV,
typename host_ordinal_type_array>
638 Node> >::scatter (KV& kokkos_new_view,
639 host_ordinal_type_array &perm_g2l,
640 host_ordinal_type_array &sendCountRows,
641 host_ordinal_type_array &sendDisplRows,
644 auto nCols = this->mv_->getNumVectors();
645 auto nRows = this->mv_->getGlobalLength();
646 auto nRows_l = this->mv_->getLocalLength();
647 auto comm = this->mv_->getMap()->getComm();
648 auto myRank = comm->getRank();
650 auto lclMV_d = this->mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
651 auto lclMV = Kokkos::create_mirror_view(lclMV_d);
652 for (
size_t j=0; j<nCols; j++) {
653 if (myRank == 0 && this->buf_.extent(0) > 0) {
654 for (global_size_t i=0; i<nRows; i++) this->buf_(i, 0) = kokkos_new_view(perm_g2l(i),j);
657 Scalar * sendbuf =
reinterpret_cast<Scalar*
> (this->buf_.extent(0) > 0 || myRank != 0 ? this->buf_.data() : &kokkos_new_view(0,j));
658 Scalar * recvbuf =
reinterpret_cast<Scalar*
> (nRows_l > 0 ? &lclMV(0,j) : lclMV.data());
659 Teuchos::scatterv<int, Scalar> (sendbuf, sendCountRows.data(), sendDisplRows.data(),
663 Kokkos::deep_copy(lclMV_d, lclMV);
668 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
674 Node> >::description()
const
676 std::ostringstream oss;
677 oss <<
"Amesos2 adapter wrapping: ";
678 oss << mv_->description();
683 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
689 Node> >::describe (Teuchos::FancyOStream& os,
690 const Teuchos::EVerbosityLevel verbLevel)
const
692 mv_->describe (os, verbLevel);
696 template <
typename Scalar,
typename LocalOrdinal,
typename GlobalOrdinal,
class Node >
697 const char* MultiVecAdapter<
701 Node> >::name =
"Amesos2 adapter for Tpetra::MultiVector";
705 #endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Copy or assign views based on memory spaces.
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.
EDistribution
Definition: Amesos2_TypeDecl.hpp:89