AbstractLinAlgPack_VectorSpaceBlocked.cpp

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 // 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 Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include <algorithm>
00032 
00033 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00034 #include "AbstractLinAlgPack_VectorMutableBlocked.hpp"
00035 #include "AbstractLinAlgPack_VectorSpaceSubSpace.hpp"
00036 #include "AbstractLinAlgPack_VectorSpaceFactory.hpp"
00037 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00038 #include "Teuchos_Workspace.hpp"
00039 #include "Teuchos_TestForException.hpp"
00040 
00041 namespace AbstractLinAlgPack {
00042 
00043 VectorSpaceBlocked::VectorSpaceBlocked(
00044   const VectorSpace::space_ptr_t         vec_spaces[]
00045   ,int                                   num_vec_spaces
00046   ,const VectorSpace::space_fcty_ptr_t   &small_vec_spc_fcty
00047   )
00048 {
00049   initialize(vec_spaces,num_vec_spaces,small_vec_spc_fcty);
00050 }
00051   
00052 void VectorSpaceBlocked::initialize(
00053   const VectorSpace::space_ptr_t         vec_spaces[]
00054   ,int                                   num_vec_spaces
00055   ,const VectorSpace::space_fcty_ptr_t   &small_vec_spc_fcty
00056   )
00057 {
00058   vector_spaces_.resize(num_vec_spaces);
00059   std::copy(vec_spaces,vec_spaces+num_vec_spaces,vector_spaces_.begin());
00060   vec_spaces_offsets_.resize(num_vec_spaces+1);
00061   vec_spaces_offsets_[0] = 0;
00062   for( int k = 1; k <= num_vec_spaces; ++k )
00063     vec_spaces_offsets_[k] = vec_spaces_offsets_[k-1] + vec_spaces[k-1]->dim();
00064   small_vec_spc_fcty_ = small_vec_spc_fcty;
00065 }
00066 
00067 void VectorSpaceBlocked::get_vector_space_position(
00068   index_type i, int* kth_vector_space, index_type* kth_global_offset ) const
00069 {
00070   // Validate the preconditions
00071 #ifdef TEUCHOS_DEBUG
00072   TEST_FOR_EXCEPTION(
00073     i < 1 || this->dim() < i, std::out_of_range
00074     ,"VectorSpaceBlocked::get_vector_space_position(...): Error, i = "
00075     << i << " is not in range [1,"<<this->dim()<<"]"
00076     );
00077 #endif
00078   *kth_vector_space  = 0;
00079   *kth_global_offset = 0;
00080   while( *kth_vector_space < vector_spaces_.size() ) {
00081     const RTOp_index_type off_kp1 = vec_spaces_offsets_[*kth_vector_space+1];
00082     if( off_kp1 + 1 > i ) {
00083       *kth_global_offset = vec_spaces_offsets_[*kth_vector_space];
00084       break;
00085     }
00086     ++(*kth_vector_space);
00087   }
00088 #ifdef TEUCHOS_DEBUG
00089   TEST_FOR_EXCEPT( !( *kth_vector_space < vector_spaces_.size() ) );
00090 #endif
00091 }
00092 
00093 // overridden from VectorSpace
00094 
00095 bool VectorSpaceBlocked::is_compatible(const VectorSpace& vec_space) const
00096 {
00097   const VectorSpaceBlocked
00098     *vec_space_comp = dynamic_cast<const VectorSpaceBlocked*>(&vec_space);
00099   if( !vec_space_comp || vec_space_comp->vector_spaces_.size() != this->vector_spaces_.size() )
00100     return false;
00101   // There is some hope that these are compatible vector spaces.
00102   for( int k = 0; k < vector_spaces_.size(); ++k ) {
00103     if( !vec_space_comp->vector_spaces_[k]->is_compatible(*this->vector_spaces_[k]) )
00104       return false;
00105   }
00106   return true; // If we get here we are compatible!
00107 }
00108   
00109 index_type VectorSpaceBlocked::dim() const
00110 {
00111   return vec_spaces_offsets_[vector_spaces_.size()];
00112 }
00113 
00114 VectorSpace::space_fcty_ptr_t
00115 VectorSpaceBlocked::small_vec_spc_fcty() const
00116 {
00117   return small_vec_spc_fcty_;
00118 }
00119 
00120 VectorSpace::vec_mut_ptr_t
00121 VectorSpaceBlocked::create_member() const
00122 {
00123   namespace rcp = MemMngPack;
00124   using Teuchos::Workspace;
00125   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00126 
00127   const int num_vec_spaces = this->num_vector_spaces();
00128 
00129   // Create the vector objects array.
00130   Workspace<VectorMutable::vec_mut_ptr_t>
00131     vecs(wss,num_vec_spaces);
00132   for( int k = 0; k < num_vec_spaces; ++k )
00133     vecs[k] = vector_spaces_[k]->create_member();
00134 
00135   return Teuchos::rcp(
00136     new VectorMutableBlocked(
00137       &vecs[0]
00138       ,Teuchos::rcp(new VectorSpaceBlocked(*this))
00139       ) );
00140 }
00141 
00142 VectorSpace::multi_vec_mut_ptr_t
00143 VectorSpaceBlocked::create_members(size_type num_vecs) const
00144 {
00145   TEST_FOR_EXCEPT(true); // ToDo: Implement using MultiVectorMutableCompositeStd!
00146   return Teuchos::null;
00147 }
00148 
00149 VectorSpace::space_ptr_t
00150 VectorSpaceBlocked::clone() const
00151 {
00152   return Teuchos::rcp(new VectorSpaceBlocked(*this)); // ToDo: Fix the behavior when needed!
00153 }
00154 
00155 VectorSpace::space_ptr_t
00156 VectorSpaceBlocked::sub_space(const Range1D& rng_in) const
00157 {
00158   namespace rcp = MemMngPack;
00159   const index_type dim = this->dim();
00160   const Range1D    rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
00161   // Validate the preconditions
00162 #ifdef TEUCHOS_DEBUG
00163   TEST_FOR_EXCEPTION(
00164     dim < rng.ubound(), std::out_of_range
00165     ,"VectorSpaceBlocked::sub_space(...): Error, rng = "
00166     << "["<<rng.lbound()<<","<<rng.ubound()<<"] is not in range [1,"<<dim<<"]" );
00167 #endif
00168   if( rng.lbound() == 1 && rng.ubound() == dim )
00169     return space_ptr_t( this, false ); // Client selected the whole composite vector space.
00170   // Get the position of the vector space object of interest
00171   int           kth_vector_space  = -1;
00172   index_type    kth_global_offset = 0;
00173   this->get_vector_space_position(rng.lbound(),&kth_vector_space,&kth_global_offset);
00174   const vector_spaces_t      &vector_spaces      = vector_spaces_;      // Need to examine in debugger!
00175   const vec_spaces_offsets_t &vec_spaces_offsets = vec_spaces_offsets_;
00176 #ifdef TEUCHOS_DEBUG
00177   TEST_FOR_EXCEPT( !(  0 <= kth_vector_space && kth_vector_space <= vector_spaces.size()  ) );
00178 #endif
00179   if( rng.lbound() == kth_global_offset + 1
00180     && rng.size() == vec_spaces_offsets[kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
00181     // The client selected a whole single constituent vector space.
00182     return vector_spaces[kth_vector_space];
00183   if( rng.ubound() <= vec_spaces_offsets[kth_vector_space+1] )
00184     // The client selected a sub-space of a single consituent vector space
00185     return vector_spaces[kth_vector_space]->sub_space(rng-vec_spaces_offsets[kth_vector_space]);
00186   // The client selected a sub-space that spans two or more constituent vector spaces
00187   // Get the position of the vector space object with the last element of interest
00188   int           end_kth_vector_space  = -1;
00189   index_type    end_kth_global_offset = 0;
00190   this->get_vector_space_position(rng.ubound(),&end_kth_vector_space,&end_kth_global_offset);
00191 #ifdef TEUCHOS_DEBUG
00192   TEST_FOR_EXCEPT( !(  0 <= end_kth_vector_space && end_kth_vector_space <= vector_spaces.size()  ) );
00193   TEST_FOR_EXCEPT( !(  end_kth_vector_space > kth_vector_space  ) );
00194 #endif
00195   // Create a VectorSpaceComposite object containing the relavant constituent vector spaces
00196   Teuchos::RCP<VectorSpaceBlocked>
00197     vec_space_comp = Teuchos::rcp(
00198       new VectorSpaceBlocked(
00199         &vector_spaces[kth_vector_space]
00200         ,end_kth_vector_space - kth_vector_space + 1 )
00201       );
00202   if( rng.lbound() == kth_global_offset + 1
00203     && rng.size() == vec_spaces_offsets[end_kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
00204     // The client selected exactly a contigous set of vector spaces
00205     return vec_space_comp;
00206   // The client selected some sub-set of elements in the contigous set of vector spaces
00207   return Teuchos::rcp(
00208     new VectorSpaceSubSpace(
00209       vec_space_comp
00210       ,Range1D( 
00211         rng.lbound()-vec_spaces_offsets[kth_vector_space]
00212         ,rng.ubound()-vec_spaces_offsets[kth_vector_space] )
00213       ) );
00214 }
00215 
00216 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:09:16 2011 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.6.3