Thyra  Version of the Day
Thyra_EpetraThyraWrappers.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
43 #include "Thyra_DefaultSpmdVectorSpace.hpp"
44 #include "Thyra_DefaultSpmdMultiVector.hpp"
45 #include "Thyra_DefaultSpmdVector.hpp"
46 #include "Thyra_DetachedVectorView.hpp"
47 #include "Thyra_DetachedMultiVectorView.hpp"
48 #include "Thyra_VectorSpaceFactoryBase.hpp"
49 #include "Thyra_ProductVectorSpaceBase.hpp"
50 #include "Teuchos_Assert.hpp"
51 #include "Teuchos_dyn_cast.hpp"
52 
53 #include "Teuchos_DefaultSerialComm.hpp"
54 #ifdef HAVE_MPI
55 #include "Teuchos_DefaultMpiComm.hpp"
56 #endif
57 
58 #include "Epetra_Map.h"
59 #include "Epetra_Comm.h"
60 #include "Epetra_SerialComm.h"
61 #ifdef HAVE_MPI
62 # include "Epetra_MpiComm.h"
63 #endif
64 #include "Epetra_MultiVector.h"
65 #include "Epetra_Vector.h"
66 
67 //
68 // Helpers
69 //
70 
71 
72 namespace {
73 
74 
75 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
76 unwrapSingleProductVectorSpace(
77  const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > &vs_in
78  )
79 {
80  using Teuchos::RCP;
81  using Teuchos::rcp_dynamic_cast;
83  const RCP<const ProductVectorSpaceBase<double> > pvs =
84  rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(vs_in);
85  if (nonnull(pvs)) {
86  TEUCHOS_ASSERT_EQUALITY( pvs->numBlocks(), 1 );
87  return pvs->getBlock(0);
88  }
89  return vs_in;
90 }
91 
92 
93 } // namespace
94 
95 
96 //
97 // Implementations of user function
98 //
99 
100 
101 Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> >
102 Thyra::create_Comm( const RCP<const Epetra_Comm> &epetraComm )
103 {
104  using Teuchos::rcp;
105  using Teuchos::rcp_dynamic_cast;
106  using Teuchos::set_extra_data;
107 
108  RCP<const Epetra_SerialComm>
109  serialEpetraComm = rcp_dynamic_cast<const Epetra_SerialComm>(epetraComm);
110  if( serialEpetraComm.get() ) {
111  RCP<const Teuchos::SerialComm<Ordinal> >
112  serialComm = rcp(new Teuchos::SerialComm<Ordinal>());
113  set_extra_data( serialEpetraComm, "serialEpetraComm", Teuchos::inOutArg(serialComm) );
114  return serialComm;
115  }
116 
117 #ifdef HAVE_MPI
118 
119  RCP<const Epetra_MpiComm>
120  mpiEpetraComm = rcp_dynamic_cast<const Epetra_MpiComm>(epetraComm);
121  if( mpiEpetraComm.get() ) {
122  RCP<const Teuchos::OpaqueWrapper<MPI_Comm> >
123  rawMpiComm = Teuchos::opaqueWrapper(mpiEpetraComm->Comm());
124  set_extra_data( mpiEpetraComm, "mpiEpetraComm", Teuchos::inOutArg(rawMpiComm) );
125  RCP<const Teuchos::MpiComm<Ordinal> >
126  mpiComm = rcp(new Teuchos::MpiComm<Ordinal>(rawMpiComm));
127  return mpiComm;
128  }
129 
130 #endif // HAVE_MPI
131 
132  // If you get here then the failed!
133  return Teuchos::null;
134 
135 }
136 
137 
138 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
140  const RCP<const Epetra_Map> &epetra_map
141  )
142 {
143  using Teuchos::as; using Teuchos::inoutArg;
144 #ifdef TEUCHOS_DEBUG
145  TEUCHOS_TEST_FOR_EXCEPTION(
146  !epetra_map.get(), std::invalid_argument,
147  "create_VectorSpace::initialize(...): Error!" );
148 #endif // TEUCHOS_DEBUG
149  RCP<const Teuchos::Comm<Ordinal> > comm =
150  create_Comm(Teuchos::rcpFromRef(epetra_map->Comm())).assert_not_null();
151  Teuchos::set_extra_data(epetra_map, "epetra_map", inoutArg(comm));
152  const Ordinal localSubDim = epetra_map->NumMyElements();
153  RCP<DefaultSpmdVectorSpace<double> > vs =
154  defaultSpmdVectorSpace<double>(
155  comm, localSubDim, epetra_map->NumGlobalElements64(),
156  !epetra_map->DistributedGlobal());
157  TEUCHOS_ASSERT_EQUALITY(vs->dim(), as<Ordinal>(epetra_map->NumGlobalElements64()));
158  // NOTE: the above assert just checks to make sure that the size of the
159  // Ordinal type can hold the size returned from NumGlobalElemenets64(). A
160  // 64 bit system will always have Ordinal=ptrdiff_t by default which will
161  // always be 64 bit so this should be fine. However, if Ordinal were
162  // defined to only be 32 bit and then this exception could trigger. Because
163  // this assert will only likely trigger in a non-debug build, we will leave
164  // the assert unguarded since it is very cheap to perform.
165  Teuchos::set_extra_data( epetra_map, "epetra_map", inoutArg(vs) );
166  return vs;
167 }
168 
169 
170 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
172  const RCP<const VectorSpaceBase<double> > &parentSpace,
173  const int dim
174  )
175 {
176  using Teuchos::rcp_dynamic_cast;
177 #ifdef TEUCHOS_DEBUG
178  TEUCHOS_TEST_FOR_EXCEPT(parentSpace.get()==NULL);
179  Teuchos::dyn_cast<const SpmdVectorSpaceBase<double> >(*parentSpace);
180  TEUCHOS_TEST_FOR_EXCEPT(dim <= 0);
181 #endif
182  return parentSpace->smallVecSpcFcty()->createVecSpc(dim);
183 }
184 
185 
186 Teuchos::RCP<Thyra::VectorBase<double> >
188  const RCP<Epetra_Vector> &epetra_v,
189  const RCP<const VectorSpaceBase<double> > &space_in
190  )
191 {
192 #ifdef TEUCHOS_DEBUG
193  TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
194 #endif
195  RCP<const SpmdVectorSpaceBase<double> >
196  space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
197  space_in,true);
198  if(!epetra_v.get())
199  return Teuchos::null;
200  // New local view of raw data
201  double *localValues;
202  epetra_v->ExtractView( &localValues );
203  // Build the Vector
204  RCP<SpmdVectorBase<double> >
205  v = Teuchos::rcp(
207  space,
208  Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
209  1
210  )
211  );
212  Teuchos::set_extra_data<RCP<Epetra_Vector> >( epetra_v, "Epetra_Vector",
213  Teuchos::inOutArg(v) );
214  return v;
215 }
216 
217 
218 Teuchos::RCP<const Thyra::VectorBase<double> >
220  const RCP<const Epetra_Vector> &epetra_v,
221  const RCP<const VectorSpaceBase<double> > &space_in
222  )
223 {
224 #ifdef TEUCHOS_DEBUG
225  TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
226 #endif
227  RCP<const SpmdVectorSpaceBase<double> >
228  space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
229  space_in,true);
230  if(!epetra_v.get())
231  return Teuchos::null;
232  // New local view of raw data
233  double *localValues;
234  epetra_v->ExtractView( &localValues );
235  // Build the Vector
236  RCP<const SpmdVectorBase<double> >
237  v = Teuchos::rcp(
239  space,
240  Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
241  1
242  )
243  );
244  Teuchos::set_extra_data<RCP<const Epetra_Vector> >( epetra_v, "Epetra_Vector",
245  Teuchos::inOutArg(v) );
246  return v;
247 }
248 
249 
250 Teuchos::RCP<Thyra::MultiVectorBase<double> >
252  const RCP<Epetra_MultiVector> &epetra_mv,
253  const RCP<const VectorSpaceBase<double> > &range_in,
254  const RCP<const VectorSpaceBase<double> > &domain_in
255  )
256 {
257  using Teuchos::rcp_dynamic_cast;
258 #ifdef TEUCHOS_DEBUG
259  TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
260 #endif
261  const RCP<const SpmdVectorSpaceBase<double> > range =
262  Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
263  unwrapSingleProductVectorSpace(range_in),
264  true
265  );
266  RCP<const ScalarProdVectorSpaceBase<double> > domain =
267  Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
268  unwrapSingleProductVectorSpace(domain_in),
269  true
270  );
271  if (!epetra_mv.get() )
272  return Teuchos::null;
273  if ( is_null(domain) ) {
274  domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
275  create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
276  );
277  }
278  // New local view of raw data
279  double *localValues; int leadingDim;
280  if( epetra_mv->ConstantStride() ) {
281  epetra_mv->ExtractView( &localValues, &leadingDim );
282  }
283  else {
284  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
285  }
286  // Build the MultiVector
287  RCP<SpmdMultiVectorBase<double> >
288  mv = Teuchos::rcp(
290  range,
291  domain,
292  Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
293  leadingDim
294  )
295  );
296  Teuchos::set_extra_data<RCP<Epetra_MultiVector> >(
297  epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
298  return mv;
299 }
300 
301 
302 Teuchos::RCP<const Thyra::MultiVectorBase<double> >
304  const RCP<const Epetra_MultiVector> &epetra_mv,
305  const RCP<const VectorSpaceBase<double> > &range_in,
306  const RCP<const VectorSpaceBase<double> > &domain_in
307  )
308 {
309  using Teuchos::rcp_dynamic_cast;
310 #ifdef TEUCHOS_DEBUG
311  TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
312 #endif
313  const RCP<const SpmdVectorSpaceBase<double> > range =
314  Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
315  unwrapSingleProductVectorSpace(range_in),
316  true
317  );
318  RCP<const ScalarProdVectorSpaceBase<double> > domain =
319  Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
320  unwrapSingleProductVectorSpace(domain_in),
321  true
322  );
323  if (!epetra_mv.get())
324  return Teuchos::null;
325  if ( is_null(domain) ) {
326  domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
327  create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
328  );
329  }
330  // New local view of raw data
331  double *localValues; int leadingDim;
332  if( epetra_mv->ConstantStride() ) {
333  epetra_mv->ExtractView( &localValues, &leadingDim );
334  }
335  else {
336  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
337  }
338  // Build the MultiVector
339  RCP<const SpmdMultiVectorBase<double> >
340  mv = Teuchos::rcp(
342  range,
343  domain,
344  Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
345  leadingDim
346  )
347  );
348  Teuchos::set_extra_data<RCP<const Epetra_MultiVector> >(
349  epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
350  return mv;
351 }
352 
353 
354 Teuchos::RCP<const Epetra_Comm>
355 Thyra::get_Epetra_Comm(const Teuchos::Comm<Ordinal>& comm_in)
356 {
357 
358  using Teuchos::rcp;
359  using Teuchos::ptrFromRef;
360  using Teuchos::ptr_dynamic_cast;
361  using Teuchos::SerialComm;
362 #ifdef HAVE_MPI
363  using Teuchos::MpiComm;
364 #endif
365 
366  const Ptr<const Teuchos::Comm<Ordinal> > comm = Teuchos::ptrFromRef(comm_in);
367 
368  const Ptr<const SerialComm<Ordinal> > serialComm =
369  ptr_dynamic_cast<const SerialComm<Ordinal> >(comm);
370 
371  RCP<const Epetra_Comm> epetraComm;
372 
373 #ifdef HAVE_MPI
374 
375  const Ptr<const MpiComm<Ordinal> > mpiComm =
376  ptr_dynamic_cast<const MpiComm<Ordinal> >(comm);
377 
378  TEUCHOS_TEST_FOR_EXCEPTION(is_null(mpiComm) && is_null(serialComm),
379  std::runtime_error,
380  "SPMD std::vector space has a communicator that is "
381  "neither a serial comm nor an MPI comm");
382 
383  if (nonnull(mpiComm)) {
384  epetraComm = rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()()));
385  }
386  else {
387  epetraComm = rcp(new Epetra_SerialComm());
388  }
389 
390 #else
391 
392  TEUCHOS_TEST_FOR_EXCEPTION(is_null(serialComm), std::runtime_error,
393  "SPMD std::vector space has a communicator that is "
394  "neither a serial comm nor an MPI comm");
395 
396  epetraComm = rcp(new Epetra_SerialComm());
397 
398 #endif
399 
400  TEUCHOS_TEST_FOR_EXCEPTION(is_null(epetraComm), std::runtime_error,
401  "null communicator created");
402 
403  return epetraComm;
404 
405 }
406 
407 
408 Teuchos::RCP<const Epetra_Map>
410  const RCP<const Epetra_Comm>& comm)
411 {
412 
413  using Teuchos::rcpFromRef;
414  using Teuchos::rcpFromPtr;
415  using Teuchos::rcp_dynamic_cast;
416  using Teuchos::ptrFromRef;
417  using Teuchos::ptr_dynamic_cast;
418 
419  const Ptr<const VectorSpaceBase<double> > vs_ptr = ptrFromRef(vs_in);
420 
421  const Ptr<const SpmdVectorSpaceBase<double> > spmd_vs =
422  ptr_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs_ptr);
423 
424  const Ptr<const ProductVectorSpaceBase<double> > &prod_vs =
425  ptr_dynamic_cast<const ProductVectorSpaceBase<double> >(vs_ptr);
426 
427  TEUCHOS_TEST_FOR_EXCEPTION( is_null(spmd_vs) && is_null(prod_vs), std::logic_error,
428  "Error, the concrete VectorSpaceBase object of type "
429  +Teuchos::demangleName(typeid(vs_in).name())+" does not support the"
430  " SpmdVectorSpaceBase or the ProductVectorSpaceBase interfaces!" );
431 
432  const int numBlocks = (nonnull(prod_vs) ? prod_vs->numBlocks() : 1);
433 
434  // Get an array of SpmdVectorBase objects for the blocks
435 
436  Array<RCP<const SpmdVectorSpaceBase<double> > > spmd_vs_blocks;
437  if (nonnull(prod_vs)) {
438  for (int block_i = 0; block_i < numBlocks; ++block_i) {
439  const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i =
440  rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
441  prod_vs->getBlock(block_i), true);
442  spmd_vs_blocks.push_back(spmd_vs_i);
443  }
444  }
445  else {
446  spmd_vs_blocks.push_back(rcpFromPtr(spmd_vs));
447  }
448 
449  // Find the number of local elements, summed over all blocks
450 
451  int myLocalElements = 0;
452  for (int block_i = 0; block_i < numBlocks; ++block_i) {
453  myLocalElements += spmd_vs_blocks[block_i]->localSubDim();
454  }
455 
456  // Find the GIDs owned by this processor, taken from all blocks
457 
458  int count=0;
459  int blockOffset = 0;
460  Array<int> myGIDs(myLocalElements);
461  for (int block_i = 0; block_i < numBlocks; ++block_i) {
462  const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i = spmd_vs_blocks[block_i];
463  const int lowGIDInBlock = spmd_vs_i->localOffset();
464  const int numLocalElementsInBlock = spmd_vs_i->localSubDim();
465  for (int i=0; i < numLocalElementsInBlock; ++i, ++count) {
466  myGIDs[count] = blockOffset + lowGIDInBlock + i;
467  }
468  blockOffset += spmd_vs_i->dim();
469  }
470 
471  const int globalDim = vs_in.dim();
472 
473  return Teuchos::rcp(
474  new Epetra_Map(globalDim, myLocalElements, myGIDs.getRawPtr(), 0, *comm));
475 
476 }
477 
478 
479 Teuchos::RCP<Epetra_Vector>
481  const Epetra_Map &map,
482  const RCP<VectorBase<double> > &v
483  )
484 {
485  using Teuchos::get_optional_extra_data;
486 #ifdef TEUCHOS_DEBUG
487  RCP<const VectorSpaceBase<double> >
488  epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
490  "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
491 #endif
492  //
493  // First, try to grab the Epetra_Vector straight out of the
494  // RCP since this is the fastest way.
495  //
496  const Ptr<const RCP<Epetra_Vector> >
497  epetra_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
498  v,"Epetra_Vector");
499  if(!is_null(epetra_v_ptr)) {
500  return *epetra_v_ptr;
501  }
502  //
503  // The assumption that we (rightly) make here is that if the vector spaces
504  // are compatible, that either the multi-vectors are both in-core or the
505  // vector spaces are both derived from SpmdVectorSpaceBase and have
506  // compatible maps.
507  //
508  const VectorSpaceBase<double> &vs = *v->range();
509  const SpmdVectorSpaceBase<double> *mpi_vs
510  = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
511  const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
512  const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
513  //
514  // Here we will extract a view of the local elements in the underlying
515  // VectorBase object. In most cases, no data will be allocated or copied
516  // and only some small objects will be created and a few virtual functions
517  // will be called so the overhead should be low and fixed.
518  //
519  // Create a *mutable* view of the local elements, this view will be set on
520  // the RCP that is returned. As a result, this view will be relased
521  // when the returned Epetra_Vector is released.
522  //
523  // Note that the input vector 'v' will be remembered through this detached
524  // view!
525  //
526  RCP<DetachedVectorView<double> >
527  emvv = Teuchos::rcp(
529  v
530  ,Range1D(localOffset,localOffset+localSubDim-1)
531  ,true // forceContiguous
532  )
533  );
534  // Create a temporary Epetra_Vector object and give it
535  // the above local view
536  RCP<Epetra_Vector>
537  epetra_v = Teuchos::rcp(
538  new Epetra_Vector(
539  ::View // CV
540  ,map // Map
541  ,const_cast<double*>(emvv->values()) // V
542  )
543  );
544  // Give the explict view object to the above Epetra_Vector smart pointer
545  // object. In this way, when the client is finished with the Epetra_Vector
546  // view the destructor from the object in emvv will automatically commit the
547  // changes to the elements in the input v VectorBase object (reguardless of
548  // its implementation). This is truly an elegant result!
549  Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_v),
550  Teuchos::PRE_DESTROY );
551  // We are done!
552  return epetra_v;
553 }
554 
555 
556 Teuchos::RCP<const Epetra_Vector>
558  const Epetra_Map &map,
559  const RCP<const VectorBase<double> > &v
560  )
561 {
562  using Teuchos::get_optional_extra_data;
563 #ifdef TEUCHOS_DEBUG
564  RCP<const VectorSpaceBase<double> >
565  epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
567  "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
568 #endif
569  //
570  // First, try to grab the Epetra_Vector straight out of the
571  // RCP since this is the fastest way.
572  //
573  const Ptr<const RCP<const Epetra_Vector> >
574  epetra_v_ptr = get_optional_extra_data<RCP<const Epetra_Vector> >(
575  v,"Epetra_Vector");
576  if(!is_null(epetra_v_ptr))
577  return *epetra_v_ptr;
578  const Ptr<const RCP<Epetra_Vector> >
579  epetra_nonconst_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
580  v,"Epetra_Vector");
581  if(!is_null(epetra_nonconst_v_ptr))
582  return *epetra_nonconst_v_ptr;
583  //
584  // Same as above function except as stated below
585  //
586  const VectorSpaceBase<double> &vs = *v->range();
587  const SpmdVectorSpaceBase<double> *mpi_vs
588  = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
589  const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
590  const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
591  // Get an explicit *non-mutable* view of all of the elements in the multi
592  // vector. Note that 'v' will be remembered by this view!
593  RCP<ConstDetachedVectorView<double> >
594  evv = Teuchos::rcp(
596  v
597  ,Range1D(localOffset,localOffset+localSubDim-1)
598  ,true // forceContiguous
599  )
600  );
601  // Create a temporary Epetra_Vector object and give it
602  // the above view
603  RCP<Epetra_Vector>
604  epetra_v = Teuchos::rcp(
605  new Epetra_Vector(
606  ::View // CV
607  ,map // Map
608  ,const_cast<double*>(evv->values()) // V
609  )
610  );
611  // This next statement will cause the destructor to free the view if
612  // needed (see above function). Since this view is non-mutable,
613  // only a releaseDetachedView(...) and not a commit will be called.
614  // This is the whole reason there is a seperate implementation for
615  // the const and non-const cases.
616  Teuchos::set_extra_data( evv, "evv", Teuchos::inOutArg(epetra_v),
617  Teuchos::PRE_DESTROY );
618  // We are done!
619  return epetra_v;
620 }
621 
622 
623 Teuchos::RCP<Epetra_MultiVector>
625  const Epetra_Map &map,
626  const RCP<MultiVectorBase<double> > &mv
627  )
628 {
629  using Teuchos::get_optional_extra_data;
630 #ifdef TEUCHOS_DEBUG
631  RCP<const VectorSpaceBase<double> >
632  epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
634  "Thyra::get_Epetra_MultiVector(map,mv)", *epetra_vs, *mv->range() );
635 #endif
636  //
637  // First, try to grab the Epetra_MultiVector straight out of the
638  // RCP since this is the fastest way.
639  //
640  const Ptr<const RCP<Epetra_MultiVector> >
641  epetra_mv_ptr = get_optional_extra_data<RCP<Epetra_MultiVector> >(
642  mv,"Epetra_MultiVector");
643  if(!is_null(epetra_mv_ptr)) {
644  return *epetra_mv_ptr;
645  }
646  //
647  // The assumption that we (rightly) make here is that if the vector spaces
648  // are compatible, that either the multi-vectors are both in-core or the
649  // vector spaces are both derived from SpmdVectorSpaceBase and have
650  // compatible maps.
651  //
652  const VectorSpaceBase<double> &vs = *mv->range();
653  const SpmdVectorSpaceBase<double> *mpi_vs
654  = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
655  const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
656  const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
657  //
658  // Here we will extract a view of the local elements in the underlying
659  // MultiVectorBase object. In most cases, no data will be allocated or
660  // copied and only some small objects will be created and a few virtual
661  // functions will be called so the overhead should be low and fixed.
662  //
663  // Create a *mutable* view of the local elements, this view will be set on
664  // the RCP that is returned. As a result, this view will be relased
665  // when the returned Epetra_MultiVector is released.
666  //
667  RCP<DetachedMultiVectorView<double> >
668  emmvv = Teuchos::rcp(
670  *mv
671  ,Range1D(localOffset,localOffset+localSubDim-1)
672  )
673  );
674  // Create a temporary Epetra_MultiVector object and give it
675  // the above local view
676  RCP<Epetra_MultiVector>
677  epetra_mv = Teuchos::rcp(
678  new Epetra_MultiVector(
679  ::View // CV
680  ,map // Map
681  ,const_cast<double*>(emmvv->values()) // A
682  ,emmvv->leadingDim() // MyLDA
683  ,emmvv->numSubCols() // NumVectors
684  )
685  );
686  // Give the explict view object to the above Epetra_MultiVector
687  // smart pointer object. In this way, when the client is finished
688  // with the Epetra_MultiVector view the destructor from the object
689  // in emmvv will automatically commit the changes to the elements in
690  // the input mv MultiVectorBase object (reguardless of its
691  // implementation). This is truly an elegant result!
692  Teuchos::set_extra_data( emmvv, "emmvv", Teuchos::inOutArg(epetra_mv),
693  Teuchos::PRE_DESTROY );
694  // Also set the mv itself as extra data just to be safe
695  Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
696  // We are done!
697  return epetra_mv;
698 }
699 
700 
701 Teuchos::RCP<const Epetra_MultiVector>
703  const Epetra_Map &map,
704  const RCP<const MultiVectorBase<double> > &mv
705  )
706 {
707  using Teuchos::get_optional_extra_data;
708 #ifdef TEUCHOS_DEBUG
709  RCP<const VectorSpaceBase<double> >
710  epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
712  "Thyra::get_Epetra_MultiVector(map,mv)",
713  *epetra_vs, *mv->range() );
714 #endif
715  //
716  // First, try to grab the Epetra_MultiVector straight out of the
717  // RCP since this is the fastest way.
718  //
719  const Ptr<const RCP<const Epetra_MultiVector> >
720  epetra_mv_ptr
721  = get_optional_extra_data<RCP<const Epetra_MultiVector> >(
722  mv,"Epetra_MultiVector" );
723  if(!is_null(epetra_mv_ptr)) {
724  return *epetra_mv_ptr;
725  }
726  //
727  // Same as above function except as stated below
728  //
729  const VectorSpaceBase<double> &vs = *mv->range();
730  const SpmdVectorSpaceBase<double> *mpi_vs
731  = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
732  const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
733  const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
734  // Get an explicit *non-mutable* view of all of the elements in
735  // the multi vector.
736  RCP<ConstDetachedMultiVectorView<double> >
737  emvv = Teuchos::rcp(
739  *mv
740  ,Range1D(localOffset,localOffset+localSubDim-1)
741  )
742  );
743  // Create a temporary Epetra_MultiVector object and give it
744  // the above view
745  RCP<Epetra_MultiVector>
746  epetra_mv = Teuchos::rcp(
747  new Epetra_MultiVector(
748  ::View // CV
749  ,map // Map
750  ,const_cast<double*>(emvv->values()) // A
751  ,emvv->leadingDim() // MyLDA
752  ,emvv->numSubCols() // NumVectors
753  )
754  );
755  // This next statement will cause the destructor to free the view if
756  // needed (see above function). Since this view is non-mutable,
757  // only a releaseDetachedView(...) and not a commit will be called.
758  // This is the whole reason there is a seperate implementation for
759  // the const and non-const cases.
760  Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_mv),
761  Teuchos::PRE_DESTROY );
762  // Also set the mv itself as extra data just to be safe
763  Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
764  // We are done!
765  return epetra_mv;
766 }
767 
768 
769 Teuchos::RCP<Epetra_MultiVector>
771  const Epetra_Map &map,
773  )
774 {
775  using Teuchos::rcpWithEmbeddedObj;
776  using Teuchos::rcpFromRef;
777  using Teuchos::ptrFromRef;
778  using Teuchos::ptr_dynamic_cast;
779  using Teuchos::outArg;
780  ArrayRCP<double> mvData;
781  Ordinal mvLeadingDim = -1;
782  Ptr<SpmdMultiVectorBase<double> > mvSpmdMv;
783  if (nonnull(mvSpmdMv = ptr_dynamic_cast<SpmdMultiVectorBase<double> >(ptrFromRef(mv)))) {
784  mvSpmdMv->getNonconstLocalData(outArg(mvData), outArg(mvLeadingDim));
785  }
786  if (nonnull(mvData)) {
787  return rcpWithEmbeddedObj(
788  new Epetra_MultiVector(
789  ::View, map, mvData.getRawPtr(), mvLeadingDim, mv.domain()->dim()
790  ),
791  mvData
792  );
793  }
794  return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
795 }
796 
797 
798 Teuchos::RCP<const Epetra_MultiVector>
800  const Epetra_Map &map,
801  const MultiVectorBase<double> &mv
802  )
803 {
804  using Teuchos::rcpWithEmbeddedObj;
805  using Teuchos::rcpFromRef;
806  using Teuchos::ptrFromRef;
807  using Teuchos::ptr_dynamic_cast;
808  using Teuchos::outArg;
809  ArrayRCP<const double> mvData;
810  Ordinal mvLeadingDim = -1;
811  Ptr<const SpmdMultiVectorBase<double> > mvSpmdMv;
812  if (nonnull(mvSpmdMv = ptr_dynamic_cast<const SpmdMultiVectorBase<double> >(ptrFromRef(mv)))) {
813  mvSpmdMv->getLocalData(outArg(mvData), outArg(mvLeadingDim));
814  }
815  if (nonnull(mvData)) {
816  return rcpWithEmbeddedObj(
817  new Epetra_MultiVector(
818  ::View,map, const_cast<double*>(mvData.getRawPtr()), mvLeadingDim, mv.domain()->dim()
819  ),
820  mvData
821  );
822  }
823  return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
824 }
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible...
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible...
Create an explicit non-mutable (const) view of a MultiVectorBase object.
virtual Teuchos::RCP< const VectorSpaceBase< Scalar > > getBlock(const int k) const =0
Returns a vector space for the kth (zero-based) block.
Base interface class for SPMD multi-vectors.
RCP< const Epetra_Map > get_Epetra_Map(const VectorSpaceBase< double > &vs, const RCP< const Epetra_Comm > &comm)
Get (or create) an Epetra_Map object given an VectorSpaceBase object an optionally an extra Epetra_Co...
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
Create an explicit non-mutable (const) view of a VectorBase object.
Abstract interface for objects that represent a space for vectors.
RCP< MultiVectorBase< double > > create_MultiVector(const RCP< Epetra_MultiVector > &epetra_mv, const RCP< const VectorSpaceBase< double > > &range, const RCP< const VectorSpaceBase< double > > &domain=Teuchos::null)
Create a non-const MultiVectorBase object from a non-const Epetra_MultiVector object.
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
Create an explicit mutable (non-const) view of a MultiVectorBase object.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Interface for a collection of column vectors called a multi-vector.
Create an explicit mutable (non-const) view of a VectorBase object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
Abstract interface for finite-dimensional dense vectors.
Efficient concrete implementation subclass for SPMD vectors.
RCP< const Epetra_Comm > get_Epetra_Comm(const Teuchos::Comm< Ordinal > &comm)
Get (or create) and Epetra_Comm given a Teuchos::Comm object.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
virtual Ordinal localOffset() const =0
Returns the offset for the local sub-vector stored on this process.
Efficient concrete implementation subclass for SPMD multi-vectors.
RCP< const VectorSpaceBase< double > > create_LocallyReplicatedVectorSpace(const RCP< const VectorSpaceBase< double > > &parentSpace, const int dim)
Create a VectorSpaceBase object that creates locally replicated vector objects.
RCP< const Teuchos::Comm< Ordinal > > create_Comm(const RCP< const Epetra_Comm > &epetraComm)
Given an Epetra_Comm object, return an equivalent Teuchos::Comm object.
Base abstract VectorSpaceBase class for all SPMD-based vector spaces.
virtual Ordinal dim() const =0
Return the dimension of the vector space.
Teuchos::Range1D Range1D
virtual Ordinal localSubDim() const =0
Returns the number of local elements stored on this process.