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