Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_TpetraMultiVector_def.hpp
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 
42 #ifndef THYRA_TPETRA_MULTIVECTOR_HPP
43 #define THYRA_TPETRA_MULTIVECTOR_HPP
44 
45 #include "Thyra_TpetraMultiVector_decl.hpp"
46 #include "Thyra_TpetraVectorSpace.hpp"
47 #include "Thyra_TpetraVector.hpp"
48 #include "Teuchos_Assert.hpp"
49 
50 
51 namespace Thyra {
52 
53 
54 // Constructors/initializers/accessors
55 
56 
57 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
59 {}
60 
61 
62 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
64  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
65  const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace,
66  const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraMultiVector
67  )
68 {
69  initializeImpl(tpetraVectorSpace, domainSpace, tpetraMultiVector);
70 }
71 
72 
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
75  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
76  const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace,
77  const RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraMultiVector
78  )
79 {
80  initializeImpl(tpetraVectorSpace, domainSpace, tpetraMultiVector);
81 }
82 
83 
84 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87 {
88  return tpetraMultiVector_.getNonconstObj();
89 }
90 
91 
92 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
95 {
96  return tpetraMultiVector_;
97 }
98 
99 
100 // Overridden public functions form MultiVectorAdapterBase
101 
102 
103 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
106 {
107  return domainSpace_;
108 }
109 
110 
111 // Overridden protected functions from MultiVectorBase
112 
113 
114 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
115 void
117 {
118  tpetraMultiVector_.getNonconstObj()->putScalar(alpha);
119 }
120 
121 
122 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125 {
126  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
127 
128  // If the cast succeeded, call Tpetra directly.
129  // Otherwise, fall back to the RTOp implementation.
130  if (nonnull(tmv)) {
131  tpetraMultiVector_.getNonconstObj()->assign(*tmv);
132  } else {
133  // This version will require/modify the host view of this vector.
134  tpetraMultiVector_.getNonconstObj()->sync_host ();
135  tpetraMultiVector_.getNonconstObj()->modify_host ();
137  }
138 }
139 
140 
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 void
144 {
145  tpetraMultiVector_.getNonconstObj()->scale(alpha);
146 }
147 
148 
149 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
151  Scalar alpha,
152  const MultiVectorBase<Scalar>& mv
153  )
154 {
155  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
156 
157  // If the cast succeeded, call Tpetra directly.
158  // Otherwise, fall back to the RTOp implementation.
159  if (nonnull(tmv)) {
161  tpetraMultiVector_.getNonconstObj()->update(alpha, *tmv, ST::one());
162  } else {
163  // This version will require/modify the host view of this vector.
164  tpetraMultiVector_.getNonconstObj()->sync_host ();
165  tpetraMultiVector_.getNonconstObj()->modify_host ();
167  }
168 }
169 
170 
171 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
173  const ArrayView<const Scalar>& alpha,
174  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
175  const Scalar& beta
176  )
177 {
178 #ifdef TEUCHOS_DEBUG
179  TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
180 #endif
181 
182  // Try to cast mv to an array of this type
183  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
184  Teuchos::Array<RCP<const TMV> > tmvs(mv.size());
185  RCP<const TMV> tmv;
186  bool allCastsSuccessful = true;
187  {
188  auto mvIter = mv.begin();
189  auto tmvIter = tmvs.begin();
190  for (; mvIter != mv.end(); ++mvIter, ++tmvIter) {
191  tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromPtr(*mvIter));
192  if (nonnull(tmv)) {
193  *tmvIter = tmv;
194  } else {
195  allCastsSuccessful = false;
196  break;
197  }
198  }
199  }
200 
201  // If casts succeeded, or input arrays are size 0, call Tpetra directly.
202  // Otherwise, fall back to the RTOp implementation.
203  auto len = tmvs.size();
204  if (len == 0) {
205  tpetraMultiVector_.getNonconstObj()->scale(beta);
206  } else if (len == 1 && allCastsSuccessful) {
207  tpetraMultiVector_.getNonconstObj()->update(alpha[0], *tmvs[0], beta);
208  } else if (len == 2 && allCastsSuccessful) {
209  tpetraMultiVector_.getNonconstObj()->update(alpha[0], *tmvs[0], alpha[1], *tmvs[1], beta);
210  } else if (allCastsSuccessful) {
212  auto tmvIter = tmvs.begin();
213  auto alphaIter = alpha.begin();
214 
215  // Check if any entry of tmvs aliases this object's wrapped vector.
216  // If so, replace that entry in the array with a copy.
217  tmv = Teuchos::null;
218  for (; tmvIter != tmvs.end(); ++tmvIter) {
219  if (tmvIter->getRawPtr() == tpetraMultiVector_.getConstObj().getRawPtr()) {
220  if (tmv.is_null()) {
221  tmv = Teuchos::rcp(new TMV(*tpetraMultiVector_.getConstObj(), Teuchos::Copy));
222  }
223  *tmvIter = tmv;
224  }
225  }
226  tmvIter = tmvs.begin();
227 
228  // We add two MVs at a time, so only scale if even num MVs,
229  // and additionally do the first addition if odd num MVs.
230  if ((tmvs.size() % 2) == 0) {
231  tpetraMultiVector_.getNonconstObj()->scale(beta);
232  } else {
233  tpetraMultiVector_.getNonconstObj()->update(*alphaIter, *(*tmvIter), beta);
234  ++tmvIter;
235  ++alphaIter;
236  }
237  for (; tmvIter != tmvs.end(); tmvIter+=2, alphaIter+=2) {
238  tpetraMultiVector_.getNonconstObj()->update(
239  *alphaIter, *(*tmvIter), *(alphaIter+1), *(*(tmvIter+1)), ST::one());
240  }
241  } else {
242  // This version will require/modify the host view of this vector.
243  tpetraMultiVector_.getNonconstObj()->sync_host ();
244  tpetraMultiVector_.getNonconstObj()->modify_host ();
246  }
247 }
248 
249 
250 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
252  const MultiVectorBase<Scalar>& mv,
253  const ArrayView<Scalar>& prods
254  ) const
255 {
256  auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
257 
258  // If the cast succeeded, call Tpetra directly.
259  // Otherwise, fall back to the RTOp implementation.
260  if (nonnull(tmv)) {
261  tpetraMultiVector_.getConstObj()->dot(*tmv, prods);
262  } else {
263  // This version will require/modify the host view of this vector.
264  tpetraMultiVector_.getNonconstObj()->sync_host ();
265  tpetraMultiVector_.getNonconstObj()->modify_host ();
267  }
268 }
269 
270 
271 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
273  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
274  ) const
275 {
276  tpetraMultiVector_.getConstObj()->norm1(norms);
277 }
278 
279 
280 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
282  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
283  ) const
284 {
285  tpetraMultiVector_.getConstObj()->norm2(norms);
286 }
287 
288 
289 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
291  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
292  ) const
293 {
294  tpetraMultiVector_.getConstObj()->normInf(norms);
295 }
296 
297 
298 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301 {
302 #ifdef TEUCHOS_DEBUG
303  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j, 0, this->domain()->dim());
304 #endif
305  return constTpetraVector<Scalar>(
306  tpetraVectorSpace_,
307  tpetraMultiVector_->getVector(j)
308  );
309 }
310 
311 
312 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
315 {
316 #ifdef TEUCHOS_DEBUG
317  TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j, 0, this->domain()->dim());
318 #endif
319  return tpetraVector<Scalar>(
320  tpetraVectorSpace_,
321  tpetraMultiVector_.getNonconstObj()->getVectorNonConst(j)
322  );
323 }
324 
325 
326 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
329  const Range1D& col_rng_in
330  ) const
331 {
332 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
333  std::cerr << "\nTpetraMultiVector::subView(Range1D) const called!\n";
334 #endif
335  const Range1D colRng = this->validateColRange(col_rng_in);
336 
338  this->getConstTpetraMultiVector()->subView(colRng);
339 
340  const RCP<const ScalarProdVectorSpaceBase<Scalar> > viewDomainSpace =
341  tpetraVectorSpace<Scalar>(
342  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
343  tpetraView->getNumVectors(),
344  tpetraView->getMap()->getComm()
345  )
346  );
347 
348  return constTpetraMultiVector(
349  tpetraVectorSpace_,
350  viewDomainSpace,
351  tpetraView
352  );
353 }
354 
355 
356 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
359  const Range1D& col_rng_in
360  )
361 {
362 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
363  std::cerr << "\nTpetraMultiVector::subView(Range1D) called!\n";
364 #endif
365  const Range1D colRng = this->validateColRange(col_rng_in);
366 
368  this->getTpetraMultiVector()->subViewNonConst(colRng);
369 
370  const RCP<const ScalarProdVectorSpaceBase<Scalar> > viewDomainSpace =
371  tpetraVectorSpace<Scalar>(
372  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
373  tpetraView->getNumVectors(),
374  tpetraView->getMap()->getComm()
375  )
376  );
377 
378  return tpetraMultiVector(
379  tpetraVectorSpace_,
380  viewDomainSpace,
381  tpetraView
382  );
383 }
384 
385 
386 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
389  const ArrayView<const int>& cols_in
390  ) const
391 {
392 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
393  std::cerr << "\nTpetraMultiVector::subView(ArrayView) const called!\n";
394 #endif
395  // Tpetra wants col indices as size_t
396  Array<std::size_t> cols(cols_in.size());
397  for (Array<std::size_t>::size_type i = 0; i < cols.size(); ++i)
398  cols[i] = static_cast<std::size_t>(cols_in[i]);
399 
401  this->getConstTpetraMultiVector()->subView(cols());
402 
403  const RCP<const ScalarProdVectorSpaceBase<Scalar> > viewDomainSpace =
404  tpetraVectorSpace<Scalar>(
405  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
406  tpetraView->getNumVectors(),
407  tpetraView->getMap()->getComm()
408  )
409  );
410 
411  return constTpetraMultiVector(
412  tpetraVectorSpace_,
413  viewDomainSpace,
414  tpetraView
415  );
416 }
417 
418 
419 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
422  const ArrayView<const int>& cols_in
423  )
424 {
425 #ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
426  std::cerr << "\nTpetraMultiVector::subView(ArrayView) called!\n";
427 #endif
428  // Tpetra wants col indices as size_t
429  Array<std::size_t> cols(cols_in.size());
430  for (Array<std::size_t>::size_type i = 0; i < cols.size(); ++i)
431  cols[i] = static_cast<std::size_t>(cols_in[i]);
432 
434  this->getTpetraMultiVector()->subViewNonConst(cols());
435 
436  const RCP<const ScalarProdVectorSpaceBase<Scalar> > viewDomainSpace =
437  tpetraVectorSpace<Scalar>(
438  Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
439  tpetraView->getNumVectors(),
440  tpetraView->getMap()->getComm()
441  )
442  );
443 
444  return tpetraMultiVector(
445  tpetraVectorSpace_,
446  viewDomainSpace,
447  tpetraView
448  );
449 }
450 
451 
452 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
455  const RTOpPack::RTOpT<Scalar> &primary_op,
456  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs,
457  const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs,
458  const ArrayView<const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
459  const Ordinal primary_global_offset
460  ) const
461 {
463 
464  // Sync any non-target Tpetra MVs to host space
465  for (auto itr = multi_vecs.begin(); itr != multi_vecs.end(); ++itr) {
466  Ptr<const TMV> tmv = Teuchos::ptr_dynamic_cast<const TMV>(*itr);
467  if (nonnull(tmv)) {
468  Teuchos::rcp_const_cast<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(
469  tmv->getConstTpetraMultiVector())-> sync_host ();
470  }
471  }
472 
473  // Sync any target Tpetra MVs and mark modified
474  for (auto itr = targ_multi_vecs.begin(); itr != targ_multi_vecs.end(); ++itr) {
475  Ptr<TMV> tmv = Teuchos::ptr_dynamic_cast<TMV>(*itr);
476  if (nonnull(tmv)) {
477  tmv->getTpetraMultiVector()->sync_host ();
478  tmv->getTpetraMultiVector()->modify_host ();
479  }
480  }
481 
483  primary_op, multi_vecs, targ_multi_vecs, reduct_objs, primary_global_offset);
484 }
485 
486 
487 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
490  const Range1D &rowRng,
491  const Range1D &colRng,
493  ) const
494 {
495  // Only viewing data, so just sync dual view to host space
496  typedef typename Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
497  Teuchos::rcp_const_cast<TMV>(
498  tpetraMultiVector_.getConstObj())->sync_host ();
499 
501  acquireDetachedMultiVectorViewImpl(rowRng, colRng, sub_mv);
502 }
503 
504 
505 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
508  const Range1D &rowRng,
509  const Range1D &colRng,
511  )
512 {
513  // Sync to host and mark as modified
514  tpetraMultiVector_.getNonconstObj()->sync_host ();
515  tpetraMultiVector_.getNonconstObj()->modify_host ();
516 
518  acquireNonconstDetachedMultiVectorViewImpl(rowRng, colRng, sub_mv);
519 }
520 
521 
522 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
526  )
527 {
530 
531  // Sync changes from host view to execution space
532  typedef typename Tpetra::MultiVector<
533  Scalar,LocalOrdinal,GlobalOrdinal,Node>::execution_space execution_space;
534  tpetraMultiVector_.getNonconstObj()->template sync<execution_space>();
535 }
536 
537 
538 /* ToDo: Implement these?
539 
540 
541 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
542 RCP<const MultiVectorBase<Scalar> >
543 TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::nonContigSubViewImpl(
544  const ArrayView<const int> &cols
545  ) const
546 {
547  THYRA_DEBUG_ASSERT_MV_COLS("nonContigSubViewImpl(cols)", cols);
548  const int numCols = cols.size();
549  const ArrayRCP<Scalar> localValuesView = createContiguousCopy(cols);
550  return defaultSpmdMultiVector<Scalar>(
551  spmdRangeSpace_,
552  createSmallScalarProdVectorSpaceBase<Scalar>(spmdRangeSpace_, numCols),
553  localValuesView
554  );
555 }
556 
557 
558 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
559 RCP<MultiVectorBase<Scalar> >
560 TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::nonconstNonContigSubViewImpl(
561  const ArrayView<const int> &cols )
562 {
563  THYRA_DEBUG_ASSERT_MV_COLS("nonContigSubViewImpl(cols)", cols);
564  const int numCols = cols.size();
565  const ArrayRCP<Scalar> localValuesView = createContiguousCopy(cols);
566  const Ordinal localSubDim = spmdRangeSpace_->localSubDim();
567  RCP<CopyBackSpmdMultiVectorEntries<Scalar> > copyBackView =
568  copyBackSpmdMultiVectorEntries<Scalar>(cols, localValuesView.getConst(),
569  localSubDim, localValues_.create_weak(), leadingDim_);
570  return Teuchos::rcpWithEmbeddedObjPreDestroy(
571  new TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(
572  spmdRangeSpace_,
573  createSmallScalarProdVectorSpaceBase<Scalar>(spmdRangeSpace_, numCols),
574  localValuesView),
575  copyBackView
576  );
577 }
578 
579 */
580 
581 
582 // Overridden protected members from SpmdMultiVectorBase
583 
584 
585 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
588 {
589  return tpetraVectorSpace_;
590 }
591 
592 
593 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
595  const Ptr<ArrayRCP<Scalar> > &localValues, const Ptr<Ordinal> &leadingDim
596  )
597 {
598  *localValues = tpetraMultiVector_.getNonconstObj()->get1dViewNonConst();
599  *leadingDim = tpetraMultiVector_->getStride();
600 }
601 
602 
603 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
605  const Ptr<ArrayRCP<const Scalar> > &localValues, const Ptr<Ordinal> &leadingDim
606  ) const
607 {
608  *localValues = tpetraMultiVector_->get1dView();
609  *leadingDim = tpetraMultiVector_->getStride();
610 }
611 
612 
613 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
615  const EOpTransp M_trans,
616  const MultiVectorBase<Scalar> &X,
617  const Ptr<MultiVectorBase<Scalar> > &Y,
618  const Scalar alpha,
619  const Scalar beta
620  ) const
621 {
622  // Try to extract Tpetra objects from X and Y
623  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
624  Teuchos::RCP<const TMV> X_tpetra = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(X));
625  Teuchos::RCP<TMV> Y_tpetra = this->getTpetraMultiVector(Teuchos::rcpFromPtr(Y));
626 
627  // If the cast succeeded, call Tpetra directly.
628  // Otherwise, fall back to the default implementation.
629  if (nonnull(X_tpetra) && nonnull(Y_tpetra)) {
630  // Sync everything to the execution space
631  typedef typename TMV::execution_space execution_space;
632  Teuchos::rcp_const_cast<TMV>(X_tpetra)->template sync<execution_space>();
633  Y_tpetra->template sync<execution_space>();
634  Teuchos::rcp_const_cast<TMV>(
635  tpetraMultiVector_.getConstObj())->template sync<execution_space>();
636 
638  TEUCHOS_TEST_FOR_EXCEPTION(ST::isComplex && (M_trans == CONJ),
639  std::logic_error,
640  "Error, conjugation without transposition is not allowed for complex scalar types!");
641 
643  switch (M_trans) {
644  case NOTRANS:
645  trans = Teuchos::NO_TRANS;
646  break;
647  case CONJ:
648  trans = Teuchos::NO_TRANS;
649  break;
650  case TRANS:
651  trans = Teuchos::TRANS;
652  break;
653  case CONJTRANS:
654  trans = Teuchos::CONJ_TRANS;
655  break;
656  }
657 
658  Y_tpetra->template modify<execution_space>();
659  Y_tpetra->multiply(trans, Teuchos::NO_TRANS, alpha, *tpetraMultiVector_.getConstObj(), *X_tpetra, beta);
660  } else {
661  Teuchos::rcp_const_cast<TMV>(
662  tpetraMultiVector_.getConstObj())->sync_host ();
663  SpmdMultiVectorDefaultBase<Scalar>::euclideanApply(M_trans, X, Y, alpha, beta);
664  }
665 
666 }
667 
668 // private
669 
670 
671 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
672 template<class TpetraMultiVector_t>
674  const RCP<const TpetraVectorSpace<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraVectorSpace,
675  const RCP<const ScalarProdVectorSpaceBase<Scalar> > &domainSpace,
676  const RCP<TpetraMultiVector_t> &tpetraMultiVector
677  )
678 {
679 #ifdef THYRA_DEBUG
680  TEUCHOS_ASSERT(nonnull(tpetraVectorSpace));
681  TEUCHOS_ASSERT(nonnull(domainSpace));
682  TEUCHOS_ASSERT(nonnull(tpetraMultiVector));
683  // ToDo: Check to make sure that tpetraMultiVector is compatible with
684  // tpetraVectorSpace.
685 #endif
686  tpetraVectorSpace_ = tpetraVectorSpace;
687  domainSpace_ = domainSpace;
688  tpetraMultiVector_.initialize(tpetraMultiVector);
689  this->updateSpmdSpace();
690 }
691 
692 
693 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
694 RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
696 getTpetraMultiVector(const RCP<MultiVectorBase<Scalar> >& mv) const
697 {
698  using Teuchos::rcp_dynamic_cast;
701 
702  RCP<TMV> tmv = rcp_dynamic_cast<TMV>(mv);
703  if (nonnull(tmv)) {
704  return tmv->getTpetraMultiVector();
705  }
706 
707  RCP<TV> tv = rcp_dynamic_cast<TV>(mv);
708  if (nonnull(tv)) {
709  return tv->getTpetraVector();
710  }
711 
712  return Teuchos::null;
713 }
714 
715 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
716 RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
718 getConstTpetraMultiVector(const RCP<const MultiVectorBase<Scalar> >& mv) const
719 {
720  using Teuchos::rcp_dynamic_cast;
723 
724  RCP<const TMV> tmv = rcp_dynamic_cast<const TMV>(mv);
725  if (nonnull(tmv)) {
726  return tmv->getConstTpetraMultiVector();
727  }
728 
729  RCP<const TV> tv = rcp_dynamic_cast<const TV>(mv);
730  if (nonnull(tv)) {
731  return tv->getConstTpetraVector();
732  }
733 
734  return Teuchos::null;
735 }
736 
737 
738 } // end namespace Thyra
739 
740 
741 #endif // THYRA_TPETRA_MULTIVECTOR_HPP
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Concrete implementation of Thyra::MultiVector in terms of Tpetra::MultiVector.
TpetraMultiVector()
Construct to uninitialized.
RCP< const VectorBase< Scalar > > colImpl(Ordinal j) const
void getNonconstLocalMultiVectorDataImpl(const Ptr< ArrayRCP< Scalar > > &localValues, const Ptr< Ordinal > &leadingDim)
Concrete implementation of an SPMD vector space for Tpetra.
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraMultiVector()
Extract the underlying non-const Tpetra::MultiVector object.
virtual void assignImpl(Scalar alpha)
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
void euclideanApply(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
Uses GEMM() and Teuchos::reduceAll() to implement.
iterator begin() const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
Default implementation of assign(MV) using RTOps.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
virtual void euclideanApply(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
size_type size() const
virtual void mvMultiReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const ArrayView< const Ptr< RTOpPack::ReductTarget > > &reduct_objs, const Ordinal primary_global_offset) const
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Default implementation of linear_combination using RTOps.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols_in) const
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
void initialize(const RCP< const TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVectorSpace, const RCP< const ScalarProdVectorSpaceBase< Scalar > > &domainSpace, const RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraMultiVector)
Initialize.
Use the transposed operator.
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
virtual void scaleImpl(Scalar alpha)
Ordinal size_type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
virtual void mvMultiReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const ArrayView< const Ptr< RTOpPack::ReductTarget > > &reduct_objs, const Ordinal primary_global_offset) const
Interface for a collection of column vectors called a multi-vector.
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols_in)
void constInitialize(const RCP< const TpetraVectorSpace< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraVectorSpace, const RCP< const ScalarProdVectorSpaceBase< Scalar > > &domainSpace, const RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraMultiVector)
Initialize.
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Concrete Thyra::SpmdVectorBase using Tpetra::Vector.
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
bool nonnull(const boost::shared_ptr< T > &p)
void getLocalMultiVectorDataImpl(const Ptr< ArrayRCP< const Scalar > > &localValues, const Ptr< Ordinal > &leadingDim) const
size_type size() const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraMultiVector() const
Extract the underlying const Tpetra::MultiVector object.
RCP< const SpmdVectorSpaceBase< Scalar > > spmdSpaceImpl() const
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
iterator begin()
RCP< const ScalarProdVectorSpaceBase< Scalar > > domainScalarProdVecSpc() const
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Default implementation of dots using RTOps.
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Default implementation of update using RTOps.