Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_TpetraMultiVecAdapter_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Amesos2: Templated Direct Sparse Solver Package
6 // Copyright 2011 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //
42 // @HEADER
43 
53 #ifndef AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
54 #define AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
55 
56 #include <type_traits>
59 
60 
61 namespace Amesos2 {
62 
63  using Tpetra::MultiVector;
64 
65  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
66  MultiVecAdapter<
67  MultiVector<Scalar,
68  LocalOrdinal,
69  GlobalOrdinal,
70  Node> >::MultiVecAdapter( const Teuchos::RCP<multivec_t>& m )
71  : mv_(m)
72  {}
73 
74  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
75  Teuchos::RCP<
76  MultiVector<Scalar,
77  LocalOrdinal,
78  GlobalOrdinal,
79  Node> >
80  MultiVecAdapter<
81  MultiVector<Scalar,
82  LocalOrdinal,
83  GlobalOrdinal,
84  Node> >::clone() const
85  {
86  using MV = MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
87  Teuchos::RCP<MV> Y (new MV (mv_->getMap(), mv_->getNumVectors(), false));
88  Y->setCopyOrView (Teuchos::View);
89  return Y;
90  }
91 
92  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
93  typename MultiVecAdapter<
94  MultiVector<Scalar,
95  LocalOrdinal,
96  GlobalOrdinal,
97  Node> >::multivec_t::impl_scalar_type *
98  MultiVecAdapter<
99  MultiVector<Scalar,
100  LocalOrdinal,
101  GlobalOrdinal,
102  Node> >::getMVPointer_impl() const
103  {
104  TEUCHOS_TEST_FOR_EXCEPTION( this->getGlobalNumVectors() != 1,
105  std::invalid_argument,
106  "Amesos2_TpetraMultiVectorAdapter: getMVPointer_impl should only be called for case with a single vector and single MPI process" );
107 
108  auto contig_local_view_2d = mv_->getLocalViewHost(Tpetra::Access::ReadWrite);
109  auto contig_local_view_1d = Kokkos::subview (contig_local_view_2d, Kokkos::ALL (), 0);
110  return contig_local_view_1d.data();
111  }
112 
113  // TODO Proper type handling:
114  // Consider a MultiVectorTraits class
115  // typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type
116  // NOTE: In this class, above already has a typedef multivec_t
117  // typedef typename multivector_type::impl_scalar_type return_scalar_type; // this is the POD type the dual_view_type is templated on
118  // Traits class needed to do this generically for the general MultiVectorAdapter interface
119 
120 
121  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
122  void
123  MultiVecAdapter<
124  MultiVector<Scalar,
125  LocalOrdinal,
126  GlobalOrdinal,
127  Node> >::get1dCopy(const Teuchos::ArrayView<scalar_t>& av,
128  size_t lda,
129  Teuchos::Ptr<
130  const Tpetra::Map<LocalOrdinal,
131  GlobalOrdinal,
132  Node> > distribution_map,
133  EDistribution distribution) const
134  {
135  using Teuchos::as;
136  using Teuchos::RCP;
137  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
138  const size_t num_vecs = getGlobalNumVectors ();
139 
140  TEUCHOS_TEST_FOR_EXCEPTION(
141  distribution_map.getRawPtr () == NULL, std::invalid_argument,
142  "Amesos2::MultiVecAdapter::get1dCopy: distribution_map argument is null.");
143  TEUCHOS_TEST_FOR_EXCEPTION(
144  mv_.is_null (), std::logic_error,
145  "Amesos2::MultiVecAdapter::get1dCopy: mv_ is null.");
146  // Check mv_ before getMap(), because the latter calls mv_->getMap().
147  TEUCHOS_TEST_FOR_EXCEPTION(
148  this->getMap ().is_null (), std::logic_error,
149  "Amesos2::MultiVecAdapter::get1dCopy: this->getMap() returns null.");
150 
151 #ifdef HAVE_AMESOS2_DEBUG
152  const size_t requested_vector_length = distribution_map->getLocalNumElements ();
153  TEUCHOS_TEST_FOR_EXCEPTION(
154  lda < requested_vector_length, std::invalid_argument,
155  "Amesos2::MultiVecAdapter::get1dCopy: On process " <<
156  distribution_map->getComm ()->getRank () << " of the distribution Map's "
157  "communicator, the given stride lda = " << lda << " is not large enough "
158  "for the local vector length " << requested_vector_length << ".");
159  TEUCHOS_TEST_FOR_EXCEPTION(
160  as<size_t> (av.size ()) < as<size_t> ((num_vecs - 1) * lda + requested_vector_length),
161  std::invalid_argument, "Amesos2::MultiVector::get1dCopy: MultiVector "
162  "storage not large enough given leading dimension and number of vectors." );
163 #endif // HAVE_AMESOS2_DEBUG
164 
165  // Special case when number vectors == 1 and single MPI process
166  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
167  mv_->get1dCopy (av, lda);
168  }
169  else {
170 
171  // (Re)compute the Export object if necessary. If not, then we
172  // don't need to clone distribution_map; we can instead just get
173  // the previously cloned target Map from the Export object.
174  RCP<const map_type> distMap;
175  if (exporter_.is_null () ||
176  ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
177  ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
178  // Since we're caching the Export object, and since the Export
179  // needs to keep the distribution Map, we have to make a copy of
180  // the latter in order to ensure that it will stick around past
181  // the scope of this function call. (Ptr is not reference
182  // counted.)
183  distMap = rcp(new map_type(*distribution_map));
184  // (Re)create the Export object.
185  exporter_ = rcp (new export_type (this->getMap (), distMap));
186  }
187  else {
188  distMap = exporter_->getTargetMap ();
189  }
190 
191  multivec_t redist_mv (distMap, num_vecs);
192 
193  // Redistribute the input (multi)vector.
194  redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
195 
196  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
197  // Do this if GIDs contiguous - existing functionality
198  // Copy the imported (multi)vector's data into the ArrayView.
199  redist_mv.get1dCopy (av, lda);
200  }
201  else {
202  // Do this if GIDs not contiguous...
203  // sync is needed for example if mv was updated on device, but will be passed through Amesos2 to solver running on host
204  //redist_mv.template sync < host_execution_space > ();
205  auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
206  if ( redist_mv.isConstantStride() ) {
207  for ( size_t j = 0; j < num_vecs; ++j) {
208  auto av_j = av(lda*j, lda);
209  for ( size_t i = 0; i < lda; ++i ) {
210  av_j[i] = contig_local_view_2d(i,j); //lda may not be correct if redist_mv is not constant stride...
211  }
212  }
213  }
214  else {
215  // ... lda should come from Teuchos::Array* allocation,
216  // not the MultiVector, since the MultiVector does NOT
217  // have constant stride in this case.
218  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
219  const size_t lclNumRows = redist_mv.getLocalLength();
220  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
221  auto av_j = av(lda*j, lclNumRows);
222  auto X_lcl_j_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
223  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
224 
225  using val_type = typename std::remove_const<typename decltype( X_lcl_j_1d )::value_type>::type;
226  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
227  Kokkos::deep_copy (umavj, X_lcl_j_1d);
228  }
229  }
230  }
231  }
232  }
233 
234  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
235  template <typename KV>
236  bool
237  MultiVecAdapter<
238  MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>
239  >::get1dCopy_kokkos_view(
240  bool bInitialize,
241  KV& kokkos_view,
242  [[maybe_unused]] size_t lda,
243  Teuchos::Ptr<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > distribution_map,
244  EDistribution distribution
245  ) const {
246  using Teuchos::as;
247  using Teuchos::RCP;
248  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
249  const size_t num_vecs = getGlobalNumVectors ();
250 
251  TEUCHOS_TEST_FOR_EXCEPTION(
252  distribution_map.getRawPtr () == NULL, std::invalid_argument,
253  "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: distribution_map argument is null.");
254  TEUCHOS_TEST_FOR_EXCEPTION(
255  mv_.is_null (), std::logic_error,
256  "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: mv_ is null.");
257  // Check mv_ before getMap(), because the latter calls mv_->getMap().
258  TEUCHOS_TEST_FOR_EXCEPTION(
259  this->getMap ().is_null (), std::logic_error,
260  "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: this->getMap() returns null.");
261 
262 #ifdef HAVE_AMESOS2_DEBUG
263  const size_t requested_vector_length = distribution_map->getLocalNumElements ();
264  TEUCHOS_TEST_FOR_EXCEPTION(
265  lda < requested_vector_length, std::invalid_argument,
266  "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: On process " <<
267  distribution_map->getComm ()->getRank () << " of the distribution Map's "
268  "communicator, the given stride lda = " << lda << " is not large enough "
269  "for the local vector length " << requested_vector_length << ".");
270 
271  // Note do not check size since deep_copy_or_assign_view below will allocate
272  // if necessary - but may just assign ptrs.
273 #endif // HAVE_AMESOS2_DEBUG
274 
275  // Special case when number vectors == 1 and single MPI process
276  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
277  if(mv_->isConstantStride()) {
278  bool bAssigned;
279  //deep_copy_or_assign_view(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
280  deep_copy_only(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
281  return bAssigned; // if bAssigned is true we are accessing the mv data directly without a copy
282  }
283  else {
284  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Resolve handling for non-constant stride.");
285  }
286  }
287  else {
288 
289  // (Re)compute the Export object if necessary. If not, then we
290  // don't need to clone distribution_map; we can instead just get
291  // the previously cloned target Map from the Export object.
292  RCP<const map_type> distMap;
293  if (exporter_.is_null () ||
294  ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
295  ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
296  // Since we're caching the Export object, and since the Export
297  // needs to keep the distribution Map, we have to make a copy of
298  // the latter in order to ensure that it will stick around past
299  // the scope of this function call. (Ptr is not reference
300  // counted.)
301  distMap = rcp(new map_type(*distribution_map));
302  // (Re)create the Export object.
303  exporter_ = rcp (new export_type (this->getMap (), distMap));
304  }
305  else {
306  distMap = exporter_->getTargetMap ();
307  }
308 
309  multivec_t redist_mv (distMap, num_vecs);
310 
311  // Redistribute the input (multi)vector.
312  redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
313 
314  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
315  // Do this if GIDs contiguous - existing functionality
316  // Copy the imported (multi)vector's data into the Kokkos View.
317  bool bAssigned;
318  deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
319  return false; // do not return bAssigned because redist_mv was already copied so always return false
320  }
321  else {
322  if(redist_mv.isConstantStride()) {
323  bool bAssigned; // deep_copy_or_assign_view sets true if assigned (no deep copy)
324  deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
325  return false; // do not return bAssigned because redist_mv was already copied so always return false
326  }
327  else {
328  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Kokkos adapter non-constant stride not imlemented.");
329  }
330  }
331  }
332  }
333 
334  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
335  Teuchos::ArrayRCP<Scalar>
336  MultiVecAdapter<
337  MultiVector<Scalar,
338  LocalOrdinal,
339  GlobalOrdinal,
340  Node> >::get1dViewNonConst (bool local)
341  {
342  // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
343  // its code was commented out, and it didn't return anything. The
344  // latter is ESPECIALLY dangerous, given that the return value is
345  // an ArrayRCP. Thus, I added the exception throw below.
346  TEUCHOS_TEST_FOR_EXCEPTION(
347  true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
348  "Not implemented.");
349 
350  // if( local ){
351  // this->localize();
352  // /* Use the global element list returned by
353  // * mv_->getMap()->getLocalElementList() to get a subCopy of mv_ which we
354  // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
355  // */
356  // if(l_l_mv_.is_null() ){
357  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
358  // = mv_->getMap()->getLocalElementList();
359  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
360 
361  // // Convert the node element to a list of size_t type objects
362  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
363  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
364  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
365  // *(it_st++) = Teuchos::as<size_t>(*it_go);
366  // }
367 
368  // // To be consistent with the globalize() function, get a view of the local mv
369  // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
370 
371  // return(l_l_mv_->get1dViewNonConst());
372  // } else {
373  // // We need to re-import values to the local, since changes may have been
374  // // made to the global structure that are not reflected in the local
375  // // view.
376  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
377  // = mv_->getMap()->getLocalElementList();
378  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
379 
380  // // Convert the node element to a list of size_t type objects
381  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
382  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
383  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
384  // *(it_st++) = Teuchos::as<size_t>(*it_go);
385  // }
386 
387  // return l_l_mv_->get1dViewNonConst();
388  // }
389  // } else {
390  // if( mv_->isDistributed() ){
391  // this->localize();
392 
393  // return l_mv_->get1dViewNonConst();
394  // } else { // not distributed, no need to import
395  // return mv_->get1dViewNonConst();
396  // }
397  // }
398  }
399 
400 
401  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
402  void
403  MultiVecAdapter<
404  MultiVector<Scalar,
405  LocalOrdinal,
406  GlobalOrdinal,
407  Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
408  size_t lda,
409  Teuchos::Ptr<
410  const Tpetra::Map<LocalOrdinal,
411  GlobalOrdinal,
412  Node> > source_map,
413  EDistribution distribution )
414  {
415  using Teuchos::RCP;
416  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
417 
418  TEUCHOS_TEST_FOR_EXCEPTION(
419  source_map.getRawPtr () == NULL, std::invalid_argument,
420  "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
421  TEUCHOS_TEST_FOR_EXCEPTION(
422  mv_.is_null (), std::logic_error,
423  "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
424  // getMap() calls mv_->getMap(), so test first whether mv_ is null.
425  TEUCHOS_TEST_FOR_EXCEPTION(
426  this->getMap ().is_null (), std::logic_error,
427  "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
428 
429  const size_t num_vecs = getGlobalNumVectors ();
430 
431  // Special case when number vectors == 1 and single MPI process
432  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
433  // num_vecs = 1; stride does not matter
434  auto mv_view_to_modify_2d = mv_->getLocalViewHost(Tpetra::Access::OverwriteAll);
435  for ( size_t i = 0; i < lda; ++i ) {
436  mv_view_to_modify_2d(i,0) = new_data[i]; // Only one vector
437  }
438  }
439  else {
440 
441  // (Re)compute the Import object if necessary. If not, then we
442  // don't need to clone source_map; we can instead just get the
443  // previously cloned source Map from the Import object.
444  RCP<const map_type> srcMap;
445  if (importer_.is_null () ||
446  ! importer_->getSourceMap ()->isSameAs (* source_map) ||
447  ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
448 
449  // Since we're caching the Import object, and since the Import
450  // needs to keep the source Map, we have to make a copy of the
451  // latter in order to ensure that it will stick around past the
452  // scope of this function call. (Ptr is not reference counted.)
453  srcMap = rcp(new map_type(*source_map));
454  importer_ = rcp (new import_type (srcMap, this->getMap ()));
455  }
456  else {
457  srcMap = importer_->getSourceMap ();
458  }
459 
460  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
461  // Do this if GIDs contiguous - existing functionality
462  // Redistribute the output (multi)vector.
463  const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
464  mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
465  }
466  else {
467  multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
468  if ( redist_mv.isConstantStride() ) {
469  auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
470  for ( size_t j = 0; j < num_vecs; ++j) {
471  auto av_j = new_data(lda*j, lda);
472  for ( size_t i = 0; i < lda; ++i ) {
473  contig_local_view_2d(i,j) = av_j[i];
474  }
475  }
476  }
477  else {
478  // ... lda should come from Teuchos::Array* allocation,
479  // not the MultiVector, since the MultiVector does NOT
480  // have constant stride in this case.
481  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
482  const size_t lclNumRows = redist_mv.getLocalLength();
483  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
484  auto av_j = new_data(lda*j, lclNumRows);
485  auto X_lcl_j_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
486  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
487 
488  using val_type = typename std::remove_const<typename decltype( X_lcl_j_1d )::value_type>::type;
489  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
490  Kokkos::deep_copy (umavj, X_lcl_j_1d);
491  }
492  }
493 
494  //typedef typename multivec_t::node_type::memory_space memory_space;
495  //redist_mv.template sync <memory_space> ();
496 
497  mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
498  }
499  }
500 
501  }
502 
503  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
504  template <typename KV>
505  void
506  MultiVecAdapter<
507  MultiVector<Scalar,
508  LocalOrdinal,
509  GlobalOrdinal,
510  Node> >::put1dData_kokkos_view(KV& kokkos_new_data,
511  size_t lda,
512  Teuchos::Ptr<
513  const Tpetra::Map<LocalOrdinal,
514  GlobalOrdinal,
515  Node> > source_map,
516  EDistribution distribution )
517  {
518  using Teuchos::RCP;
519  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
520 
521  TEUCHOS_TEST_FOR_EXCEPTION(
522  source_map.getRawPtr () == NULL, std::invalid_argument,
523  "Amesos2::MultiVecAdapter::put1dData_kokkos_view: source_map argument is null.");
524  TEUCHOS_TEST_FOR_EXCEPTION(
525  mv_.is_null (), std::logic_error,
526  "Amesos2::MultiVecAdapter::put1dData_kokkos_view: the internal MultiVector mv_ is null.");
527  // getMap() calls mv_->getMap(), so test first whether mv_ is null.
528  TEUCHOS_TEST_FOR_EXCEPTION(
529  this->getMap ().is_null (), std::logic_error,
530  "Amesos2::MultiVecAdapter::put1dData_kokkos_view: this->getMap() returns null.");
531 
532  const size_t num_vecs = getGlobalNumVectors ();
533 
534  // Special case when number vectors == 1 and single MPI process
535  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
536 
537  // num_vecs = 1; stride does not matter
538 
539  // If this is the optimized path then kokkos_new_data will be the dst
540  auto mv_view_to_modify_2d = mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
541  //deep_copy_or_assign_view(mv_view_to_modify_2d, kokkos_new_data);
542  deep_copy_only(mv_view_to_modify_2d, kokkos_new_data);
543  }
544  else {
545 
546  // (Re)compute the Import object if necessary. If not, then we
547  // don't need to clone source_map; we can instead just get the
548  // previously cloned source Map from the Import object.
549  RCP<const map_type> srcMap;
550  if (importer_.is_null () ||
551  ! importer_->getSourceMap ()->isSameAs (* source_map) ||
552  ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
553 
554  // Since we're caching the Import object, and since the Import
555  // needs to keep the source Map, we have to make a copy of the
556  // latter in order to ensure that it will stick around past the
557  // scope of this function call. (Ptr is not reference counted.)
558  srcMap = rcp(new map_type(*source_map));
559  importer_ = rcp (new import_type (srcMap, this->getMap ()));
560  }
561  else {
562  srcMap = importer_->getSourceMap ();
563  }
564 
565  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
566  // Use View scalar type, not MV Scalar because we want Kokkos::complex, not
567  // std::complex to avoid a Kokkos::complex<double> to std::complex<float>
568  // conversion which would require a double copy and fail here. Then we'll be
569  // setup to safely reinterpret_cast complex to std if necessary.
570  typedef typename multivec_t::dual_view_type::t_host::value_type tpetra_mv_view_type;
571  Kokkos::View<tpetra_mv_view_type**,typename KV::array_layout,
572  Kokkos::HostSpace> convert_kokkos_new_data;
573  deep_copy_or_assign_view(convert_kokkos_new_data, kokkos_new_data);
574 #ifdef HAVE_TEUCHOS_COMPLEX
575  // convert_kokkos_new_data may be Kokkos::complex and Scalar could be std::complex
576  auto pData = reinterpret_cast<Scalar*>(convert_kokkos_new_data.data());
577 #else
578  auto pData = convert_kokkos_new_data.data();
579 #endif
580 
581  const multivec_t source_mv (srcMap, Teuchos::ArrayView<const scalar_t>(
582  pData, kokkos_new_data.size()), lda, num_vecs);
583  mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
584  }
585  else {
586  multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
587  // Cuda solvers won't currently use this path since they are just serial
588  // right now, so this mirror should be harmless (and not strictly necessary).
589  // Adding it for future possibilities though we may then refactor this
590  // for better efficiency if the kokkos_new_data view is on device.
591  auto host_kokkos_new_data = Kokkos::create_mirror_view(kokkos_new_data);
592  Kokkos::deep_copy(host_kokkos_new_data, kokkos_new_data);
593  if ( redist_mv.isConstantStride() ) {
594  auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
595  for ( size_t j = 0; j < num_vecs; ++j) {
596  auto av_j = Kokkos::subview(host_kokkos_new_data, Kokkos::ALL, j);
597  for ( size_t i = 0; i < lda; ++i ) {
598  contig_local_view_2d(i,j) = av_j(i);
599  }
600  }
601  }
602  else {
603  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Kokkos adapter "
604  "CONTIGUOUS_AND_ROOTED not implemented for put1dData_kokkos_view "
605  "with non constant stride.");
606  }
607 
608  //typedef typename multivec_t::node_type::memory_space memory_space;
609  //redist_mv.template sync <memory_space> ();
610 
611  mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
612  }
613  }
614 
615  }
616 
617  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
618  std::string
619  MultiVecAdapter<
620  MultiVector<Scalar,
621  LocalOrdinal,
622  GlobalOrdinal,
623  Node> >::description() const
624  {
625  std::ostringstream oss;
626  oss << "Amesos2 adapter wrapping: ";
627  oss << mv_->description();
628  return oss.str();
629  }
630 
631 
632  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
633  void
634  MultiVecAdapter<
635  MultiVector<Scalar,
636  LocalOrdinal,
637  GlobalOrdinal,
638  Node> >::describe (Teuchos::FancyOStream& os,
639  const Teuchos::EVerbosityLevel verbLevel) const
640  {
641  mv_->describe (os, verbLevel);
642  }
643 
644 
645  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
646  const char* MultiVecAdapter<
647  MultiVector<Scalar,
648  LocalOrdinal,
649  GlobalOrdinal,
650  Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
651 
652 } // end namespace Amesos2
653 
654 #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:123