Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultProductMultiVector_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_MULTI_VECTOR_DEF_HPP
43#define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
44
45#include "Thyra_DefaultProductMultiVector_decl.hpp"
46#include "Thyra_DefaultProductVectorSpace.hpp"
47#include "Thyra_DefaultProductVector.hpp"
48#include "Thyra_AssertOp.hpp"
49
50
51namespace Thyra {
52
53
54// Constructors/initializers/accessors
55
56
57template<class Scalar>
61
62
63template<class Scalar>
65 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
66 const int numMembers
67 )
68{
69#ifdef TEUCHOS_DEBUG
70 TEUCHOS_TEST_FOR_EXCEPT( is_null(productSpace_in) );
71 TEUCHOS_TEST_FOR_EXCEPT( numMembers <= 0 );
72#endif
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);
78}
79
80
81template<class Scalar>
83 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
84 const ArrayView<const RCP<MultiVectorBase<Scalar> > > &multiVecs
85 )
86{
87 initializeImpl(productSpace_in,multiVecs);
88}
89
90
91template<class Scalar>
93 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
94 const ArrayView<const RCP<const MultiVectorBase<Scalar> > > &multiVecs
95 )
96{
97 initializeImpl(productSpace_in,multiVecs);
98}
99
100
101template<class Scalar>
103{
104 productSpace_ = Teuchos::null;
105 multiVecs_.resize(0);
106 numBlocks_ = 0;
107}
108
109
110// Overridden public functions from Teuchos::Describable
111
112
113template<class Scalar>
115{
116 std::ostringstream oss;
117 oss
119 << "{"
120 << "rangeDim="<<this->range()->dim()
121 << ",domainDim="<<this->domain()->dim()
122 << ",numBlocks = "<<numBlocks_
123 << "}";
124 return oss.str();
125}
126
127
128template<class Scalar>
130 FancyOStream &out_arg,
131 const Teuchos::EVerbosityLevel verbLevel
132 ) const
133{
134 using Teuchos::OSTab;
135 using Teuchos::describe;
136 if (verbLevel == Teuchos::VERB_NONE)
137 return;
138 RCP<FancyOStream> out = rcp(&out_arg,false);
139 OSTab tab(out);
140 switch(verbLevel) {
142 break;
143 case Teuchos::VERB_DEFAULT: // fall through
144 case Teuchos::VERB_LOW: // fall through
145 *out << this->description() << std::endl;
146 break;
147 case Teuchos::VERB_MEDIUM: // fall through
148 case Teuchos::VERB_HIGH: // fall through
150 {
151 *out
153 << "rangeDim="<<this->range()->dim()
154 << ",domainDim="<<this->domain()->dim()
155 << "}\n";
156 OSTab tab2(out);
157 *out
158 << "numBlocks="<< numBlocks_ << std::endl
159 << "Constituent multi-vector objects V[0], V[1], ... V[numBlocks-1]:\n";
160 OSTab tab3(out);
161 for( int k = 0; k < numBlocks_; ++k ) {
162 *out << "V["<<k<<"] = "
163 << Teuchos::describe(*multiVecs_[k].getConstObj(),verbLevel);
164 }
165 break;
166 }
167 default:
168 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
169 }
170}
171
172
173// Overridden public functions from ProductMultiVectorBase
174
175
176template<class Scalar>
179{
180 return productSpace_;
181}
182
183
184template<class Scalar>
186{
187 return multiVecs_[k].isConst();
188}
189
190
191template<class Scalar>
194{
195 return multiVecs_[k].getNonconstObj();
196}
197
198
199template<class Scalar>
202{
203 return multiVecs_[k].getConstObj();
204}
205
206
207// Overridden public functions from MultiVectorBase
208
209
210template<class Scalar>
213{
214 assertInitialized();
216 for ( int k = 0; k < numBlocks_; ++k )
217 blocks.push_back(multiVecs_[k].getConstObj()->clone_mv());
218 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
219}
220
221
222// Overriden public functions from LinearOpBase
223
224
225template<class Scalar>
228{
229 return productSpace_;
230}
231
232
233template<class Scalar>
236{
237 if (is_null(productSpace_))
238 return Teuchos::null;
239 return multiVecs_[0].getConstObj()->domain();
240}
241
242
243// protected
244
245
246// Overriden protected functions from MultiVectorBase
247
248
249template<class Scalar>
251{
252 for ( int k = 0; k < numBlocks_; ++k ) {
253 multiVecs_[k].getNonconstObj()->assign(alpha);
254 }
255}
256
257
258template<class Scalar>
261 )
262{
263#ifdef TEUCHOS_DEBUG
264 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), mv.domain()->dim());
265 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
266#endif
267
268 // Apply operation block-by-block if mv is also a ProductMultiVector
271 if (nonnull(pv)) {
272 for (int k = 0; k < numBlocks_; ++k) {
273 multiVecs_[k].getNonconstObj()->assign(*pv->getMultiVectorBlock(k));
274 }
275 } else {
277 }
278}
279
280
281template<class Scalar>
283{
284 for (int k = 0; k < numBlocks_; ++k) {
285 multiVecs_[k].getNonconstObj()->scale(alpha);
286 }
287}
288
289
290template<class Scalar>
292 Scalar alpha,
294 )
295{
296#ifdef TEUCHOS_DEBUG
297 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), mv.domain()->dim());
298 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
299#endif
300
301 // Apply operation block-by-block if mv is also a ProductMultiVector
304 if (nonnull(pv)) {
305 for (int k = 0; k < numBlocks_; ++k) {
306 multiVecs_[k].getNonconstObj()->update(alpha, *pv->getMultiVectorBlock(k));
307 }
308 } else {
310 }
311}
312
313
314template<class Scalar>
316 const ArrayView<const Scalar>& alpha,
317 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > >& mv,
318 const Scalar& beta
319 )
320{
321#ifdef TEUCHOS_DEBUG
322 TEUCHOS_ASSERT_EQUALITY(alpha.size(), mv.size());
323 for (Ordinal i = 0; i < mv.size(); ++i) {
324 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), mv[i]->domain()->dim());
325 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv[i]->range()));
326 }
327#endif
328
329 // Apply operation block-by-block if each element of mv is also a ProductMultiVector
330 bool allCastsSuccessful = true;
332 for (Ordinal i = 0; i < mv.size(); ++i) {
334 if (pvs[i].is_null()) {
335 allCastsSuccessful = false;
336 break;
337 }
338 }
339
340 if (allCastsSuccessful) {
341 Array<RCP<const MultiVectorBase<Scalar> > > blocks_rcp(pvs.size());
343 for (int k = 0; k < numBlocks_; ++k) {
344 for (Ordinal i = 0; i < pvs.size(); ++i) {
345 blocks_rcp[i] = pvs[i]->getMultiVectorBlock(k);
346 blocks[i] = blocks_rcp[i].ptr();
347 }
348 multiVecs_[k].getNonconstObj()->linear_combination(alpha, blocks(), beta);
349 }
350 } else {
352 }
353}
354
355
356template<class Scalar>
358 const MultiVectorBase<Scalar>& mv,
359 const ArrayView<Scalar>& prods
360 ) const
361{
362#ifdef TEUCHOS_DEBUG
363 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), mv.domain()->dim());
364 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), prods.size());
365 TEUCHOS_ASSERT(productSpace_->isCompatible(*mv.range()));
366#endif
367 // Apply operation block-by-block if mv is also a ProductMultiVector
370 if (nonnull(pv)) {
371 for (Ordinal i = 0; i < prods.size(); ++i)
372 prods[i] = ScalarTraits<Scalar>::zero();
373
374 Array<Scalar> subProds(prods.size());
375 for (int k = 0; k < numBlocks_; ++k) {
376 multiVecs_[k].getConstObj()->dots(*pv->getMultiVectorBlock(k), subProds());
377 for (Ordinal i = 0; i < prods.size(); ++i) {
378 prods[i] += subProds[i];
379 }
380 }
381 } else {
383 }
384}
385
386
387template<class Scalar>
390 ) const
391{
392#ifdef TEUCHOS_DEBUG
393 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), norms.size());
394#endif
395 typedef ScalarTraits<Scalar> ST;
396 for (Ordinal i = 0; i < norms.size(); ++i)
397 norms[i] = ST::magnitude(ST::zero());
398
399 Array<typename ST::magnitudeType> subNorms(norms.size());
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];
404 }
405 }
406}
407
408
409template<class Scalar>
412 ) const
413{
414#ifdef TEUCHOS_DEBUG
415 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), norms.size());
416#endif
417 typedef ScalarTraits<Scalar> ST;
418 for (Ordinal i = 0; i < norms.size(); ++i)
419 norms[i] = ST::magnitude(ST::zero());
420
421 Array<typename ST::magnitudeType> subNorms(norms.size());
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];
426 }
427 }
428
429 for (Ordinal i = 0; i < norms.size(); ++i) {
430 norms[i] = ST::magnitude(ST::squareroot(norms[i]));
431 }
432}
433
434
435template<class Scalar>
438 ) const
439{
440#ifdef TEUCHOS_DEBUG
441 TEUCHOS_ASSERT_EQUALITY(this->domain()->dim(), norms.size());
442#endif
443 typedef ScalarTraits<Scalar> ST;
444 for (Ordinal i = 0; i < norms.size(); ++i)
445 norms[i] = ST::magnitude(ST::zero());
446
447 Array<typename ST::magnitudeType> subNorms(norms.size());
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];
453 }
454 }
455}
456
457
458template<class Scalar>
461{
462 validateColIndex(j);
464 for ( int k = 0; k < numBlocks_; ++k )
465 cols_.push_back(multiVecs_[k].getConstObj()->col(j));
466 return defaultProductVector<Scalar>(productSpace_, cols_());
467}
468
469
470template<class Scalar>
473{
474 validateColIndex(j);
476 for ( int k = 0; k < numBlocks_; ++k )
477 cols_.push_back(multiVecs_[k].getNonconstObj()->col(j));
478 return defaultProductVector<Scalar>(productSpace_, cols_());
479}
480
481
482template<class Scalar>
485{
486 assertInitialized();
488 for ( int k = 0; k < numBlocks_; ++k )
489 blocks.push_back(multiVecs_[k].getConstObj()->subView(colRng));
490 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
491}
492
493
494template<class Scalar>
497{
498 assertInitialized();
500 for ( int k = 0; k < numBlocks_; ++k )
501 blocks.push_back(multiVecs_[k].getNonconstObj()->subView(colRng));
502 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
503}
504
505
506template<class Scalar>
509 const ArrayView<const int> &cols
510 ) const
511{
512 assertInitialized();
514 for ( int k = 0; k < numBlocks_; ++k )
515 blocks.push_back(multiVecs_[k].getConstObj()->subView(cols));
516 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
517}
518
519
520template<class Scalar>
523 const ArrayView<const int> &cols
524 )
525{
526 assertInitialized();
528 for ( int k = 0; k < numBlocks_; ++k )
529 blocks.push_back(multiVecs_[k].getNonconstObj()->subView(cols));
530 return defaultProductMultiVector<Scalar>(productSpace_, blocks());
531}
532
533
534template<class Scalar>
536 const RTOpPack::RTOpT<Scalar> &primary_op,
537 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs_in,
538 const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs_inout,
539 const ArrayView<const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
540 const Ordinal primary_global_offset_in
541 ) const
542{
543
544 using Teuchos::ptr_dynamic_cast;
545 using Thyra::applyOp;
546
547 assertInitialized();
548
549#ifdef TEUCHOS_DEBUG
550 for ( int j = 0; j < multi_vecs_in.size(); ++j ) {
552 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
553 *this->range(), *multi_vecs_in[j]->range()
554 );
556 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
557 *this->domain(), *multi_vecs_in[j]->domain()
558 );
559 }
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()
564 );
566 "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
567 *this->domain(), *targ_multi_vecs_inout[j]->domain()
568 );
569 }
570#endif // TEUCHOS_DEBUG
571
572 //
573 // Try to dynamic cast all of the multi-vector objects to the
574 // ProductMultiVectoBase interface.
575 //
576
577 bool allProductMultiVectors = true;
578
580 for ( int j = 0; j < multi_vecs_in.size() && allProductMultiVectors; ++j ) {
581#ifdef TEUCHOS_DEBUG
582 TEUCHOS_TEST_FOR_EXCEPT( is_null(multi_vecs_in[j]) );
583#endif
585 multi_vecs_j = ptr_dynamic_cast<const ProductMultiVectorBase<Scalar> >(
586 multi_vecs_in[j]
587 );
588 if ( !is_null(multi_vecs_j) ) {
589 multi_vecs.push_back(multi_vecs_j);
590 }
591 else {
592 allProductMultiVectors = false;
593 }
594 }
595
597 for ( int j = 0; j < targ_multi_vecs_inout.size() && allProductMultiVectors; ++j )
598 {
599#ifdef TEUCHOS_DEBUG
600 TEUCHOS_TEST_FOR_EXCEPT( is_null(targ_multi_vecs_inout[j]) );
601#endif
603 targ_multi_vecs_j = ptr_dynamic_cast<ProductMultiVectorBase<Scalar> >(
604 targ_multi_vecs_inout[j]
605 );
606 if (!is_null(targ_multi_vecs_j)) {
607 targ_multi_vecs.push_back(targ_multi_vecs_j);
608 }
609 else {
610 allProductMultiVectors = false;
611 }
612 }
613
614 //
615 // Do the reduction operations
616 //
617
618 if ( allProductMultiVectors ) {
619
620 // All of the multi-vector objects support the ProductMultiVectorBase
621 // interface so we can do the reductions block by block. Note, this is
622 // not the most efficient implementation in an SPMD program but this is
623 // easy to code up and use!
624
625 // We must set up temporary arrays for the pointers to the MultiVectorBase
626 // blocks for each block of objects! What a mess!
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());
635
636 Ordinal g_off = primary_global_offset_in;
637
638 for ( int k = 0; k < numBlocks_; ++k ) {
639
640 const Ordinal dim_k = productSpace_->getBlock(k)->dim();
641
642 // Fill the MultiVector array objects for this block
643
644 for ( int j = 0; j < multi_vecs_in.size(); ++j ) {
645 RCP<const MultiVectorBase<Scalar> > multi_vecs_rcp_block_k_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();
649 }
650
651 for ( int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
652 RCP<MultiVectorBase<Scalar> > targ_multi_vecs_rcp_block_k_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();
656 }
657
658 // Apply the RTOp object to the MultiVectors for this block
659
660 Thyra::applyOp<Scalar>(
661 primary_op, multi_vecs_block_k(), targ_multi_vecs_block_k(),
662 reduct_objs,
663 g_off
664 );
665
666 g_off += dim_k;
667 }
668
669 }
670 else {
671
672 // All of the multi-vector objects do not support the
673 // ProductMultiVectorBase interface but if we got here (in debug mode)
674 // then the spaces said that they are compatible so fall back on the
675 // column-by-column implementation that will work correctly in serial.
676
678 primary_op, multi_vecs_in(), targ_multi_vecs_inout(),
679 reduct_objs, primary_global_offset_in);
680
681 }
682
683}
684
685
686template<class Scalar>
688 const Range1D &rowRng,
689 const Range1D &colRng,
691 ) const
692{
694 rowRng, colRng, sub_mv );
695 // ToDo: Override this implementation if needed!
696}
697
698
699template<class Scalar>
702 ) const
703{
705 sub_mv );
706 // ToDo: Override this implementation if needed!
707}
708
709
710template<class Scalar>
712 const Range1D &rowRng,
713 const Range1D &colRng,
715 )
716{
718 rowRng,colRng,sub_mv );
719 // ToDo: Override this implementation if needed!
720}
721
722
723template<class Scalar>
731
732
733// Overridden protected functions from LinearOpBase
734
735
736template<class Scalar>
738{
739 return true; // We can do it all!
740}
741
742
743template<class Scalar>
745 const EOpTransp M_trans,
746 const MultiVectorBase<Scalar> &X_in,
747 const Ptr<MultiVectorBase<Scalar> > &Y_inout,
748 const Scalar alpha,
749 const Scalar beta
750 ) const
751{
752
754 using Teuchos::dyn_cast;
755 using Thyra::apply;
756
757#ifdef TEUCHOS_DEBUG
759 "DefaultProductMultiVector<Scalar>::apply(...)",
760 *this, M_trans, X_in, &*Y_inout );
761#endif
762
763 if ( real_trans(M_trans) == NOTRANS ) {
764 //
765 // Y = b*Y + a*M*X
766 //
767 // =>
768 //
769 // Y[k] = b*Y[k] + a*M[k]*X, k = 0...numBlocks-1
770 //
772 &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
773 for ( int k = 0; k < numBlocks_; ++k ) {
774 Thyra::apply(
775 *multiVecs_[k].getConstObj(), M_trans,
776 X_in, Y.getNonconstMultiVectorBlock(k).ptr(),
777 alpha, beta );
778 }
779 }
780 else {
781 //
782 // Y = b*Y + a*trans(M)*X
783 //
784 // =>
785 //
786 // Y = b*Y + sum( a * trans(M[k]) * X[k], k=0...numBlocks-1 )
787 //
789 &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
790 for ( int k = 0; k < numBlocks_; ++k ) {
792 M_k = multiVecs_[k].getConstObj(),
793 X_k = X.getMultiVectorBlock(k);
794 if ( 0 == k ) {
795 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, beta );
796 }
797 else {
798 Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, ST::one() );
799 }
800 }
801 }
802}
803
804
805// private
806
807
808template<class Scalar>
809template<class MultiVectorType>
811 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
812 const ArrayView<const RCP<MultiVectorType> > &multiVecs
813 )
814{
815 // This function provides the "strong" guarantee (i.e. if an exception is
816 // thrown, then *this will be left in the original state as before the
817 // function was called)!
818#ifdef TEUCHOS_DEBUG
819 TEUCHOS_ASSERT(nonnull(productSpace_in));
820 TEUCHOS_ASSERT_EQUALITY(multiVecs.size(), productSpace_in->numBlocks());
821#endif // TEUCHOS_DEBUG
823 theDomain = multiVecs[0]->domain();
824 const int numBlocks = productSpace_in->numBlocks();
825#ifdef TEUCHOS_DEBUG
826 for ( int k = 0; k < numBlocks; ++k ) {
829 *theDomain, *multiVecs[k]->domain()
830 );
831 }
832#endif
833 productSpace_ = productSpace_in;
834 numBlocks_ = numBlocks;
835 multiVecs_.assign(multiVecs.begin(),multiVecs.end());
836}
837
838
839#ifdef TEUCHOS_DEBUG
840
841
842template<class Scalar>
843void DefaultProductMultiVector<Scalar>::assertInitialized() const
844{
846 is_null(productSpace_), std::logic_error,
847 "Error, this DefaultProductMultiVector object is not intialized!"
848 );
849}
850
851
852template<class Scalar>
853void DefaultProductMultiVector<Scalar>::validateColIndex(const int j) const
854{
855 assertInitialized();
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<<"]!"
860 );
861}
862
863
864#endif // TEUCHOS_DEBUG
865
866
867} // namespace Thyra
868
869
870template<class Scalar>
872Thyra::defaultProductMultiVector()
873{
874 return Teuchos::rcp(new DefaultProductMultiVector<Scalar>);
875}
876
877
878template<class Scalar>
880Thyra::defaultProductMultiVector(
881 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
882 const int numMembers
883 )
884{
885 RCP<DefaultProductMultiVector<Scalar> > pmv = defaultProductMultiVector<Scalar>();
886 pmv->initialize(productSpace, numMembers);
887 return pmv;
888}
889
890
891template<class Scalar>
893Thyra::defaultProductMultiVector(
894 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
895 const ArrayView<const RCP<MultiVectorBase<Scalar> > > &multiVecs
896 )
897{
898 const RCP<DefaultProductMultiVector<Scalar> > pmv =
899 defaultProductMultiVector<Scalar>();
900 pmv->initialize(productSpace, multiVecs);
901 return pmv;
902}
903
904
905template<class Scalar>
907Thyra::defaultProductMultiVector(
908 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
909 const ArrayView<const RCP<const MultiVectorBase<Scalar> > > &multiVecs
910 )
911{
912 const RCP<DefaultProductMultiVector<Scalar> > pmv =
913 defaultProductMultiVector<Scalar>();
914 pmv->initialize(productSpace, multiVecs);
915 return pmv;
916}
917
918
919template<class Scalar>
921Thyra::castOrCreateSingleBlockProductMultiVector(
922 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
923 const RCP<const MultiVectorBase<Scalar> > &mv
924 )
925{
926 const RCP<const ProductMultiVectorBase<Scalar> > pmv =
928 if (nonnull(pmv))
929 return pmv;
930 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
931}
932
933
934template<class Scalar>
936Thyra::nonconstCastOrCreateSingleBlockProductMultiVector(
937 const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
938 const RCP<MultiVectorBase<Scalar> > &mv
939 )
940{
941 const RCP<ProductMultiVectorBase<Scalar> > pmv =
943 if (nonnull(pmv))
944 return pmv;
945 return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
946}
947
948
949//
950// Explicit instantiation macro
951//
952// Must be expanded from within the Thyra namespace!
953//
954
955
956#define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_INSTANT(SCALAR) \
957 \
958 template class DefaultProductMultiVector<SCALAR >; \
959 \
960 template RCP<DefaultProductMultiVector<SCALAR > > \
961 defaultProductMultiVector<SCALAR >(); \
962 \
963 \
964 template RCP<DefaultProductMultiVector<SCALAR > > \
965 defaultProductMultiVector( \
966 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
967 const int numMembers \
968 ); \
969 \
970 \
971 template RCP<DefaultProductMultiVector<SCALAR > > \
972 defaultProductMultiVector( \
973 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
974 const ArrayView<const RCP<MultiVectorBase<SCALAR > > > &multiVecs \
975 ); \
976 \
977 \
978 template RCP<DefaultProductMultiVector<SCALAR > > \
979 defaultProductMultiVector( \
980 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
981 const ArrayView<const RCP<const MultiVectorBase<SCALAR > > > &multiVecs \
982 ); \
983 \
984 \
985 template RCP<const ProductMultiVectorBase<SCALAR > > \
986 castOrCreateSingleBlockProductMultiVector( \
987 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
988 const RCP<const MultiVectorBase<SCALAR > > &mv \
989 ); \
990 \
991 \
992 template RCP<ProductMultiVectorBase<SCALAR > > \
993 nonconstCastOrCreateSingleBlockProductMultiVector( \
994 const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
995 const RCP<MultiVectorBase<SCALAR > > &mv \
996 );
997
998
999#endif // THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
size_type size() const
size_type size() const
void push_back(const value_type &x)
virtual std::string description() const
Ptr< T > ptr() const
Concrete implementation of a product multi-vector.
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const int numMembers)
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
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)
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)
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
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)
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)