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