Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_TpetraMultiVecAdapter_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Amesos2: Templated Direct Sparse Solver Package
4//
5// Copyright 2011 NTESS and the Amesos2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
19#ifndef AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
20#define AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
21
22#include <type_traits>
23#include "Amesos2_TpetraMultiVecAdapter_decl.hpp"
24#include "Amesos2_Kokkos_View_Copy_Assign.hpp"
25#include "Teuchos_CompilerCodeTweakMacros.hpp"
26
27
28namespace Amesos2 {
29
30 using Tpetra::MultiVector;
31
32 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
33 MultiVecAdapter<
34 MultiVector<Scalar,
35 LocalOrdinal,
36 GlobalOrdinal,
37 Node> >::MultiVecAdapter( const Teuchos::RCP<multivec_t>& m )
38 : mv_(m)
39 {}
40
41 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
42 Teuchos::RCP<
43 MultiVector<Scalar,
44 LocalOrdinal,
45 GlobalOrdinal,
46 Node> >
47 MultiVecAdapter<
48 MultiVector<Scalar,
49 LocalOrdinal,
50 GlobalOrdinal,
51 Node> >::clone() const
52 {
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);
56 return Y;
57 }
58
59 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
60 typename MultiVecAdapter<
61 MultiVector<Scalar,
62 LocalOrdinal,
63 GlobalOrdinal,
64 Node> >::multivec_t::impl_scalar_type *
65 MultiVecAdapter<
66 MultiVector<Scalar,
67 LocalOrdinal,
68 GlobalOrdinal,
69 Node> >::getMVPointer_impl() const
70 {
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" );
74
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();
78 }
79
80 // TODO Proper type handling:
81 // Consider a MultiVectorTraits class
82 // typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type
83 // NOTE: In this class, above already has a typedef multivec_t
84 // typedef typename multivector_type::impl_scalar_type return_scalar_type; // this is the POD type the dual_view_type is templated on
85 // Traits class needed to do this generically for the general MultiVectorAdapter interface
86
87
88 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
89 void
90 MultiVecAdapter<
91 MultiVector<Scalar,
92 LocalOrdinal,
93 GlobalOrdinal,
94 Node> >::get1dCopy(const Teuchos::ArrayView<scalar_t>& av,
95 size_t lda,
96 Teuchos::Ptr<
97 const Tpetra::Map<LocalOrdinal,
98 GlobalOrdinal,
99 Node> > distribution_map,
100 EDistribution distribution) const
101 {
102 using Teuchos::as;
103 using Teuchos::RCP;
104 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
105 const size_t num_vecs = getGlobalNumVectors ();
106
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.");
113 // Check mv_ before getMap(), because the latter calls mv_->getMap().
114 TEUCHOS_TEST_FOR_EXCEPTION(
115 this->getMap ().is_null (), std::logic_error,
116 "Amesos2::MultiVecAdapter::get1dCopy: this->getMap() returns null.");
117
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
131
132 // Special case when number vectors == 1 and single MPI process
133 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
134 mv_->get1dCopy (av, lda);
135 }
136 else {
137
138 // (Re)compute the Export object if necessary. If not, then we
139 // don't need to clone distribution_map; we can instead just get
140 // the previously cloned target Map from the Export object.
141 RCP<const map_type> distMap;
142 if (exporter_.is_null () ||
143 ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
144 ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
145 // Since we're caching the Export object, and since the Export
146 // needs to keep the distribution Map, we have to make a copy of
147 // the latter in order to ensure that it will stick around past
148 // the scope of this function call. (Ptr is not reference
149 // counted.)
150 distMap = rcp(new map_type(*distribution_map));
151 // (Re)create the Export object.
152 exporter_ = rcp (new export_type (this->getMap (), distMap));
153 }
154 else {
155 distMap = exporter_->getTargetMap ();
156 }
157
158 multivec_t redist_mv (distMap, num_vecs);
159
160 // Redistribute the input (multi)vector.
161 redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
162
163 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
164 // Do this if GIDs contiguous - existing functionality
165 // Copy the imported (multi)vector's data into the ArrayView.
166 redist_mv.get1dCopy (av, lda);
167 }
168 else {
169 // Do this if GIDs not contiguous...
170 // sync is needed for example if mv was updated on device, but will be passed through Amesos2 to solver running on host
171 //redist_mv.template sync < host_execution_space > ();
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); //lda may not be correct if redist_mv is not constant stride...
178 }
179 }
180 }
181 else {
182 // ... lda should come from Teuchos::Array* allocation,
183 // not the MultiVector, since the MultiVector does NOT
184 // have constant stride in this case.
185 // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
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);
191
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);
195 }
196 }
197 }
198 }
199 }
200
201 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
202 template <typename KV>
203 bool
204 MultiVecAdapter<
205 MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>
206 >::get1dCopy_kokkos_view(
207 bool bInitialize,
208 KV& kokkos_view,
209 [[maybe_unused]] size_t lda,
210 Teuchos::Ptr<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > distribution_map,
211 EDistribution distribution
212 ) const {
213 using Teuchos::as;
214 using Teuchos::RCP;
215 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
216 const size_t num_vecs = getGlobalNumVectors ();
217
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.");
224 // Check mv_ before getMap(), because the latter calls mv_->getMap().
225 TEUCHOS_TEST_FOR_EXCEPTION(
226 this->getMap ().is_null (), std::logic_error,
227 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: this->getMap() returns null.");
228
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 << ".");
237
238 // Note do not check size since deep_copy_or_assign_view below will allocate
239 // if necessary - but may just assign ptrs.
240#endif // HAVE_AMESOS2_DEBUG
241
242 // Special case when number vectors == 1 and single MPI process
243 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
244 if(mv_->isConstantStride()) {
245 bool bAssigned;
246 //deep_copy_or_assign_view(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
247 deep_copy_only(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
248 return bAssigned; // if bAssigned is true we are accessing the mv data directly without a copy
249 }
250 else {
251 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Resolve handling for non-constant stride.");
252 TEUCHOS_UNREACHABLE_RETURN(false);
253 }
254 }
255 else {
256
257 // (Re)compute the Export object if necessary. If not, then we
258 // don't need to clone distribution_map; we can instead just get
259 // the previously cloned target Map from the Export object.
260 RCP<const map_type> distMap;
261 if (exporter_.is_null () ||
262 ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
263 ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
264 // Since we're caching the Export object, and since the Export
265 // needs to keep the distribution Map, we have to make a copy of
266 // the latter in order to ensure that it will stick around past
267 // the scope of this function call. (Ptr is not reference
268 // counted.)
269 distMap = rcp(new map_type(*distribution_map));
270 // (Re)create the Export object.
271 exporter_ = rcp (new export_type (this->getMap (), distMap));
272 }
273 else {
274 distMap = exporter_->getTargetMap ();
275 }
276
277 multivec_t redist_mv (distMap, num_vecs);
278
279 // Redistribute the input (multi)vector.
280 redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
281
282 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
283 // Do this if GIDs contiguous - existing functionality
284 // Copy the imported (multi)vector's data into the Kokkos View.
285 bool bAssigned;
286 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
287 return false; // do not return bAssigned because redist_mv was already copied so always return false
288 }
289 else {
290 if(redist_mv.isConstantStride()) {
291 bool bAssigned; // deep_copy_or_assign_view sets true if assigned (no deep copy)
292 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
293 return false; // do not return bAssigned because redist_mv was already copied so always return false
294 }
295 else {
296 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Kokkos adapter non-constant stride not imlemented.");
297 TEUCHOS_UNREACHABLE_RETURN(false);
298 }
299 }
300 }
301 }
302
303 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
304 Teuchos::ArrayRCP<Scalar>
305 MultiVecAdapter<
306 MultiVector<Scalar,
307 LocalOrdinal,
308 GlobalOrdinal,
309 Node> >::get1dViewNonConst (bool local)
310 {
311 // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
312 // its code was commented out, and it didn't return anything. The
313 // latter is ESPECIALLY dangerous, given that the return value is
314 // an ArrayRCP. Thus, I added the exception throw below.
315 TEUCHOS_TEST_FOR_EXCEPTION(
316 true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
317 "Not implemented.");
318
319 // if( local ){
320 // this->localize();
321 // /* Use the global element list returned by
322 // * mv_->getMap()->getLocalElementList() to get a subCopy of mv_ which we
323 // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
324 // */
325 // if(l_l_mv_.is_null() ){
326 // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
327 // = mv_->getMap()->getLocalElementList();
328 // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
329
330 // // Convert the node element to a list of size_t type objects
331 // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
332 // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
333 // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
334 // *(it_st++) = Teuchos::as<size_t>(*it_go);
335 // }
336
337 // // To be consistent with the globalize() function, get a view of the local mv
338 // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
339
340 // return(l_l_mv_->get1dViewNonConst());
341 // } else {
342 // // We need to re-import values to the local, since changes may have been
343 // // made to the global structure that are not reflected in the local
344 // // view.
345 // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
346 // = mv_->getMap()->getLocalElementList();
347 // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
348
349 // // Convert the node element to a list of size_t type objects
350 // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
351 // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
352 // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
353 // *(it_st++) = Teuchos::as<size_t>(*it_go);
354 // }
355
356 // return l_l_mv_->get1dViewNonConst();
357 // }
358 // } else {
359 // if( mv_->isDistributed() ){
360 // this->localize();
361
362 // return l_mv_->get1dViewNonConst();
363 // } else { // not distributed, no need to import
364 // return mv_->get1dViewNonConst();
365 // }
366 // }
367 }
368
369
370 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
371 void
372 MultiVecAdapter<
373 MultiVector<Scalar,
374 LocalOrdinal,
375 GlobalOrdinal,
376 Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
377 size_t lda,
378 Teuchos::Ptr<
379 const Tpetra::Map<LocalOrdinal,
380 GlobalOrdinal,
381 Node> > source_map,
382 EDistribution distribution )
383 {
384 using Teuchos::RCP;
385 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
386
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.");
393 // getMap() calls mv_->getMap(), so test first whether mv_ is null.
394 TEUCHOS_TEST_FOR_EXCEPTION(
395 this->getMap ().is_null (), std::logic_error,
396 "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
397
398 const size_t num_vecs = getGlobalNumVectors ();
399
400 // Special case when number vectors == 1 and single MPI process
401 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
402 // num_vecs = 1; stride does not matter
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]; // Only one vector
406 }
407 }
408 else {
409
410 // (Re)compute the Import object if necessary. If not, then we
411 // don't need to clone source_map; we can instead just get the
412 // previously cloned source Map from the Import object.
413 RCP<const map_type> srcMap;
414 if (importer_.is_null () ||
415 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
416 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
417
418 // Since we're caching the Import object, and since the Import
419 // needs to keep the source Map, we have to make a copy of the
420 // latter in order to ensure that it will stick around past the
421 // scope of this function call. (Ptr is not reference counted.)
422 srcMap = rcp(new map_type(*source_map));
423 importer_ = rcp (new import_type (srcMap, this->getMap ()));
424 }
425 else {
426 srcMap = importer_->getSourceMap ();
427 }
428
429 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
430 // Do this if GIDs contiguous - existing functionality
431 // Redistribute the output (multi)vector.
432 const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
433 mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
434 }
435 else {
436 multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
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];
443 }
444 }
445 }
446 else {
447 // ... lda should come from Teuchos::Array* allocation,
448 // not the MultiVector, since the MultiVector does NOT
449 // have constant stride in this case.
450 // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
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);
456
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);
460 }
461 }
462
463 //typedef typename multivec_t::node_type::memory_space memory_space;
464 //redist_mv.template sync <memory_space> ();
465
466 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
467 }
468 }
469
470 }
471
472 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
473 template <typename KV>
474 void
475 MultiVecAdapter<
476 MultiVector<Scalar,
477 LocalOrdinal,
478 GlobalOrdinal,
479 Node> >::put1dData_kokkos_view(KV& kokkos_new_data,
480 size_t lda,
481 Teuchos::Ptr<
482 const Tpetra::Map<LocalOrdinal,
483 GlobalOrdinal,
484 Node> > source_map,
485 EDistribution distribution )
486 {
487 using Teuchos::RCP;
488 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
489
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.");
496 // getMap() calls mv_->getMap(), so test first whether 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.");
500
501 const size_t num_vecs = getGlobalNumVectors ();
502
503 // Special case when number vectors == 1 and single MPI process
504 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
505
506 // num_vecs = 1; stride does not matter
507
508 // If this is the optimized path then kokkos_new_data will be the dst
509 auto mv_view_to_modify_2d = mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
510 //deep_copy_or_assign_view(mv_view_to_modify_2d, kokkos_new_data);
511 deep_copy_only(mv_view_to_modify_2d, kokkos_new_data);
512 }
513 else {
514
515 // (Re)compute the Import object if necessary. If not, then we
516 // don't need to clone source_map; we can instead just get the
517 // previously cloned source Map from the Import object.
518 RCP<const map_type> srcMap;
519 if (importer_.is_null () ||
520 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
521 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
522
523 // Since we're caching the Import object, and since the Import
524 // needs to keep the source Map, we have to make a copy of the
525 // latter in order to ensure that it will stick around past the
526 // scope of this function call. (Ptr is not reference counted.)
527 srcMap = rcp(new map_type(*source_map));
528 importer_ = rcp (new import_type (srcMap, this->getMap ()));
529 }
530 else {
531 srcMap = importer_->getSourceMap ();
532 }
533
534 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
535 // Use View scalar type, not MV Scalar because we want Kokkos::complex, not
536 // std::complex to avoid a Kokkos::complex<double> to std::complex<float>
537 // conversion which would require a double copy and fail here. Then we'll be
538 // setup to safely reinterpret_cast complex to std if necessary.
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
544 // convert_kokkos_new_data may be Kokkos::complex and Scalar could be std::complex
545 auto pData = reinterpret_cast<Scalar*>(convert_kokkos_new_data.data());
546#else
547 auto pData = convert_kokkos_new_data.data();
548#endif
549
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);
553 }
554 else {
555 multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
556 // Cuda solvers won't currently use this path since they are just serial
557 // right now, so this mirror should be harmless (and not strictly necessary).
558 // Adding it for future possibilities though we may then refactor this
559 // for better efficiency if the kokkos_new_data view is on device.
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);
568 }
569 }
570 }
571 else {
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.");
575 }
576
577 //typedef typename multivec_t::node_type::memory_space memory_space;
578 //redist_mv.template sync <memory_space> ();
579
580 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
581 }
582 }
583
584 }
585
586 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
587 template<typename KV, typename host_ordinal_type_array>
588 LocalOrdinal
589 MultiVecAdapter<
590 MultiVector<Scalar,
591 LocalOrdinal,
592 GlobalOrdinal,
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,
597 EDistribution distribution) const
598 {
599 using Teuchos::as;
600 using mv_scalar_t = typename KV::non_const_value_type;
601 size_t nCols = this->mv_->getNumVectors();
602 size_t nRows = this->mv_->getGlobalLength();
603 size_t nRows_l = this->mv_->getLocalLength();
604 auto comm = this->mv_->getMap()->getComm();
605 auto myRank = comm->getRank();
606
607 bool mixed_precision = !(std::is_same<Scalar, mv_scalar_t>::value);
608 if (myRank == 0) {
609 if (kokkos_new_view.extent(0) != nRows || kokkos_new_view.extent(1) != nCols) {
610 Kokkos::resize(kokkos_new_view, nRows, nCols);
611 }
612 if (perm_g2l.extent(0) == nRows || mixed_precision) {
613 if (this->buf_.extent(0) < nRows) Kokkos::resize(this->buf_, nRows, 1);
614 } else {
615 Kokkos::resize(this->buf_, 0, 1);
616 }
617 }
618 {
619 auto lclMV_d = this->mv_->getLocalViewDevice(Tpetra::Access::ReadOnly);
620 auto lclMV = Kokkos::create_mirror_view(lclMV_d);
621 Kokkos::deep_copy(lclMV, lclMV_d);
622 for (size_t j=0; j<nCols; j++) {
623 // lclMV with ReadOnly = const sendbuf
624 const Scalar * sendbuf = reinterpret_cast<const Scalar*> (nRows_l > 0 ? &lclMV(0,j) : lclMV.data());
625 void * recvbuf = (this->buf_.extent(0) > 0 || myRank != 0 ? (void*)(this->buf_.data()) : (void*)(&kokkos_new_view(0,j)));
626 Teuchos::gatherv<int, Scalar> (sendbuf, nRows_l,
627 (Scalar*)recvbuf, recvCountRows.data(), recvDisplRows.data(),
628 0, *comm);
629 if (myRank == 0 && this->buf_.extent(0) > 0) {
630 typedef Kokkos::DefaultHostExecutionSpace HostExecSpaceType;
631 if (perm_g2l.extent(0) == nRows) {
632 // permute, and type-casting if needed
633 Kokkos::parallel_for("Amesos2::TpetraMultiVecAdapter::gather", Kokkos::RangePolicy<HostExecSpaceType>(0, nRows),
634 [&](const int i) { kokkos_new_view(perm_g2l(i),j) = as<mv_scalar_t>(this->buf_(i,0)); });
635 } else {
636 // type-casting
637 Kokkos::parallel_for("Amesos2::TpetraMultiVecAdapter::gather", Kokkos::RangePolicy<HostExecSpaceType>(0, nRows),
638 [&](const int i) { kokkos_new_view(i,j) = as<mv_scalar_t>(this->buf_(i,0)); });
639 }
640 }
641 }
642 }
643 return 0;
644 }
645
646 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
647 template<typename KV, typename host_ordinal_type_array>
648 LocalOrdinal
649 MultiVecAdapter<
650 MultiVector<Scalar,
651 LocalOrdinal,
652 GlobalOrdinal,
653 Node> >::scatter (KV& kokkos_new_view,
654 host_ordinal_type_array &perm_g2l,
655 host_ordinal_type_array &sendCountRows,
656 host_ordinal_type_array &sendDisplRows,
657 EDistribution distribution) const
658 {
659 using Teuchos::as;
660 using mv_scalar_t = typename KV::non_const_value_type;
661 size_t nCols = this->mv_->getNumVectors();
662 size_t nRows = this->mv_->getGlobalLength();
663 size_t nRows_l = this->mv_->getLocalLength();
664 auto comm = this->mv_->getMap()->getComm();
665 auto myRank = comm->getRank();
666
667 {
668 auto lclMV_d = this->mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
669 auto lclMV = Kokkos::create_mirror_view(lclMV_d);
670 for (size_t j=0; j<nCols; j++) {
671 if (myRank == 0 && this->buf_.extent(0) > 0) {
672 typedef Kokkos::DefaultHostExecutionSpace HostExecSpaceType;
673 if (perm_g2l.extent(0) == nRows) {
674 // permute, and type-casting if needed
675 Kokkos::parallel_for("Amesos2::TpetraMultiVecAdapter::scatter", Kokkos::RangePolicy<HostExecSpaceType>(0, nRows),
676 [&](const int i) { this->buf_(i, 0) = as<mv_scalar_t>(kokkos_new_view(perm_g2l(i),j)); });
677 } else {
678 // type-casting
679 Kokkos::parallel_for("Amesos2::TpetraMultiVecAdapter::scatter", Kokkos::RangePolicy<HostExecSpaceType>(0, nRows),
680 [&](const int i) { this->buf_(i, 0) = as<mv_scalar_t>(kokkos_new_view(i,j)); });
681 }
682 }
683 // lclMV with OverwriteAll
684 void * sendbuf = (this->buf_.extent(0) > 0 || myRank != 0 ? (void*)(this->buf_.data()) : (void*)(&kokkos_new_view(0,j)));
685 Scalar * recvbuf = reinterpret_cast<Scalar*> (nRows_l > 0 ? &lclMV(0,j) : lclMV.data());
686 Teuchos::scatterv<int, Scalar> ((Scalar*)sendbuf, sendCountRows.data(), sendDisplRows.data(),
687 recvbuf, nRows_l,
688 0, *comm);
689 }
690 Kokkos::deep_copy(lclMV_d, lclMV);
691 }
692 return 0;
693 }
694
695 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
696 std::string
697 MultiVecAdapter<
698 MultiVector<Scalar,
699 LocalOrdinal,
700 GlobalOrdinal,
701 Node> >::description() const
702 {
703 std::ostringstream oss;
704 oss << "Amesos2 adapter wrapping: ";
705 oss << mv_->description();
706 return oss.str();
707 }
708
709
710 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
711 void
712 MultiVecAdapter<
713 MultiVector<Scalar,
714 LocalOrdinal,
715 GlobalOrdinal,
716 Node> >::describe (Teuchos::FancyOStream& os,
717 const Teuchos::EVerbosityLevel verbLevel) const
718 {
719 mv_->describe (os, verbLevel);
720 }
721
722
723 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
724 const char* MultiVecAdapter<
725 MultiVector<Scalar,
726 LocalOrdinal,
727 GlobalOrdinal,
728 Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
729
730} // end namespace Amesos2
731
732#endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Amesos2_Kokkos_View_Copy_Assign.hpp
Copy or assign views based on memory spaces.
Amesos2_TpetraMultiVecAdapter_decl.hpp
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.
Generated by doxygen 1.9.8