42 #ifndef THYRA_DEFAULT_PRODUCT_VECTOR_DEF_HPP 43 #define THYRA_DEFAULT_PRODUCT_VECTOR_DEF_HPP 46 #include "Thyra_DefaultProductVector_decl.hpp" 47 #include "Thyra_DefaultProductVectorSpace.hpp" 48 #include "Teuchos_Workspace.hpp" 57 template <
class Scalar>
65 template <
class Scalar>
75 template <
class Scalar>
81 numBlocks_ = productSpace_in->numBlocks();
82 productSpace_ = productSpace_in;
83 vecs_.resize(numBlocks_);
84 for(
int k = 0; k < numBlocks_; ++k )
85 vecs_[k].initialize(createMember(productSpace_in->getBlock(k)));
89 template <
class Scalar>
97 TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(productSpace_in->numBlocks()),
98 as<Ordinal>(vecs.size()) );
100 numBlocks_ = productSpace_in->numBlocks();
101 productSpace_ = productSpace_in;
102 vecs_.resize(numBlocks_);
103 for(
int k = 0; k < numBlocks_; ++k )
104 vecs_[k].initialize(vecs[k]);
108 template <
class Scalar>
116 TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(productSpace_in->numBlocks()),
117 as<Ordinal>(vecs.size()) );
119 numBlocks_ = productSpace_in->numBlocks();
120 productSpace_ = productSpace_in;
121 vecs_.resize(numBlocks_);
122 for(
int k = 0; k < numBlocks_; ++k )
123 vecs_[k].initialize(vecs[k]);
127 template <
class Scalar>
130 productSpace_ = Teuchos::null;
139 template<
class Scalar>
142 const RCP<const VectorSpaceBase<Scalar> > vs = this->space();
143 std::ostringstream oss;
145 << Teuchos::Describable::description()
147 <<
"dim="<<(nonnull(vs) ? vs->dim() : 0)
148 <<
", numBlocks = "<<numBlocks_
154 template<
class Scalar>
156 Teuchos::FancyOStream &out_arg,
157 const Teuchos::EVerbosityLevel verbLevel
160 using Teuchos::FancyOStream;
161 using Teuchos::OSTab;
162 using Teuchos::describe;
163 RCP<FancyOStream> out = rcp(&out_arg,
false);
166 case Teuchos::VERB_NONE:
168 case Teuchos::VERB_DEFAULT:
169 case Teuchos::VERB_LOW:
170 *out << this->description() << std::endl;
172 case Teuchos::VERB_MEDIUM:
173 case Teuchos::VERB_HIGH:
174 case Teuchos::VERB_EXTREME:
177 << Teuchos::Describable::description() <<
"{" 178 <<
"dim=" << this->space()->dim()
182 <<
"numBlocks="<< numBlocks_ << std::endl
183 <<
"Constituent vector objects v[0], v[1], ... v[numBlocks-1]:\n";
185 for(
int k = 0; k < numBlocks_; ++k ) {
186 *out <<
"v["<<k<<
"] = " << describe(*vecs_[k].getConstObj(),verbLevel);
191 TEUCHOS_TEST_FOR_EXCEPT(
true);
199 template <
class Scalar>
205 TEUCHOS_TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
206 TEUCHOS_TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
212 template <
class Scalar>
218 TEUCHOS_TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
219 TEUCHOS_TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
228 template <
class Scalar>
229 RCP<VectorBase<Scalar> >
233 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
235 return vecs_[k].getNonconstObj();
239 template <
class Scalar>
240 RCP<const VectorBase<Scalar> >
244 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
246 return vecs_[k].getConstObj();
253 template <
class Scalar>
254 RCP<const ProductVectorSpaceBase<Scalar> >
257 return productSpace_;
261 template <
class Scalar>
265 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
267 return vecs_[k].isConst();
271 template <
class Scalar>
272 RCP<MultiVectorBase<Scalar> >
275 return getNonconstVectorBlock(k);
279 template <
class Scalar>
280 RCP<const MultiVectorBase<Scalar> >
283 return getVectorBlock(k);
290 template <
class Scalar>
291 RCP< const VectorSpaceBase<Scalar> >
294 return productSpace_;
298 template <
class Scalar>
303 const Ptr<RTOpPack::ReductTarget> &reduct_obj,
315 using Teuchos::ptr_dynamic_cast;
316 using Teuchos::describe;
321 const Ordinal n = productSpace_->dim();
322 const int num_vecs = vecs.size();
323 const int num_targ_vecs = targ_vecs.size();
328 for(
int k = 0; k < num_vecs; ++k) {
329 test_failed = !this->space()->isCompatible(*vecs[k]->space());
330 TEUCHOS_TEST_FOR_EXCEPTION(
332 ,
"DefaultProductVector::applyOp(...): Error vecs["<<k<<
"]->space() = " 333 <<vecs[k]->space()->description()<<
"\' is not compatible with this " 334 <<
"vector space = "<<this->space()->description()<<
"!" 337 for(
int k = 0; k < num_targ_vecs; ++k) {
338 test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
339 TEUCHOS_TEST_FOR_EXCEPTION(
341 ,
"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<
"]->space() = " 342 <<targ_vecs[k]->space()->description()<<
"\' is not compatible with this " 343 <<
"vector space = "<<this->space()->description()<<
"!" 355 Array<RCP<const ProductVectorBase<Scalar> > > vecs_args_store(num_vecs);
356 Array<Ptr<const ProductVectorBase<Scalar> > > vecs_args(num_vecs);
357 for(
int k = 0; k < num_vecs; ++k) {
359 castOrCreateProductVectorBase<Scalar>(rcpFromPtr(vecs[k]));
360 vecs_args[k] = vecs_args_store[k].ptr();
363 Array<RCP<ProductVectorBase<Scalar> > > targ_vecs_args_store(num_targ_vecs);
364 Array<Ptr<ProductVectorBase<Scalar> > > targ_vecs_args(num_targ_vecs);
365 for(
int k = 0; k < num_targ_vecs; ++k) {
366 targ_vecs_args_store[k] =
367 castOrCreateNonconstProductVectorBase<Scalar>(rcpFromPtr(targ_vecs[k]));
368 targ_vecs_args[k] = targ_vecs_args_store[k].ptr();
376 Ordinal num_elements_remaining = dim;
377 const int numBlocks = productSpace_->numBlocks();
378 Array<RCP<const VectorBase<Scalar> > >
379 sub_vecs_rcps(num_vecs);
380 Array<Ptr<const VectorBase<Scalar> > >
382 Array<RCP<VectorBase<Scalar> > >
383 sub_targ_vecs_rcps(num_targ_vecs);
384 Array<Ptr<VectorBase<Scalar> > >
385 sub_targ_vecs(num_targ_vecs);
387 for(
int k = 0; k < numBlocks; ++k) {
388 const Ordinal dim_k = productSpace_->getBlock(k)->dim();
390 for(
int i = 0; i < num_vecs; ++i ) {
391 sub_vecs_rcps[i] = vecs_args[i]->getVectorBlock(k);
392 sub_vecs[i] = sub_vecs_rcps[i].ptr();
395 for(
int j = 0; j < num_targ_vecs; ++j ) {
396 sub_targ_vecs_rcps[j] = targ_vecs_args[j]->getNonconstVectorBlock(k);
397 sub_targ_vecs[j] = sub_targ_vecs_rcps[j].ptr();
399 Thyra::applyOp<Scalar>(
400 op, sub_vecs(), sub_targ_vecs(),
402 global_offset_in + g_off
405 num_elements_remaining -= dim_k;
407 TEUCHOS_TEST_FOR_EXCEPT(!(num_elements_remaining==0));
418 template <
class Scalar>
424 rng = rng_in.full_range() ?
Range1D(0,productSpace_->dim()-1) : rng_in;
425 int kth_vector_space = -1;
427 productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
429 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
432 rng.lbound() + rng.size()
433 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
438 &*vecs_[kth_vector_space].getConstObj()
439 )->acquireDetachedView( rng - kth_global_offset, sub_vec );
453 template <
class Scalar>
458 if( sub_vec->
values().get() == NULL )
return;
459 int kth_vector_space = -1;
461 productSpace_->getVecSpcPoss(sub_vec->
globalOffset(),&kth_vector_space,&kth_global_offset);
463 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
467 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
472 vecs_[kth_vector_space].getConstObj()->releaseDetachedView(sub_vec);
481 template <
class Scalar>
487 rng = rng_in.full_range() ?
Range1D(0,productSpace_->dim()-1) : rng_in;
488 int kth_vector_space = -1;
490 productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
492 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
495 rng.lbound() + rng.size()
496 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
500 vecs_[kth_vector_space].getConstObj()->acquireDetachedView(
501 rng - kth_global_offset, sub_vec
516 template <
class Scalar>
521 if( sub_vec->
values().get() == NULL )
return;
522 int kth_vector_space = -1;
524 productSpace_->getVecSpcPoss(sub_vec->
globalOffset(),&kth_vector_space,&kth_global_offset);
526 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
530 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
535 vecs_[kth_vector_space].getNonconstObj()->commitDetachedView(sub_vec);
544 template <
class Scalar>
549 int kth_vector_space = -1;
551 productSpace_->getVecSpcPoss(sub_vec.
globalOffset(),&kth_vector_space,&kth_global_offset);
553 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
557 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
563 vecs_[kth_vector_space].getNonconstObj()->setSubVector(sub_vec_g);
577 template <
class Scalar>
580 const int num_vecs = vecs_.size();
581 for(
int k = 0; k < num_vecs; ++k) {
582 vecs_[k].getNonconstObj()->assign(alpha);
590 template<
class Scalar>
591 Teuchos::RCP<Thyra::ProductVectorBase<Scalar> >
592 Thyra::castOrCreateNonconstProductVectorBase(
const RCP<VectorBase<Scalar> > v)
594 using Teuchos::rcp_dynamic_cast;
595 using Teuchos::tuple;
596 const RCP<ProductVectorBase<Scalar> > prod_v =
597 rcp_dynamic_cast<ProductVectorBase<Scalar> >(v);
598 if (nonnull(prod_v)) {
601 return defaultProductVector<Scalar>(
602 productVectorSpace<Scalar>(tuple(v->space())()),
608 template<
class Scalar>
609 Teuchos::RCP<const Thyra::ProductVectorBase<Scalar> >
610 Thyra::castOrCreateProductVectorBase(
const RCP<
const VectorBase<Scalar> > v)
612 using Teuchos::rcp_dynamic_cast;
613 using Teuchos::tuple;
614 const RCP<const ProductVectorBase<Scalar> > prod_v =
615 rcp_dynamic_cast<
const ProductVectorBase<Scalar> >(v);
616 if (nonnull(prod_v)) {
619 return defaultProductVector<Scalar>(
620 productVectorSpace<Scalar>(tuple(v->space())()),
630 #define THYRA_DEFAULT_PRODUCT_VECTOR_INSTANT(SCALAR) \ 632 template class DefaultProductVector<SCALAR >; \ 634 template RCP<ProductVectorBase<SCALAR > > \ 635 castOrCreateNonconstProductVectorBase(const RCP<VectorBase<SCALAR > > v); \ 637 template RCP<const ProductVectorBase<SCALAR > > \ 638 castOrCreateProductVectorBase(const RCP<const VectorBase<SCALAR > > v); \ 642 #endif // THYRA_DEFAULT_PRODUCT_VECTOR_DEF_HPP void setSubVector(const RTOpPack::SparseSubVectorT< Scalar > &sub_vec)
Calls setSubVectorImpl().
void setSubVectorImpl(const RTOpPack::SparseSubVectorT< Scalar > &sub_vec)
void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)
Thrown if vector spaces are incompatible.
void applyOpImpl(const RTOpPack::RTOpT< Scalar > &op, const ArrayView< const Ptr< const VectorBase< Scalar > > > &vecs, const ArrayView< const Ptr< VectorBase< Scalar > > > &targ_vecs, const Ptr< RTOpPack::ReductTarget > &reduct_obj, const Ordinal global_offset) const
std::string description() const
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
void uninitialize()
Uninitialize.
Ordinal globalOffset() const
Teuchos_Ordinal subDim() const
const ArrayRCP< const Scalar > values() const
virtual void releaseDetachedVectorViewImpl(RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
void setNonconstBlock(int i, const RCP< VectorBase< Scalar > > &b)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
DefaultProductVector()
Construct to uninitialized.
void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
bool blockIsConst(const int k) const
Abstract interface for finite-dimensional dense vectors.
Teuchos_Ordinal globalOffset() const
virtual void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const
void releaseDetachedVectorViewImpl(RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
RCP< VectorBase< Scalar > > getNonconstVectorBlock(const int k)
void setGlobalOffset(Ordinal globalOffset_in)
void setGlobalOffset(Teuchos_Ordinal globalOffset_in)
RCP< const VectorSpaceBase< Scalar > > space() const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
const ArrayRCP< Scalar > values() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
void assignImpl(Scalar alpha)
virtual void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
void setBlock(int i, const RCP< const VectorBase< Scalar > > &b)
void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace)
Initialize.
virtual void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)