MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_VectorSpaceBlocked.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include <algorithm>
45 
51 #include "Teuchos_Workspace.hpp"
52 #include "Teuchos_Assert.hpp"
53 
54 namespace AbstractLinAlgPack {
55 
57  const VectorSpace::space_ptr_t vec_spaces[]
58  ,int num_vec_spaces
59  ,const VectorSpace::space_fcty_ptr_t &small_vec_spc_fcty
60  )
61 {
62  initialize(vec_spaces,num_vec_spaces,small_vec_spc_fcty);
63 }
64 
66  const VectorSpace::space_ptr_t vec_spaces[]
67  ,int num_vec_spaces
68  ,const VectorSpace::space_fcty_ptr_t &small_vec_spc_fcty
69  )
70 {
71  vector_spaces_.resize(num_vec_spaces);
72  std::copy(vec_spaces,vec_spaces+num_vec_spaces,vector_spaces_.begin());
73  vec_spaces_offsets_.resize(num_vec_spaces+1);
74  vec_spaces_offsets_[0] = 0;
75  for( int k = 1; k <= num_vec_spaces; ++k )
76  vec_spaces_offsets_[k] = vec_spaces_offsets_[k-1] + vec_spaces[k-1]->dim();
77  small_vec_spc_fcty_ = small_vec_spc_fcty;
78 }
79 
81  index_type i, int* kth_vector_space, index_type* kth_global_offset ) const
82 {
83  // Validate the preconditions
84 #ifdef TEUCHOS_DEBUG
86  i < 1 || this->dim() < i, std::out_of_range
87  ,"VectorSpaceBlocked::get_vector_space_position(...): Error, i = "
88  << i << " is not in range [1,"<<this->dim()<<"]"
89  );
90 #endif
91  *kth_vector_space = 0;
92  *kth_global_offset = 0;
93  while( *kth_vector_space < vector_spaces_.size() ) {
94  const RTOp_index_type off_kp1 = vec_spaces_offsets_[*kth_vector_space+1];
95  if( off_kp1 + 1 > i ) {
96  *kth_global_offset = vec_spaces_offsets_[*kth_vector_space];
97  break;
98  }
99  ++(*kth_vector_space);
100  }
101 #ifdef TEUCHOS_DEBUG
102  TEUCHOS_TEST_FOR_EXCEPT( !( *kth_vector_space < vector_spaces_.size() ) );
103 #endif
104 }
105 
106 // overridden from VectorSpace
107 
108 bool VectorSpaceBlocked::is_compatible(const VectorSpace& vec_space) const
109 {
110  const VectorSpaceBlocked
111  *vec_space_comp = dynamic_cast<const VectorSpaceBlocked*>(&vec_space);
112  if( !vec_space_comp || vec_space_comp->vector_spaces_.size() != this->vector_spaces_.size() )
113  return false;
114  // There is some hope that these are compatible vector spaces.
115  for( int k = 0; k < vector_spaces_.size(); ++k ) {
116  if( !vec_space_comp->vector_spaces_[k]->is_compatible(*this->vector_spaces_[k]) )
117  return false;
118  }
119  return true; // If we get here we are compatible!
120 }
121 
123 {
124  return vec_spaces_offsets_[vector_spaces_.size()];
125 }
126 
129 {
130  return small_vec_spc_fcty_;
131 }
132 
135 {
136  namespace rcp = MemMngPack;
137  using Teuchos::Workspace;
139 
140  const int num_vec_spaces = this->num_vector_spaces();
141 
142  // Create the vector objects array.
143  Workspace<VectorMutable::vec_mut_ptr_t>
144  vecs(wss,num_vec_spaces);
145  for( int k = 0; k < num_vec_spaces; ++k )
146  vecs[k] = vector_spaces_[k]->create_member();
147 
148  return Teuchos::rcp(
150  &vecs[0]
151  ,Teuchos::rcp(new VectorSpaceBlocked(*this))
152  ) );
153 }
154 
157 {
158  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement using MultiVectorMutableCompositeStd!
159  return Teuchos::null;
160 }
161 
164 {
165  return Teuchos::rcp(new VectorSpaceBlocked(*this)); // ToDo: Fix the behavior when needed!
166 }
167 
170 {
171  namespace rcp = MemMngPack;
172  const index_type dim = this->dim();
173  const Range1D rng = rng_in.full_range() ? Range1D(1,dim) : rng_in;
174  // Validate the preconditions
175 #ifdef TEUCHOS_DEBUG
177  dim < rng.ubound(), std::out_of_range
178  ,"VectorSpaceBlocked::sub_space(...): Error, rng = "
179  << "["<<rng.lbound()<<","<<rng.ubound()<<"] is not in range [1,"<<dim<<"]" );
180 #endif
181  if( rng.lbound() == 1 && rng.ubound() == dim )
182  return space_ptr_t( this, false ); // Client selected the whole composite vector space.
183  // Get the position of the vector space object of interest
184  int kth_vector_space = -1;
185  index_type kth_global_offset = 0;
186  this->get_vector_space_position(rng.lbound(),&kth_vector_space,&kth_global_offset);
187  const vector_spaces_t &vector_spaces = vector_spaces_; // Need to examine in debugger!
188  const vec_spaces_offsets_t &vec_spaces_offsets = vec_spaces_offsets_;
189 #ifdef TEUCHOS_DEBUG
190  TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= kth_vector_space && kth_vector_space <= vector_spaces.size() ) );
191 #endif
192  if( rng.lbound() == kth_global_offset + 1
193  && rng.size() == vec_spaces_offsets[kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
194  // The client selected a whole single constituent vector space.
195  return vector_spaces[kth_vector_space];
196  if( rng.ubound() <= vec_spaces_offsets[kth_vector_space+1] )
197  // The client selected a sub-space of a single consituent vector space
198  return vector_spaces[kth_vector_space]->sub_space(rng-vec_spaces_offsets[kth_vector_space]);
199  // The client selected a sub-space that spans two or more constituent vector spaces
200  // Get the position of the vector space object with the last element of interest
201  int end_kth_vector_space = -1;
202  index_type end_kth_global_offset = 0;
203  this->get_vector_space_position(rng.ubound(),&end_kth_vector_space,&end_kth_global_offset);
204 #ifdef TEUCHOS_DEBUG
205  TEUCHOS_TEST_FOR_EXCEPT( !( 0 <= end_kth_vector_space && end_kth_vector_space <= vector_spaces.size() ) );
206  TEUCHOS_TEST_FOR_EXCEPT( !( end_kth_vector_space > kth_vector_space ) );
207 #endif
208  // Create a VectorSpaceComposite object containing the relavant constituent vector spaces
210  vec_space_comp = Teuchos::rcp(
211  new VectorSpaceBlocked(
212  &vector_spaces[kth_vector_space]
213  ,end_kth_vector_space - kth_vector_space + 1 )
214  );
215  if( rng.lbound() == kth_global_offset + 1
216  && rng.size() == vec_spaces_offsets[end_kth_vector_space+1] - vec_spaces_offsets[kth_vector_space] )
217  // The client selected exactly a contigous set of vector spaces
218  return vec_space_comp;
219  // The client selected some sub-set of elements in the contigous set of vector spaces
220  return Teuchos::rcp(
222  vec_space_comp
223  ,Range1D(
224  rng.lbound()-vec_spaces_offsets[kth_vector_space]
225  ,rng.ubound()-vec_spaces_offsets[kth_vector_space] )
226  ) );
227 }
228 
229 } // end namespace AbstractLinAlgPack
Concrete subclass for a default sub-space of a vector.
void initialize(const VectorSpace::space_ptr_t vec_spaces[], int num_vector_spaces, const VectorSpace::space_fcty_ptr_t &small_vec_spc_fcty=Teuchos::null)
Initialize with a set of vector space objects.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Index size() const
Return the size of the range (ubound() - lbound() + 1)
Index ubound() const
Return upper bound of the range.
multi_vec_mut_ptr_t create_members(size_type num_vecs) const
space_fcty_ptr_t small_vec_spc_fcty() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
. One-based subregion index range class.
Teuchos::RCP< const VectorSpace > space_ptr_t
Abstract interface for objects that represent a space for mutable coordinate vectors.
void copy(const f_int &N, const f_dbl_prec *X, const f_int &INCX, f_dbl_prec *Y, const f_int &INCY)
VectorSpace subclass for the composite of one or more VectorSpace objects.
bool full_range() const
Returns true if the range represents the entire region (constructed from Range1D()) ...
const VectorSpace::space_ptr_t * vector_spaces() const
std::vector< VectorSpace::space_ptr_t > vector_spaces_t
Index lbound() const
Return lower bound of the range.
int num_vector_spaces() const
Return the value of num_vec_spaces passed to this->initialize().
bool is_compatible(const VectorSpace &vec_space) const
Returns true if same type and has compatible vector spaces.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void get_vector_space_position(index_type i, int *kth_vector_space, index_type *kth_global_offset) const
Get the position of the vector space object and its offset into the composite vector of the vector sp...
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
Teuchos_Ordinal RTOp_index_type
Definition: RTOp.h:68