MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_VectorMutableBlocked.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) 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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <typeinfo>
00043 #include <algorithm>
00044 
00045 #include "AbstractLinAlgPack_VectorMutableBlocked.hpp"
00046 #include "AbstractLinAlgPack_VectorMutableSubView.hpp"
00047 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00048 #include "Teuchos_Workspace.hpp"
00049 #include "Teuchos_Assert.hpp"
00050 #include "Teuchos_FancyOStream.hpp"
00051 
00052 // Uncomment to ignore cache of reduction data
00053 //#define ALAP_VECTOR_MUTABLE_BLOCKED_IGNORE_CACHE_DATA
00054 
00055 namespace {
00056 template< class T >
00057 inline
00058 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00059 template< class T >
00060 inline
00061 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00062 } // end namespace
00063 
00064 namespace AbstractLinAlgPack {
00065 
00066 VectorMutableBlocked::VectorMutableBlocked(
00067   VectorMutable::vec_mut_ptr_t*  vecs
00068   ,const vec_space_comp_ptr_t&         vec_space
00069   )
00070   // The members will be initialized in this->has_changed()
00071 {
00072   initialize(vecs,vec_space);
00073 }
00074 
00075 void VectorMutableBlocked::initialize(
00076   VectorMutable::vec_mut_ptr_t*  vecs
00077   ,const vec_space_comp_ptr_t&         vec_space
00078   )
00079 {
00080   TEUCHOS_TEST_FOR_EXCEPTION(
00081     vec_space.get() == NULL, std::logic_error
00082     ,"VectorMutableBlocked::initialize(...): Error, Must be constructed from "
00083     "a non-null block vector space!" );
00084   const int num_vec_spaces = vec_space->num_vector_spaces();
00085   vecs_.resize(num_vec_spaces);
00086   std::copy(vecs,vecs+num_vec_spaces,vecs_.begin());
00087   vec_space_ = vec_space;
00088   this->has_changed();
00089 }
00090   
00091 // overriddend form Vector
00092 
00093 index_type VectorMutableBlocked::dim() const
00094 {
00095   return vec_space_->dim();
00096 }
00097 
00098 const VectorSpace& VectorMutableBlocked::space() const
00099 {
00100   return *vec_space_;
00101 }
00102 
00103 void VectorMutableBlocked::apply_op(
00104   const RTOpPack::RTOp& op
00105   ,const size_t num_vecs, const Vector* vecs[]
00106   ,const size_t num_targ_vecs, VectorMutable* targ_vecs[]
00107   ,RTOpPack::ReductTarget *reduct_obj
00108   ,const index_type first_ele_in, const index_type sub_dim_in, const index_type global_offset_in
00109   ) const
00110 {
00111   using Teuchos::Workspace;
00112   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00113   
00114   const index_type
00115     n = this->dim();
00116 
00117   // Validate the compatibility of the vectors!
00118 #ifdef TEUCHOS_DEBUG
00119   TEUCHOS_TEST_FOR_EXCEPTION(
00120     !(1 <= first_ele_in && first_ele_in <= n), std::out_of_range
00121     ,"VectorMutableBlocked::apply_op(...): Error, "
00122     "first_ele_in = " << first_ele_in << " is not in range [1,"<<n<<"]" );
00123   TEUCHOS_TEST_FOR_EXCEPTION(
00124     sub_dim_in < 0, std::invalid_argument
00125     ,"VectorMutableBlocked::apply_op(...): Error, "
00126     "sub_dim_in = " << sub_dim_in << " is not acceptable" );
00127   TEUCHOS_TEST_FOR_EXCEPTION(
00128     global_offset_in < 0, std::invalid_argument
00129     ,"VectorMutableBlocked::apply_op(...): Error, "
00130     "global_offset_in = " << global_offset_in << " is not acceptable" );
00131   TEUCHOS_TEST_FOR_EXCEPTION(
00132     sub_dim_in > 0 && sub_dim_in - (first_ele_in - 1) > n, std::length_error
00133     ,"VectorMutableBlocked::apply_op(...): Error, "
00134     "global_offset_in = " << global_offset_in << ", sub_dim_in = " << sub_dim_in
00135     << "first_ele_in = " << first_ele_in << " and n = " << n
00136     << " are not compatible" );
00137   bool test_failed;
00138   {for(int k = 0; k < num_vecs; ++k) {
00139     test_failed = !this->space().is_compatible(vecs[k]->space());
00140     TEUCHOS_TEST_FOR_EXCEPTION(
00141       test_failed, VectorSpace::IncompatibleVectorSpaces
00142       ,"VectorMutableBlocked::apply_op(...): Error vecs["<<k<<"]->space() "
00143       <<"of type \'"<<typeName(vecs[k]->space())<<"\' is not compatible with this "
00144       <<"\'VectorSpaceBlocked\' vector space!"
00145       );
00146   }}
00147   {for(int k = 0; k < num_targ_vecs; ++k) {
00148     test_failed = !this->space().is_compatible(targ_vecs[k]->space());
00149     TEUCHOS_TEST_FOR_EXCEPTION(
00150       test_failed, VectorSpace::IncompatibleVectorSpaces
00151       ,"VectorMutableBlocked::apply_op(...): Error targ_vecs["<<k<<"]->space() "
00152       <<"of type \'"<<typeName(vecs[k]->space())<<"\' is not compatible with this "
00153       <<"\'VectorSpaceBlocked\' vector space!"
00154       );
00155   }}
00156 #endif
00157   
00158   // Dynamic cast the pointers for the vector arguments
00159   Workspace<const VectorMutableBlocked*>
00160     vecs_args(wss,num_vecs);
00161   {for(int k = 0; k < num_vecs; ++k) {
00162     vecs_args[k] = dynamic_cast<const VectorMutableBlocked*>(vecs[k]);
00163 #ifdef TEUCHOS_DEBUG
00164     TEUCHOS_TEST_FOR_EXCEPTION(
00165       vecs_args[k] == NULL, VectorSpace::IncompatibleVectorSpaces
00166       ,"VectorMutableBlocked::apply_op(...): Error vecs["<<k<<"] "
00167       <<"of type \'"<<typeName(*vecs[k])<<"\' does not support the "
00168       <<"\'VectorMutableBlocked\' interface!"
00169       );
00170 #endif
00171   }}
00172   Workspace<VectorMutableBlocked*>
00173     targ_vecs_args(wss,num_targ_vecs);
00174   {for(int k = 0; k < num_targ_vecs; ++k) {
00175     targ_vecs_args[k] = dynamic_cast<VectorMutableBlocked*>(targ_vecs[k]);
00176 #ifdef TEUCHOS_DEBUG
00177     TEUCHOS_TEST_FOR_EXCEPTION(
00178       targ_vecs_args[k] == NULL, VectorSpace::IncompatibleVectorSpaces
00179       ,"VectorMutableBlocked::apply_op(...): Error targ_vecs["<<k<<"] "
00180       <<"of type \'"<<typeName(*targ_vecs[k])<<"\' does not support the "
00181       <<"\'VectorMutableBlocked\' interface!"
00182       );
00183 #endif
00184   }}
00185 
00186   // Perform the reduction on each vector segment at a time.
00187   // ToDo: There could be a parallel version of this method!
00188   const index_type this_dim = this->dim();
00189   const index_type sub_dim  = ( sub_dim_in == 0
00190                   ? this_dim - (first_ele_in - 1)
00191                   : sub_dim_in );
00192   index_type num_elements_remaining = sub_dim;
00193   const int  num_vec_spaces = vec_space_->num_vector_spaces();
00194   Workspace<const Vector*>
00195     sub_vecs(wss,num_vecs);
00196   Workspace<VectorMutable*>
00197     sub_targ_vecs(wss,num_targ_vecs);
00198   index_type g_off = -first_ele_in + 1;
00199   for(int k = 0; k < num_vec_spaces; ++k) {
00200     const index_type local_dim = vecs_[k]->dim();
00201     if( g_off < 0 && -g_off > local_dim ) {
00202       g_off += local_dim;
00203       continue;
00204     }
00205     const index_type
00206       local_sub_dim = ( g_off >= 0
00207                 ? my_min( local_dim, num_elements_remaining )
00208                 : my_min( local_dim + g_off, num_elements_remaining ) );
00209     if( local_sub_dim <= 0 )
00210       break;
00211     for( int i = 0; i < num_vecs; ++i ) // Fill constituent vectors for segment k
00212       sub_vecs[i] = vecs_args[i]->vecs_[k].get();
00213     for( int j = 0; j < num_targ_vecs; ++j ) // Fill constituent target vectors for segment k
00214       sub_targ_vecs[j] = targ_vecs_args[j]->vecs_[k].get();
00215     AbstractLinAlgPack::apply_op( 
00216       op
00217       ,num_vecs,num_vecs?&sub_vecs[0]:NULL
00218       ,num_targ_vecs,num_targ_vecs?&sub_targ_vecs[0]:NULL
00219       ,reduct_obj
00220       ,g_off < 0 ? -g_off + 1 : 1                                // first_ele
00221       ,local_sub_dim                                             // sub_dim
00222       ,g_off < 0 ? global_offset_in : global_offset_in + g_off   // global_offset
00223       );
00224     g_off += local_dim;
00225     num_elements_remaining -= local_sub_dim;
00226   }
00227   TEUCHOS_TEST_FOR_EXCEPT( !( num_elements_remaining == 0 ) );
00228   // Must allert all of the block vectors that they may have changed!
00229   {for(int k = 0; k < num_targ_vecs; ++k) {
00230     targ_vecs[k]->has_changed();
00231   }}
00232 }
00233 
00234 index_type VectorMutableBlocked::nz() const
00235 {
00236   if( nz_ < 0.0 ) {
00237     const int num_vec_spaces = vec_space_->num_vector_spaces();
00238     nz_ = 0;
00239     for( int k = 0; k < num_vec_spaces; ++k )
00240       nz_ += vecs_[k]->nz();
00241   }
00242   return nz_;
00243 }
00244 
00245 std::ostream& VectorMutableBlocked::output(
00246   std::ostream& out_arg, bool print_dim, bool newline
00247   ,index_type global_offset
00248   ) const
00249 {
00250   Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00251   Teuchos::OSTab tab(out);
00252   if(print_dim)
00253     *out << this->dim() << std::endl;
00254   size_type off = global_offset;
00255   const int num_vec_spaces = vec_space_->num_vector_spaces();
00256   for( int k = 0; k < num_vec_spaces; ++k ) {
00257     vecs_[k]->output(*out,false,false,off);
00258     off += vecs_[k]->dim();
00259   }
00260   if(newline)
00261     *out << std::endl;
00262   return out_arg;
00263 }
00264 
00265 value_type VectorMutableBlocked::get_ele(index_type i) const
00266 {
00267   int         kth_vector_space  = -1;
00268   index_type  kth_global_offset = 0;
00269   vec_space_->get_vector_space_position(i,&kth_vector_space,&kth_global_offset);
00270 #ifdef TEUCHOS_DEBUG
00271   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs_.size()  ) );
00272 #endif
00273   return vecs_[kth_vector_space]->get_ele( i - kth_global_offset );
00274 }
00275 
00276 value_type VectorMutableBlocked::norm_1() const
00277 {
00278 #ifdef ALAP_VECTOR_MUTABLE_BLOCKED_IGNORE_CACHE_DATA
00279   if( true ) {
00280 #else
00281   if( norm_1_ < 0.0 ) {
00282 #endif
00283     const int num_vec_spaces = vec_space_->num_vector_spaces();
00284     norm_1_ = 0.0;
00285     for( int k = 0; k < num_vec_spaces; ++k )
00286       norm_1_ += vecs_[k]->norm_1();
00287   }
00288   return norm_1_;
00289 }
00290 
00291 value_type VectorMutableBlocked::norm_inf() const
00292 {
00293 #ifdef ALAP_VECTOR_MUTABLE_BLOCKED_IGNORE_CACHE_DATA
00294   if( true ) {
00295 #else
00296   if( norm_inf_ < 0.0 ) {
00297 #endif
00298     const int num_vec_spaces = vec_space_->num_vector_spaces();
00299     norm_inf_ = 0.0;
00300     for( int k = 0; k < num_vec_spaces; ++k )
00301       norm_inf_  = my_max( norm_inf_, vecs_[k]->norm_inf() );
00302   }
00303   return norm_inf_;
00304 }
00305 
00306 value_type VectorMutableBlocked::inner_product( const Vector& v ) const
00307 {
00308   return Vector::inner_product(v); // ToDo: Specialize? (why bother?)
00309 }
00310 
00311 void VectorMutableBlocked::get_sub_vector( const Range1D& rng_in, RTOpPack::SubVector* sub_vec ) const
00312 {
00313   const Range1D
00314     rng = rng_in.full_range() ? Range1D( 1, this->dim()) : rng_in;
00315   int         kth_vector_space  = -1;
00316   index_type  kth_global_offset = 0;
00317   vec_space_->get_vector_space_position(rng.lbound(),&kth_vector_space,&kth_global_offset);
00318 #ifdef TEUCHOS_DEBUG
00319   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs_.size()  ) );
00320 #endif
00321   if( rng.lbound() + rng.size() <= kth_global_offset + 1 + vecs_[kth_vector_space]->dim() ) {
00322     // This involves only one sub-vector so just return it.
00323     static_cast<const Vector*>(vecs_[kth_vector_space].get())->get_sub_vector(
00324       rng - kth_global_offset, sub_vec );
00325     sub_vec->setGlobalOffset( sub_vec->globalOffset() + kth_global_offset );
00326   }
00327   else {
00328     // Just let the default implementation handle this.  ToDo: In the futrue
00329     // we could manually construct an explicit sub-vector that spanned
00330     // two or more consitituent vectors but this would be a lot of work.
00331     // However, this would require the use of temporary memory but
00332     // so what.
00333     Vector::get_sub_vector(rng_in,sub_vec);
00334   }
00335 }
00336 
00337 void VectorMutableBlocked::free_sub_vector( RTOpPack::SubVector* sub_vec ) const
00338 {
00339   if( sub_vec->values() == NULL ) return;
00340   int         kth_vector_space  = -1;
00341   index_type  kth_global_offset = 0;
00342   vec_space_->get_vector_space_position(
00343     sub_vec->globalOffset()+1,&kth_vector_space,&kth_global_offset);
00344 #ifdef TEUCHOS_DEBUG
00345   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs_.size()  ) );
00346 #endif
00347   if( sub_vec->globalOffset() + sub_vec->subDim() <= kth_global_offset +  vecs_[kth_vector_space]->dim() ) {
00348     // This sub_vec was extracted from a single constituent vector
00349     sub_vec->setGlobalOffset( sub_vec->globalOffset() - kth_global_offset );
00350     vecs_[kth_vector_space].get()->free_sub_vector(sub_vec);
00351   }
00352   else {
00353     // This sub_vec was created by the default implementation!
00354     Vector::free_sub_vector(sub_vec);
00355   }
00356 }
00357 
00358 void VectorMutableBlocked::has_changed() const
00359 {
00360   nz_ = -1;                   // set to uninitalized!
00361   norm_1_ = norm_inf_ = -1.0;
00362   Vector::has_changed();
00363 }
00364 
00365 // overridden from VectorMutable
00366 
00367 VectorMutable::vec_mut_ptr_t
00368 VectorMutableBlocked::sub_view( const Range1D& rng_in )
00369 {
00370   namespace rcp = MemMngPack;
00371   const index_type dim = this->dim();
00372   const Range1D    rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
00373   // Validate the preconditions
00374 #ifdef TEUCHOS_DEBUG
00375   TEUCHOS_TEST_FOR_EXCEPTION(
00376     dim < rng.ubound(), std::out_of_range
00377     ,"VectorMutableBlocked::sub_view(...): Error, rng = "
00378     << "["<<rng.lbound()<<","<<rng.ubound()<<"] is not in range [1,"<<dim<<"]" );
00379 #endif
00380   vecs_t &vecs = this->vecs_; // Need to examine in debugger!
00381   // See if the entire vector is being returned or just NULL
00382   if( rng.lbound() == 1 && rng.ubound() == dim ) {
00383     if( vecs.size() == 1 )   return vecs[0]->sub_view(Range1D());
00384     else                     return Teuchos::rcp(this,false);
00385   }
00386   // From here on out we will return a view that could change the
00387   // elements of one or more of the constituent vectors so we had
00388   // better wipe out the cache
00389   this->has_changed();
00390   // Get the position of the vector space object of interest
00391   int           kth_vector_space  = -1;
00392   index_type    kth_global_offset = 0;
00393   vec_space_->get_vector_space_position(rng.lbound(),&kth_vector_space,&kth_global_offset);
00394   const VectorSpace::space_ptr_t*  vector_spaces      = vec_space_->vector_spaces();
00395   const index_type*                vec_spaces_offsets = vec_space_->vector_spaces_offsets();
00396 #ifdef TEUCHOS_DEBUG
00397   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs.size()  ) );
00398 #endif
00399   if( rng.lbound() == kth_global_offset + 1
00400     && rng.size() == vec_spaces_offsets[kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
00401     // The client selected a whole single constituent vector.
00402     return vecs[kth_vector_space]->sub_view(Range1D());
00403   if( rng.ubound() <= vec_spaces_offsets[kth_vector_space+1] )
00404     // The client selected a sub-vector of a single consituent vector
00405     return vecs[kth_vector_space]->sub_view( rng - vec_spaces_offsets[kth_vector_space] );
00406   // The client selected a sub-vector that spans two or more constituent vectors.
00407   // Get the position of the vector space object with the last element of interest
00408   int           end_kth_vector_space  = -1;
00409   index_type    end_kth_global_offset = 0;
00410   vec_space_->get_vector_space_position(rng.ubound(),&end_kth_vector_space,&end_kth_global_offset);
00411 #ifdef TEUCHOS_DEBUG
00412   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= end_kth_vector_space && end_kth_vector_space <= vecs.size()  ) );
00413   TEUCHOS_TEST_FOR_EXCEPT( !(  end_kth_vector_space > kth_vector_space  ) );
00414 #endif
00415   // Create a VectorWithOpMutableCompsiteStd object containing the relavant constituent vectors
00416   Teuchos::RCP<VectorMutableBlocked>
00417     vec_comp = Teuchos::rcp(new VectorMutableBlocked(
00418       &vecs[kth_vector_space]
00419       ,Teuchos::rcp(new VectorSpaceBlocked(
00420         &vector_spaces[kth_vector_space]
00421         ,end_kth_vector_space - kth_vector_space + 1 ))
00422       ));
00423   if( rng.lbound() == kth_global_offset + 1
00424     && rng.size() == vec_spaces_offsets[end_kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
00425     // The client selected exactly a contigous set of vectors
00426     return vec_comp;
00427   // The client selected some sub-set of elements in the contigous set of vectors
00428   return Teuchos::rcp(
00429     new VectorMutableSubView(
00430       vec_comp
00431       ,Range1D( 
00432         rng.lbound()-vec_spaces_offsets[kth_vector_space]
00433         ,rng.ubound()-vec_spaces_offsets[kth_vector_space] )
00434       ) );
00435 }
00436 
00437 void VectorMutableBlocked::axpy( value_type alpha, const Vector& x )
00438 {
00439   VectorMutable::axpy(alpha,x); // ToDo: Specialize? (why bother?)
00440   this->has_changed();
00441 }
00442 
00443 VectorMutable& VectorMutableBlocked::operator=(value_type alpha)
00444 {
00445   const int num_vec_spaces = vec_space_->num_vector_spaces();
00446   for( int k = 0; k < num_vec_spaces; ++k )
00447     *vecs_[k] = alpha;
00448   this->has_changed();
00449   return *this;
00450 }
00451 
00452 VectorMutable& VectorMutableBlocked::operator=(const Vector& v)
00453 {
00454   VectorMutable::operator=(v); // ToDo: Specialize (why bother?)
00455   this->has_changed();
00456   return *this;
00457 }
00458 
00459 void VectorMutableBlocked::set_ele( index_type i, value_type val )
00460 {
00461   int         kth_vector_space  = -1;
00462   index_type  kth_global_offset = 0;
00463   vec_space_->get_vector_space_position(i,&kth_vector_space,&kth_global_offset);
00464 #ifdef TEUCHOS_DEBUG
00465   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs_.size()  ) );
00466 #endif
00467   vecs_[kth_vector_space]->set_ele( i - kth_global_offset, val );
00468   this->has_changed();
00469 }
00470 
00471 void VectorMutableBlocked::set_sub_vector( const RTOpPack::SparseSubVector& sub_vec )
00472 {
00473   int         kth_vector_space  = -1;
00474   index_type  kth_global_offset = 0;
00475   vec_space_->get_vector_space_position(
00476     sub_vec.globalOffset()+1,&kth_vector_space,&kth_global_offset);
00477 #ifdef TEUCHOS_DEBUG
00478   TEUCHOS_TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vecs_.size()  ) );
00479 #endif
00480   if( sub_vec.globalOffset() + sub_vec.subDim() <= kth_global_offset +  vecs_[kth_vector_space]->dim() ) {
00481     // This sub-vector fits into a single constituent vector
00482     RTOpPack::SparseSubVector sub_vec_g = sub_vec;
00483     sub_vec_g.setGlobalOffset( sub_vec_g.globalOffset() - kth_global_offset );
00484     vecs_[kth_vector_space]->set_sub_vector(sub_vec_g);
00485   }
00486   else {
00487     // Let the default implementation take care of this.  ToDo: In the futrue
00488     // it would be possible to manualy set the relavent constituent
00489     // vectors with no temp memory allocations.
00490     VectorMutable::set_sub_vector(sub_vec);
00491   }
00492   this->has_changed();
00493 }
00494 
00495 void VectorMutableBlocked::assert_in_range(int k) const
00496 {
00497   assert_initialized();
00498   TEUCHOS_TEST_FOR_EXCEPTION(
00499     k >= vec_space_->num_vector_spaces(), std::logic_error
00500     ,"VectorMutableBlocked::assert_in_range(k) : Error, k = " << k << " is not "
00501     "in range [0," << vec_space_->num_vector_spaces() << "]" );
00502 }
00503 
00504 void VectorMutableBlocked::assert_initialized() const
00505 {
00506   TEUCHOS_TEST_FOR_EXCEPTION(!vec_space_.get(),std::logic_error,"Error, this is not initalized!");
00507 }
00508 
00509 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines