Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultProductVector_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_DEFAULT_PRODUCT_VECTOR_DEF_HPP
43#define THYRA_DEFAULT_PRODUCT_VECTOR_DEF_HPP
44
45
46#include "Thyra_DefaultProductVector_decl.hpp"
47#include "Thyra_DefaultProductVectorSpace.hpp"
48#include "Teuchos_Workspace.hpp"
49
50
51namespace Thyra {
52
53
54// Constructors/initializers/accessors
55
56
57template <class Scalar>
63
64
65template <class Scalar>
67 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in
68 )
69 : numBlocks_(0)
70{
71 initialize(productSpace_in);
72}
73
74
75template <class Scalar>
77 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in
78 )
79{
80 // ToDo: Validate input!
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)));
86}
87
88
89template <class Scalar>
91 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
92 const ArrayView<const RCP<VectorBase<Scalar> > > &vecs
93 )
94{
95 using Teuchos::as;
96#ifdef TEUCHOS_DEBUG
97 TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(productSpace_in->numBlocks()),
98 as<Ordinal>(vecs.size()) );
99#endif
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]);
105}
106
107
108template <class Scalar>
110 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
111 const ArrayView<const RCP<const VectorBase<Scalar> > > &vecs
112 )
113{
114 using Teuchos::as;
115#ifdef TEUCHOS_DEBUG
116 TEUCHOS_ASSERT_EQUALITY( as<Ordinal>(productSpace_in->numBlocks()),
117 as<Ordinal>(vecs.size()) );
118#endif
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]);
124}
125
126
127template <class Scalar>
129{
130 productSpace_ = Teuchos::null;
131 vecs_.resize(0);
132 numBlocks_ = 0;
133}
134
135
136// Overridden from Teuchos::Describable
137
138
139template<class Scalar>
141{
142 const RCP<const VectorSpaceBase<Scalar> > vs = this->space();
143 std::ostringstream oss;
144 oss
146 << "{"
147 << "dim="<<(nonnull(vs) ? vs->dim() : 0)
148 << ", numBlocks = "<<numBlocks_
149 << "}";
150 return oss.str();
151}
152
153
154template<class Scalar>
156 Teuchos::FancyOStream &out_arg,
157 const Teuchos::EVerbosityLevel verbLevel
158 ) const
159{
161 using Teuchos::OSTab;
162 using Teuchos::describe;
163 RCP<FancyOStream> out = rcp(&out_arg,false);
164 OSTab tab(out);
165 switch(verbLevel) {
167 break;
170 *out << this->description() << std::endl;
171 break;
175 {
176 *out
178 << "dim=" << this->space()->dim()
179 << "}\n";
180 OSTab tab2(out);
181 *out
182 << "numBlocks="<< numBlocks_ << std::endl
183 << "Constituent vector objects v[0], v[1], ... v[numBlocks-1]:\n";
184 OSTab tab3(out);
185 for( int k = 0; k < numBlocks_; ++k ) {
186 *out << "v["<<k<<"] = " << describe(*vecs_[k].getConstObj(),verbLevel);
187 }
188 break;
189 }
190 default:
191 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
192 }
193}
194
195
196// Extensions to ProductVectorBase suitable for physically-blocked vectors
197
198
199template <class Scalar>
201 int i, const RCP<const VectorBase<Scalar> >& b
202 )
203{
204#ifdef TEUCHOS_DEBUG
205 TEUCHOS_TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
206 TEUCHOS_TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
207#endif
208 vecs_[i] = b;
209}
210
211
212template <class Scalar>
214 int i, const RCP<VectorBase<Scalar> >& b
215 )
216{
217#ifdef TEUCHOS_DEBUG
218 TEUCHOS_TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
219 TEUCHOS_TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
220#endif
221 vecs_[i] = b;
222}
223
224
225// Overridden from ProductVectorBase
226
227
228template <class Scalar>
231{
232#ifdef TEUCHOS_DEBUG
233 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
234#endif
235 return vecs_[k].getNonconstObj();
236}
237
238
239template <class Scalar>
242{
243#ifdef TEUCHOS_DEBUG
244 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
245#endif
246 return vecs_[k].getConstObj();
247}
248
249
250// Overridden from ProductMultiVectorBase
251
252
253template <class Scalar>
256{
257 return productSpace_;
258}
259
260
261template <class Scalar>
263{
264#ifdef TEUCHOS_DEBUG
265 TEUCHOS_TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
266#endif
267 return vecs_[k].isConst();
268}
269
270
271template <class Scalar>
274{
275 return getNonconstVectorBlock(k);
276}
277
278
279template <class Scalar>
282{
283 return getVectorBlock(k);
284}
285
286
287// Overridden from VectorBase
288
289
290template <class Scalar>
293{
294 return productSpace_;
295}
296
297
298/*template <class Scalar>
299void DefaultProductVector<Scalar>::randomizeImpl(
300 Scalar l,
301 Scalar u
302 )
303{
304 for(int k = 0; k < numBlocks_; ++k) {
305 vecs_[k].getNonconstObj()->randomize(l, u);
306 }
307}*/
308
309
310template <class Scalar>
312 const VectorBase<Scalar>& x
313 )
314{
315 // Apply operation block-by-block if mv is also a ProductVector
318 if (nonnull(pv)) {
319#ifdef TEUCHOS_DEBUG
320 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
321#endif
322 for(int k = 0; k < numBlocks_; ++k) {
323 vecs_[k].getNonconstObj()->abs(*pv->getVectorBlock(k));
324 }
325 } else {
327 }
328}
329
330
331template <class Scalar>
333 const VectorBase<Scalar>& x
334 )
335{
336 // Apply operation block-by-block if mv is also a ProductVector
339 if (nonnull(pv)) {
340#ifdef TEUCHOS_DEBUG
341 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
342#endif
343 for(int k = 0; k < numBlocks_; ++k) {
344 vecs_[k].getNonconstObj()->reciprocal(*pv->getVectorBlock(k));
345 }
346 } else {
348 }
349}
350
351
352template <class Scalar>
354 const VectorBase<Scalar>& x
355 )
356{
357 // Apply operation block-by-block if mv is also a ProductVector
360 if (nonnull(pv)) {
361#ifdef TEUCHOS_DEBUG
362 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
363#endif
364 for(int k = 0; k < numBlocks_; ++k) {
365 vecs_[k].getNonconstObj()->ele_wise_scale(*pv->getVectorBlock(k));
366 }
367 } else {
369 }
370}
371
372
373template <class Scalar>
376 const VectorBase<Scalar>& x
377 ) const
378{
379 // Apply operation block-by-block if mv is also a ProductVector
380 typedef ScalarTraits<Scalar> ST;
383 if (nonnull(pv)) {
384#ifdef TEUCHOS_DEBUG
385 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
386#endif
387 typename ST::magnitudeType norm = ST::magnitude(ST::zero());
388 for(int k = 0; k < numBlocks_; ++k) {
389 typename ST::magnitudeType subNorm
390 = vecs_[k].getConstObj()->norm_2(*pv->getVectorBlock(k));
391 norm += subNorm*subNorm;
392 }
393 return ST::magnitude(ST::squareroot(norm));
394 } else {
396 }
397}
398
399
400template <class Scalar>
402 const RTOpPack::RTOpT<Scalar> &op,
403 const ArrayView<const Ptr<const VectorBase<Scalar> > > &vecs,
404 const ArrayView<const Ptr<VectorBase<Scalar> > > &targ_vecs,
405 const Ptr<RTOpPack::ReductTarget> &reduct_obj,
406 const Ordinal global_offset_in
407 ) const
408{
409
410 // 2008/02/20: rabartl: ToDo: Upgrade Teuchos::Workspace<T> to implicitly
411 // convert to Teuchos::ArrayView<T>. This will allow the calls to
412 // applyOp(...) with sub_vecs and sub_targ_vecs to work without trouble!
413 // For now, I just want to get this done. It is likely that this function
414 // is going to change in major ways soon anyway!
415
416 //using Teuchos::Workspace;
417 using Teuchos::ptr_dynamic_cast;
418 using Teuchos::describe;
419 using Teuchos::null;
420
421 //Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
422
423 const Ordinal n = productSpace_->dim();
424 const int num_vecs = vecs.size();
425 const int num_targ_vecs = targ_vecs.size();
426
427 // Validate the compatibility of the vectors!
428#ifdef TEUCHOS_DEBUG
429 bool test_failed;
430 for(int k = 0; k < num_vecs; ++k) {
431 test_failed = !this->space()->isCompatible(*vecs[k]->space());
434 ,"DefaultProductVector::applyOp(...): Error vecs["<<k<<"]->space() = "
435 <<vecs[k]->space()->description()<<"\' is not compatible with this "
436 <<"vector space = "<<this->space()->description()<<"!"
437 );
438 }
439 for(int k = 0; k < num_targ_vecs; ++k) {
440 test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
443 ,"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<"]->space() = "
444 <<targ_vecs[k]->space()->description()<<"\' is not compatible with this "
445 <<"vector space = "<<this->space()->description()<<"!"
446 );
447 }
448#endif
449
450 //
451 // Dynamic cast each of the vector arguments to the ProductVectorBase interface
452 //
453 // NOTE: If the constituent vector is not a product vector, then a product
454 // vector of one component is created.
455 //
456
457 Array<RCP<const ProductVectorBase<Scalar> > > vecs_args_store(num_vecs);
458 Array<Ptr<const ProductVectorBase<Scalar> > > vecs_args(num_vecs);
459 for(int k = 0; k < num_vecs; ++k) {
460 vecs_args_store[k] =
461 castOrCreateProductVectorBase<Scalar>(rcpFromPtr(vecs[k]));
462 vecs_args[k] = vecs_args_store[k].ptr();
463 }
464
465 Array<RCP<ProductVectorBase<Scalar> > > targ_vecs_args_store(num_targ_vecs);
466 Array<Ptr<ProductVectorBase<Scalar> > > targ_vecs_args(num_targ_vecs);
467 for(int k = 0; k < num_targ_vecs; ++k) {
468 targ_vecs_args_store[k] =
469 castOrCreateNonconstProductVectorBase<Scalar>(rcpFromPtr(targ_vecs[k]));
470 targ_vecs_args[k] = targ_vecs_args_store[k].ptr();
471 }
472
473 //
474 // If we get here, then we will implement the applyOpImpl(...) one vector
475 // block at a time.
476 //
477 const Ordinal dim = n;
478 Ordinal num_elements_remaining = dim;
479 const int numBlocks = productSpace_->numBlocks();
481 sub_vecs_rcps(num_vecs);
483 sub_vecs(num_vecs);
485 sub_targ_vecs_rcps(num_targ_vecs);
487 sub_targ_vecs(num_targ_vecs);
488 Ordinal g_off = 0;
489 for(int k = 0; k < numBlocks; ++k) {
490 const Ordinal dim_k = productSpace_->getBlock(k)->dim();
491 // Fill constituent vectors for block k
492 for( int i = 0; i < num_vecs; ++i ) {
493 sub_vecs_rcps[i] = vecs_args[i]->getVectorBlock(k);
494 sub_vecs[i] = sub_vecs_rcps[i].ptr();
495 }
496 // Fill constituent target vectors for block k
497 for( int j = 0; j < num_targ_vecs; ++j ) {
498 sub_targ_vecs_rcps[j] = targ_vecs_args[j]->getNonconstVectorBlock(k);
499 sub_targ_vecs[j] = sub_targ_vecs_rcps[j].ptr();
500 }
501 Thyra::applyOp<Scalar>(
502 op, sub_vecs(), sub_targ_vecs(),
503 reduct_obj,
504 global_offset_in + g_off
505 );
506 g_off += dim_k;
507 num_elements_remaining -= dim_k;
508 }
509 TEUCHOS_TEST_FOR_EXCEPT(!(num_elements_remaining==0));
510
511}
512
513
514// protected
515
516
517// Overridden protected functions from VectorBase
518
519
520template <class Scalar>
522 const Range1D& rng_in, RTOpPack::ConstSubVectorView<Scalar>* sub_vec
523 ) const
524{
525 const Range1D
526 rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
527 int kth_vector_space = -1;
528 Ordinal kth_global_offset = 0;
529 productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
530#ifdef TEUCHOS_DEBUG
531 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
532#endif
533 if(
534 rng.lbound() + rng.size()
535 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
536 )
537 {
538 // This involves only one sub-vector so just return it.
539 const_cast<const VectorBase<Scalar>*>(
540 &*vecs_[kth_vector_space].getConstObj()
541 )->acquireDetachedView( rng - kth_global_offset, sub_vec );
542 sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
543 }
544 else {
545 // Just let the default implementation handle this. ToDo: In the future
546 // we could manually construct an explicit sub-vector that spanned
547 // two or more constituent vectors but this would be a lot of work.
548 // However, this would require the use of temporary memory but
549 // so what.
551 }
552}
553
554
555template <class Scalar>
558 ) const
559{
560 if( sub_vec->values().get() == NULL ) return;
561 int kth_vector_space = -1;
562 Ordinal kth_global_offset = 0;
563 productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
564#ifdef TEUCHOS_DEBUG
565 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
566#endif
567 if(
568 sub_vec->globalOffset() + sub_vec->subDim()
569 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
570 )
571 {
572 // This sub_vec was extracted from a single constituent vector
573 sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
574 vecs_[kth_vector_space].getConstObj()->releaseDetachedView(sub_vec);
575 }
576 else {
577 // This sub_vec was created by the default implementation!
579 }
580}
581
582
583template <class Scalar>
585 const Range1D& rng_in, RTOpPack::SubVectorView<Scalar>* sub_vec
586 )
587{
588 const Range1D
589 rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
590 int kth_vector_space = -1;
591 Ordinal kth_global_offset = 0;
592 productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
593#ifdef TEUCHOS_DEBUG
594 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
595#endif
596 if(
597 rng.lbound() + rng.size()
598 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
599 )
600 {
601 // This involves only one sub-vector so just return it.
602 vecs_[kth_vector_space].getConstObj()->acquireDetachedView(
603 rng - kth_global_offset, sub_vec
604 );
605 sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
606 }
607 else {
608 // Just let the default implementation handle this. ToDo: In the future
609 // we could manually construct an explicit sub-vector that spanned
610 // two or more constituent vectors but this would be a lot of work.
611 // However, this would require the use of temporary memory but
612 // so what.
614 }
615}
616
617
618template <class Scalar>
621 )
622{
623 if( sub_vec->values().get() == NULL ) return;
624 int kth_vector_space = -1;
625 Ordinal kth_global_offset = 0;
626 productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
627#ifdef TEUCHOS_DEBUG
628 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
629#endif
630 if(
631 sub_vec->globalOffset() + sub_vec->subDim()
632 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
633 )
634 {
635 // This sub_vec was extracted from a single constituent vector
636 sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
637 vecs_[kth_vector_space].getNonconstObj()->commitDetachedView(sub_vec);
638 }
639 else {
640 // This sub_vec was created by the default implementation!
642 }
643}
644
645
646template <class Scalar>
649 )
650{
651 int kth_vector_space = -1;
652 Ordinal kth_global_offset = 0;
653 productSpace_->getVecSpcPoss(sub_vec.globalOffset(),&kth_vector_space,&kth_global_offset);
654#ifdef TEUCHOS_DEBUG
655 TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
656#endif
657 if(
658 sub_vec.globalOffset() + sub_vec.subDim()
659 <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
660 )
661 {
662 // This sub-vector fits into a single constituent vector
663 RTOpPack::SparseSubVectorT<Scalar> sub_vec_g = sub_vec;
664 sub_vec_g.setGlobalOffset( sub_vec_g.globalOffset() - kth_global_offset );
665 vecs_[kth_vector_space].getNonconstObj()->setSubVector(sub_vec_g);
666 }
667 else {
668 // Let the default implementation take care of this. ToDo: In the future
669 // it would be possible to manually set the relevant constituent
670 // vectors with no temp memory allocations.
672 }
673}
674
675
676// Overridden protected functions from MultiVectorBase
677
678
679template <class Scalar>
681{
682 for(int k = 0; k < numBlocks_; ++k) {
683 vecs_[k].getNonconstObj()->assign(alpha);
684 }
685}
686
687
688template <class Scalar>
691 )
692{
693#ifdef TEUCHOS_DEBUG
694 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
695#endif
698 if (nonnull(pv)) {
699#ifdef TEUCHOS_DEBUG
700 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
701#endif
702 for(int k = 0; k < numBlocks_; ++k) {
703 vecs_[k].getNonconstObj()->assign(*pv->getMultiVectorBlock(k));
704 }
705 } else {
707 }
708}
709
710
711template <class Scalar>
713{
714 for(int k = 0; k < numBlocks_; ++k) {
715 vecs_[k].getNonconstObj()->scale(alpha);
716 }
717}
718
719
720template <class Scalar>
722 Scalar alpha,
724 )
725{
726#ifdef TEUCHOS_DEBUG
727 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
728#endif
731 if (nonnull(pv)) {
732#ifdef TEUCHOS_DEBUG
733 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
734#endif
735 for(int k = 0; k < numBlocks_; ++k) {
736 vecs_[k].getNonconstObj()->update(alpha, *pv->getMultiVectorBlock(k));
737 }
738 } else {
740 }
741}
742
743
744template <class Scalar>
746 const ArrayView<const Scalar>& alpha,
747 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
748 const Scalar& beta
749 )
750{
751#ifdef TEUCHOS_DEBUG
752 TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
753#endif
754
755 // Apply operation block-by-block if each element of mv is also a ProductMultiVector
756 bool allCastsSuccessful = true;
758 for (Ordinal i = 0; i < mv.size(); ++i) {
760 if (pvs[i].is_null()) {
761 allCastsSuccessful = false;
762 break;
763 }
764#ifdef TEUCHOS_DEBUG
765 TEUCHOS_ASSERT_EQUALITY(pvs[i]->domain()->dim(), 1);
766 TEUCHOS_ASSERT(productSpace_->isCompatible(*pvs[i]->productSpace()));
767#endif
768 }
769
770 if (allCastsSuccessful) {
771 Array<RCP<const MultiVectorBase<Scalar> > > blocks_rcp(pvs.size());
773 for (int k = 0; k < numBlocks_; ++k) {
774 for (Ordinal i = 0; i < pvs.size(); ++i) {
775 blocks_rcp[i] = pvs[i]->getMultiVectorBlock(k);
776 blocks[i] = blocks_rcp[i].ptr();
777 }
778 vecs_[k].getNonconstObj()->linear_combination(alpha, blocks(), beta);
779 }
780 } else {
782 }
783}
784
785
786template <class Scalar>
788 const MultiVectorBase<Scalar>& mv,
789 const ArrayView<Scalar>& prods
790 ) const
791{
792#ifdef TEUCHOS_DEBUG
793 TEUCHOS_ASSERT_EQUALITY(prods.size(), 1);
794 TEUCHOS_ASSERT_EQUALITY(mv.domain()->dim(), 1);
795#endif
798 if (nonnull(pv)) {
799#ifdef TEUCHOS_DEBUG
800 TEUCHOS_ASSERT(productSpace_->isCompatible(*pv->productSpace()));
801#endif
802 prods[0] = ScalarTraits<Scalar>::zero();
803 for(int k = 0; k < numBlocks_; ++k) {
804 Scalar prod;
805 vecs_[k].getConstObj()->dots(*pv->getMultiVectorBlock(k), Teuchos::arrayView(&prod, 1));
806 prods[0] += prod;
807 }
808 } else {
810 }
811}
812
813
814template <class Scalar>
817 ) const
818{
819#ifdef TEUCHOS_DEBUG
820 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
821#endif
822 typedef ScalarTraits<Scalar> ST;
823 norms[0] = ST::magnitude(ST::zero());
824 for(int k = 0; k < numBlocks_; ++k) {
825 norms[0] += vecs_[k].getConstObj()->norm_1();
826 }
827}
828
829
830template <class Scalar>
833 ) const
834{
835#ifdef TEUCHOS_DEBUG
836 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
837#endif
838 typedef ScalarTraits<Scalar> ST;
839 norms[0] = ST::magnitude(ST::zero());
840 for(int k = 0; k < numBlocks_; ++k) {
841 typename ST::magnitudeType subNorm
842 = vecs_[k].getConstObj()->norm_2();
843 norms[0] += subNorm*subNorm;
844 }
845 norms[0] = ST::magnitude(ST::squareroot(norms[0]));
846}
847
848
849template <class Scalar>
852 ) const
853{
854#ifdef TEUCHOS_DEBUG
855 TEUCHOS_ASSERT_EQUALITY(norms.size(), 1);
856#endif
857 typedef ScalarTraits<Scalar> ST;
858 norms[0] = ST::magnitude(ST::zero());
859 for(int k = 0; k < numBlocks_; ++k) {
860 typename ST::magnitudeType subNorm = vecs_[k].getConstObj()->norm_inf();
861 if (subNorm > norms[0]) {
862 norms[0] = subNorm;
863 }
864 }
865}
866
867
868} // namespace Thyra
869
870
871template<class Scalar>
873Thyra::castOrCreateNonconstProductVectorBase(const RCP<VectorBase<Scalar> > v)
874{
875 using Teuchos::rcp_dynamic_cast;
876 using Teuchos::tuple;
877 const RCP<ProductVectorBase<Scalar> > prod_v =
878 rcp_dynamic_cast<ProductVectorBase<Scalar> >(v);
879 if (nonnull(prod_v)) {
880 return prod_v;
881 }
882 return defaultProductVector<Scalar>(
883 productVectorSpace<Scalar>(tuple(v->space())()),
884 tuple(v)()
885 );
886}
887
888
889template<class Scalar>
891Thyra::castOrCreateProductVectorBase(const RCP<const VectorBase<Scalar> > v)
892{
893 using Teuchos::rcp_dynamic_cast;
894 using Teuchos::tuple;
895 const RCP<const ProductVectorBase<Scalar> > prod_v =
896 rcp_dynamic_cast<const ProductVectorBase<Scalar> >(v);
897 if (nonnull(prod_v)) {
898 return prod_v;
899 }
900 return defaultProductVector<Scalar>(
901 productVectorSpace<Scalar>(tuple(v->space())()),
902 tuple(v)()
903 );
904}
905
906
907//
908// Explicit instant macro
909//
910
911#define THYRA_DEFAULT_PRODUCT_VECTOR_INSTANT(SCALAR) \
912 \
913 template class DefaultProductVector<SCALAR >; \
914 \
915 template RCP<ProductVectorBase<SCALAR > > \
916 castOrCreateNonconstProductVectorBase(const RCP<VectorBase<SCALAR > > v); \
917 \
918 template RCP<const ProductVectorBase<SCALAR > > \
919 castOrCreateProductVectorBase(const RCP<const VectorBase<SCALAR > > v); \
920
921
922
923#endif // THYRA_DEFAULT_PRODUCT_VECTOR_DEF_HPP
Ordinal globalOffset() const
void setGlobalOffset(Ordinal globalOffset_in)
const ArrayRCP< const Scalar > values() const
Teuchos_Ordinal globalOffset() const
Teuchos_Ordinal subDim() const
void setGlobalOffset(Teuchos_Ordinal globalOffset_in)
const ArrayRCP< Scalar > values() const
T * get() const
size_type size() const
size_type size() const
virtual std::string description() const
bool full_range() const
Ordinal size() const
Ordinal lbound() const
void setSubVectorImpl(const RTOpPack::SparseSubVectorT< Scalar > &sub_vec)
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace)
Initialize.
DefaultProductVector()
Construct to uninitialized.
void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const VectorBase< Scalar > > getVectorBlock(const int k) const
void releaseDetachedVectorViewImpl(RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
virtual Teuchos::ScalarTraits< Scalar >::magnitudeType norm2WeightedImpl(const VectorBase< Scalar > &x) const
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
void setNonconstBlock(int i, const RCP< VectorBase< Scalar > > &b)
void setBlock(int i, const RCP< const VectorBase< Scalar > > &b)
virtual void absImpl(const VectorBase< Scalar > &x)
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
virtual void reciprocalImpl(const VectorBase< Scalar > &x)
RCP< VectorBase< Scalar > > getNonconstVectorBlock(const int k)
RCP< const VectorSpaceBase< Scalar > > space() const
virtual void eleWiseScaleImpl(const VectorBase< Scalar > &x)
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
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
void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Thrown if vector spaces are incompatible.
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 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.
Abstract interface for finite-dimensional dense vectors.
void setSubVector(const RTOpPack::SparseSubVectorT< Scalar > &sub_vec)
Calls setSubVectorImpl().
virtual void acquireNonconstDetachedVectorViewImpl(const Range1D &rng, RTOpPack::SubVectorView< Scalar > *sub_vec)
virtual Teuchos::ScalarTraits< Scalar >::magnitudeType norm2WeightedImpl(const VectorBase< Scalar > &x) const
Default implementation of norm_2 (weighted) using RTOps.
virtual void acquireDetachedVectorViewImpl(const Range1D &rng, RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
virtual void reciprocalImpl(const VectorBase< Scalar > &x)
Default implementation of reciprocal using RTOps.
virtual void absImpl(const VectorBase< Scalar > &x)
Default implementation of abs using RTOps.
virtual void eleWiseScaleImpl(const VectorBase< Scalar > &x)
Default implementation of ele_wise_scale using RTOps.
virtual void releaseDetachedVectorViewImpl(RTOpPack::ConstSubVectorView< Scalar > *sub_vec) const
virtual void commitNonconstDetachedVectorViewImpl(RTOpPack::SubVectorView< Scalar > *sub_vec)
#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 nonnull(const std::shared_ptr< T > &p)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
Teuchos::Range1D Range1D
TypeTo as(const TypeFrom &t)
T_To & dyn_cast(T_From &from)