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
35namespace Amesos2 {
36
37MultiVecAdapter<Epetra_MultiVector>::MultiVecAdapter(const MultiVecAdapter<multivec_t>& adapter)
38 : mv_(adapter.mv_)
39 , mv_map_(adapter.mv_map_)
40{ }
41
42MultiVecAdapter<Epetra_MultiVector>::MultiVecAdapter(const Teuchos::RCP<multivec_t>& m)
43 : mv_(m)
44{
45 mv_map_ = Teuchos::rcpFromRef(mv_->Map());
46}
47
48
49Teuchos::RCP<Epetra_MultiVector>
50MultiVecAdapter<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
58bool MultiVecAdapter<Epetra_MultiVector>::isLocallyIndexed() const
59{
60 return !mv_->DistributedGlobal();
61}
62
63bool MultiVecAdapter<Epetra_MultiVector>::isGloballyIndexed() const
64{
65 return mv_->DistributedGlobal();
66}
67
68
69const Teuchos::RCP<const Teuchos::Comm<int> >
70MultiVecAdapter<Epetra_MultiVector>::getComm() const
71{
72 return Util::to_teuchos_comm(Teuchos::rcpFromRef(mv_->Comm()));
73}
74
75
76size_t MultiVecAdapter<Epetra_MultiVector>::getLocalLength() const
77{
78 return Teuchos::as<size_t>(mv_->MyLength());
79}
80
81
82size_t MultiVecAdapter<Epetra_MultiVector>::getLocalNumVectors() const
83{
84 return Teuchos::as<size_t>(mv_->NumVectors());
85}
86
87
88MultiVecAdapter<Epetra_MultiVector>::global_size_t
89MultiVecAdapter<Epetra_MultiVector>::getGlobalLength() const
90{
91 return Teuchos::as<global_size_t>(mv_->GlobalLength());
92}
93
94
95size_t MultiVecAdapter<Epetra_MultiVector>::getGlobalNumVectors() const
96{
97 return Teuchos::as<size_t>(mv_->NumVectors());
98}
99
100
101size_t MultiVecAdapter<Epetra_MultiVector>::getStride() const
102{
103 return Teuchos::as<size_t>(mv_->Stride());
104}
105
106
107bool MultiVecAdapter<Epetra_MultiVector>::isConstantStride() const
108{
109 return mv_->ConstantStride();
110}
111
112
113Teuchos::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> >
117MultiVecAdapter<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()
145Teuchos::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> >
149MultiVecAdapter<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
176double * 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
187void 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
233void 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
282Teuchos::ArrayRCP<MultiVecAdapter<Epetra_MultiVector>::scalar_t>
283MultiVecAdapter<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
338void
339MultiVecAdapter<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
377void
378MultiVecAdapter<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
411std::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
419void 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
431Teuchos::RCP<const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
432 MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
433 MultiVecAdapter<Epetra_MultiVector>::node_t> >
434MultiVecAdapter<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
440const 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
Amesos2_EpetraMultiVecAdapter_decl.hpp
Amesos2::MultiVecAdapter specialization for the Epetra_MultiVector class.
Amesos2::EDistribution
EDistribution
Definition Amesos2_TypeDecl.hpp:89
Amesos2_Util.hpp
Utility functions for Amesos2.
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:743
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.
size
const int size
Definition klu2_simple.cpp:50
Amesos2::MultiVecAdapter
A templated MultiVector class adapter for Amesos2.
Definition Amesos2_MultiVecAdapter_decl.hpp:142
Generated by doxygen 1.9.8