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 
164  // Since we're caching the Export object, and since the Export
165  // needs to keep the distribution Map, we have to make a copy of
166  // the latter in order to ensure that it will stick around past
167  // the scope of this function call. (Ptr is not reference
168  // counted.) Map's clone() method suffices, even though it only
169  // makes a shallow copy of some of Map's data, because Map is
170  // immutable and those data are reference-counted (e.g.,
171  // ArrayRCP or RCP).
172  distMap = distribution_map->template clone<Node> (distribution_map->getNode ());
173 
174  // (Re)create the Export object.
175  exporter_ = rcp (new export_type (this->getMap (), distMap));
176  }
177  else {
178  distMap = exporter_->getTargetMap ();
179  }
180 
181  multivec_t redist_mv (distMap, num_vecs);
182 
183  // Redistribute the input (multi)vector.
184  redist_mv.doExport (*mv_, *exporter_, Tpetra::REPLACE);
185 
186  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
187  // Do this if GIDs contiguous - existing functionality
188  // Copy the imported (multi)vector's data into the ArrayView.
189  redist_mv.get1dCopy (av, lda);
190  }
191  else {
192  // Do this if GIDs not contiguous...
193  // sync is needed for example if mv was updated on device, but will be passed through Amesos2 to solver running on host
194  typedef typename multivec_t::dual_view_type dual_view_type;
195  typedef typename dual_view_type::host_mirror_space host_execution_space;
196  redist_mv.template sync < host_execution_space > ();
197 
198  auto contig_local_view_2d = redist_mv.template getLocalView<host_execution_space>();
199  if ( redist_mv.isConstantStride() ) {
200  for ( size_t j = 0; j < num_vecs; ++j) {
201  auto av_j = av(lda*j, lda);
202  for ( size_t i = 0; i < lda; ++i ) {
203  av_j[i] = contig_local_view_2d(i,j); //lda may not be correct if redist_mv is not constant stride...
204  }
205  }
206  }
207  else {
208  // ... lda should come from Teuchos::Array* allocation,
209  // not the MultiVector, since the MultiVector does NOT
210  // have constant stride in this case.
211  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
212  const size_t lclNumRows = redist_mv.getLocalLength();
213  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
214  auto av_j = av(lda*j, lclNumRows);
215  auto X_lcl_j_2d = redist_mv.template getLocalView<host_execution_space> ();
216  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
217 
218  using val_type = typename decltype( X_lcl_j_1d )::value_type;
219  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
220  Kokkos::deep_copy (umavj, X_lcl_j_1d);
221  }
222  }
223  }
224  }
225  }
226 
227  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
228  Teuchos::ArrayRCP<Scalar>
229  MultiVecAdapter<
230  MultiVector<Scalar,
231  LocalOrdinal,
232  GlobalOrdinal,
233  Node> >::get1dViewNonConst (bool local)
234  {
235  // FIXME (mfh 22 Jan 2014) When I first found this routine, all of
236  // its code was commented out, and it didn't return anything. The
237  // latter is ESPECIALLY dangerous, given that the return value is
238  // an ArrayRCP. Thus, I added the exception throw below.
239  TEUCHOS_TEST_FOR_EXCEPTION(
240  true, std::logic_error, "Amesos2::MultiVecAdapter::get1dViewNonConst: "
241  "Not implemented.");
242 
243  // if( local ){
244  // this->localize();
245  // /* Use the global element list returned by
246  // * mv_->getMap()->getNodeElementList() to get a subCopy of mv_ which we
247  // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
248  // */
249  // if(l_l_mv_.is_null() ){
250  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
251  // = mv_->getMap()->getNodeElementList();
252  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
253 
254  // // Convert the node element to a list of size_t type objects
255  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
256  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
257  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
258  // *(it_st++) = Teuchos::as<size_t>(*it_go);
259  // }
260 
261  // // To be consistent with the globalize() function, get a view of the local mv
262  // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
263 
264  // return(l_l_mv_->get1dViewNonConst());
265  // } else {
266  // // We need to re-import values to the local, since changes may have been
267  // // made to the global structure that are not reflected in the local
268  // // view.
269  // Teuchos::ArrayView<const GlobalOrdinal> nodeElements_go
270  // = mv_->getMap()->getNodeElementList();
271  // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
272 
273  // // Convert the node element to a list of size_t type objects
274  // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
275  // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
276  // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
277  // *(it_st++) = Teuchos::as<size_t>(*it_go);
278  // }
279 
280  // return l_l_mv_->get1dViewNonConst();
281  // }
282  // } else {
283  // if( mv_->isDistributed() ){
284  // this->localize();
285 
286  // return l_mv_->get1dViewNonConst();
287  // } else { // not distributed, no need to import
288  // return mv_->get1dViewNonConst();
289  // }
290  // }
291  }
292 
293 
294  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node>
295  void
296  MultiVecAdapter<
297  MultiVector<Scalar,
298  LocalOrdinal,
299  GlobalOrdinal,
300  Node> >::put1dData(const Teuchos::ArrayView<const scalar_t>& new_data,
301  size_t lda,
302  Teuchos::Ptr<
303  const Tpetra::Map<LocalOrdinal,
304  GlobalOrdinal,
305  Node> > source_map,
306  EDistribution distribution )
307  {
308  using Teuchos::RCP;
309  typedef Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> map_type;
310 
311  TEUCHOS_TEST_FOR_EXCEPTION(
312  source_map.getRawPtr () == NULL, std::invalid_argument,
313  "Amesos2::MultiVecAdapter::put1dData: source_map argument is null.");
314  TEUCHOS_TEST_FOR_EXCEPTION(
315  mv_.is_null (), std::logic_error,
316  "Amesos2::MultiVecAdapter::put1dData: the internal MultiVector mv_ is null.");
317  // getMap() calls mv_->getMap(), so test first whether mv_ is null.
318  TEUCHOS_TEST_FOR_EXCEPTION(
319  this->getMap ().is_null (), std::logic_error,
320  "Amesos2::MultiVecAdapter::put1dData: this->getMap() returns null.");
321 
322  const size_t num_vecs = getGlobalNumVectors ();
323 
324  // Special case when number vectors == 1 and single MPI process
325  if ( num_vecs == 1 && this->getComm()->getRank() == 0 && this->getComm()->getSize() == 1 ) {
326  typedef typename multivec_t::dual_view_type::host_mirror_space host_execution_space;
327  // num_vecs = 1; stride does not matter
328  auto mv_view_to_modify_2d = mv_->template getLocalView<host_execution_space>();
329  for ( size_t i = 0; i < lda; ++i ) {
330  mv_view_to_modify_2d(i,0) = new_data[i]; // Only one vector
331  }
332  }
333  else {
334 
335  // (Re)compute the Import object if necessary. If not, then we
336  // don't need to clone source_map; we can instead just get the
337  // previously cloned source Map from the Import object.
338  RCP<const map_type> srcMap;
339  if (importer_.is_null () ||
340  ! importer_->getSourceMap ()->isSameAs (* source_map) ||
341  ! importer_->getTargetMap ()->isSameAs (* (this->getMap ()))) {
342 
343  // Since we're caching the Import object, and since the Import
344  // needs to keep the source Map, we have to make a copy of the
345  // latter in order to ensure that it will stick around past the
346  // scope of this function call. (Ptr is not reference counted.)
347  // Map's clone() method suffices, even though it only makes a
348  // shallow copy of some of Map's data, because Map is immutable
349  // and those data are reference-counted (e.g., ArrayRCP or RCP).
350  srcMap = source_map->template clone<Node> (source_map->getNode ());
351  importer_ = rcp (new import_type (srcMap, this->getMap ()));
352  }
353  else {
354  srcMap = importer_->getSourceMap ();
355  }
356 
357  if ( distribution != CONTIGUOUS_AND_ROOTED ) {
358  // Do this if GIDs contiguous - existing functionality
359  // Redistribute the output (multi)vector.
360  const multivec_t source_mv (srcMap, new_data, lda, num_vecs);
361  mv_->doImport (source_mv, *importer_, Tpetra::REPLACE);
362  }
363  else {
364  multivec_t redist_mv (srcMap, num_vecs); // unused for ROOTED case
365  typedef typename multivec_t::dual_view_type dual_view_type;
366  typedef typename dual_view_type::host_mirror_space host_execution_space;
367  redist_mv.template modify< host_execution_space > ();
368 
369  if ( redist_mv.isConstantStride() ) {
370  auto contig_local_view_2d = redist_mv.template getLocalView<host_execution_space>();
371  for ( size_t j = 0; j < num_vecs; ++j) {
372  auto av_j = new_data(lda*j, lda);
373  for ( size_t i = 0; i < lda; ++i ) {
374  contig_local_view_2d(i,j) = av_j[i];
375  }
376  }
377  }
378  else {
379  // ... lda should come from Teuchos::Array* allocation,
380  // not the MultiVector, since the MultiVector does NOT
381  // have constant stride in this case.
382  // TODO lda comes from X->getGlobalLength() in solve_impl - should this be changed???
383  const size_t lclNumRows = redist_mv.getLocalLength();
384  for (size_t j = 0; j < redist_mv.getNumVectors(); ++j) {
385  auto av_j = new_data(lda*j, lclNumRows);
386  auto X_lcl_j_2d = redist_mv.template getLocalView<host_execution_space> ();
387  auto X_lcl_j_1d = Kokkos::subview (X_lcl_j_2d, Kokkos::ALL (), j);
388 
389  using val_type = typename decltype( X_lcl_j_1d )::value_type;
390  Kokkos::View<val_type*, Kokkos::HostSpace> umavj ( const_cast< val_type* > ( reinterpret_cast<const val_type*> ( av_j.getRawPtr () ) ), av_j.size () );
391  Kokkos::deep_copy (umavj, X_lcl_j_1d);
392  }
393  }
394 
395  typedef typename multivec_t::node_type::memory_space memory_space;
396  redist_mv.template sync <memory_space> ();
397 
398  mv_->doImport (redist_mv, *importer_, Tpetra::REPLACE);
399  }
400  }
401 
402  }
403 
404 
405  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
406  std::string
407  MultiVecAdapter<
408  MultiVector<Scalar,
409  LocalOrdinal,
410  GlobalOrdinal,
411  Node> >::description() const
412  {
413  std::ostringstream oss;
414  oss << "Amesos2 adapter wrapping: ";
415  oss << mv_->description();
416  return oss.str();
417  }
418 
419 
420  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
421  void
422  MultiVecAdapter<
423  MultiVector<Scalar,
424  LocalOrdinal,
425  GlobalOrdinal,
426  Node> >::describe (Teuchos::FancyOStream& os,
427  const Teuchos::EVerbosityLevel verbLevel) const
428  {
429  mv_->describe (os, verbLevel);
430  }
431 
432 
433  template <typename Scalar, typename LocalOrdinal, typename GlobalOrdinal, class Node >
434  const char* MultiVecAdapter<
435  MultiVector<Scalar,
436  LocalOrdinal,
437  GlobalOrdinal,
438  Node> >::name = "Amesos2 adapter for Tpetra::MultiVector";
439 
440 } // end namespace Amesos2
441 
442 #endif // AMESOS2_TPETRA_MULTIVEC_ADAPTER_DEF_HPP
Definition: Amesos2_AbstractConcreteMatrixAdapter.hpp:48
Amesos2::MultiVecAdapter specialization for the Tpetra::MultiVector class.
EDistribution
Definition: Amesos2_TypeDecl.hpp:123