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