Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
Amesos2_EpetraMultiVecAdapter_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_EPETRA_MULTIVEC_ADAPTER_DEF_HPP
54 #define AMESOS2_EPETRA_MULTIVEC_ADAPTER_DEF_HPP
55 
56 #include <Teuchos_as.hpp>
57 
58 #include <Epetra_SerialComm.h>
59 #ifdef HAVE_MPI
60 #include <Epetra_MpiComm.h>
61 #endif
62 #include <Epetra_LocalMap.h>
63 #include <Epetra_Import.h>
64 #include <Epetra_Export.h>
65 
67 #include "Amesos2_Util.hpp"
68 
69 namespace Amesos2 {
70 
72  : mv_(adapter.mv_)
73  , mv_map_(adapter.mv_map_)
74 { }
75 
76 MultiVecAdapter<Epetra_MultiVector>::MultiVecAdapter(const Teuchos::RCP<multivec_t>& m)
77  : mv_(m)
78 {
79  mv_map_ = Teuchos::rcpFromRef(mv_->Map());
80 }
81 
82 
84 {
85  return !mv_->DistributedGlobal();
86 }
87 
89 {
90  return mv_->DistributedGlobal();
91 }
92 
93 
94 const Teuchos::RCP<const Teuchos::Comm<int> >
96 {
97  return Util::to_teuchos_comm(Teuchos::rcpFromRef(mv_->Comm()));
98 }
99 
100 
102 {
103  return Teuchos::as<size_t>(mv_->MyLength());
104 }
105 
106 
108 {
109  return Teuchos::as<size_t>(mv_->NumVectors());
110 }
111 
112 
115 {
116  return Teuchos::as<global_size_t>(mv_->GlobalLength());
117 }
118 
119 
121 {
122  return Teuchos::as<size_t>(mv_->NumVectors());
123 }
124 
125 
127 {
128  return Teuchos::as<size_t>(mv_->Stride());
129 }
130 
131 
133 {
134  return mv_->ConstantStride();
135 }
136 
137 
138 Teuchos::RCP<const Tpetra::Vector<MultiVecAdapter<Epetra_MultiVector>::scalar_t,
143 {
144  using Teuchos::RCP;
145  using Teuchos::rcp;
146  using Teuchos::ArrayRCP;
147  using Tpetra::MultiVector;
148 
149  typedef scalar_t st;
150  typedef local_ordinal_t lot;
151  typedef global_ordinal_t got;
152  typedef node_t nt;
153 
154  RCP<MultiVector<st,lot,got,nt> > vector = rcp(new MultiVector<st,lot,got,nt>(this->getMap(),1));
155 
156  // Copy vector contents into Tpetra multi-vector
157  ArrayRCP<st> it = vector->getDataNonConst(0);
158  double* vector_data = mv_->operator[](Teuchos::as<int>(j)); // values from j^th vector
159  Tpetra::global_size_t size = vector->getGlobalLength();
160 
161  for( Tpetra::global_size_t i = 0; i < size; ++i ){
162  *it = vector_data[i];
163  }
164 
165  return vector->getVector(j);
166 }
167 
168 
169 // Implementation is essentially the same as getVector()
170 Teuchos::RCP<Tpetra::Vector<MultiVecAdapter<Epetra_MultiVector>::scalar_t,
175 {
176  using Teuchos::RCP;
177  using Teuchos::rcp;
178  using Teuchos::ArrayRCP;
179  using Tpetra::MultiVector;
180 
181  typedef scalar_t st;
182  typedef local_ordinal_t lot;
183  typedef global_ordinal_t got;
184  typedef node_t nt;
185 
186  RCP<MultiVector<st,lot,got,nt> > vector = rcp(new MultiVector<st,lot,got,nt>(this->getMap(),1));
187 
188  // Copy vector contents into Tpetra multi-vector
189  ArrayRCP<st> it = vector->getDataNonConst(0);
190  double* vector_data = mv_->operator[](Teuchos::as<int>(j)); // values from j^th vector
191  Tpetra::global_size_t size = vector->getGlobalLength();
192 
193  for( Tpetra::global_size_t i = 0; i < size; ++i ){
194  *it = vector_data[i];
195  }
196 
197  return vector->getVectorNonConst(j);
198 }
199 
200 
202 {
203  TEUCHOS_TEST_FOR_EXCEPTION( this->getGlobalNumVectors() != 1,
204  std::invalid_argument,
205  "Amesos2_EpetraMultiVectorAdapter: getMVPointer_impl should only be called for case with a single vector and single MPI process" );
206 
207  double* vector_data = mv_->operator[](Teuchos::as<int>(0)); // raw pointer to data from 0^th vector
208  return vector_data;
209 }
210 
211 
213  const Teuchos::ArrayView<MultiVecAdapter<Epetra_MultiVector>::scalar_t>& av,
214  size_t lda,
215  Teuchos::Ptr<
219  EDistribution /* distribution */) const
220 {
221  using Teuchos::rcpFromPtr;
222  using Teuchos::as;
223 
224  const size_t num_vecs = getGlobalNumVectors();
225 
226 #ifdef HAVE_AMESOS2_DEBUG
227  const size_t requested_vector_length = distribution_map->getNodeNumElements();
228  TEUCHOS_TEST_FOR_EXCEPTION( lda < requested_vector_length,
229  std::invalid_argument,
230  "Given stride is not large enough for local vector length" );
231  TEUCHOS_TEST_FOR_EXCEPTION( as<size_t>(av.size()) < (num_vecs-1) * lda + requested_vector_length,
232  std::invalid_argument,
233  "MultiVector storage not large enough given leading dimension "
234  "and number of vectors" );
235 #endif
236 
237  // Optimization for ROOTED and single MPI process
238  if ( num_vecs == 1 && mv_->Comm().MyPID() == 0 && mv_->Comm().NumProc() == 1 ) {
239  mv_->ExtractCopy(av.getRawPtr(), lda);
240  }
241  else {
242  Epetra_Map e_dist_map
243  = *Util::tpetra_map_to_epetra_map<local_ordinal_t,
244  global_ordinal_t,
245  global_size_t,
246  node_t>(*distribution_map);
247 
248  multivec_t redist_mv(e_dist_map, as<int>(num_vecs));
249  const Epetra_Import importer(e_dist_map, *mv_map_); // Note, target/source order is reversed in Tpetra
250  redist_mv.Import(*mv_, importer, Insert);
251 
252  // Finally, do copy
253  redist_mv.ExtractCopy(av.getRawPtr(), lda);
254  }
255 
256 }
257 
259  Kokkos::View<scalar_t**, Kokkos::LayoutLeft, Kokkos::HostSpace> & host_view,
260  size_t lda,
261  Teuchos::Ptr<
265  EDistribution /* distribution */) const
266 {
267  using Teuchos::rcpFromPtr;
268  using Teuchos::as;
269 
270  const size_t num_vecs = getGlobalNumVectors();
271 
272  #ifdef HAVE_AMESOS2_DEBUG
273  const size_t requested_vector_length = distribution_map->getNodeNumElements();
274  TEUCHOS_TEST_FOR_EXCEPTION( lda < requested_vector_length,
275  std::invalid_argument,
276  "Given stride is not large enough for local vector length" );
277  #endif
278 
279  // First make a host view
280  host_view = Kokkos::View<scalar_t**, Kokkos::LayoutLeft, Kokkos::HostSpace>(
281  Kokkos::ViewAllocateWithoutInitializing("get1dCopy_kokkos_view"),
282  distribution_map->getNodeNumElements(), num_vecs);
283 
284  // Optimization for ROOTED and single MPI process
285  if ( num_vecs == 1 && this->mv_->Comm().MyPID() == 0 && this->mv_->Comm().NumProc() == 1 ) {
286  mv_->ExtractCopy(host_view.data(), lda);
287  }
288  else {
289  Epetra_Map e_dist_map
290  = *Util::tpetra_map_to_epetra_map<local_ordinal_t,
291  global_ordinal_t,
292  global_size_t,
293  node_t>(*distribution_map);
294 
295  multivec_t redist_mv(e_dist_map, as<int>(num_vecs));
296  const Epetra_Import importer(e_dist_map, *mv_map_); // Note, target/source order is reversed in Tpetra
297  redist_mv.Import(*mv_, importer, Insert);
298 
299  // Finally, access data
300 
301  // Can we consider direct ptr usage with ExtractView?
302  // For now I will just copy - this was discussed as low priority for now.
303  redist_mv.ExtractCopy(host_view.data(), lda);
304  }
305 }
306 
307 Teuchos::ArrayRCP<MultiVecAdapter<Epetra_MultiVector>::scalar_t>
309 {
310  ((void) local);
311  // TEUCHOS_TEST_FOR_EXCEPTION( !this->isConstantStride(),
312  // std::logic_error,
313  // "get1dViewNonConst() : can only get 1d view if stride is constant");
314 
315  // if( local ){
316  // TEUCHOS_TEST_FOR_EXCEPTION(
317  // true,
318  // std::logic_error,
319  // "Amesos2::MultiVecAdapter<Epetra_MultiVector> : 1D views not yet supported for local-local Epetra multi-vectors");
320 
321  // // localize();
322  // // /* Use the global element list returned by
323  // // * mv_->getMap()->getNodeElementList() to get a subCopy of mv_ which we
324  // // * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
325  // // */
326  // // l_l_mv_ = Teuchos::null;
327 
328  // // Teuchos::Array<GlobalOrdinal> nodeElements_go(mv_->Map().NumMyElements());
329  // // mv_->Map().MyGlobalElements(nodeElements_go.getRawPtr());
330  // // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
331 
332  // // // Convert the node element to a list of size_t type objects
333  // // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
334  // // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
335  // // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
336  // // *(it_st++) = Teuchos::as<size_t>(*it_go);
337  // // }
338 
339  // // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
340 
341  // // return(l_l_mv_->get1dViewNonConst());
342  // } else {
343  // scalar_t* values;
344  // int lda;
345 
346  // if( !isLocal() ){
347  // this->localize();
348  // l_mv_->ExtractView(&values, &lda);
349  // } else {
350  // mv_->ExtractView(&values, &lda);
351  // }
352 
353  // TEUCHOS_TEST_FOR_EXCEPTION( lda != Teuchos::as<int>(this->getStride()),
354  // std::logic_error,
355  // "Stride reported during extraction not consistent with what multivector reports");
356 
357  // return Teuchos::arcp(values,0,lda*this->getGlobalNumVectors(),false);
358  // }
359  return Teuchos::null;
360 }
361 
362 
363 void
365  const Teuchos::ArrayView<const MultiVecAdapter<Epetra_MultiVector>::scalar_t>& new_data,
366  size_t lda,
367  Teuchos::Ptr<
371  EDistribution /* distribution */)
372 {
373  using Teuchos::rcpFromPtr;
374  using Teuchos::as;
375 
376  const size_t num_vecs = getGlobalNumVectors();
377  // TODO: check that the following const_cast is safe
378  double* data_ptr = const_cast<double*>(new_data.getRawPtr());
379 
380  // Optimization for ROOTED and single MPI process
381  if ( num_vecs == 1 && mv_->Comm().MyPID() == 0 && mv_->Comm().NumProc() == 1 ) {
382  // First, functioning impl
383  //const multivec_t source_mv(Copy, *mv_map_, data_ptr, as<int>(lda), as<int>(num_vecs));
384  //const Epetra_Import importer(*mv_map_, *mv_map_); //trivial - map does not change
385  //mv_->Import(source_mv, importer, Insert);
386  // Element-wise copy rather than using importer
387  auto vector = mv_->Pointers();
388  for ( size_t i = 0; i < lda; ++i ) {
389  vector[0][i] = data_ptr[i];
390  }
391  }
392  else {
393  const Epetra_BlockMap e_source_map
394  = *Util::tpetra_map_to_epetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*source_map);
395  const multivec_t source_mv(Copy, e_source_map, data_ptr, as<int>(lda), as<int>(num_vecs));
396  const Epetra_Import importer(*mv_map_, e_source_map);
397 
398  mv_->Import(source_mv, importer, Insert);
399  }
400 }
401 
402 void
404  Kokkos::View<scalar_t**, Kokkos::LayoutLeft, Kokkos::HostSpace> & host_new_data,
405  size_t lda,
406  Teuchos::Ptr<
410  EDistribution /* distribution */)
411 {
412  using Teuchos::rcpFromPtr;
413  using Teuchos::as;
414 
415  const size_t num_vecs = getGlobalNumVectors();
416 
417  double* data_ptr = host_new_data.data();
418 
419  // Optimization for ROOTED and single MPI process
420  if ( num_vecs == 1 && mv_->Comm().MyPID() == 0 && mv_->Comm().NumProc() == 1 ) {
421  auto vector = mv_->Pointers();
422  for ( size_t i = 0; i < lda; ++i ) {
423  vector[0][i] = data_ptr[i];
424  }
425  }
426  else {
427  const Epetra_BlockMap e_source_map
428  = *Util::tpetra_map_to_epetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*source_map);
429  const multivec_t source_mv(Copy, e_source_map, data_ptr, as<int>(lda), as<int>(num_vecs));
430  const Epetra_Import importer(*mv_map_, e_source_map);
431 
432  mv_->Import(source_mv, importer, Insert);
433  }
434 }
435 
437 {
438  std::ostringstream oss;
439  oss << "Amesos2 adapter wrapping: Epetra_MultiVector";
440  return oss.str();
441 }
442 
443 
445  Teuchos::FancyOStream& os,
446  const Teuchos::EVerbosityLevel verbLevel) const
447 {
448  // TODO: implement!
449  if(verbLevel != Teuchos::VERB_NONE)
450  {
451  os << "TODO: implement! ";
452  }
453 }
454 
455 
456 Teuchos::RCP<const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
460 {
461  return Util::epetra_map_to_tpetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*mv_map_);
462 }
463 
464 
466 = "Amesos2 adapter for Epetra_MultiVector";
467 
468 
469 } // end namespace Amesos2
470 
471 #endif // AMESOS2_EPETRA_MULTIVEC_ADAPTER_DEF_HPP
Utility functions for Amesos2.
const RCP< const Teuchos::Comm< int > > to_teuchos_comm(RCP< const Epetra_Comm > c)
Transform an Epetra_Comm object into a Teuchos::Comm object.
Amesos2::MultiVecAdapter specialization for the Epetra_MultiVector class.
RCP< Epetra_Map > tpetra_map_to_epetra_map(const Tpetra::Map< LO, GO, Node > &map)
Transform a Tpetra::Map object into an Epetra_Map.
Definition: Amesos2_Util.hpp:884
EDistribution
Definition: Amesos2_TypeDecl.hpp:123
A templated MultiVector class adapter for Amesos2.
Definition: Amesos2_MultiVecAdapter_decl.hpp:176