AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_SparseVectorClassDef.hpp
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 #ifndef SPARSE_VECTOR_CLASS_DEF_H
43 #define SPARSE_VECTOR_CLASS_DEF_H
44 
45 #include <algorithm>
46 
47 #include "AbstractLinAlgPack_SparseVectorClassDecl.hpp"
48 #include "AbstractLinAlgPack_compare_element_indexes.hpp"
49 #include "Teuchos_Assert.hpp"
50 
51 namespace AbstractLinAlgPack {
52 
53 // Non-member non-public utility functions
54 
55 // /////////////////////////////////////////////////////////////////////////////////////
56 // Definitions of non-member functions
57 
58 template<class T_Element>
59 SparseVectorSlice<T_Element>
60 create_slice(const SparseVectorUtilityPack::SpVecIndexLookup<T_Element>& index_lookup
61  , size_type size, Range1D rng)
62 {
63  // Check preconditions
64  if(rng.full_range()) {
65  rng = Range1D(1,size);
66  }
67  else {
68 #ifdef TEUCHOS_DEBUG
69  TEUCHOS_TEST_FOR_EXCEPT( !( rng.ubound() <= size ) );
70 #endif
71  }
72 
73  // If there are no elements then any subregion will also not have any elements.
74  if(!index_lookup.nz())
75  return SparseVectorSlice<T_Element>(0, 0, 0, rng.ubound() - rng.lbound() + 1, true);
76 
77  // Create the slice (assumed sorted oviously).
78  typedef SparseVectorUtilityPack::SpVecIndexLookup<T_Element> SpVecIndexLookup;
79  index_lookup.validate_state();
80  typename SpVecIndexLookup::poss_type
81  lower_poss = index_lookup.find_poss(rng.lbound(), SpVecIndexLookup::LOWER_ELE),
82  upper_poss = index_lookup.find_poss(rng.ubound(), SpVecIndexLookup::UPPER_ELE);
83  if( lower_poss.poss == upper_poss.poss
84  && ( lower_poss.rel == SpVecIndexLookup::AFTER_ELE
85  || upper_poss.rel == SpVecIndexLookup::BEFORE_ELE )
86  )
87  {
88  // The requested subvector does not contain any elements.
89  return SparseVectorSlice<T_Element>(0, 0, 0, rng.ubound() - rng.lbound() + 1, true);
90  }
91  else {
92  // There are nonzero elements
93  return SparseVectorSlice<T_Element>(index_lookup.ele() + lower_poss.poss
94  , upper_poss.poss - lower_poss.poss + 1, index_lookup.offset() - rng.lbound() + 1
95  , rng.ubound() - rng.lbound() + 1, true);
96  }
97 }
98 
99 // /////////////////////////////////////////////////////////////////////////////////////
100 // Definitions of members for SparseVector<>
101 
102 template <class T_Element, class T_Alloc>
104  const SparseVector<T_Element,T_Alloc>& sp_vec )
105  :
106  alloc_(sp_vec.alloc_), size_(sp_vec.size_), max_nz_(sp_vec.max_nz_)
107  , assume_sorted_(sp_vec.assume_sorted_)
108  , know_is_sorted_(sp_vec.know_is_sorted_)
109 {
110  // Allocate the memory for the elements and set the memory of the sparse vector.
111  index_lookup_.set_sp_vec(
112 #ifdef _PG_CXX
113  new element_type[max_nz_]
114 #else
115  alloc_.allocate(max_nz_,NULL)
116 #endif
117  ,sp_vec.nz(),sp_vec.offset());
118  // Perform an uninitialized copy of the elements
119  iterator ele_to_itr = index_lookup_.ele();
120  const_iterator ele_from_itr = sp_vec.begin();
121  while(ele_from_itr != sp_vec.end()) {
122 #ifdef _PG_CXX
123  new (ele_to_itr++) element_type(*ele_from_itr++);
124 #else
125  alloc_.construct(ele_to_itr++,*ele_from_itr++);
126 #endif
127  }
128 }
129 
130 template <class T_Element, class T_Alloc>
133  , const allocator_type& alloc )
134  :
135  alloc_(alloc), size_(sp_vec_slc.dim()), max_nz_(sp_vec_slc.nz())
136  , assume_sorted_(sp_vec_slc.is_sorted())
137  , know_is_sorted_(false)
138 {
139  // Allocate the memory for the elements and set the memory of the sparse vector.
140  index_lookup_.set_sp_vec(
141 #ifdef _PG_CXX
142  new element_type[max_nz_]
143 #else
144  alloc_.allocate(max_nz_,NULL)
145 #endif
146  , sp_vec_slc.nz()
147  ,sp_vec_slc.offset() );
148  // Perform an uninitialized copy of the elements
149  iterator ele_to_itr = index_lookup_.ele();
150  const_iterator ele_from_itr = sp_vec_slc.begin();
151  while(ele_from_itr != sp_vec_slc.end()) {
152 #ifdef _PG_CXX
153  new (ele_to_itr++) element_type(*ele_from_itr++);
154 #else
155  alloc_.construct(ele_to_itr++,*ele_from_itr++);
156 #endif
157  }
158 }
159 
160 template <class T_Element, class T_Alloc>
163  const SparseVector<T_Element,T_Alloc>& sp_vec)
164 {
165  if(this == &sp_vec) return *this; // assignment to self
166 
167  know_is_sorted_ = sp_vec.know_is_sorted_;
168  assume_sorted_ = sp_vec.assume_sorted_;
169 
170  if( max_nz() < sp_vec.nz() ) {
171  // need to allocate more storage
172  resize(0,0); // free current storage
173  resize(sp_vec.dim(),sp_vec.nz(),sp_vec.offset());
174  }
175  else if( nz() ) {
176  // Don't allocate new memory, just call distructors on current elements
177  // and reset to uninitialized.
178  for(iterator ele_itr = begin(); ele_itr != end();) {
179 #ifdef _PG_CXX
180  (ele_itr++)->~element_type();
181 #else
182  alloc_.destroy(ele_itr++);
183 #endif
184  }
185  // Set the other data
186  size_ = sp_vec.size_;
187  }
188 
189  // set nz and offset
190  index_lookup_.set_sp_vec(index_lookup_.ele(),sp_vec.nz(),sp_vec.offset());
191 
192  if( sp_vec.nz() ) {
193  // Perform an uninitialized copy of the elements
194  iterator ele_to_itr = index_lookup_.ele();
195  const_iterator ele_from_itr = sp_vec.begin();
196  while(ele_from_itr != sp_vec.end()) {
197 #ifdef _PG_CXX
198  new (ele_to_itr++) element_type(*ele_from_itr++);
199 #else
200  alloc_.construct(ele_to_itr++,*ele_from_itr++);
201 #endif
202  }
203  }
204 
205  return *this;
206 }
207 
208 template <class T_Element, class T_Alloc>
211  const SparseVectorSlice<T_Element>& sp_vec_slc )
212 {
213  know_is_sorted_ = false;
214  assume_sorted_ = sp_vec_slc.is_sorted();
215 
216  if(max_nz() < sp_vec_slc.nz()) {
217  // need to allocate more storage
218  resize(0,0); // free current storage
219  resize(sp_vec_slc.dim(),sp_vec_slc.nz(),sp_vec_slc.offset());
220  }
221  else if( nz() ) {
222  // Don't allocate new memory, just call distructors on current elements
223  // and reset to uninitialized.
224  for(iterator ele_itr = begin(); ele_itr != end();) {
225 #ifdef _PG_CXX
226  (ele_itr++)->~element_type();
227 #else
228  alloc_.destroy(ele_itr++);
229 #endif
230  }
231  // Set the other data
232  size_ = sp_vec_slc.dim();
233  }
234 
235  // set nz and offset
236  index_lookup_.set_sp_vec(index_lookup_.ele(),sp_vec_slc.nz()
237  ,sp_vec_slc.offset());
238 
239  if( sp_vec_slc.nz() ) {
240  // Perform an uninitialized copy of the elements
241  iterator ele_to_itr = index_lookup_.ele();
242  const_iterator ele_from_itr = sp_vec_slc.begin();
243  while(ele_from_itr != sp_vec_slc.end()) {
244 #ifdef _PG_CXX
245  new (ele_to_itr++) element_type(*ele_from_itr++);
246 #else
247  alloc_.construct(ele_to_itr++,*ele_from_itr++);
248 #endif
249  }
250  }
251 
252  return *this;
253 }
254 
255 template <class T_Element, class T_Alloc>
257 {
258  if(!sv.dim()) return DenseLinAlgPack::NO_OVERLAP;
259 
260  const_iterator this_begin = begin();
261  typename SparseVectorSlice<T_Element>::const_iterator sv_begin = sv.begin();
262 
263  if( this_begin == sv_begin && end() == sv.end() )
264  {
265  return DenseLinAlgPack::SAME_MEM;
266  }
267 
268  if( ( this_begin < sv_begin && end() < sv_begin )
269  || ( sv_begin < this_begin && sv.end() < this_begin ) )
270  {
271  return DenseLinAlgPack::NO_OVERLAP;
272  }
273 
274  return DenseLinAlgPack::SOME_OVERLAP;
275 }
276 
277 template <class T_Element, class T_Alloc>
279  , difference_type offset)
280 {
281  // free existing storage
282  if(index_lookup_.ele()) {
283  for(element_type* p = index_lookup_.ele(); p < index_lookup_.ele() + index_lookup_.nz(); ++p) {
284 #ifdef _PG_CXX
285  p->~element_type();
286 #else
287  alloc_.destroy(p);
288 #endif
289  }
290 #ifdef _PG_CXX
291  delete [] index_lookup_.ele();
292 #else
293  alloc_.deallocate(index_lookup_.ele(), max_nz_);
294 #endif
295  }
296 
297  // reinitialize
298  max_nz_ = 0;
299  know_is_sorted_ = false;
300 
301  if(max_nz) {
302  // reallocate storage
303 #ifdef _PG_CXX
304  index_lookup_.set_sp_vec( new element_type[max_nz_ = max_nz], 0, offset );
305 #else
306  index_lookup_.set_sp_vec( alloc_.allocate( max_nz_ = max_nz, 0 ), 0, offset );
307 #endif
308  size_ = size;
309  }
310  else {
311  // reinitialize to no storage
312  index_lookup_.set_sp_vec( 0, 0, offset );
313  size_ = size; // allow size to be nonzero with nz = 0
314  }
315 }
316 
317 template <class T_Element, class T_Alloc>
319  , difference_type offset)
320 {
321 #ifdef TEUCHOS_DEBUG
323  nz > max_nz, std::length_error
324  ,"SparseVector<...>::uninitialized_resize(...) : nz can not be greater"
325  " than max_nz" );
326 #endif
327  resize(size,max_nz,offset);
328  index_lookup_.set_sp_vec(index_lookup_.ele(), nz, index_lookup_.offset());
329 }
330 
331 template <class T_Element, class T_Alloc>
333 {
334  assert_space(1);
335  assert_is_sorted();
336  // Find the insertion point
337  if( nz() ) {
339  typedef typename SpVecIndexLookup::poss_type poss_type;
340  index_lookup_.validate_state();
341  poss_type poss
342  = ( nz()
343  ? index_lookup_.find_poss(ele.index(), SpVecIndexLookup::LOWER_ELE)
344  : poss_type(0,SpVecIndexLookup::BEFORE_ELE)
345  );
346  // Make sure this element does not already exist!
347 #ifdef TEUCHOS_DEBUG
349  nz() && poss.rel == SpVecIndexLookup::EQUAL_TO_ELE, std::length_error
350  ,"SparseVector<...>::insert_element(...) : Error, this index"
351  " all ready exists!" );
352 #endif
353  const size_type
354  insert_poss = (poss.rel == SpVecIndexLookup::BEFORE_ELE ? poss.poss : poss.poss+1);
355  // Copy elements out of the way to make room for inserted element
356  std::copy_backward( // This assumes element_type supports assignment!
357  index_lookup_.ele() + insert_poss, index_lookup_.ele() + index_lookup_.nz()
358  , index_lookup_.ele() + index_lookup_.nz() + 1 );
359  index_lookup_.ele()[insert_poss] = ele;
360  index_lookup_.incr_nz();
361  }
362  else { // The first element we are adding!
363  index_lookup_.ele()[0] = ele;
364  index_lookup_.incr_nz();
365  }
366 }
367 
368 template <class T_Element, class T_Alloc>
370  if( index_lookup_.nz() > 0 )
371  std::stable_sort(begin(),end(),compare_element_indexes_less<element_type>());
372  know_is_sorted_ = true;
373 }
374 
375 template <class T_Element, class T_Alloc>
377 {
378 
379  if(!index_lookup_.nz()) return; // An empty sparse vector is certainly valid
380 
381  // Loop through the elements. If they are sorted then duplicate
382  // elements will be adjacent to each other so they will be easy to
383  // find.
384  typename T_Element::index_type last_index;
385  for(T_Element* p = index_lookup_.ele();
386  p < index_lookup_.ele() + index_lookup_.nz(); ++p)
387  {
388  typename T_Element::index_type curr_index = p->index() + offset();
389 #ifdef TEUCHOS_DEBUG
391  (1 > curr_index) || (curr_index > dim()), std::out_of_range
392  ,"SparseVector<...>::assert_valid_and_sorted():"
393  << " Error, not in range: element (0-based) " << p - index_lookup_.ele() - 1
394  << " with index + offset = " << curr_index
395  << " is not in range" );
396 #endif
397  if(p == index_lookup_.ele()) { // skip these tests for the first element
398  last_index = curr_index;
399  continue;
400  }
401 #ifdef TEUCHOS_DEBUG
403  curr_index < last_index, NotSortedException
404  ,"SparseVector<...>::assert_valid_and_sorted():"
405  << " Error, not sorted: element (0-based) " << p - index_lookup_.ele() - 1
406  << " and " << p - index_lookup_.ele() << " are not in assending order" );
408  curr_index == last_index, DuplicateIndexesException
409  ,"SparseVector<...>::assert_valid_and_sorted():"
410  << " Error, duplicate indexes: element (0-based) " << p - index_lookup_.ele() - 1
411  << " and " << p - index_lookup_.ele() << " have duplicate indexes" );
412 #endif
413  last_index = curr_index;
414  }
415 }
416 
417 
418 // /////////////////////////////////////////////////////////////////////////////////////
419 // Definitions of members for SparseVectorSlice<>
420 
421 template <class T_Element>
423 {
424  if(!sv.dim()) return DenseLinAlgPack::NO_OVERLAP;
425 
426  const_iterator this_begin = begin(),
427  sv_begin = sv.begin();
428 
429  if( this_begin == sv_begin && end() == sv.end() )
430  {
431  return DenseLinAlgPack::SAME_MEM;
432  }
433 
434  if( ( this_begin < sv_begin && end() < sv_begin )
435  || ( sv_begin < this_begin && sv.end() < this_begin ) )
436  {
437  return DenseLinAlgPack::NO_OVERLAP;
438  }
439 
440  return DenseLinAlgPack::SOME_OVERLAP;
441 }
442 
443 } // end namespace AbstractLinAlgPack
444 
445 #endif // SPARSE_VECTOR_CLASS_DEF_H
SparseVector< T_Element, T_Alloc > & operator=(const SparseVector< T_Element, T_Alloc > &sp_vec)
Assignment operator.
difference_type offset() const
Return the offset for the indexes (ith real index = begin()[i-1]->index() + offset()# ...
void resize(size_type size, size_type max_nz, difference_type offset=0)
Resize to #size# with a maximum of max_nz# non-zero elements.
size_type dim() const
Return the number of elements in the full vector.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void set_sp_vec(element_type *ele, size_type nz, difference_type offset)
Set a new sparse vector.
bool is_sorted() const
Return true if the sequence is assumed sorted.
void assert_valid_and_sorted() const
Assert that sparse vector is sorted.
size_type dim() const
Return the number of elements in the full vector.
size_type nz() const
Return the number of non-zero elements.
size_t size_type
void uninitialized_resize(size_type size, size_type nz, size_type max_nz, difference_type offset=0)
Resize to #size# with a max_nz# uninitialized non-zero elements.
EOverLap overlap(const SparseVectorSlice< T_Element > &sv) const
difference_type offset() const
Return the offset for the indexes (ith real index = begin()[i-1]->index() + offset()# ...
void sort()
Sort the elements into assending order by index.
SparseVector(const allocator_type &alloc=allocator_type())
Constructs a sparse vector with no elements (nz() == dim() == 0#) and assumes the elements are not so...
EOverLap overlap(const SparseVectorSlice< T_Element > &sv) const
iterator begin()
Returns iterator that iterates forward through the nonzero elements.
size_type nz() const
Return the number of non-zero elements.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void insert_element(element_type ele)
Add an element into a sorted sequence.
void validate_state() const
Called by client to ensure that the internal state is valid.