42#ifndef THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
43#define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
45#include "Thyra_DefaultProductMultiVector_decl.hpp"
46#include "Thyra_DefaultProductVectorSpace.hpp"
47#include "Thyra_DefaultProductVector.hpp"
48#include "Thyra_AssertOp.hpp"
74 const int numBlocks = productSpace_in->numBlocks();
75 for (
int k = 0; k < numBlocks; ++k )
76 multiVecs.
push_back(createMembers(productSpace_in->getBlock(k),numMembers));
77 initialize(productSpace_in,multiVecs);
87 initializeImpl(productSpace_in,multiVecs);
97 initializeImpl(productSpace_in,multiVecs);
101template<
class Scalar>
104 productSpace_ = Teuchos::null;
105 multiVecs_.resize(0);
113template<
class Scalar>
116 std::ostringstream oss;
120 <<
"rangeDim="<<this->range()->dim()
121 <<
",domainDim="<<this->domain()->dim()
122 <<
",numBlocks = "<<numBlocks_
128template<
class Scalar>
130 FancyOStream &out_arg,
135 using Teuchos::describe;
145 *out << this->description() << std::endl;
153 <<
"rangeDim="<<this->range()->dim()
154 <<
",domainDim="<<this->domain()->dim()
158 <<
"numBlocks="<< numBlocks_ << std::endl
159 <<
"Constituent multi-vector objects V[0], V[1], ... V[numBlocks-1]:\n";
161 for(
int k = 0; k < numBlocks_; ++k ) {
162 *out <<
"V["<<k<<
"] = "
163 << Teuchos::describe(*multiVecs_[k].getConstObj(),verbLevel);
176template<
class Scalar>
180 return productSpace_;
184template<
class Scalar>
187 return multiVecs_[k].isConst();
191template<
class Scalar>
195 return multiVecs_[k].getNonconstObj();
199template<
class Scalar>
203 return multiVecs_[k].getConstObj();
210template<
class Scalar>
216 for (
int k = 0; k < numBlocks_; ++k )
217 blocks.
push_back(multiVecs_[k].getConstObj()->clone_mv());
218 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
225template<
class Scalar>
229 return productSpace_;
233template<
class Scalar>
237 if (is_null(productSpace_))
238 return Teuchos::null;
239 return multiVecs_[0].getConstObj()->domain();
249template<
class Scalar>
252 for (
int k = 0; k < numBlocks_; ++k ) {
253 multiVecs_[k].getNonconstObj()->assign(alpha);
258template<
class Scalar>
270 = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
272 for (
int k = 0; k < numBlocks_; ++k) {
273 multiVecs_[k].getNonconstObj()->assign(*pv->getMultiVectorBlock(k));
281template<
class Scalar>
284 for (
int k = 0; k < numBlocks_; ++k) {
285 multiVecs_[k].getNonconstObj()->scale(alpha);
290template<
class Scalar>
303 = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
305 for (
int k = 0; k < numBlocks_; ++k) {
306 multiVecs_[k].getNonconstObj()->update(alpha, *pv->getMultiVectorBlock(k));
314template<
class Scalar>
323 for (
Ordinal i = 0; i < mv.size(); ++i) {
330 bool allCastsSuccessful =
true;
332 for (
Ordinal i = 0; i < mv.size(); ++i) {
333 pvs[i] = Teuchos::ptr_dynamic_cast<const ProductMultiVectorBase<Scalar> >(mv[i]);
334 if (pvs[i].is_null()) {
335 allCastsSuccessful =
false;
340 if (allCastsSuccessful) {
343 for (
int k = 0; k < numBlocks_; ++k) {
345 blocks_rcp[i] = pvs[i]->getMultiVectorBlock(k);
346 blocks[i] = blocks_rcp[i].ptr();
348 multiVecs_[k].getNonconstObj()->linear_combination(alpha, blocks(), beta);
356template<
class Scalar>
369 = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
375 for (
int k = 0; k < numBlocks_; ++k) {
376 multiVecs_[k].getConstObj()->dots(*pv->getMultiVectorBlock(k), subProds());
378 prods[i] += subProds[i];
387template<
class Scalar>
396 for (
Ordinal i = 0; i < norms.size(); ++i)
397 norms[i] = ST::magnitude(ST::zero());
400 for (
int k = 0; k < numBlocks_; ++k) {
401 multiVecs_[k].getConstObj()->norms_1(subNorms());
402 for (
Ordinal i = 0; i < norms.size(); ++i) {
403 norms[i] += subNorms[i];
409template<
class Scalar>
418 for (
Ordinal i = 0; i < norms.size(); ++i)
419 norms[i] = ST::magnitude(ST::zero());
422 for (
int k = 0; k < numBlocks_; ++k) {
423 multiVecs_[k].getConstObj()->norms_2(subNorms());
424 for (
Ordinal i = 0; i < norms.size(); ++i) {
425 norms[i] += subNorms[i]*subNorms[i];
429 for (
Ordinal i = 0; i < norms.size(); ++i) {
430 norms[i] = ST::magnitude(ST::squareroot(norms[i]));
435template<
class Scalar>
444 for (
Ordinal i = 0; i < norms.size(); ++i)
445 norms[i] = ST::magnitude(ST::zero());
448 for (
int k = 0; k < numBlocks_; ++k) {
449 multiVecs_[k].getConstObj()->norms_inf(subNorms());
450 for (
Ordinal i = 0; i < norms.size(); ++i) {
451 if (subNorms[i] > norms[i])
452 norms[i] = subNorms[i];
458template<
class Scalar>
464 for (
int k = 0; k < numBlocks_; ++k )
465 cols_.
push_back(multiVecs_[k].getConstObj()->col(j));
466 return defaultProductVector<Scalar>(productSpace_, cols_());
470template<
class Scalar>
476 for (
int k = 0; k < numBlocks_; ++k )
477 cols_.
push_back(multiVecs_[k].getNonconstObj()->col(j));
478 return defaultProductVector<Scalar>(productSpace_, cols_());
482template<
class Scalar>
488 for (
int k = 0; k < numBlocks_; ++k )
489 blocks.
push_back(multiVecs_[k].getConstObj()->subView(colRng));
490 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
494template<
class Scalar>
500 for (
int k = 0; k < numBlocks_; ++k )
501 blocks.
push_back(multiVecs_[k].getNonconstObj()->subView(colRng));
502 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
506template<
class Scalar>
514 for (
int k = 0; k < numBlocks_; ++k )
515 blocks.
push_back(multiVecs_[k].getConstObj()->subView(cols));
516 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
520template<
class Scalar>
528 for (
int k = 0; k < numBlocks_; ++k )
529 blocks.
push_back(multiVecs_[k].getNonconstObj()->subView(cols));
530 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
534template<
class Scalar>
540 const Ordinal primary_global_offset_in
544 using Teuchos::ptr_dynamic_cast;
545 using Thyra::applyOp;
550 for (
int j = 0; j < multi_vecs_in.size(); ++j ) {
552 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
553 *this->range(), *multi_vecs_in[j]->range()
556 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
557 *this->domain(), *multi_vecs_in[j]->domain()
560 for (
int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
562 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
563 *this->range(), *targ_multi_vecs_inout[j]->range()
566 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
567 *this->domain(), *targ_multi_vecs_inout[j]->domain()
577 bool allProductMultiVectors =
true;
580 for (
int j = 0; j < multi_vecs_in.size() && allProductMultiVectors; ++j ) {
585 multi_vecs_j = ptr_dynamic_cast<const ProductMultiVectorBase<Scalar> >(
588 if ( !is_null(multi_vecs_j) ) {
592 allProductMultiVectors =
false;
597 for (
int j = 0; j < targ_multi_vecs_inout.size() && allProductMultiVectors; ++j )
603 targ_multi_vecs_j = ptr_dynamic_cast<ProductMultiVectorBase<Scalar> >(
604 targ_multi_vecs_inout[j]
606 if (!is_null(targ_multi_vecs_j)) {
607 targ_multi_vecs.
push_back(targ_multi_vecs_j);
610 allProductMultiVectors =
false;
618 if ( allProductMultiVectors ) {
628 multi_vecs_rcp_block_k(multi_vecs_in.size());
630 multi_vecs_block_k(multi_vecs_in.size());
632 targ_multi_vecs_rcp_block_k(targ_multi_vecs_inout.size());
634 targ_multi_vecs_block_k(targ_multi_vecs_inout.size());
636 Ordinal g_off = primary_global_offset_in;
638 for (
int k = 0; k < numBlocks_; ++k ) {
640 const Ordinal dim_k = productSpace_->getBlock(k)->dim();
644 for (
int j = 0; j < multi_vecs_in.size(); ++j ) {
646 multi_vecs[j]->getMultiVectorBlock(k);
647 multi_vecs_rcp_block_k[j] = multi_vecs_rcp_block_k_j;
648 multi_vecs_block_k[j] = multi_vecs_rcp_block_k_j.
ptr();
651 for (
int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
653 targ_multi_vecs[j]->getNonconstMultiVectorBlock(k);
654 targ_multi_vecs_rcp_block_k[j] = targ_multi_vecs_rcp_block_k_j;
655 targ_multi_vecs_block_k[j] = targ_multi_vecs_rcp_block_k_j.
ptr();
660 Thyra::applyOp<Scalar>(
661 primary_op, multi_vecs_block_k(), targ_multi_vecs_block_k(),
678 primary_op, multi_vecs_in(), targ_multi_vecs_inout(),
679 reduct_objs, primary_global_offset_in);
686template<
class Scalar>
694 rowRng, colRng, sub_mv );
699template<
class Scalar>
710template<
class Scalar>
718 rowRng,colRng,sub_mv );
723template<
class Scalar>
736template<
class Scalar>
743template<
class Scalar>
759 "DefaultProductMultiVector<Scalar>::apply(...)",
760 *
this, M_trans, X_in, &*Y_inout );
772 &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
773 for (
int k = 0; k < numBlocks_; ++k ) {
775 *multiVecs_[k].getConstObj(), M_trans,
789 &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
790 for (
int k = 0; k < numBlocks_; ++k ) {
792 M_k = multiVecs_[k].getConstObj(),
795 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.
ptr(), alpha, beta );
798 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.
ptr(), alpha, ST::one() );
808template<
class Scalar>
809template<
class MultiVectorType>
823 theDomain = multiVecs[0]->domain();
824 const int numBlocks = productSpace_in->numBlocks();
826 for (
int k = 0; k < numBlocks; ++k ) {
829 *theDomain, *multiVecs[k]->domain()
833 productSpace_ = productSpace_in;
834 numBlocks_ = numBlocks;
835 multiVecs_.assign(multiVecs.begin(),multiVecs.end());
842template<
class Scalar>
843void DefaultProductMultiVector<Scalar>::assertInitialized()
const
846 is_null(productSpace_), std::logic_error,
847 "Error, this DefaultProductMultiVector object is not intialized!"
852template<
class Scalar>
853void DefaultProductMultiVector<Scalar>::validateColIndex(
const int j)
const
856 const int domainDim = multiVecs_[0].getConstObj()->domain()->dim();
858 ! ( 0 <= j && j < domainDim ), std::logic_error,
859 "Error, the column index j = " << j <<
" does not fall in the range [0,"<<domainDim<<
"]!"
870template<
class Scalar>
872Thyra::defaultProductMultiVector()
874 return Teuchos::rcp(
new DefaultProductMultiVector<Scalar>);
878template<
class Scalar>
880Thyra::defaultProductMultiVector(
881 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
885 RCP<DefaultProductMultiVector<Scalar> > pmv = defaultProductMultiVector<Scalar>();
886 pmv->initialize(productSpace, numMembers);
891template<
class Scalar>
893Thyra::defaultProductMultiVector(
894 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
895 const ArrayView<
const RCP<MultiVectorBase<Scalar> > > &multiVecs
898 const RCP<DefaultProductMultiVector<Scalar> > pmv =
899 defaultProductMultiVector<Scalar>();
900 pmv->initialize(productSpace, multiVecs);
905template<
class Scalar>
907Thyra::defaultProductMultiVector(
908 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
909 const ArrayView<
const RCP<
const MultiVectorBase<Scalar> > > &multiVecs
912 const RCP<DefaultProductMultiVector<Scalar> > pmv =
913 defaultProductMultiVector<Scalar>();
914 pmv->initialize(productSpace, multiVecs);
919template<
class Scalar>
921Thyra::castOrCreateSingleBlockProductMultiVector(
922 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
923 const RCP<
const MultiVectorBase<Scalar> > &mv
926 const RCP<const ProductMultiVectorBase<Scalar> > pmv =
927 Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(mv);
930 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
934template<
class Scalar>
936Thyra::nonconstCastOrCreateSingleBlockProductMultiVector(
937 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
938 const RCP<MultiVectorBase<Scalar> > &mv
941 const RCP<ProductMultiVectorBase<Scalar> > pmv =
942 Teuchos::rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(mv);
945 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
956#define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_INSTANT(SCALAR) \
958 template class DefaultProductMultiVector<SCALAR >; \
960 template RCP<DefaultProductMultiVector<SCALAR > > \
961 defaultProductMultiVector<SCALAR >(); \
964 template RCP<DefaultProductMultiVector<SCALAR > > \
965 defaultProductMultiVector( \
966 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
967 const int numMembers \
971 template RCP<DefaultProductMultiVector<SCALAR > > \
972 defaultProductMultiVector( \
973 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
974 const ArrayView<const RCP<MultiVectorBase<SCALAR > > > &multiVecs \
978 template RCP<DefaultProductMultiVector<SCALAR > > \
979 defaultProductMultiVector( \
980 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
981 const ArrayView<const RCP<const MultiVectorBase<SCALAR > > > &multiVecs \
985 template RCP<const ProductMultiVectorBase<SCALAR > > \
986 castOrCreateSingleBlockProductMultiVector( \
987 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
988 const RCP<const MultiVectorBase<SCALAR > > &mv \
992 template RCP<ProductMultiVectorBase<SCALAR > > \
993 nonconstCastOrCreateSingleBlockProductMultiVector( \
994 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
995 const RCP<MultiVectorBase<SCALAR > > &mv \
void push_back(const value_type &x)
virtual std::string description() const
Concrete implementation of a product multi-vector.
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const int numMembers)
DefaultProductMultiVector()
Construct to uninitialized.
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
virtual void scaleImpl(Scalar alpha)
void uninitialize()
Uninitialize.
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
std::string description() const
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
RCP< const VectorSpaceBase< Scalar > > range() const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
virtual void assignImpl(Scalar alpha)
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)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< MultiVectorBase< Scalar > > clone_mv() const
bool opSupportedImpl(EOpTransp M_trans) const
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
RCP< const VectorSpaceBase< Scalar > > domain() const
RCP< const VectorBase< Scalar > > colImpl(Ordinal j) const
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
bool blockIsConst(const int k) const
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
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 releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
virtual void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
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.
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 commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
Base interface for product multi-vectors.
virtual Teuchos::RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)=0
Returns a non-persisting non-const view of the zero-based kth block multi-vector.
virtual Teuchos::RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const =0
Returns a non-persisting const view of the (zero-based) kth block multi-vector.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
bool is_null(const std::shared_ptr< T > &p)
bool nonnull(const std::shared_ptr< T > &p)
#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.
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
@ NOTRANS
Use the non-transposed operator.
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)