42#ifndef THYRA_TPETRA_MULTIVECTOR_HPP
43#define THYRA_TPETRA_MULTIVECTOR_HPP
45#include "Thyra_TpetraMultiVector_decl.hpp"
46#include "Thyra_TpetraVectorSpace.hpp"
47#include "Thyra_TpetraVector.hpp"
48#include "Teuchos_Assert.hpp"
57template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
62template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 const RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraMultiVector
69 initializeImpl(tpetraVectorSpace, domainSpace, tpetraMultiVector);
73template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 const RCP<
const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraMultiVector
80 initializeImpl(tpetraVectorSpace, domainSpace, tpetraMultiVector);
84template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 return tpetraMultiVector_.getNonconstObj();
92template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 return tpetraMultiVector_;
103template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
114template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 tpetraMultiVector_.getNonconstObj()->putScalar(alpha);
122template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
131 tpetraMultiVector_.getNonconstObj()->assign(*tmv);
138template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
142 tpetraMultiVector_.getNonconstObj()->scale(alpha);
146template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
152 auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
158 tpetraMultiVector_.getNonconstObj()->update(alpha, *tmv, ST::one());
165template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
180 bool allCastsSuccessful =
true;
182 auto mvIter = mv.begin();
183 auto tmvIter = tmvs.
begin();
184 for (; mvIter != mv.end(); ++mvIter, ++tmvIter) {
185 tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromPtr(*mvIter));
189 allCastsSuccessful =
false;
197 auto len = tmvs.
size();
199 tpetraMultiVector_.getNonconstObj()->scale(beta);
200 }
else if (len == 1 && allCastsSuccessful) {
201 tpetraMultiVector_.getNonconstObj()->update(alpha[0], *tmvs[0], beta);
202 }
else if (len == 2 && allCastsSuccessful) {
203 tpetraMultiVector_.getNonconstObj()->update(alpha[0], *tmvs[0], alpha[1], *tmvs[1], beta);
204 }
else if (allCastsSuccessful) {
206 auto tmvIter = tmvs.
begin();
207 auto alphaIter = alpha.
begin();
212 for (; tmvIter != tmvs.
end(); ++tmvIter) {
213 if (tmvIter->getRawPtr() == tpetraMultiVector_.getConstObj().getRawPtr()) {
220 tmvIter = tmvs.
begin();
224 if ((tmvs.
size() % 2) == 0) {
225 tpetraMultiVector_.getNonconstObj()->scale(beta);
227 tpetraMultiVector_.getNonconstObj()->update(*alphaIter, *(*tmvIter), beta);
231 for (; tmvIter != tmvs.
end(); tmvIter+=2, alphaIter+=2) {
232 tpetraMultiVector_.getNonconstObj()->update(
233 *alphaIter, *(*tmvIter), *(alphaIter+1), *(*(tmvIter+1)), ST::one());
241template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
247 auto tmv = this->getConstTpetraMultiVector(Teuchos::rcpFromRef(mv));
252 tpetraMultiVector_.getConstObj()->dot(*tmv, prods);
259template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
264 tpetraMultiVector_.getConstObj()->norm1(norms);
268template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
273 tpetraMultiVector_.getConstObj()->norm2(norms);
277template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
282 tpetraMultiVector_.getConstObj()->normInf(norms);
286template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
293 return constTpetraVector<Scalar>(
295 tpetraMultiVector_->getVector(j)
300template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
307 return tpetraVector<Scalar>(
309 tpetraMultiVector_.getNonconstObj()->getVectorNonConst(j)
314template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
320#ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
321 std::cerr <<
"\nTpetraMultiVector::subView(Range1D) const called!\n";
323 const Range1D colRng = this->validateColRange(col_rng_in);
326 this->getConstTpetraMultiVector()->subView(colRng);
329 tpetraVectorSpace<Scalar>(
330 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
331 tpetraView->getNumVectors(),
332 tpetraView->getMap()->getComm()
336 return constTpetraMultiVector(
344template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
350#ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
351 std::cerr <<
"\nTpetraMultiVector::subView(Range1D) called!\n";
353 const Range1D colRng = this->validateColRange(col_rng_in);
356 this->getTpetraMultiVector()->subViewNonConst(colRng);
359 tpetraVectorSpace<Scalar>(
360 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
361 tpetraView->getNumVectors(),
362 tpetraView->getMap()->getComm()
366 return tpetraMultiVector(
374template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
380#ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
381 std::cerr <<
"\nTpetraMultiVector::subView(ArrayView) const called!\n";
386 cols[i] =
static_cast<std::size_t
>(cols_in[i]);
389 this->getConstTpetraMultiVector()->subView(cols());
392 tpetraVectorSpace<Scalar>(
393 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
394 tpetraView->getNumVectors(),
395 tpetraView->getMap()->getComm()
399 return constTpetraMultiVector(
407template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
413#ifdef THYRA_DEFAULT_SPMD_MULTI_VECTOR_VERBOSE_TO_ERROR_OUT
414 std::cerr <<
"\nTpetraMultiVector::subView(ArrayView) called!\n";
419 cols[i] =
static_cast<std::size_t
>(cols_in[i]);
422 this->getTpetraMultiVector()->subViewNonConst(cols());
425 tpetraVectorSpace<Scalar>(
426 Tpetra::createLocalMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
427 tpetraView->getNumVectors(),
428 tpetraView->getMap()->getComm()
432 return tpetraMultiVector(
440template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
447 const Ordinal primary_global_offset
452 primary_op, multi_vecs, targ_multi_vecs, reduct_objs, primary_global_offset);
456template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
469template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
482template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
541template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
545 return tpetraVectorSpace_;
549template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
554 *localValues = tpetraMultiVector_.getNonconstObj()->get1dViewNonConst();
555 *leadingDim = tpetraMultiVector_->getStride();
559template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
564 *localValues = tpetraMultiVector_->get1dView();
565 *leadingDim = tpetraMultiVector_->getStride();
569template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
579 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> TMV;
585 if (nonnull(X_tpetra) && nonnull(Y_tpetra)) {
589 "Error, conjugation without transposition is not allowed for complex scalar types!");
607 Y_tpetra->multiply(trans,
Teuchos::NO_TRANS, alpha, *tpetraMultiVector_.getConstObj(), *X_tpetra, beta);
618template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
619template<
class TpetraMultiVector_t>
633 tpetraVectorSpace_ = tpetraVectorSpace;
634 domainSpace_ = domainSpace;
635 tpetraMultiVector_.initialize(tpetraMultiVector);
636 this->updateSpmdSpace();
640template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
641RCP<Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
645 using Teuchos::rcp_dynamic_cast;
649 RCP<TMV> tmv = rcp_dynamic_cast<TMV>(mv);
651 return tmv->getTpetraMultiVector();
654 RCP<TV> tv = rcp_dynamic_cast<TV>(mv);
656 return tv->getTpetraVector();
659 return Teuchos::null;
662template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
663RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
667 using Teuchos::rcp_dynamic_cast;
671 RCP<const TMV> tmv = rcp_dynamic_cast<const TMV>(mv);
673 return tmv->getConstTpetraMultiVector();
676 RCP<const TV> tv = rcp_dynamic_cast<const TV>(mv);
678 return tv->getConstTpetraVector();
681 return Teuchos::null;
Interface for a collection of column vectors called a multi-vector.
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 =0
Apply a reduction/transformation operator column by column and return an array of the reduction objec...
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Default implementation of dots using RTOps.
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.
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
Default implementation of assign(MV) using RTOps.
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Default implementation of update using RTOps.
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
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.
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
Concrete implementation of Thyra::MultiVector in terms of Tpetra::MultiVector.
RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraMultiVector()
Extract the underlying non-const Tpetra::MultiVector object.
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.
RCP< const SpmdVectorSpaceBase< Scalar > > spmdSpaceImpl() const
RCP< const ScalarProdVectorSpaceBase< Scalar > > domainScalarProdVecSpc() const
RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraMultiVector() const
Extract the underlying const Tpetra::MultiVector object.
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
RCP< const VectorBase< Scalar > > colImpl(Ordinal j) const
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
void getLocalMultiVectorDataImpl(const Ptr< ArrayRCP< const Scalar > > &localValues, const Ptr< Ordinal > &leadingDim) const
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols_in)
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols_in) const
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
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 norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
TpetraMultiVector()
Construct to uninitialized.
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
virtual void euclideanApply(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
virtual void scaleImpl(Scalar alpha)
void getNonconstLocalMultiVectorDataImpl(const Ptr< ArrayRCP< Scalar > > &localValues, const Ptr< Ordinal > &leadingDim)
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 assignImpl(Scalar alpha)
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
Concrete implementation of an SPMD vector space for Tpetra.
Concrete Thyra::SpmdVectorBase using Tpetra::Vector.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
bool nonnull(const std::shared_ptr< T > &p)
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
@ TRANS
Use the transposed operator.
@ NOTRANS
Use the non-transposed operator.
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
@ CONJ
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)