Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultProductMultiVector_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
00030 #define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
00031 
00032 #include "Thyra_DefaultProductMultiVector_decl.hpp"
00033 #include "Thyra_DefaultProductVectorSpace.hpp"
00034 #include "Thyra_DefaultProductVector.hpp"
00035 #include "Thyra_AssertOp.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00041 // Constructors/initializers/accessors
00042 
00043 
00044 template<class Scalar>
00045 DefaultProductMultiVector<Scalar>::DefaultProductMultiVector()
00046   :numBlocks_(0)
00047 {}
00048 
00049 
00050 template<class Scalar>
00051 void DefaultProductMultiVector<Scalar>::initialize(
00052   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
00053   const int numMembers
00054   )
00055 {
00056 #ifdef TEUCHOS_DEBUG
00057   TEST_FOR_EXCEPT( is_null(productSpace_in) );
00058   TEST_FOR_EXCEPT( numMembers <= 0 );
00059 #endif
00060   Array<RCP<MultiVectorBase<Scalar> > > multiVecs;
00061   const int numBlocks = productSpace_in->numBlocks();
00062   for ( int k = 0; k < numBlocks; ++k )
00063     multiVecs.push_back(createMembers(productSpace_in->getBlock(k),numMembers));
00064   initialize(productSpace_in,multiVecs);
00065 }
00066 
00067 
00068 template<class Scalar>
00069 void DefaultProductMultiVector<Scalar>::initialize(
00070   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
00071   const ArrayView<const RCP<MultiVectorBase<Scalar> > > &multiVecs
00072   )
00073 {
00074   initializeImpl(productSpace_in,multiVecs);
00075 }
00076 
00077 
00078 template<class Scalar>
00079 void DefaultProductMultiVector<Scalar>::initialize(
00080   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
00081   const ArrayView<const RCP<const MultiVectorBase<Scalar> > > &multiVecs
00082   )
00083 {
00084   initializeImpl(productSpace_in,multiVecs);
00085 }
00086 
00087 
00088 template<class Scalar>
00089 void DefaultProductMultiVector<Scalar>::uninitialize()
00090 {
00091   productSpace_ = Teuchos::null;
00092   multiVecs_.resize(0);
00093   numBlocks_ = 0;
00094 }
00095 
00096 
00097 // Overridden public functions from Teuchos::Describable
00098 
00099                                                 
00100 template<class Scalar>
00101 std::string DefaultProductMultiVector<Scalar>::description() const
00102 {
00103   std::ostringstream oss;
00104   oss
00105     << Teuchos::Describable::description()
00106     << "{"
00107     << "rangeDim="<<this->range()->dim()
00108     << ",domainDim="<<this->domain()->dim()
00109     << ",numBlocks = "<<numBlocks_
00110     << "}";
00111   return oss.str();
00112 }
00113 
00114 
00115 template<class Scalar>
00116 void DefaultProductMultiVector<Scalar>::describe(
00117   FancyOStream &out_arg,
00118   const Teuchos::EVerbosityLevel verbLevel
00119   ) const
00120 {
00121   typedef Teuchos::ScalarTraits<Scalar>  ST;
00122   using Teuchos::OSTab;
00123   using Teuchos::describe;
00124   if (verbLevel == Teuchos::VERB_NONE)
00125     return;
00126   RCP<FancyOStream> out = rcp(&out_arg,false);
00127   OSTab tab(out);
00128   switch(verbLevel) {
00129     case Teuchos::VERB_NONE:
00130       break;
00131     case Teuchos::VERB_DEFAULT: // fall through
00132     case Teuchos::VERB_LOW: // fall through
00133       *out << this->description() << std::endl;
00134       break;
00135     case Teuchos::VERB_MEDIUM: // fall through
00136     case Teuchos::VERB_HIGH: // fall through
00137     case Teuchos::VERB_EXTREME:
00138     {
00139       *out
00140         << Teuchos::Describable::description() << "{"
00141         << "rangeDim="<<this->range()->dim()
00142         << ",domainDim="<<this->domain()->dim()
00143         << "}\n";
00144       OSTab tab2(out);
00145       *out
00146         <<  "numBlocks="<< numBlocks_ << std::endl
00147         <<  "Constituent multi-vector objects V[0], V[1], ... V[numBlocks-1]:\n";
00148       OSTab tab3(out);
00149       for( int k = 0; k < numBlocks_; ++k ) {
00150         *out << "V["<<k<<"] = "
00151              << Teuchos::describe(*multiVecs_[k].getConstObj(),verbLevel);
00152       }
00153       break;
00154     }
00155     default:
00156       TEST_FOR_EXCEPT(true); // Should never get here!
00157   }
00158 }
00159 
00160 
00161 // Overridden public functions from ProductMultiVectorBase
00162 
00163 
00164 template<class Scalar>
00165 RCP<const ProductVectorSpaceBase<Scalar> >
00166 DefaultProductMultiVector<Scalar>::productSpace() const
00167 {
00168   return productSpace_;
00169 }
00170 
00171 
00172 template<class Scalar>
00173 bool DefaultProductMultiVector<Scalar>::blockIsConst(const int k) const
00174 {
00175   return multiVecs_[k].isConst();
00176 }
00177 
00178 
00179 template<class Scalar>
00180 RCP<MultiVectorBase<Scalar> >
00181 DefaultProductMultiVector<Scalar>::getNonconstMultiVectorBlock(const int k)
00182 {
00183   return multiVecs_[k].getNonconstObj();
00184 }
00185 
00186 
00187 template<class Scalar>
00188 RCP<const MultiVectorBase<Scalar> >
00189 DefaultProductMultiVector<Scalar>::getMultiVectorBlock(const int k) const
00190 {
00191   return multiVecs_[k].getConstObj();
00192 }
00193 
00194 
00195 // Overridden public functions from MultiVectorBase
00196 
00197 
00198 template<class Scalar>
00199 RCP<MultiVectorBase<Scalar> >
00200 DefaultProductMultiVector<Scalar>::clone_mv() const
00201 {
00202   assertInitialized();
00203   Array<RCP<MultiVectorBase<Scalar> > > blocks;
00204   for ( int k = 0; k < numBlocks_; ++k )
00205     blocks.push_back(multiVecs_[k].getConstObj()->clone_mv());
00206   return defaultProductMultiVector<Scalar>(productSpace_, blocks());
00207 }
00208 
00209 
00210 // Overriden public functions from LinearOpBase
00211 
00212 
00213 template<class Scalar>
00214 RCP< const VectorSpaceBase<Scalar> >
00215 DefaultProductMultiVector<Scalar>::range() const
00216 {
00217   return productSpace_;
00218 }
00219 
00220 
00221 template<class Scalar>
00222 RCP< const VectorSpaceBase<Scalar> >
00223 DefaultProductMultiVector<Scalar>::domain() const
00224 {
00225   if (is_null(productSpace_))
00226     return Teuchos::null;
00227   return multiVecs_[0].getConstObj()->domain();
00228 }
00229 
00230 
00231 // protected
00232 
00233 
00234 // Overriden protected functions from MultiVectorBase
00235 
00236 
00237 template<class Scalar>
00238 RCP<const VectorBase<Scalar> >
00239 DefaultProductMultiVector<Scalar>::colImpl(Ordinal j) const
00240 {
00241   validateColIndex(j);
00242   Array<RCP<const VectorBase<Scalar> > > cols_;
00243   for ( int k = 0; k < numBlocks_; ++k )
00244     cols_.push_back(multiVecs_[k].getConstObj()->col(j));
00245   return defaultProductVector<Scalar>(productSpace_,&cols_[0]);
00246 }
00247 
00248 
00249 template<class Scalar>
00250 RCP<VectorBase<Scalar> >
00251 DefaultProductMultiVector<Scalar>::nonconstColImpl(Ordinal j)
00252 {
00253   validateColIndex(j);
00254   Array<RCP<VectorBase<Scalar> > > cols_;
00255   for ( int k = 0; k < numBlocks_; ++k )
00256     cols_.push_back(multiVecs_[k].getNonconstObj()->col(j));
00257   return defaultProductVector<Scalar>(productSpace_,&cols_[0]);
00258 }
00259 
00260 
00261 template<class Scalar>
00262 RCP<const MultiVectorBase<Scalar> >
00263 DefaultProductMultiVector<Scalar>::contigSubViewImpl( const Range1D& colRng ) const
00264 {
00265   assertInitialized();
00266   Array<RCP<const MultiVectorBase<Scalar> > > blocks;
00267   for ( int k = 0; k < numBlocks_; ++k )
00268     blocks.push_back(multiVecs_[k].getConstObj()->subView(colRng));
00269   return defaultProductMultiVector<Scalar>(productSpace_, blocks());
00270 }
00271 
00272 
00273 template<class Scalar>
00274 RCP<MultiVectorBase<Scalar> >
00275 DefaultProductMultiVector<Scalar>::nonconstContigSubViewImpl( const Range1D& colRng )
00276 {
00277   assertInitialized();
00278   Array<RCP<MultiVectorBase<Scalar> > > blocks;
00279   for ( int k = 0; k < numBlocks_; ++k )
00280     blocks.push_back(multiVecs_[k].getNonconstObj()->subView(colRng));
00281   return defaultProductMultiVector<Scalar>(productSpace_, blocks());
00282 }
00283 
00284 
00285 template<class Scalar>
00286 RCP<const MultiVectorBase<Scalar> >
00287 DefaultProductMultiVector<Scalar>::nonContigSubViewImpl(
00288   const ArrayView<const int> &cols
00289   ) const
00290 {
00291   assertInitialized();
00292   Array<RCP<const MultiVectorBase<Scalar> > > blocks;
00293   for ( int k = 0; k < numBlocks_; ++k )
00294     blocks.push_back(multiVecs_[k].getConstObj()->subView(cols));
00295   return defaultProductMultiVector<Scalar>(productSpace_, blocks());
00296 }
00297 
00298 
00299 template<class Scalar>
00300 RCP<MultiVectorBase<Scalar> >
00301 DefaultProductMultiVector<Scalar>::nonconstNonContigSubViewImpl(
00302   const ArrayView<const int> &cols
00303   )
00304 {
00305   assertInitialized();
00306   Array<RCP<MultiVectorBase<Scalar> > > blocks;
00307   for ( int k = 0; k < numBlocks_; ++k )
00308     blocks.push_back(multiVecs_[k].getNonconstObj()->subView(cols));
00309   return defaultProductMultiVector<Scalar>(productSpace_, blocks());
00310 }
00311 
00312 
00313 template<class Scalar>
00314 void DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(
00315   const RTOpPack::RTOpT<Scalar> &primary_op,
00316   const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &multi_vecs_in,
00317   const ArrayView<const Ptr<MultiVectorBase<Scalar> > > &targ_multi_vecs_inout,
00318   const ArrayView<const Ptr<RTOpPack::ReductTarget> > &reduct_objs,
00319   const Ordinal primary_global_offset_in
00320   ) const
00321 {
00322 
00323   using Teuchos::ptr_dynamic_cast;
00324   using Thyra::applyOp;
00325 
00326   assertInitialized();
00327 
00328 #ifdef TEUCHOS_DEBUG
00329   for ( int j = 0; j < multi_vecs_in.size(); ++j ) {
00330     THYRA_ASSERT_VEC_SPACES(
00331       "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
00332       *this->range(), *multi_vecs_in[j]->range()
00333       );
00334     THYRA_ASSERT_VEC_SPACES(
00335       "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
00336       *this->domain(), *multi_vecs_in[j]->domain()
00337       );
00338   }
00339   for ( int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
00340     THYRA_ASSERT_VEC_SPACES(
00341       "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
00342       *this->range(), *targ_multi_vecs_inout[j]->range()
00343       );
00344     THYRA_ASSERT_VEC_SPACES(
00345       "DefaultProductMultiVector<Scalar>::mvMultiReductApplyOpImpl(...)",
00346       *this->domain(), *targ_multi_vecs_inout[j]->domain()
00347       );
00348   }
00349 #endif //  TEUCHOS_DEBUG
00350 
00351   //
00352   // Try to dynamic cast all of the multi-vector objects to the
00353   // ProductMultiVectoBase interface.
00354   //
00355 
00356   bool allProductMultiVectors = true;
00357 
00358   Array<Ptr<const ProductMultiVectorBase<Scalar> > > multi_vecs;
00359   for ( int j = 0; j < multi_vecs_in.size() && allProductMultiVectors; ++j ) {
00360 #ifdef TEUCHOS_DEBUG
00361     TEST_FOR_EXCEPT( is_null(multi_vecs_in[j]) );
00362 #endif
00363     const Ptr<const ProductMultiVectorBase<Scalar> >
00364       multi_vecs_j = ptr_dynamic_cast<const ProductMultiVectorBase<Scalar> >(
00365         multi_vecs_in[j]
00366         );
00367     if ( !is_null(multi_vecs_j) ) {
00368       multi_vecs.push_back(multi_vecs_j);
00369     }
00370     else {
00371       allProductMultiVectors = false;
00372     }
00373   }
00374 
00375   Array<Ptr<ProductMultiVectorBase<Scalar> > > targ_multi_vecs;
00376   for ( int j = 0; j < targ_multi_vecs_inout.size() && allProductMultiVectors; ++j )
00377   {
00378 #ifdef TEUCHOS_DEBUG
00379     TEST_FOR_EXCEPT( is_null(targ_multi_vecs_inout[j]) );
00380 #endif
00381     Ptr<ProductMultiVectorBase<Scalar> >
00382       targ_multi_vecs_j = ptr_dynamic_cast<ProductMultiVectorBase<Scalar> >(
00383         targ_multi_vecs_inout[j]
00384         );
00385     if (!is_null(targ_multi_vecs_j)) {
00386       targ_multi_vecs.push_back(targ_multi_vecs_j);
00387     }
00388     else {
00389       allProductMultiVectors = false;
00390     }
00391   }
00392 
00393   //
00394   // Do the reduction operations
00395   //
00396   
00397   if ( allProductMultiVectors ) {
00398 
00399     // All of the multi-vector objects support the ProductMultiVectorBase
00400     // interface so we can do the reductions block by block.  Note, this is
00401     // not the most efficient implementation in an SPMD program but this is
00402     // easy to code up and use!
00403 
00404     // We must set up temporary arrays for the pointers to the MultiVectorBase
00405     // blocks for each block of objects!  What a mess!
00406     Array<RCP<const MultiVectorBase<Scalar> > >
00407       multi_vecs_rcp_block_k(multi_vecs_in.size());
00408     Array<Ptr<const MultiVectorBase<Scalar> > >
00409       multi_vecs_block_k(multi_vecs_in.size());
00410     Array<RCP<MultiVectorBase<Scalar> > >
00411       targ_multi_vecs_rcp_block_k(targ_multi_vecs_inout.size());
00412     Array<Ptr<MultiVectorBase<Scalar> > >
00413       targ_multi_vecs_block_k(targ_multi_vecs_inout.size());
00414 
00415     Ordinal g_off = primary_global_offset_in;
00416 
00417     for ( int k = 0; k < numBlocks_; ++k ) {
00418 
00419       const Ordinal dim_k = productSpace_->getBlock(k)->dim();
00420 
00421       // Fill the MultiVector array objects for this block
00422 
00423       for ( int j = 0; j < multi_vecs_in.size(); ++j ) {
00424         RCP<const MultiVectorBase<Scalar> > multi_vecs_rcp_block_k_j =
00425           multi_vecs[j]->getMultiVectorBlock(k);
00426         multi_vecs_rcp_block_k[j] = multi_vecs_rcp_block_k_j;
00427         multi_vecs_block_k[j] = multi_vecs_rcp_block_k_j.ptr();
00428       }
00429 
00430       for ( int j = 0; j < targ_multi_vecs_inout.size(); ++j ) {
00431         RCP<MultiVectorBase<Scalar> > targ_multi_vecs_rcp_block_k_j =
00432           targ_multi_vecs[j]->getNonconstMultiVectorBlock(k);
00433         targ_multi_vecs_rcp_block_k[j] = targ_multi_vecs_rcp_block_k_j;
00434         targ_multi_vecs_block_k[j] = targ_multi_vecs_rcp_block_k_j.ptr();
00435       }
00436 
00437       // Apply the RTOp object to the MultiVectors for this block
00438 
00439       Thyra::applyOp<Scalar>(
00440         primary_op, multi_vecs_block_k(), targ_multi_vecs_block_k(),
00441         reduct_objs,
00442         g_off
00443         );
00444 
00445       g_off += dim_k;
00446     }
00447 
00448   }
00449   else {
00450 
00451     // All of the multi-vector objects do not support the
00452     // ProductMultiVectorBase interface but if we got here (in debug mode)
00453     // then the spaces said that they are compatible so fall back on the
00454     // column-by-column implementation that will work correctly in serial.
00455 
00456     MultiVectorDefaultBase<Scalar>::mvMultiReductApplyOpImpl(
00457       primary_op, multi_vecs_in(), targ_multi_vecs_inout(),
00458       reduct_objs, primary_global_offset_in);
00459 
00460   }
00461 
00462 }
00463 
00464 
00465 template<class Scalar>
00466 void DefaultProductMultiVector<Scalar>::acquireDetachedMultiVectorViewImpl(
00467   const Range1D &rowRng,
00468   const Range1D &colRng,
00469   RTOpPack::ConstSubMultiVectorView<Scalar> *sub_mv
00470   ) const
00471 {
00472   return MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(
00473     rowRng, colRng, sub_mv );
00474   // ToDo: Override this implementation if needed!
00475 }
00476 
00477 
00478 template<class Scalar>
00479 void DefaultProductMultiVector<Scalar>::releaseDetachedMultiVectorViewImpl(
00480   RTOpPack::ConstSubMultiVectorView<Scalar>* sub_mv
00481   ) const
00482 {
00483   return MultiVectorDefaultBase<Scalar>::releaseDetachedMultiVectorViewImpl(
00484     sub_mv );
00485   // ToDo: Override this implementation if needed!
00486 }
00487 
00488 
00489 template<class Scalar>
00490 void DefaultProductMultiVector<Scalar>::acquireNonconstDetachedMultiVectorViewImpl(
00491   const Range1D &rowRng,
00492   const Range1D &colRng,
00493   RTOpPack::SubMultiVectorView<Scalar> *sub_mv
00494   )
00495 {
00496   return MultiVectorDefaultBase<Scalar>::acquireNonconstDetachedMultiVectorViewImpl(
00497     rowRng,colRng,sub_mv );
00498   // ToDo: Override this implementation if needed!
00499 }
00500 
00501 
00502 template<class Scalar>
00503 void DefaultProductMultiVector<Scalar>::commitNonconstDetachedMultiVectorViewImpl(
00504   RTOpPack::SubMultiVectorView<Scalar>* sub_mv
00505   )
00506 {
00507   return MultiVectorDefaultBase<Scalar>::commitNonconstDetachedMultiVectorViewImpl(sub_mv);
00508   // ToDo: Override this implementation if needed!
00509 }
00510 
00511 
00512 // Overridden protected functions from LinearOpBase
00513 
00514 
00515 template<class Scalar>
00516 bool DefaultProductMultiVector<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00517 {
00518   return true; // We can do it all!
00519 }
00520 
00521 
00522 template<class Scalar>
00523 void DefaultProductMultiVector<Scalar>::applyImpl(
00524   const EOpTransp M_trans,
00525   const MultiVectorBase<Scalar> &X_in,
00526   const Ptr<MultiVectorBase<Scalar> > &Y_inout,
00527   const Scalar alpha,
00528   const Scalar beta
00529   ) const
00530 {
00531 
00532   typedef Teuchos::ScalarTraits<Scalar> ST;
00533   using Teuchos::dyn_cast;
00534   using Thyra::apply;
00535 
00536 #ifdef TEUCHOS_DEBUG
00537   THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(
00538     "DefaultProductMultiVector<Scalar>::apply(...)",
00539     *this, M_trans, X_in, &*Y_inout );
00540 #endif
00541 
00542   if ( real_trans(M_trans) == NOTRANS ) {
00543     //
00544     // Y = b*Y + a*M*X
00545     //
00546     //    =>
00547     //
00548     // Y[k] = b*Y[k] + a*M[k]*X, k = 0...numBlocks-1
00549     //
00550     ProductMultiVectorBase<Scalar>
00551       &Y = dyn_cast<ProductMultiVectorBase<Scalar> >(*Y_inout);
00552     for ( int k = 0; k < numBlocks_; ++k ) {
00553       Thyra::apply(
00554         *multiVecs_[k].getConstObj(), M_trans,
00555         X_in, Y.getNonconstMultiVectorBlock(k).ptr(),
00556         alpha, beta ); 
00557     }
00558   }
00559   else {
00560     //
00561     // Y = b*Y + a*trans(M)*X
00562     //
00563     //    =>
00564     //
00565     // Y = b*Y + sum( a * trans(M[k]) * X[k], k=0...numBlocks-1 )
00566     //
00567     const ProductMultiVectorBase<Scalar>
00568       &X = dyn_cast<const ProductMultiVectorBase<Scalar> >(X_in);
00569     for ( int k = 0; k < numBlocks_; ++k ) {
00570       RCP<const MultiVectorBase<Scalar> >
00571         M_k = multiVecs_[k].getConstObj(),
00572         X_k = X.getMultiVectorBlock(k);
00573       if ( 0 == k ) {
00574         Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, beta ); 
00575       }
00576       else {
00577         Thyra::apply( *M_k, M_trans, *X_k, Y_inout.ptr(), alpha, ST::one() ); 
00578       }
00579     }
00580   }
00581 }
00582 
00583 
00584 // private
00585 
00586 
00587 template<class Scalar>
00588 template<class MultiVectorType>
00589 void DefaultProductMultiVector<Scalar>::initializeImpl(
00590   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace_in,
00591   const ArrayView<const RCP<MultiVectorType> > &multiVecs
00592   )
00593 {
00594   // This function provides the "strong" guarantee (i.e. if an exception is
00595   // thrown, then *this will be left in the original state as before the
00596   // function was called)!
00597 #ifdef TEUCHOS_DEBUG
00598   TEUCHOS_ASSERT(nonnull(productSpace_in));
00599   TEUCHOS_ASSERT_EQUALITY(multiVecs.size(), productSpace_in->numBlocks());
00600 #endif // TEUCHOS_DEBUG
00601   const RCP<const VectorSpaceBase<Scalar> >
00602     theDomain = multiVecs[0]->domain();
00603   const int numBlocks = productSpace_in->numBlocks();
00604 #ifdef TEUCHOS_DEBUG
00605   for ( int k = 0; k < numBlocks; ++k ) {
00606     THYRA_ASSERT_VEC_SPACES(
00607       Teuchos::TypeNameTraits<DefaultProductMultiVector<Scalar> >::name(),
00608       *theDomain, *multiVecs[k]->domain()
00609       );
00610   }
00611 #endif
00612   productSpace_ = productSpace_in;
00613   numBlocks_ = numBlocks;
00614   multiVecs_.assign(multiVecs.begin(),multiVecs.end());
00615 }
00616 
00617 
00618 #ifdef TEUCHOS_DEBUG
00619 
00620 
00621 template<class Scalar>
00622 void DefaultProductMultiVector<Scalar>::assertInitialized() const
00623 {
00624   TEST_FOR_EXCEPTION(
00625     is_null(productSpace_), std::logic_error,
00626     "Error, this DefaultProductMultiVector object is not intialized!"
00627     );
00628 }
00629 
00630 
00631 template<class Scalar>
00632 void DefaultProductMultiVector<Scalar>::validateColIndex(const int j) const
00633 {
00634   assertInitialized();
00635   const int domainDim = multiVecs_[0].getConstObj()->domain()->dim();
00636   TEST_FOR_EXCEPTION(
00637     ! ( 0 <= j && j < domainDim ), std::logic_error,
00638     "Error, the column index j = " << j << " does not fall in the range [0,"<<domainDim<<"]!"
00639     );
00640 }
00641 
00642 
00643 #endif // TEUCHOS_DEBUG
00644 
00645 
00646 } // namespace Thyra
00647 
00648 
00649 template<class Scalar>
00650 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
00651 Thyra::defaultProductMultiVector()
00652 {
00653   return Teuchos::rcp(new DefaultProductMultiVector<Scalar>);
00654 }
00655 
00656 
00657 template<class Scalar>
00658 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
00659 Thyra::defaultProductMultiVector(
00660   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
00661   const int numMembers
00662   )
00663 {
00664   RCP<DefaultProductMultiVector<Scalar> > pmv = defaultProductMultiVector<Scalar>();
00665   pmv->initialize(productSpace, numMembers);
00666   return pmv;
00667 }
00668 
00669 
00670 template<class Scalar>
00671 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
00672 Thyra::defaultProductMultiVector(
00673   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
00674   const ArrayView<const RCP<MultiVectorBase<Scalar> > > &multiVecs
00675   )
00676 {
00677   const RCP<DefaultProductMultiVector<Scalar> > pmv =
00678     defaultProductMultiVector<Scalar>();
00679   pmv->initialize(productSpace, multiVecs);
00680   return pmv;
00681 }
00682 
00683 
00684 template<class Scalar>
00685 Teuchos::RCP<Thyra::DefaultProductMultiVector<Scalar> >
00686 Thyra::defaultProductMultiVector(
00687   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
00688   const ArrayView<const RCP<const MultiVectorBase<Scalar> > > &multiVecs
00689   )
00690 {
00691   const RCP<DefaultProductMultiVector<Scalar> > pmv =
00692     defaultProductMultiVector<Scalar>();
00693   pmv->initialize(productSpace, multiVecs);
00694   return pmv;
00695 }
00696 
00697 
00698 template<class Scalar>
00699 Teuchos::RCP<const Thyra::ProductMultiVectorBase<Scalar> >
00700 Thyra::castOrCreateSingleBlockProductMultiVector(
00701   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
00702   const RCP<const MultiVectorBase<Scalar> > &mv
00703   )
00704 {
00705   const RCP<const ProductMultiVectorBase<Scalar> > pmv =
00706     Teuchos::rcp_dynamic_cast<const ProductMultiVectorBase<Scalar> >(mv);
00707   if (nonnull(pmv))
00708     return pmv;
00709   return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
00710 }
00711 
00712 
00713 template<class Scalar>
00714 Teuchos::RCP<Thyra::ProductMultiVectorBase<Scalar> >
00715 Thyra::nonconstCastOrCreateSingleBlockProductMultiVector(
00716   const RCP<const DefaultProductVectorSpace<Scalar> > &productSpace,
00717   const RCP<MultiVectorBase<Scalar> > &mv
00718   )
00719 {
00720   const RCP<ProductMultiVectorBase<Scalar> > pmv =
00721     Teuchos::rcp_dynamic_cast<ProductMultiVectorBase<Scalar> >(mv);
00722   if (nonnull(pmv))
00723     return pmv;
00724   return defaultProductMultiVector<Scalar>(productSpace, Teuchos::tuple(mv)());
00725 }
00726 
00727 
00728 //
00729 // Explicit instantiation macro
00730 //
00731 // Must be expanded from within the Thyra namespace!
00732 //
00733 
00734 
00735 #define THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_INSTANT(SCALAR) \
00736   \
00737   template class DefaultProductMultiVector<SCALAR >; \
00738   \
00739   template RCP<DefaultProductMultiVector<SCALAR > >  \
00740   defaultProductMultiVector<SCALAR >();  \
00741     \
00742     \
00743   template RCP<DefaultProductMultiVector<SCALAR > >  \
00744   defaultProductMultiVector(  \
00745     const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace,  \
00746     const int numMembers  \
00747     );  \
00748     \
00749     \
00750   template RCP<DefaultProductMultiVector<SCALAR > >  \
00751   defaultProductMultiVector(  \
00752     const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace,  \
00753     const ArrayView<const RCP<MultiVectorBase<SCALAR > > > &multiVecs  \
00754     );  \
00755     \
00756     \
00757   template RCP<DefaultProductMultiVector<SCALAR > >  \
00758   defaultProductMultiVector(  \
00759     const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace,  \
00760     const ArrayView<const RCP<const MultiVectorBase<SCALAR > > > &multiVecs  \
00761     );  \
00762     \
00763     \
00764   template RCP<const ProductMultiVectorBase<SCALAR > >  \
00765   castOrCreateSingleBlockProductMultiVector(  \
00766     const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
00767     const RCP<const MultiVectorBase<SCALAR > > &mv  \
00768     );  \
00769     \
00770     \
00771   template RCP<ProductMultiVectorBase<SCALAR > >  \
00772   nonconstCastOrCreateSingleBlockProductMultiVector(  \
00773     const RCP<const DefaultProductVectorSpace<SCALAR > > &productSpace, \
00774     const RCP<MultiVectorBase<SCALAR > > &mv  \
00775     );
00776 
00777 
00778 #endif // THYRA_DEFAULT_PRODUCT_MULTI_VECTOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines