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