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" 57 template<
class Scalar>
63 template<
class Scalar>
70 TEUCHOS_TEST_FOR_EXCEPT( is_null(productSpace_in) );
71 TEUCHOS_TEST_FOR_EXCEPT( numMembers <= 0 );
73 Array<RCP<MultiVectorBase<Scalar> > > multiVecs;
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);
81 template<
class Scalar>
87 initializeImpl(productSpace_in,multiVecs);
91 template<
class Scalar>
97 initializeImpl(productSpace_in,multiVecs);
101 template<
class Scalar>
104 productSpace_ = Teuchos::null;
105 multiVecs_.resize(0);
113 template<
class Scalar>
116 std::ostringstream oss;
118 << Teuchos::Describable::description()
120 <<
"rangeDim="<<this->range()->dim()
121 <<
",domainDim="<<this->domain()->dim()
122 <<
",numBlocks = "<<numBlocks_
128 template<
class Scalar>
130 FancyOStream &out_arg,
131 const Teuchos::EVerbosityLevel verbLevel
134 using Teuchos::OSTab;
135 using Teuchos::describe;
136 if (verbLevel == Teuchos::VERB_NONE)
138 RCP<FancyOStream> out = rcp(&out_arg,
false);
141 case Teuchos::VERB_NONE:
143 case Teuchos::VERB_DEFAULT:
144 case Teuchos::VERB_LOW:
145 *out << this->description() << std::endl;
147 case Teuchos::VERB_MEDIUM:
148 case Teuchos::VERB_HIGH:
149 case Teuchos::VERB_EXTREME:
152 << Teuchos::Describable::description() <<
"{" 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);
168 TEUCHOS_TEST_FOR_EXCEPT(
true);
176 template<
class Scalar>
177 RCP<const ProductVectorSpaceBase<Scalar> >
180 return productSpace_;
184 template<
class Scalar>
187 return multiVecs_[k].isConst();
191 template<
class Scalar>
192 RCP<MultiVectorBase<Scalar> >
195 return multiVecs_[k].getNonconstObj();
199 template<
class Scalar>
200 RCP<const MultiVectorBase<Scalar> >
203 return multiVecs_[k].getConstObj();
210 template<
class Scalar>
211 RCP<MultiVectorBase<Scalar> >
215 Array<RCP<MultiVectorBase<Scalar> > > blocks;
216 for (
int k = 0; k < numBlocks_; ++k )
217 blocks.push_back(multiVecs_[k].getConstObj()->clone_mv());
218 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
225 template<
class Scalar>
226 RCP< const VectorSpaceBase<Scalar> >
229 return productSpace_;
233 template<
class Scalar>
234 RCP< const VectorSpaceBase<Scalar> >
237 if (is_null(productSpace_))
238 return Teuchos::null;
239 return multiVecs_[0].getConstObj()->domain();
249 template<
class Scalar>
252 for (
int k = 0; k < numBlocks_; ++k ) {
253 multiVecs_[k].getNonconstObj()->assign(alpha);
258 template<
class Scalar>
259 RCP<const VectorBase<Scalar> >
263 Array<RCP<const VectorBase<Scalar> > > cols_;
264 for (
int k = 0; k < numBlocks_; ++k )
265 cols_.push_back(multiVecs_[k].getConstObj()->col(j));
266 return defaultProductVector<Scalar>(productSpace_, cols_());
270 template<
class Scalar>
271 RCP<VectorBase<Scalar> >
275 Array<RCP<VectorBase<Scalar> > > cols_;
276 for (
int k = 0; k < numBlocks_; ++k )
277 cols_.push_back(multiVecs_[k].getNonconstObj()->col(j));
278 return defaultProductVector<Scalar>(productSpace_, cols_());
282 template<
class Scalar>
283 RCP<const MultiVectorBase<Scalar> >
287 Array<RCP<const MultiVectorBase<Scalar> > > blocks;
288 for (
int k = 0; k < numBlocks_; ++k )
289 blocks.push_back(multiVecs_[k].getConstObj()->subView(colRng));
290 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
294 template<
class Scalar>
295 RCP<MultiVectorBase<Scalar> >
299 Array<RCP<MultiVectorBase<Scalar> > > blocks;
300 for (
int k = 0; k < numBlocks_; ++k )
301 blocks.push_back(multiVecs_[k].getNonconstObj()->subView(colRng));
302 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
306 template<
class Scalar>
307 RCP<const MultiVectorBase<Scalar> >
309 const ArrayView<const int> &cols
313 Array<RCP<const MultiVectorBase<Scalar> > > blocks;
314 for (
int k = 0; k < numBlocks_; ++k )
315 blocks.push_back(multiVecs_[k].getConstObj()->subView(cols));
316 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
320 template<
class Scalar>
321 RCP<MultiVectorBase<Scalar> >
323 const ArrayView<const int> &cols
327 Array<RCP<MultiVectorBase<Scalar> > > blocks;
328 for (
int k = 0; k < numBlocks_; ++k )
329 blocks.push_back(multiVecs_[k].getNonconstObj()->subView(cols));
330 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
334 template<
class Scalar>
339 const ArrayView<
const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
340 const Ordinal primary_global_offset_in
344 using Teuchos::ptr_dynamic_cast;
345 using Thyra::applyOp;
350 for (
int j = 0; j < multi_vecs_in.size(); ++j ) {
352 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
353 *this->range(), *multi_vecs_in[j]->range()
356 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
357 *this->domain(), *multi_vecs_in[j]->domain()
360 for (
int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
362 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
363 *this->range(), *targ_multi_vecs_inout[j]->range()
366 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
367 *this->domain(), *targ_multi_vecs_inout[j]->domain()
370 #endif // TEUCHOS_DEBUG 377 bool allProductMultiVectors =
true;
379 Array<Ptr<const ProductMultiVectorBase<Scalar> > > multi_vecs;
380 for (
int j = 0; j < multi_vecs_in.size() && allProductMultiVectors; ++j ) {
382 TEUCHOS_TEST_FOR_EXCEPT( is_null(multi_vecs_in[j]) );
384 const Ptr<const ProductMultiVectorBase<Scalar> >
388 if ( !is_null(multi_vecs_j) ) {
389 multi_vecs.push_back(multi_vecs_j);
392 allProductMultiVectors =
false;
396 Array<Ptr<ProductMultiVectorBase<Scalar> > > targ_multi_vecs;
397 for (
int j = 0; j < targ_multi_vecs_inout.size() && allProductMultiVectors; ++j )
400 TEUCHOS_TEST_FOR_EXCEPT( is_null(targ_multi_vecs_inout[j]) );
402 Ptr<ProductMultiVectorBase<Scalar> >
404 targ_multi_vecs_inout[j]
406 if (!is_null(targ_multi_vecs_j)) {
407 targ_multi_vecs.push_back(targ_multi_vecs_j);
410 allProductMultiVectors =
false;
418 if ( allProductMultiVectors ) {
427 Array<RCP<const MultiVectorBase<Scalar> > >
428 multi_vecs_rcp_block_k(multi_vecs_in.size());
429 Array<Ptr<const MultiVectorBase<Scalar> > >
430 multi_vecs_block_k(multi_vecs_in.size());
431 Array<RCP<MultiVectorBase<Scalar> > >
432 targ_multi_vecs_rcp_block_k(targ_multi_vecs_inout.size());
433 Array<Ptr<MultiVectorBase<Scalar> > >
434 targ_multi_vecs_block_k(targ_multi_vecs_inout.size());
436 Ordinal g_off = primary_global_offset_in;
438 for (
int k = 0; k < numBlocks_; ++k ) {
440 const Ordinal dim_k = productSpace_->getBlock(k)->dim();
444 for (
int j = 0; j < multi_vecs_in.size(); ++j ) {
445 RCP<const MultiVectorBase<Scalar> > multi_vecs_rcp_block_k_j =
447 multi_vecs_rcp_block_k[j] = multi_vecs_rcp_block_k_j;
448 multi_vecs_block_k[j] = multi_vecs_rcp_block_k_j.ptr();
451 for (
int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
452 RCP<MultiVectorBase<Scalar> > targ_multi_vecs_rcp_block_k_j =
453 targ_multi_vecs[j]->getNonconstMultiVectorBlock(k);
454 targ_multi_vecs_rcp_block_k[j] = targ_multi_vecs_rcp_block_k_j;
455 targ_multi_vecs_block_k[j] = targ_multi_vecs_rcp_block_k_j.ptr();
460 Thyra::applyOp<Scalar>(
461 primary_op, multi_vecs_block_k(), targ_multi_vecs_block_k(),
478 primary_op, multi_vecs_in(), targ_multi_vecs_inout(),
479 reduct_objs, primary_global_offset_in);
486 template<
class Scalar>
494 rowRng, colRng, sub_mv );
499 template<
class Scalar>
510 template<
class Scalar>
518 rowRng,colRng,sub_mv );
523 template<
class Scalar>
536 template<
class Scalar>
543 template<
class Scalar>
553 typedef Teuchos::ScalarTraits<Scalar> ST;
554 using Teuchos::dyn_cast;
559 "DefaultProductMultiVector<Scalar>::apply(...)",
560 *
this, M_trans, X_in, &*Y_inout );
573 for (
int k = 0; k < numBlocks_; ++k ) {
575 *multiVecs_[k].getConstObj(), M_trans,
590 for (
int k = 0; k < numBlocks_; ++k ) {
591 RCP<const MultiVectorBase<Scalar> >
592 M_k = multiVecs_[k].getConstObj(),
595 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, beta );
598 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, ST::one() );
608 template<
class Scalar>
609 template<
class MultiVectorType>
612 const ArrayView<
const RCP<MultiVectorType> > &multiVecs
619 TEUCHOS_ASSERT(nonnull(productSpace_in));
620 TEUCHOS_ASSERT_EQUALITY(multiVecs.size(), productSpace_in->numBlocks());
621 #endif // TEUCHOS_DEBUG 622 const RCP<const VectorSpaceBase<Scalar> >
623 theDomain = multiVecs[0]->domain();
624 const int numBlocks = productSpace_in->numBlocks();
626 for (
int k = 0; k < numBlocks; ++k ) {
629 *theDomain, *multiVecs[k]->domain()
633 productSpace_ = productSpace_in;
634 numBlocks_ = numBlocks;
635 multiVecs_.assign(multiVecs.begin(),multiVecs.end());
642 template<
class Scalar>
643 void DefaultProductMultiVector<Scalar>::assertInitialized()
const 645 TEUCHOS_TEST_FOR_EXCEPTION(
646 is_null(productSpace_), std::logic_error,
647 "Error, this DefaultProductMultiVector object is not intialized!" 652 template<
class Scalar>
653 void DefaultProductMultiVector<Scalar>::validateColIndex(
const int j)
const 656 const int domainDim = multiVecs_[0].getConstObj()->domain()->dim();
657 TEUCHOS_TEST_FOR_EXCEPTION(
658 ! ( 0 <= j && j < domainDim ), std::logic_error,
659 "Error, the column index j = " << j <<
" does not fall in the range [0,"<<domainDim<<
"]!" 664 #endif // TEUCHOS_DEBUG 670 template<
class Scalar>
671 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
672 Thyra::defaultProductMultiVector()
674 return Teuchos::rcp(
new DefaultProductMultiVector<Scalar>);
678 template<
class Scalar>
679 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
680 Thyra::defaultProductMultiVector(
681 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
685 RCP<DefaultProductMultiVector<Scalar> > pmv = defaultProductMultiVector<Scalar>();
686 pmv->initialize(productSpace, numMembers);
691 template<
class Scalar>
692 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
693 Thyra::defaultProductMultiVector(
694 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
695 const ArrayView<
const RCP<MultiVectorBase<Scalar> > > &multiVecs
698 const RCP<DefaultProductMultiVector<Scalar> > pmv =
699 defaultProductMultiVector<Scalar>();
700 pmv->initialize(productSpace, multiVecs);
705 template<
class Scalar>
706 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
707 Thyra::defaultProductMultiVector(
708 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
709 const ArrayView<
const RCP<
const MultiVectorBase<Scalar> > > &multiVecs
712 const RCP<DefaultProductMultiVector<Scalar> > pmv =
713 defaultProductMultiVector<Scalar>();
714 pmv->initialize(productSpace, multiVecs);
719 template<
class Scalar>
720 Teuchos::RCP<const Thyra::ProductMultiVectorBase<Scalar> >
721 Thyra::castOrCreateSingleBlockProductMultiVector(
722 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
723 const RCP<
const MultiVectorBase<Scalar> > &mv
726 const RCP<const ProductMultiVectorBase<Scalar> > pmv =
727 Teuchos::rcp_dynamic_cast<
const ProductMultiVectorBase<Scalar> >(mv);
730 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
734 template<
class Scalar>
735 Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> >
736 Thyra::nonconstCastOrCreateSingleBlockProductMultiVector(
737 const RCP<
const DefaultProductVectorSpace<Scalar> > &productSpace,
738 const RCP<MultiVectorBase<Scalar> > &mv
741 const RCP<ProductMultiVectorBase<Scalar> > pmv =
742 Teuchos::rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(mv);
745 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
756 #define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_INSTANT(SCALAR) \ 758 template class DefaultProductMultiVector<SCALAR >; \ 760 template RCP<DefaultProductMultiVector<SCALAR > > \ 761 defaultProductMultiVector<SCALAR >(); \ 764 template RCP<DefaultProductMultiVector<SCALAR > > \ 765 defaultProductMultiVector( \ 766 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \ 767 const int numMembers \ 771 template RCP<DefaultProductMultiVector<SCALAR > > \ 772 defaultProductMultiVector( \ 773 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \ 774 const ArrayView<const RCP<MultiVectorBase<SCALAR > > > &multiVecs \ 778 template RCP<DefaultProductMultiVector<SCALAR > > \ 779 defaultProductMultiVector( \ 780 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \ 781 const ArrayView<const RCP<const MultiVectorBase<SCALAR > > > &multiVecs \ 785 template RCP<const ProductMultiVectorBase<SCALAR > > \ 786 castOrCreateSingleBlockProductMultiVector( \ 787 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \ 788 const RCP<const MultiVectorBase<SCALAR > > &mv \ 792 template RCP<ProductMultiVectorBase<SCALAR > > \ 793 nonconstCastOrCreateSingleBlockProductMultiVector( \ 794 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \ 795 const RCP<MultiVectorBase<SCALAR > > &mv \ 799 #endif // THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP Base interface for product multi-vectors.
#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...
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
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
#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...
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
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
virtual void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
RCP< const VectorBase< Scalar > > colImpl(Ordinal j) const
Use the non-transposed operator.
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
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. ...
EOpTransp real_trans(EOpTransp transp)
Return NOTRANS or TRANS for real scalar valued operators and this also is used for determining struct...
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
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
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.
std::string description() const
virtual void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
Interface for a collection of column vectors called a multi-vector.
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
virtual void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
Concrete implementation of a product multi-vector.
bool opSupportedImpl(EOpTransp M_trans) const
RCP< MultiVectorBase< Scalar > > clone_mv() const
virtual void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< const VectorSpaceBase< Scalar > > range() const
bool blockIsConst(const int k) const
void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const int numMembers)
RCP< const VectorSpaceBase< Scalar > > domain() const
void uninitialize()
Uninitialize.
void assignImpl(Scalar alpha)
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
DefaultProductMultiVector()
Construct to uninitialized.