Thyra_DefaultProductVector.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_PRODUCT_VECTOR_HPP
00030 #define THYRA_PRODUCT_VECTOR_HPP
00031 
00032 
00033 #include "Thyra_DefaultProductVectorDecl.hpp"
00034 #include "Thyra_DefaultProductVectorSpace.hpp"
00035 #include "Teuchos_Workspace.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00041 // Constructors/initializers/accessors
00042 
00043 
00044 template <class Scalar>
00045 DefaultProductVector<Scalar>::DefaultProductVector()
00046 {
00047   uninitialize();
00048 }
00049 
00050 
00051 template <class Scalar>
00052 DefaultProductVector<Scalar>::DefaultProductVector(
00053   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00054   )
00055 {
00056   initialize(productSpace);
00057 }
00058 
00059 
00060 template <class Scalar>
00061 DefaultProductVector<Scalar>::DefaultProductVector(
00062   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00063   ,const Teuchos::RCP<VectorBase<Scalar> >                      vecs[]
00064   )
00065 {
00066   initialize(productSpace,vecs);
00067 }
00068 
00069 
00070 template <class Scalar>
00071 DefaultProductVector<Scalar>::DefaultProductVector(
00072   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00073   ,const Teuchos::RCP<const VectorBase<Scalar> >                vecs[]
00074   )
00075 {
00076   initialize(productSpace,vecs);
00077 }
00078 
00079 
00080 template <class Scalar>
00081 void DefaultProductVector<Scalar>::initialize(
00082   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00083   )
00084 {
00085   // ToDo: Validate input!
00086   numBlocks_ = productSpace->numBlocks();
00087   productSpace_ = productSpace;
00088   vecs_.resize(numBlocks_);
00089   for( int k = 0; k < numBlocks_; ++k )
00090     vecs_[k].initialize(createMember(productSpace->getBlock(k)));
00091 }
00092 
00093 
00094 template <class Scalar>
00095 void DefaultProductVector<Scalar>::initialize(
00096   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00097   ,const Teuchos::RCP<VectorBase<Scalar> >                      vecs[]
00098   )
00099 {
00100   // ToDo: Validate input!
00101   numBlocks_ = productSpace->numBlocks();
00102   productSpace_ = productSpace;
00103   vecs_.resize(numBlocks_);
00104   for( int k = 0; k < numBlocks_; ++k )
00105     vecs_[k].initialize(vecs[k]);
00106 }
00107 
00108 
00109 template <class Scalar>
00110 void DefaultProductVector<Scalar>::initialize(
00111   const Teuchos::RCP<const DefaultProductVectorSpace<Scalar> >  &productSpace
00112   ,const Teuchos::RCP<const VectorBase<Scalar> >                vecs[]
00113   )
00114 {
00115   // ToDo: Validate input!
00116   numBlocks_ = productSpace->numBlocks();
00117   productSpace_ = productSpace;
00118   vecs_.resize(numBlocks_);
00119   for( int k = 0; k < numBlocks_; ++k )
00120     vecs_[k].initialize(vecs[k]);
00121 }
00122 
00123 
00124 template <class Scalar>
00125 void DefaultProductVector<Scalar>::uninitialize()
00126 {
00127   productSpace_ = Teuchos::null;
00128   vecs_.resize(0);
00129   numBlocks_ = 0;
00130 }
00131 
00132 
00133 // Overridden from Teuchos::Describable
00134 
00135                                                 
00136 template<class Scalar>
00137 std::string DefaultProductVector<Scalar>::description() const
00138 {
00139   std::ostringstream oss;
00140   oss
00141     << Teuchos::Describable::description()
00142     << "{"
00143     << "dim="<<this->space()->dim()
00144     << ",numBlocks = "<<numBlocks_
00145     << "}";
00146   return oss.str();
00147 }
00148 
00149 
00150 template<class Scalar>
00151 void DefaultProductVector<Scalar>::describe(
00152   Teuchos::FancyOStream &out_arg,
00153   const Teuchos::EVerbosityLevel verbLevel
00154   ) const
00155 {
00156   typedef Teuchos::ScalarTraits<Scalar>  ST;
00157   using Teuchos::RCP;
00158   using Teuchos::FancyOStream;
00159   using Teuchos::OSTab;
00160   using Teuchos::describe;
00161   RCP<FancyOStream> out = rcp(&out_arg,false);
00162   OSTab tab(out);
00163   switch(verbLevel) {
00164     case Teuchos::VERB_DEFAULT:
00165     case Teuchos::VERB_LOW:
00166       *out << this->description() << std::endl;
00167       break;
00168     case Teuchos::VERB_MEDIUM:
00169     case Teuchos::VERB_HIGH:
00170     case Teuchos::VERB_EXTREME:
00171     {
00172       *out
00173         << Teuchos::Describable::description() << "{"
00174         << "dim=" << this->space()->dim()
00175         << "}\n";
00176       OSTab tab(out);
00177       *out
00178         <<  "numBlocks="<< numBlocks_ << std::endl
00179         <<  "Constituent vector objects v[0], v[1], ... v[numBlocks-1]:\n";
00180       tab.incrTab();
00181       for( int k = 0; k < numBlocks_; ++k ) {
00182         *out << "v["<<k<<"] = " << describe(*vecs_[k].getConstObj(),verbLevel);
00183       }
00184       break;
00185     }
00186     default:
00187       TEST_FOR_EXCEPT(true); // Should never get here!
00188   }
00189 }
00190 
00191 
00192 // Overridden from ProductVectorBase
00193 
00194 
00195 template <class Scalar>
00196 Teuchos::RCP<VectorBase<Scalar> >
00197 DefaultProductVector<Scalar>::getNonconstVectorBlock(const int k)
00198 {
00199   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00200   return vecs_[k].getNonconstObj();
00201 }
00202 
00203 
00204 template <class Scalar>
00205 Teuchos::RCP<const VectorBase<Scalar> >
00206 DefaultProductVector<Scalar>::getVectorBlock(const int k) const
00207 {
00208   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00209   return vecs_[k].getConstObj();
00210 }
00211 
00212 
00213 // Overridden from ProductMultiVectorBase
00214 
00215 
00216 template <class Scalar>
00217 Teuchos::RCP<const ProductVectorSpaceBase<Scalar> >
00218 DefaultProductVector<Scalar>::productSpace() const
00219 {
00220   return productSpace_;
00221 }
00222 
00223 
00224 template <class Scalar>
00225 bool DefaultProductVector<Scalar>::blockIsConst(const int k) const
00226 {
00227   TEST_FOR_EXCEPT( k < 0 || numBlocks_-1 < k);
00228   return vecs_[k].isConst();
00229 }
00230 
00231 
00232 template <class Scalar>
00233 Teuchos::RCP<MultiVectorBase<Scalar> >
00234 DefaultProductVector<Scalar>::getNonconstMultiVectorBlock(const int k)
00235 {
00236   return getNonconstVectorBlock(k);
00237 }
00238 
00239 
00240 template <class Scalar>
00241 void DefaultProductVector<Scalar>
00242 ::setBlock(int i, 
00243            const Teuchos::RCP<const VectorBase<Scalar> >& b)
00244 {
00245   TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
00246   TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
00247   vecs_[i] = b;
00248 }
00249 
00250 
00251 template <class Scalar>
00252 void DefaultProductVector<Scalar>
00253 ::setNonconstBlock(int i, 
00254                    const Teuchos::RCP<VectorBase<Scalar> >& b)
00255 {
00256   TEST_FOR_EXCEPT(i < 0 || i >= numBlocks_);
00257   TEST_FOR_EXCEPT(!productSpace_->getBlock(i)->isCompatible(*(b->space())));
00258   vecs_[i] = b;
00259 }
00260 
00261 
00262 template <class Scalar>
00263 Teuchos::RCP<const MultiVectorBase<Scalar> >
00264 DefaultProductVector<Scalar>::getMultiVectorBlock(const int k) const
00265 {
00266   return getVectorBlock(k);
00267 }
00268 
00269 
00270 // Overridden from VectorBase
00271 
00272 
00273 template <class Scalar>
00274 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00275 DefaultProductVector<Scalar>::space() const
00276 {
00277   return productSpace_;
00278 }
00279 
00280 
00281 template <class Scalar>
00282 void DefaultProductVector<Scalar>::applyOp(
00283   const RTOpPack::RTOpT<Scalar>    &op
00284   ,const int                       num_vecs
00285   ,const VectorBase<Scalar>*const  vecs[]
00286   ,const int                       num_targ_vecs
00287   ,VectorBase<Scalar>*const        targ_vecs[]
00288   ,RTOpPack::ReductTarget          *reduct_obj
00289   ,const Index                     first_ele_offset_in
00290   ,const Index                     sub_dim_in
00291   ,const Index                     global_offset_in
00292   ) const
00293 {
00294   using Teuchos::Workspace;
00295   using Teuchos::RCP;
00296   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00297   //
00298   const Index n = productSpace_->dim();
00299   // Validate the compatibility of the vectors!
00300 #ifdef TEUCHOS_DEBUG
00301   TEST_FOR_EXCEPTION(
00302     !(0 <= first_ele_offset_in && first_ele_offset_in < n), std::out_of_range
00303     ,"DefaultProductVector::applyOp(...): Error, "
00304     "first_ele_offset_in = " << first_ele_offset_in << " is not in range [0,"<<(n-1)<<"]" );
00305   TEST_FOR_EXCEPTION(
00306     global_offset_in < 0, std::invalid_argument
00307     ,"DefaultProductVector::applyOp(...): Error, "
00308     "global_offset_in = " << global_offset_in << " is not acceptable" );
00309   TEST_FOR_EXCEPTION(
00310     sub_dim_in >= 0 && sub_dim_in - first_ele_offset_in > n, std::length_error
00311     ,"DefaultProductVector::applyOp(...): Error, "
00312     "global_offset_in = " << global_offset_in << ", sub_dim_in = " << sub_dim_in
00313     << "first_ele_offset_in = " << first_ele_offset_in << " and n = " << n
00314     << " are not compatible" );
00315   bool test_failed;
00316   for(int k = 0; k < num_vecs; ++k) {
00317     test_failed = !this->space()->isCompatible(*vecs[k]->space());
00318     TEST_FOR_EXCEPTION(
00319       test_failed, Exceptions::IncompatibleVectorSpaces
00320       ,"DefaultProductVector::applyOp(...): Error vecs["<<k<<"]->space() "
00321       <<"of type \'"<<typeName(*vecs[k]->space())<<"\' is not compatible with this "
00322       <<"\'VectorSpaceBlocked\' vector space!"
00323       );
00324   }
00325   for(int k = 0; k < num_targ_vecs; ++k) {
00326     test_failed = !this->space()->isCompatible(*targ_vecs[k]->space());
00327     TEST_FOR_EXCEPTION(
00328       test_failed, Exceptions::IncompatibleVectorSpaces
00329       ,"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<"]->space() "
00330       <<"of type \'"<<typeName(*vecs[k]->space())<<"\' is not compatible with this "
00331       <<"\'VectorSpaceBlocked\' vector space!"
00332       );
00333   }
00334 #endif
00335   //
00336   // The first thing that we do is to see if all of the vectors involved are
00337   // incore, serial vectors.  In this case, the vectors should be compatible.
00338   // To accomplish this we will pick the continguous vector to implement
00339   // the applyOp(...) function.
00340   //
00341   // Get the index of an incore-only input vector and input/output vector?
00342   const bool this_isInCore = productSpace_->hasInCoreView(
00343     Range1D(),VIEW_TYPE_DETACHED,STRIDE_TYPE_NONUNIT
00344     );
00345   int incore_vec_k = -1, incore_targ_vec_k = -1;
00346   // Dynamic cast the pointers for the vector arguments
00347   Workspace<const ProductVectorBase<Scalar>*>
00348     vecs_args(wss,num_vecs,false);
00349   for(int k = 0; k < num_vecs; ++k) {
00350     vecs_args[k] = dynamic_cast<const ProductVectorBase<Scalar>*>(vecs[k]);
00351     if( vecs_args[k] == NULL ) {
00352       const bool isInCore_k = vecs[k]->space()->hasInCoreView();
00353       if( this_isInCore && isInCore_k ) {
00354         incore_vec_k = k;
00355         break;
00356       }
00357       TEST_FOR_EXCEPTION(
00358         !this_isInCore || (this_isInCore && !isInCore_k), Exceptions::IncompatibleVectorSpaces
00359         ,"DefaultProductVector::applyOp(...): Error vecs["<<k<<"] "
00360         <<"of type \'"<<typeName(*vecs[k])<<"\' does not support the "
00361         <<"\'DefaultProductVector<Scalar>\' interface and is not an incore vector or this is not an incore vector!"
00362         );
00363     }
00364   }
00365   Workspace<ProductVectorBase<Scalar>*>
00366     targ_vecs_args(wss,num_targ_vecs,false);
00367   for(int k = 0; k < num_targ_vecs; ++k) {
00368     targ_vecs_args[k] = dynamic_cast<ProductVectorBase<Scalar>*>(targ_vecs[k]);
00369     if( targ_vecs_args[k] == NULL ) {
00370       const bool isInCore_k = targ_vecs[k]->space()->hasInCoreView();
00371       if( this_isInCore && isInCore_k ) {
00372         incore_targ_vec_k = k;
00373         break;
00374       }
00375       TEST_FOR_EXCEPTION(
00376         !this_isInCore || (this_isInCore && !isInCore_k), Exceptions::IncompatibleVectorSpaces
00377         ,"DefaultProductVector::applyOp(...): Error targ_vecs["<<k<<"] "
00378         <<"of type \'"<<typeName(*targ_vecs[k])<<"\' does not support the "
00379         <<"\'DefaultProductVector<Scalar>\' interface and is not an incore vector or this is not an incore vector!"
00380         );
00381     }
00382   }
00383   // Let a incore-only vector with a contiguous view handle this through
00384   // explicit vector access?
00385   if( incore_vec_k >= 0 ) {
00386     vecs[incore_vec_k]->applyOp(
00387       op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj
00388       ,first_ele_offset_in,sub_dim_in,global_offset_in
00389       );
00390     return;
00391   }
00392   else if ( incore_targ_vec_k >= 0 ) {
00393     targ_vecs[incore_targ_vec_k]->applyOp(
00394       op,num_vecs,vecs,num_targ_vecs,targ_vecs,reduct_obj
00395       ,first_ele_offset_in,sub_dim_in,global_offset_in
00396       );
00397     return;
00398   }
00399   //
00400   // If we get here, then we will implement the applyOp(...) one vector block
00401   // at a time.
00402   //
00403   const Index this_dim = n;
00404   const Index sub_dim
00405     = (
00406       sub_dim_in < 0
00407       ? this_dim - first_ele_offset_in
00408       : sub_dim_in
00409       );
00410   Index num_elements_remaining = sub_dim;
00411   const int numBlocks = productSpace_->numBlocks();
00412   Workspace<RCP<const VectorBase<Scalar> > >
00413     sub_vecs_rcps(wss,num_vecs);
00414   Workspace<const VectorBase<Scalar>*>
00415     sub_vecs(wss,num_vecs,false);
00416   Workspace<RCP<VectorBase<Scalar> > >
00417     sub_targ_vecs_rcps(wss,num_targ_vecs);
00418   Workspace<VectorBase<Scalar>*>
00419     sub_targ_vecs(wss,num_targ_vecs,false);
00420   Index g_off = -first_ele_offset_in;
00421   for(int k = 0; k < numBlocks; ++k) {
00422     const Index local_dim = productSpace_->getBlock(k)->dim();
00423     if( g_off < 0 && -g_off+1 > local_dim ) {
00424       g_off += local_dim;
00425       continue;
00426     }
00427     const Index
00428       local_sub_dim
00429       = ( g_off >= 0
00430           ? std::min( local_dim, num_elements_remaining )
00431           : std::min( local_dim + g_off, num_elements_remaining ) );
00432     if( local_sub_dim <= 0 )
00433       break;
00434     // Fill constituent vectors for block k
00435     for( int i = 0; i < num_vecs; ++i ) {
00436       sub_vecs_rcps[i] = vecs_args[i]->getVectorBlock(k);
00437       sub_vecs[i] = &*sub_vecs_rcps[i];
00438     }
00439     // Fill constituent target vectors for block k
00440     for( int j = 0; j < num_targ_vecs; ++j ) {
00441       sub_targ_vecs_rcps[j] = targ_vecs_args[j]->getNonconstVectorBlock(k);
00442       sub_targ_vecs[j] = &*sub_targ_vecs_rcps[j];
00443     }
00444     Thyra::applyOp( 
00445       op
00446       ,num_vecs,num_vecs?&sub_vecs[0]:NULL
00447       ,num_targ_vecs,num_targ_vecs?&sub_targ_vecs[0]:NULL
00448       ,reduct_obj
00449       ,g_off < 0 ? -g_off : 0                                    // first_ele_offset
00450       ,local_sub_dim                                             // sub_dim
00451       ,g_off < 0 ? global_offset_in : global_offset_in + g_off   // global_offset
00452       );
00453     g_off += local_dim;
00454     num_elements_remaining -= local_sub_dim;
00455   }
00456   TEST_FOR_EXCEPT(!(num_elements_remaining==0));
00457 }
00458 
00459 
00460 template <class Scalar>
00461 void DefaultProductVector<Scalar>::acquireDetachedVectorViewImpl(
00462   const Range1D& rng_in, RTOpPack::ConstSubVectorView<Scalar>* sub_vec
00463   ) const
00464 {
00465   const Range1D
00466     rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
00467   int    kth_vector_space  = -1;
00468   Index  kth_global_offset = 0;
00469   productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00470 #ifdef TEUCHOS_DEBUG
00471   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00472 #endif
00473   if(
00474     rng.lbound() + rng.size()
00475     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00476     )
00477   {
00478     // This involves only one sub-vector so just return it.
00479     const_cast<const VectorBase<Scalar>*>(
00480       &*vecs_[kth_vector_space].getConstObj()
00481       )->acquireDetachedView( rng - kth_global_offset, sub_vec );
00482     sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
00483   }
00484   else {
00485     // Just let the default implementation handle this.  ToDo: In the future
00486     // we could manually construct an explicit sub-vector that spanned
00487     // two or more constituent vectors but this would be a lot of work.
00488     // However, this would require the use of temporary memory but
00489     // so what.
00490     VectorDefaultBase<Scalar>::acquireDetachedVectorViewImpl(rng_in,sub_vec);
00491   }
00492 }
00493 
00494 
00495 template <class Scalar>
00496 void DefaultProductVector<Scalar>::releaseDetachedVectorViewImpl(
00497   RTOpPack::ConstSubVectorView<Scalar>* sub_vec
00498   ) const
00499 {
00500   if( sub_vec->values() == NULL ) return;
00501   int    kth_vector_space  = -1;
00502   Index  kth_global_offset = 0;
00503   productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
00504 #ifdef TEUCHOS_DEBUG
00505   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00506 #endif
00507   if(
00508     sub_vec->globalOffset() + sub_vec->subDim()
00509     <= kth_global_offset +  vecs_[kth_vector_space].getConstObj()->space()->dim()
00510     )
00511   {
00512     // This sub_vec was extracted from a single constituent vector
00513     sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
00514     vecs_[kth_vector_space].getConstObj()->releaseDetachedView(sub_vec);
00515   }
00516   else {
00517     // This sub_vec was created by the default implementation!
00518     VectorDefaultBase<Scalar>::releaseDetachedVectorViewImpl(sub_vec);
00519   }
00520 }
00521 
00522 
00523 template <class Scalar>
00524 void DefaultProductVector<Scalar>::acquireNonconstDetachedVectorViewImpl(
00525   const Range1D& rng_in, RTOpPack::SubVectorView<Scalar>* sub_vec
00526   )
00527 {
00528   const Range1D
00529     rng = rng_in.full_range() ? Range1D(0,productSpace_->dim()-1) : rng_in;
00530   int    kth_vector_space  = -1;
00531   Index  kth_global_offset = 0;
00532   productSpace_->getVecSpcPoss(rng.lbound(),&kth_vector_space,&kth_global_offset);
00533 #ifdef TEUCHOS_DEBUG
00534   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00535 #endif
00536   if(
00537     rng.lbound() + rng.size()
00538     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00539     )
00540   {
00541     // This involves only one sub-vector so just return it.
00542     vecs_[kth_vector_space].getConstObj()->acquireDetachedView(
00543       rng - kth_global_offset, sub_vec
00544       );
00545     sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
00546   }
00547   else {
00548     // Just let the default implementation handle this.  ToDo: In the future
00549     // we could manually construct an explicit sub-vector that spanned
00550     // two or more constituent vectors but this would be a lot of work.
00551     // However, this would require the use of temporary memory but
00552     // so what.
00553     VectorDefaultBase<Scalar>::acquireNonconstDetachedVectorViewImpl(rng_in,sub_vec);
00554   }
00555 }
00556 
00557 
00558 template <class Scalar>
00559 void DefaultProductVector<Scalar>::commitNonconstDetachedVectorViewImpl(
00560   RTOpPack::SubVectorView<Scalar>* sub_vec
00561   )
00562 {
00563   if( sub_vec->values() == NULL ) return;
00564   int    kth_vector_space  = -1;
00565   Index  kth_global_offset = 0;
00566   productSpace_->getVecSpcPoss(sub_vec->globalOffset(),&kth_vector_space,&kth_global_offset);
00567 #ifdef TEUCHOS_DEBUG
00568   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00569 #endif
00570   if(
00571     sub_vec->globalOffset() + sub_vec->subDim()
00572     <= kth_global_offset +  vecs_[kth_vector_space].getConstObj()->space()->dim()
00573     )
00574   {
00575     // This sub_vec was extracted from a single constituent vector
00576     sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
00577     vecs_[kth_vector_space].getNonconstObj()->commitDetachedView(sub_vec);
00578   }
00579   else {
00580     // This sub_vec was created by the default implementation!
00581     VectorDefaultBase<Scalar>::commitNonconstDetachedVectorViewImpl(sub_vec);
00582   }
00583 }
00584 
00585 
00586 template <class Scalar>
00587 void DefaultProductVector<Scalar>::setSubVector(
00588   const RTOpPack::SparseSubVectorT<Scalar>& sub_vec
00589   )
00590 {
00591   int    kth_vector_space  = -1;
00592   Index  kth_global_offset = 0;
00593   productSpace_->getVecSpcPoss(sub_vec.globalOffset(),&kth_vector_space,&kth_global_offset);
00594 #ifdef TEUCHOS_DEBUG
00595   TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= numBlocks_ ) );
00596 #endif
00597   if(
00598     sub_vec.globalOffset() + sub_vec.subDim()
00599     <= kth_global_offset + vecs_[kth_vector_space].getConstObj()->space()->dim()
00600     )
00601   {
00602     // This sub-vector fits into a single constituent vector
00603     RTOpPack::SparseSubVectorT<Scalar> sub_vec_g = sub_vec;
00604     sub_vec_g.setGlobalOffset( sub_vec_g.globalOffset() - kth_global_offset );
00605     vecs_[kth_vector_space].getNonconstObj()->setSubVector(sub_vec_g);
00606   }
00607   else {
00608     // Let the default implementation take care of this.  ToDo: In the future
00609     // it would be possible to manually set the relevant constituent
00610     // vectors with no temp memory allocations.
00611     VectorDefaultBase<Scalar>::setSubVector(sub_vec);
00612   }
00613 }
00614 
00615 
00616 } // namespace Thyra
00617 
00618 
00619 #endif // THYRA_PRODUCT_VECTOR_HPP

Generated on Tue Oct 20 12:47:25 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7