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
61namespace 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,
239 LocalOrdinal,
240 GlobalOrdinal,
241 Node> >::get1dCopy_kokkos_view(
242 bool bInitialize,
243 KV& kokkos_view,
244 size_t lda,
245 Teuchos::Ptr<
246 const Tpetra::Map<LocalOrdinal,
247 GlobalOrdinal,
248 Node> > distribution_map,
249 EDistribution distribution) const
250 {
251 using Teuchos::as;
252 using Teuchos::RCP;
253 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
254 const size_t num_vecs = getGlobalNumVectors ();
255
256 TEUCHOS_TEST_FOR_EXCEPTION(
257 distribution_map.getRawPtr () == NULL, std::invalid_argument,
258 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: distribution_map argument is null.");
259 TEUCHOS_TEST_FOR_EXCEPTION(
260 mv_.is_null (), std::logic_error,
261 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: mv_ is null.");
262 // Check mv_ before getMap(), because the latter calls mv_->getMap().
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 this->getMap ().is_null (), std::logic_error,
265 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: this->getMap() returns null.");
266
267#ifdef HAVE_AMESOS2_DEBUG
268 const size_t requested_vector_length = distribution_map->getLocalNumElements ();
269 TEUCHOS_TEST_FOR_EXCEPTION(
270 lda < requested_vector_length, std::invalid_argument,
271 "Amesos2::MultiVecAdapter::get1dCopy_kokkos_view: On process " <<
272 distribution_map->getComm ()->getRank () << " of the distribution Map's "
273 "communicator, the given stride lda = " << lda << " is not large enough "
274 "for the local vector length " << requested_vector_length << ".");
275
276 // Note do not check size since deep_copy_or_assign_view below will allocate
277 // if necessary - but may just assign ptrs.
278#endif // HAVE_AMESOS2_DEBUG
279
280 // Special case when number vectors == 1 and single MPI process
281 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
282 if(mv_->isConstantStride()) {
283 bool bAssigned;
284 //deep_copy_or_assign_view(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
285 deep_copy_only(bInitialize, kokkos_view, mv_->getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
286 return bAssigned; // if bAssigned is true we are accessing the mv data directly without a copy
287 }
288 else {
289 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Resolve handling for non-constant stride.");
290 }
291 }
292 else {
293
294 // (Re)compute the Export object if necessary. If not, then we
295 // don't need to clone distribution_map; we can instead just get
296 // the previously cloned target Map from the Export object.
297 RCP<const map_type> distMap;
298 if (exporter_.is_null () ||
299 ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
300 ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
301 // Since we're caching the Export object, and since the Export
302 // needs to keep the distribution Map, we have to make a copy of
303 // the latter in order to ensure that it will stick around past
304 // the scope of this function call. (Ptr is not reference
305 // counted.)
306 distMap = rcp(new map_type(*distribution_map));
307 // (Re)create the Export object.
308 exporter_ = rcp (new export_type (this->getMap (), distMap));
309 }
310 else {
311 distMap = exporter_->getTargetMap ();
312 }
313
314 multivec_t redist_mv (distMap, num_vecs);
315
316 // Redistribute the input (multi)vector.
317 redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
318
319 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
320 // Do this if GIDs contiguous - existing functionality
321 // Copy the imported (multi)vector's data into the Kokkos View.
322 bool bAssigned;
323 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
324 return false; // do not return bAssigned because redist_mv was already copied so always return false
325 }
326 else {
327 if(redist_mv.isConstantStride()) {
328 bool bAssigned; // deep_copy_or_assign_view sets true if assigned (no deep copy)
329 deep_copy_or_assign_view(bInitialize, kokkos_view, redist_mv.getLocalViewDevice(Tpetra::Access::ReadOnly), bAssigned);
330 return false; // do not return bAssigned because redist_mv was already copied so always return false
331 }
332 else {
333 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Kokkos adapter non-constant stride not imlemented.");
334 }
335 }
336 }
337 }
338
339 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
340 Teuchos::ArrayRCP<Scalar>
341 MultiVecAdapter<
342 MultiVector<Scalar,
343 LocalOrdinal,
344 GlobalOrdinal,
345 Node> >::get1dViewNonConst (bool local)
346 {
347 // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
348 // its code was commented out, and it didn't return anything. The
349 // latter is ESPECIALLY dangerous, given that the return value is
350 // an ArrayRCP. Thus, I added the exception throw below.
351 TEUCHOS_TEST_FOR_EXCEPTION(
352 true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
353 "Not implemented.");
354
355 // if( local ){
356 // this->localize();
357 // /* Use the global element list returned by
358 // * mv_->getMap()->getLocalElementList() to get a subCopy of mv_ which we
359 // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
360 // */
361 // if(l_l_mv_.is_null() ){
362 // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
363 // = mv_->getMap()->getLocalElementList();
364 // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
365
366 // // Convert the node element to a list of size_t type objects
367 // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
368 // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
369 // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
370 // *(it_st++) = Teuchos::as<size_t>(*it_go);
371 // }
372
373 // // To be consistent with the globalize() function, get a view of the local mv
374 // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
375
376 // return(l_l_mv_->get1dViewNonConst());
377 // } else {
378 // // We need to re-import values to the local, since changes may have been
379 // // made to the global structure that are not reflected in the local
380 // // view.
381 // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
382 // = mv_->getMap()->getLocalElementList();
383 // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
384
385 // // Convert the node element to a list of size_t type objects
386 // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
387 // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
388 // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
389 // *(it_st++) = Teuchos::as<size_t>(*it_go);
390 // }
391
392 // return l_l_mv_->get1dViewNonConst();
393 // }
394 // } else {
395 // if( mv_->isDistributed() ){
396 // this->localize();
397
398 // return l_mv_->get1dViewNonConst();
399 // } else { // not distributed, no need to import
400 // return mv_->get1dViewNonConst();
401 // }
402 // }
403 }
404
405
406 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
407 void
408 MultiVecAdapter<
409 MultiVector<Scalar,
410 LocalOrdinal,
411 GlobalOrdinal,
412 Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
413 size_t lda,
414 Teuchos::Ptr<
415 const Tpetra::Map<LocalOrdinal,
416 GlobalOrdinal,
417 Node> > source_map,
418 EDistribution distribution )
419 {
420 using Teuchos::RCP;
421 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
422
423 TEUCHOS_TEST_FOR_EXCEPTION(
424 source_map.getRawPtr () == NULL, std::invalid_argument,
425 "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
426 TEUCHOS_TEST_FOR_EXCEPTION(
427 mv_.is_null (), std::logic_error,
428 "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
429 // getMap() calls mv_->getMap(), so test first whether mv_ is null.
430 TEUCHOS_TEST_FOR_EXCEPTION(
431 this->getMap ().is_null (), std::logic_error,
432 "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
433
434 const size_t num_vecs = getGlobalNumVectors ();
435
436 // Special case when number vectors == 1 and single MPI process
437 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
438 // num_vecs = 1; stride does not matter
439 auto mv_view_to_modify_2d = mv_->getLocalViewHost(Tpetra::Access::OverwriteAll);
440 for ( size_t i = 0; i < lda; ++i ) {
441 mv_view_to_modify_2d(i,0) = new_data[i]; // Only one vector
442 }
443 }
444 else {
445
446 // (Re)compute the Import object if necessary. If not, then we
447 // don't need to clone source_map; we can instead just get the
448 // previously cloned source Map from the Import object.
449 RCP<const map_type> srcMap;
450 if (importer_.is_null () ||
451 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
452 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
453
454 // Since we're caching the Import object, and since the Import
455 // needs to keep the source Map, we have to make a copy of the
456 // latter in order to ensure that it will stick around past the
457 // scope of this function call. (Ptr is not reference counted.)
458 srcMap = rcp(new map_type(*source_map));
459 importer_ = rcp (new import_type (srcMap, this->getMap ()));
460 }
461 else {
462 srcMap = importer_->getSourceMap ();
463 }
464
465 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
466 // Do this if GIDs contiguous - existing functionality
467 // Redistribute the output (multi)vector.
468 const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
469 mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
470 }
471 else {
472 multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
473 if ( redist_mv.isConstantStride() ) {
474 auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
475 for ( size_t j = 0; j < num_vecs; ++j) {
476 auto av_j = new_data(lda*j, lda);
477 for ( size_t i = 0; i < lda; ++i ) {
478 contig_local_view_2d(i,j) = av_j[i];
479 }
480 }
481 }
482 else {
483 // ... lda should come from Teuchos::Array* allocation,
484 // not the MultiVector, since the MultiVector does NOT
485 // have constant stride in this case.
486 // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
487 const size_t lclNumRows = redist_mv.getLocalLength();
488 for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
489 auto av_j = new_data(lda*j, lclNumRows);
490 auto X_lcl_j_2d = redist_mv.getLocalViewHost(Tpetra::Access::ReadOnly);
491 auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
492
493 using val_type = typename std::remove_const<typename decltype( X_lcl_j_1d )::value_type>::type;
494 Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
495 Kokkos::deep_copy (umavj, X_lcl_j_1d);
496 }
497 }
498
499 //typedef typename multivec_t::node_type::memory_space memory_space;
500 //redist_mv.template sync <memory_space> ();
501
502 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
503 }
504 }
505
506 }
507
508 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
509 template <typename KV>
510 void
511 MultiVecAdapter<
512 MultiVector<Scalar,
513 LocalOrdinal,
514 GlobalOrdinal,
515 Node> >::put1dData_kokkos_view(KV& kokkos_new_data,
516 size_t lda,
517 Teuchos::Ptr<
518 const Tpetra::Map<LocalOrdinal,
519 GlobalOrdinal,
520 Node> > source_map,
521 EDistribution distribution )
522 {
523 using Teuchos::RCP;
524 typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
525
526 TEUCHOS_TEST_FOR_EXCEPTION(
527 source_map.getRawPtr () == NULL, std::invalid_argument,
528 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: source_map argument is null.");
529 TEUCHOS_TEST_FOR_EXCEPTION(
530 mv_.is_null (), std::logic_error,
531 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: the internal MultiVector mv_ is null.");
532 // getMap() calls mv_->getMap(), so test first whether mv_ is null.
533 TEUCHOS_TEST_FOR_EXCEPTION(
534 this->getMap ().is_null (), std::logic_error,
535 "Amesos2::MultiVecAdapter::put1dData_kokkos_view: this->getMap() returns null.");
536
537 const size_t num_vecs = getGlobalNumVectors ();
538
539 // Special case when number vectors == 1 and single MPI process
540 if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
541
542 // num_vecs = 1; stride does not matter
543
544 // If this is the optimized path then kokkos_new_data will be the dst
545 auto mv_view_to_modify_2d = mv_->getLocalViewDevice(Tpetra::Access::OverwriteAll);
546 //deep_copy_or_assign_view(mv_view_to_modify_2d, kokkos_new_data);
547 deep_copy_only(mv_view_to_modify_2d, kokkos_new_data);
548 }
549 else {
550
551 // (Re)compute the Import object if necessary. If not, then we
552 // don't need to clone source_map; we can instead just get the
553 // previously cloned source Map from the Import object.
554 RCP<const map_type> srcMap;
555 if (importer_.is_null () ||
556 ! importer_->getSourceMap ()->isSameAs (* source_map) ||
557 ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
558
559 // Since we're caching the Import object, and since the Import
560 // needs to keep the source Map, we have to make a copy of the
561 // latter in order to ensure that it will stick around past the
562 // scope of this function call. (Ptr is not reference counted.)
563 srcMap = rcp(new map_type(*source_map));
564 importer_ = rcp (new import_type (srcMap, this->getMap ()));
565 }
566 else {
567 srcMap = importer_->getSourceMap ();
568 }
569
570 if ( distribution != CONTIGUOUS_AND_ROOTED ) {
571 // Use View scalar type, not MV Scalar because we want Kokkos::complex, not
572 // std::complex to avoid a Kokkos::complex<double> to std::complex<float>
573 // conversion which would require a double copy and fail here. Then we'll be
574 // setup to safely reinterpret_cast complex to std if necessary.
575 typedef typename multivec_t::dual_view_type::t_host::value_type tpetra_mv_view_type;
576 Kokkos::View<tpetra_mv_view_type**,typename KV::array_layout,
577 Kokkos::HostSpace> convert_kokkos_new_data;
578 deep_copy_or_assign_view(convert_kokkos_new_data, kokkos_new_data);
579#ifdef HAVE_TEUCHOS_COMPLEX
580 // convert_kokkos_new_data may be Kokkos::complex and Scalar could be std::complex
581 auto pData = reinterpret_cast<Scalar*>(convert_kokkos_new_data.data());
582#else
583 auto pData = convert_kokkos_new_data.data();
584#endif
585
586 const multivec_t source_mv (srcMap, Teuchos::ArrayView<const scalar_t>(
587 pData, kokkos_new_data.size()), lda, num_vecs);
588 mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
589 }
590 else {
591 multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
592 // Cuda solvers won't currently use this path since they are just serial
593 // right now, so this mirror should be harmless (and not strictly necessary).
594 // Adding it for future possibilities though we may then refactor this
595 // for better efficiency if the kokkos_new_data view is on device.
596 auto host_kokkos_new_data = Kokkos::create_mirror_view(kokkos_new_data);
597 Kokkos::deep_copy(host_kokkos_new_data, kokkos_new_data);
598 if ( redist_mv.isConstantStride() ) {
599 auto contig_local_view_2d = redist_mv.getLocalViewHost(Tpetra::Access::OverwriteAll);
600 for ( size_t j = 0; j < num_vecs; ++j) {
601 auto av_j = Kokkos::subview(host_kokkos_new_data, Kokkos::ALL, j);
602 for ( size_t i = 0; i < lda; ++i ) {
603 contig_local_view_2d(i,j) = av_j(i);
604 }
605 }
606 }
607 else {
608 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Kokkos adapter "
609 "CONTIGUOUS_AND_ROOTED not implemented for put1dData_kokkos_view "
610 "with non constant stride.");
611 }
612
613 //typedef typename multivec_t::node_type::memory_space memory_space;
614 //redist_mv.template sync <memory_space> ();
615
616 mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
617 }
618 }
619
620 }
621
622 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
623 std::string
624 MultiVecAdapter<
625 MultiVector<Scalar,
626 LocalOrdinal,
627 GlobalOrdinal,
628 Node> >::description() const
629 {
630 std::ostringstream oss;
631 oss << "Amesos2 adapter wrapping: ";
632 oss << mv_->description();
633 return oss.str();
634 }
635
636
637 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
638 void
639 MultiVecAdapter<
640 MultiVector<Scalar,
641 LocalOrdinal,
642 GlobalOrdinal,
643 Node> >::describe (Teuchos::FancyOStream& os,
644 const Teuchos::EVerbosityLevel verbLevel) const
645 {
646 mv_->describe (os, verbLevel);
647 }
648
649
650 template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
651 const char* MultiVecAdapter<
652 MultiVector<Scalar,
653 LocalOrdinal,
654 GlobalOrdinal,
655 Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
656
657} // end namespace Amesos2
658
659#endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Copy or assign views based on memory spaces.
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.