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>
58 
59 
60 namespace Amesos2 {
61 
62  using Tpetra::MultiVector;
63 
64  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
65  MultiVecAdapter<
66  MultiVector<Scalar,
67  LocalOrdinal,
68  GlobalOrdinal,
69  Node> >::MultiVecAdapter( const Teuchos::RCP<multivec_t>& m )
70  : mv_(m)
71  {}
72 
73 
74  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
75  typename MultiVecAdapter<
76  MultiVector<Scalar,
77  LocalOrdinal,
78  GlobalOrdinal,
79  Node> >::multivec_t::impl_scalar_type *
80  MultiVecAdapter<
81  MultiVector<Scalar,
82  LocalOrdinal,
83  GlobalOrdinal,
84  Node> >::getMVPointer_impl() const
85  {
86  TEUCHOS_TEST_FOR_EXCEPTION( this->getGlobalNumVectors() != 1,
87  std::invalid_argument,
88  "Amesos2_TpetraMultiVectorAdapter: getMVPointer_impl should only be called for case with a single vector and single MPI process" );
89 
90  typedef typename multivec_t::dual_view_type dual_view_type;
91  typedef typename dual_view_type::host_mirror_space host_execution_space;
92  mv_->template sync<host_execution_space> ();
93  auto contig_local_view_2d = mv_->template getLocalView<host_execution_space>();
94  auto contig_local_view_1d = Kokkos::subview (contig_local_view_2d, Kokkos::ALL (), 0);
95  return contig_local_view_1d.data();
96  }
97 
98  // TODO Proper type handling:
99  // Consider a MultiVectorTraits class
100  // typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> multivector_type
101  // NOTE: In this class, above already has a typedef multivec_t
102  // typedef typename multivector_type::impl_scalar_type return_scalar_type; // this is the POD type the dual_view_type is templated on
103  // Traits class needed to do this generically for the general MultiVectorAdapter interface
104 
105 
106  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
107  void
108  MultiVecAdapter<
109  MultiVector<Scalar,
110  LocalOrdinal,
111  GlobalOrdinal,
112  Node> >::get1dCopy(const Teuchos::ArrayView<scalar_t>& av,
113  size_t lda,
114  Teuchos::Ptr<
115  const Tpetra::Map<LocalOrdinal,
116  GlobalOrdinal,
117  Node> > distribution_map,
118  EDistribution distribution) const
119  {
120  using Teuchos::as;
121  using Teuchos::RCP;
122  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
123  const size_t num_vecs = getGlobalNumVectors ();
124 
125  TEUCHOS_TEST_FOR_EXCEPTION(
126  distribution_map.getRawPtr () == NULL, std::invalid_argument,
127  "Amesos2::MultiVecAdapter::get1dCopy: distribution_map argument is null.");
128  TEUCHOS_TEST_FOR_EXCEPTION(
129  mv_.is_null (), std::logic_error,
130  "Amesos2::MultiVecAdapter::get1dCopy: mv_ is null.");
131  // Check mv_ before getMap(), because the latter calls mv_->getMap().
132  TEUCHOS_TEST_FOR_EXCEPTION(
133  this->getMap ().is_null (), std::logic_error,
134  "Amesos2::MultiVecAdapter::get1dCopy: this->getMap() returns null.");
135 
136 #ifdef HAVE_AMESOS2_DEBUG
137  const size_t requested_vector_length = distribution_map->getNodeNumElements ();
138  TEUCHOS_TEST_FOR_EXCEPTION(
139  lda < requested_vector_length, std::invalid_argument,
140  "Amesos2::MultiVecAdapter::get1dCopy: On process " <<
141  distribution_map->getComm ()->getRank () << " of the distribution Map's "
142  "communicator, the given stride lda = " << lda << " is not large enough "
143  "for the local vector length " << requested_vector_length << ".");
144  TEUCHOS_TEST_FOR_EXCEPTION(
145  as<size_t> (av.size ()) < as<size_t> ((num_vecs - 1) * lda + requested_vector_length),
146  std::invalid_argument, "Amesos2::MultiVector::get1dCopy: MultiVector "
147  "storage not large enough given leading dimension and number of vectors." );
148 #endif // HAVE_AMESOS2_DEBUG
149 
150  // Special case when number vectors == 1 and single MPI process
151  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
152  mv_->get1dCopy (av, lda);
153  }
154  else {
155 
156  // (Re)compute the Export object if necessary. If not, then we
157  // don't need to clone distribution_map; we can instead just get
158  // the previously cloned target Map from the Export object.
159  RCP<const map_type> distMap;
160  if (exporter_.is_null () ||
161  ! exporter_->getSourceMap ()->isSameAs (* (this->getMap ())) ||
162  ! exporter_->getTargetMap ()->isSameAs (* distribution_map)) {
163  // Since we're caching the Export object, and since the Export
164  // needs to keep the distribution Map, we have to make a copy of
165  // the latter in order to ensure that it will stick around past
166  // the scope of this function call. (Ptr is not reference
167  // counted.)
168  distMap = rcp(new map_type(*distribution_map));
169  // (Re)create the Export object.
170  exporter_ = rcp (new export_type (this->getMap (), distMap));
171  }
172  else {
173  distMap = exporter_->getTargetMap ();
174  }
175 
176  multivec_t redist_mv (distMap, num_vecs);
177 
178  // Redistribute the input (multi)vector.
179  redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
180 
181  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
182  // Do this if GIDs contiguous - existing functionality
183  // Copy the imported (multi)vector's data into the ArrayView.
184  redist_mv.get1dCopy (av, lda);
185  }
186  else {
187  // Do this if GIDs not contiguous...
188  // sync is needed for example if mv was updated on device, but will be passed through Amesos2 to solver running on host
189  typedef typename multivec_t::dual_view_type dual_view_type;
190  typedef typename dual_view_type::host_mirror_space host_execution_space;
191  redist_mv.template sync < host_execution_space > ();
192 
193  auto contig_local_view_2d = redist_mv.template getLocalView<host_execution_space>();
194  if ( redist_mv.isConstantStride() ) {
195  for ( size_t j = 0; j < num_vecs; ++j) {
196  auto av_j = av(lda*j, lda);
197  for ( size_t i = 0; i < lda; ++i ) {
198  av_j[i] = contig_local_view_2d(i,j); //lda may not be correct if redist_mv is not constant stride...
199  }
200  }
201  }
202  else {
203  // ... lda should come from Teuchos::Array* allocation,
204  // not the MultiVector, since the MultiVector does NOT
205  // have constant stride in this case.
206  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
207  const size_t lclNumRows = redist_mv.getLocalLength();
208  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
209  auto av_j = av(lda*j, lclNumRows);
210  auto X_lcl_j_2d = redist_mv.template getLocalView<host_execution_space> ();
211  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
212 
213  using val_type = typename decltype( X_lcl_j_1d )::value_type;
214  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
215  Kokkos::deep_copy (umavj, X_lcl_j_1d);
216  }
217  }
218  }
219  }
220  }
221 
222  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
223  Teuchos::ArrayRCP<Scalar>
224  MultiVecAdapter<
225  MultiVector<Scalar,
226  LocalOrdinal,
227  GlobalOrdinal,
228  Node> >::get1dViewNonConst (bool local)
229  {
230  // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
231  // its code was commented out, and it didn't return anything. The
232  // latter is ESPECIALLY dangerous, given that the return value is
233  // an ArrayRCP. Thus, I added the exception throw below.
234  TEUCHOS_TEST_FOR_EXCEPTION(
235  true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
236  "Not implemented.");
237 
238  // if( local ){
239  // this->localize();
240  // /* Use the global element list returned by
241  // * mv_->getMap()->getNodeElementList() to get a subCopy of mv_ which we
242  // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
243  // */
244  // if(l_l_mv_.is_null() ){
245  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
246  // = mv_->getMap()->getNodeElementList();
247  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
248 
249  // // Convert the node element to a list of size_t type objects
250  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
251  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
252  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
253  // *(it_st++) = Teuchos::as<size_t>(*it_go);
254  // }
255 
256  // // To be consistent with the globalize() function, get a view of the local mv
257  // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
258 
259  // return(l_l_mv_->get1dViewNonConst());
260  // } else {
261  // // We need to re-import values to the local, since changes may have been
262  // // made to the global structure that are not reflected in the local
263  // // view.
264  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
265  // = mv_->getMap()->getNodeElementList();
266  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
267 
268  // // Convert the node element to a list of size_t type objects
269  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
270  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
271  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
272  // *(it_st++) = Teuchos::as<size_t>(*it_go);
273  // }
274 
275  // return l_l_mv_->get1dViewNonConst();
276  // }
277  // } else {
278  // if( mv_->isDistributed() ){
279  // this->localize();
280 
281  // return l_mv_->get1dViewNonConst();
282  // } else { // not distributed, no need to import
283  // return mv_->get1dViewNonConst();
284  // }
285  // }
286  }
287 
288 
289  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
290  void
291  MultiVecAdapter<
292  MultiVector<Scalar,
293  LocalOrdinal,
294  GlobalOrdinal,
295  Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
296  size_t lda,
297  Teuchos::Ptr<
298  const Tpetra::Map<LocalOrdinal,
299  GlobalOrdinal,
300  Node> > source_map,
301  EDistribution distribution )
302  {
303  using Teuchos::RCP;
304  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
305 
306  TEUCHOS_TEST_FOR_EXCEPTION(
307  source_map.getRawPtr () == NULL, std::invalid_argument,
308  "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
309  TEUCHOS_TEST_FOR_EXCEPTION(
310  mv_.is_null (), std::logic_error,
311  "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
312  // getMap() calls mv_->getMap(), so test first whether mv_ is null.
313  TEUCHOS_TEST_FOR_EXCEPTION(
314  this->getMap ().is_null (), std::logic_error,
315  "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
316 
317  const size_t num_vecs = getGlobalNumVectors ();
318 
319  // Special case when number vectors == 1 and single MPI process
320  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
321  typedef typename multivec_t::dual_view_type::host_mirror_space host_execution_space;
322  // num_vecs = 1; stride does not matter
323  auto mv_view_to_modify_2d = mv_->template getLocalView<host_execution_space>();
324  for ( size_t i = 0; i < lda; ++i ) {
325  mv_view_to_modify_2d(i,0) = new_data[i]; // Only one vector
326  }
327  }
328  else {
329 
330  // (Re)compute the Import object if necessary. If not, then we
331  // don't need to clone source_map; we can instead just get the
332  // previously cloned source Map from the Import object.
333  RCP<const map_type> srcMap;
334  if (importer_.is_null () ||
335  ! importer_->getSourceMap ()->isSameAs (* source_map) ||
336  ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
337 
338  // Since we're caching the Import object, and since the Import
339  // needs to keep the source Map, we have to make a copy of the
340  // latter in order to ensure that it will stick around past the
341  // scope of this function call. (Ptr is not reference counted.)
342  srcMap = rcp(new map_type(*source_map));
343  importer_ = rcp (new import_type (srcMap, this->getMap ()));
344  }
345  else {
346  srcMap = importer_->getSourceMap ();
347  }
348 
349  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
350  // Do this if GIDs contiguous - existing functionality
351  // Redistribute the output (multi)vector.
352  const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
353  mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
354  }
355  else {
356  multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
357  typedef typename multivec_t::dual_view_type dual_view_type;
358  typedef typename dual_view_type::host_mirror_space host_execution_space;
359  redist_mv.template modify< host_execution_space > ();
360 
361  if ( redist_mv.isConstantStride() ) {
362  auto contig_local_view_2d = redist_mv.template getLocalView<host_execution_space>();
363  for ( size_t j = 0; j < num_vecs; ++j) {
364  auto av_j = new_data(lda*j, lda);
365  for ( size_t i = 0; i < lda; ++i ) {
366  contig_local_view_2d(i,j) = av_j[i];
367  }
368  }
369  }
370  else {
371  // ... lda should come from Teuchos::Array* allocation,
372  // not the MultiVector, since the MultiVector does NOT
373  // have constant stride in this case.
374  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
375  const size_t lclNumRows = redist_mv.getLocalLength();
376  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
377  auto av_j = new_data(lda*j, lclNumRows);
378  auto X_lcl_j_2d = redist_mv.template getLocalView<host_execution_space> ();
379  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
380 
381  using val_type = typename decltype( X_lcl_j_1d )::value_type;
382  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
383  Kokkos::deep_copy (umavj, X_lcl_j_1d);
384  }
385  }
386 
387  typedef typename multivec_t::node_type::memory_space memory_space;
388  redist_mv.template sync <memory_space> ();
389 
390  mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
391  }
392  }
393 
394  }
395 
396 
397  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
398  std::string
399  MultiVecAdapter<
400  MultiVector<Scalar,
401  LocalOrdinal,
402  GlobalOrdinal,
403  Node> >::description() const
404  {
405  std::ostringstream oss;
406  oss << "Amesos2 adapter wrapping: ";
407  oss << mv_->description();
408  return oss.str();
409  }
410 
411 
412  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
413  void
414  MultiVecAdapter<
415  MultiVector<Scalar,
416  LocalOrdinal,
417  GlobalOrdinal,
418  Node> >::describe (Teuchos::FancyOStream& os,
419  const Teuchos::EVerbosityLevel verbLevel) const
420  {
421  mv_->describe (os, verbLevel);
422  }
423 
424 
425  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
426  const char* MultiVecAdapter<
427  MultiVector<Scalar,
428  LocalOrdinal,
429  GlobalOrdinal,
430  Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
431 
432 } // end namespace Amesos2
433 
434 #endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.
EDistribution
Definition: Amesos2_TypeDecl.hpp:123