Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
51 namespace Thyra {
52 
53 
54 // Constructors/initializers/accessors
55 
56 
57 template<class Scalar>
59  :numBlocks_(0)
60 {}
61 
62 
63 template<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 
81 template<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 
91 template<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 
101 template<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 
113 template<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 
128 template<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) {
141  case Teuchos::VERB_NONE:
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 
176 template<class Scalar>
179 {
180  return productSpace_;
181 }
182 
183 
184 template<class Scalar>
186 {
187  return multiVecs_[k].isConst();
188 }
189 
190 
191 template<class Scalar>
194 {
195  return multiVecs_[k].getNonconstObj();
196 }
197 
198 
199 template<class Scalar>
202 {
203  return multiVecs_[k].getConstObj();
204 }
205 
206 
207 // Overridden public functions from MultiVectorBase
208 
209 
210 template<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 
225 template<class Scalar>
228 {
229  return productSpace_;
230 }
231 
232 
233 template<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 
249 template<class Scalar>
251 {
252  for ( int k = 0; k < numBlocks_; ++k ) {
253  multiVecs_[k].getNonconstObj()->assign(alpha);
254  }
255 }
256 
257 
258 template<class Scalar>
260  const MultiVectorBase<Scalar>& mv
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
270  = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
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 
281 template<class Scalar>
283 {
284  for (int k = 0; k < numBlocks_; ++k) {
285  multiVecs_[k].getNonconstObj()->scale(alpha);
286  }
287 }
288 
289 
290 template<class Scalar>
292  Scalar alpha,
293  const MultiVectorBase<Scalar>& mv
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
303  = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
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 
314 template<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) {
333  pvs[i] = Teuchos::ptr_dynamic_cast<const ProductMultiVectorBase<Scalar> >(mv[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());
342  Array<Ptr<const MultiVectorBase<Scalar> > > blocks(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 
356 template<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
369  = Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(Teuchos::rcpFromRef(mv));
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 
387 template<class Scalar>
389  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
409 template<class Scalar>
411  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
435 template<class Scalar>
437  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType>& norms
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 
458 template<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 
470 template<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 
482 template<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 
494 template<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 
506 template<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 
520 template<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 
534 template<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 
596  Array<Ptr<ProductMultiVectorBase<Scalar> > > targ_multi_vecs;
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 
686 template<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 
699 template<class Scalar>
702  ) const
703 {
705  sub_mv );
706  // ToDo: Override this implementation if needed!
707 }
708 
709 
710 template<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 
723 template<class Scalar>
726  )
727 {
729  // ToDo: Override this implementation if needed!
730 }
731 
732 
733 // Overridden protected functions from LinearOpBase
734 
735 
736 template<class Scalar>
738 {
739  return true; // We can do it all!
740 }
741 
742 
743 template<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 
808 template<class Scalar>
809 template<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 
842 template<class Scalar>
843 void DefaultProductMultiVector<Scalar>::assertInitialized() const
844 {
846  is_null(productSpace_), std::logic_error,
847  "Error, this DefaultProductMultiVector object is not intialized!"
848  );
849 }
850 
851 
852 template<class Scalar>
853 void 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 
870 template<class Scalar>
872 Thyra::defaultProductMultiVector()
873 {
874  return Teuchos::rcp(new DefaultProductMultiVector<Scalar>);
875 }
876 
877 
878 template<class Scalar>
880 Thyra::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 
891 template<class Scalar>
893 Thyra::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 
905 template<class Scalar>
907 Thyra::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 
919 template<class Scalar>
921 Thyra::castOrCreateSingleBlockProductMultiVector(
922  const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
923  const RCP<const MultiVectorBase<Scalar> > &mv
924  )
925 {
926  const RCP<const ProductMultiVectorBase<Scalar> > pmv =
927  Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(mv);
928  if (nonnull(pmv))
929  return pmv;
930  return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
931 }
932 
933 
934 template<class Scalar>
936 Thyra::nonconstCastOrCreateSingleBlockProductMultiVector(
937  const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
938  const RCP<MultiVectorBase<Scalar> > &mv
939  )
940 {
941  const RCP<ProductMultiVectorBase<Scalar> > pmv =
942  Teuchos::rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(mv);
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
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 updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
bool is_null(const boost::shared_ptr< T > &p)
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
virtual void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
basic_OSTab< char > OSTab
#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 acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, 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.
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Use the non-transposed operator.
RCP< VectorBase< Scalar > > nonconstColImpl(Ordinal j)
T_To & dyn_cast(T_From &from)
size_type size() const
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
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 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 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...
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) 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.
RCP< const VectorSpaceBase< Scalar > > range() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols) const
void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
RCP< MultiVectorBase< Scalar > > getNonconstMultiVectorBlock(const int k)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
Interface for a collection of column vectors called a multi-vector.
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
Ptr< T > ptr() const
Concrete implementation of a product multi-vector.
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
void update(Scalar alpha, const MultiVectorBase< Scalar > &mv)
virtual std::string description() const
RCP< MultiVectorBase< Scalar > > clone_mv() const
virtual void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
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
void assign(Scalar alpha)
V = alpha.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
bool nonnull(const boost::shared_ptr< T > &p)
void push_back(const value_type &x)
RCP< const VectorSpaceBase< Scalar > > domain() const
void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
void initialize(const RCP< const DefaultProductVectorSpace< Scalar > > &productSpace, const int numMembers)
RCP< const MultiVectorBase< Scalar > > getMultiVectorBlock(const int k) const
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
RCP< const ProductVectorSpaceBase< Scalar > > productSpace() const
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
#define TEUCHOS_ASSERT(assertion_test)
virtual void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
RCP< const VectorBase< Scalar > > colImpl(Ordinal j) const
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Default implementation of dots using RTOps.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Default implementation of update using RTOps.