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>
25 #include "Teuchos_CompilerCodeTweakMacros.hpp"
26 
27 
28 namespace 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  auto nCols = this->mv_->getNumVectors();
600  auto nRows = this->mv_->getGlobalLength();
601  auto comm = this->mv_->getMap()->getComm();
602  auto myRank = comm->getRank();
603  if (myRank == 0) {
604  Kokkos::resize(kokkos_new_view, nRows, nCols);
605  if (perm_g2l.extent(0) == nRows) {
606  Kokkos::resize(this->buf_, nRows, 1);
607  } else {
608  Kokkos::resize(this->buf_, 0, 1);
609  }
610  }
611  {
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++) {
617  // lclMV with ReadOnly = const sendbuf
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(),
622  0, *comm);
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);
625  }
626  }
627  }
628  return 0;
629  }
630 
631  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
632  template<typename KV, typename host_ordinal_type_array>
633  LocalOrdinal
634  MultiVecAdapter<
635  MultiVector<Scalar,
636  LocalOrdinal,
637  GlobalOrdinal,
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,
642  EDistribution distribution) const
643  {
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();
649  {
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);
655  }
656  // lclMV with OverwriteAll
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(),
660  recvbuf, nRows_l,
661  0, *comm);
662  }
663  Kokkos::deep_copy(lclMV_d, lclMV);
664  }
665  return 0;
666  }
667 
668  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
669  std::string
670  MultiVecAdapter<
671  MultiVector<Scalar,
672  LocalOrdinal,
673  GlobalOrdinal,
674  Node> >::description() const
675  {
676  std::ostringstream oss;
677  oss << "Amesos2 adapter wrapping: ";
678  oss << mv_->description();
679  return oss.str();
680  }
681 
682 
683  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
684  void
685  MultiVecAdapter<
686  MultiVector<Scalar,
687  LocalOrdinal,
688  GlobalOrdinal,
689  Node> >::describe (Teuchos::FancyOStream& os,
690  const Teuchos::EVerbosityLevel verbLevel) const
691  {
692  mv_->describe (os, verbLevel);
693  }
694 
695 
696  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
697  const char* MultiVecAdapter<
698  MultiVector<Scalar,
699  LocalOrdinal,
700  GlobalOrdinal,
701  Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
702 
703 } // end namespace Amesos2
704 
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