AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_SparseVectorClassDef.hpp
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 #ifndef SPARSE_VECTOR_CLASS_DEF_H
00043 #define SPARSE_VECTOR_CLASS_DEF_H
00044 
00045 #include <algorithm>
00046 
00047 #include "AbstractLinAlgPack_SparseVectorClassDecl.hpp"
00048 #include "AbstractLinAlgPack_compare_element_indexes.hpp"
00049 #include "Teuchos_Assert.hpp"
00050 
00051 namespace AbstractLinAlgPack {
00052 
00053 // Non-member non-public utility functions
00054 
00055 // /////////////////////////////////////////////////////////////////////////////////////
00056 // Definitions of non-member functions
00057 
00058 template<class T_Element>
00059 SparseVectorSlice<T_Element>
00060 create_slice(const SparseVectorUtilityPack::SpVecIndexLookup<T_Element>& index_lookup
00061   , size_type size, Range1D rng)
00062 {
00063   // Check preconditions
00064   if(rng.full_range()) {
00065     rng = Range1D(1,size);
00066   }
00067   else {
00068 #ifdef TEUCHOS_DEBUG
00069     TEUCHOS_TEST_FOR_EXCEPT( !(  rng.ubound() <= size  ) );
00070 #endif
00071   }
00072 
00073   // If there are no elements then any subregion will also not have any elements.
00074   if(!index_lookup.nz())
00075     return SparseVectorSlice<T_Element>(0, 0, 0, rng.ubound() - rng.lbound() + 1, true);
00076 
00077   // Create the slice (assumed sorted oviously).
00078   typedef SparseVectorUtilityPack::SpVecIndexLookup<T_Element> SpVecIndexLookup;
00079   index_lookup.validate_state();
00080   typename SpVecIndexLookup::poss_type
00081     lower_poss = index_lookup.find_poss(rng.lbound(), SpVecIndexLookup::LOWER_ELE),
00082     upper_poss = index_lookup.find_poss(rng.ubound(), SpVecIndexLookup::UPPER_ELE);
00083   if( lower_poss.poss == upper_poss.poss
00084     && (  lower_poss.rel == SpVecIndexLookup::AFTER_ELE
00085       ||  upper_poss.rel == SpVecIndexLookup::BEFORE_ELE )
00086     )
00087   {
00088     // The requested subvector does not contain any elements.
00089     return SparseVectorSlice<T_Element>(0, 0, 0, rng.ubound() - rng.lbound() + 1, true);
00090   }
00091   else {
00092     // There are nonzero elements
00093     return SparseVectorSlice<T_Element>(index_lookup.ele() + lower_poss.poss
00094       , upper_poss.poss - lower_poss.poss + 1, index_lookup.offset() - rng.lbound() + 1
00095       , rng.ubound() - rng.lbound() + 1, true);
00096   }
00097 }
00098 
00099 // /////////////////////////////////////////////////////////////////////////////////////
00100 // Definitions of members for SparseVector<>
00101 
00102 template <class T_Element, class T_Alloc>
00103 SparseVector<T_Element,T_Alloc>::SparseVector(
00104     const SparseVector<T_Element,T_Alloc>& sp_vec )
00105   :
00106       alloc_(sp_vec.alloc_), size_(sp_vec.size_), max_nz_(sp_vec.max_nz_)
00107     , assume_sorted_(sp_vec.assume_sorted_)
00108     , know_is_sorted_(sp_vec.know_is_sorted_)
00109 {
00110   // Allocate the memory for the elements and set the memory of the sparse vector.
00111   index_lookup_.set_sp_vec(
00112 #ifdef _PG_CXX
00113     new element_type[max_nz_]
00114 #else
00115     alloc_.allocate(max_nz_,NULL)
00116 #endif
00117     ,sp_vec.nz(),sp_vec.offset());
00118   // Perform an uninitialized copy of the elements
00119   iterator    ele_to_itr    = index_lookup_.ele();
00120   const_iterator  ele_from_itr  = sp_vec.begin();
00121   while(ele_from_itr != sp_vec.end()) {
00122 #ifdef _PG_CXX
00123     new (ele_to_itr++) element_type(*ele_from_itr++);
00124 #else
00125     alloc_.construct(ele_to_itr++,*ele_from_itr++);
00126 #endif
00127   }
00128 }
00129 
00130 template <class T_Element, class T_Alloc>
00131 SparseVector<T_Element,T_Alloc>::SparseVector(
00132       SparseVectorSlice<T_Element> sp_vec_slc
00133     , const allocator_type& alloc )
00134   :
00135     alloc_(alloc), size_(sp_vec_slc.dim()), max_nz_(sp_vec_slc.nz())
00136     , assume_sorted_(sp_vec_slc.is_sorted())
00137     , know_is_sorted_(false)
00138 {
00139   // Allocate the memory for the elements and set the memory of the sparse vector.
00140   index_lookup_.set_sp_vec(
00141 #ifdef _PG_CXX
00142     new element_type[max_nz_]
00143 #else
00144     alloc_.allocate(max_nz_,NULL)
00145 #endif
00146     , sp_vec_slc.nz()
00147     ,sp_vec_slc.offset() );
00148   // Perform an uninitialized copy of the elements
00149   iterator    ele_to_itr    = index_lookup_.ele();
00150   const_iterator  ele_from_itr  = sp_vec_slc.begin();
00151   while(ele_from_itr != sp_vec_slc.end()) {
00152 #ifdef _PG_CXX
00153     new (ele_to_itr++) element_type(*ele_from_itr++);
00154 #else
00155     alloc_.construct(ele_to_itr++,*ele_from_itr++);
00156 #endif
00157   }
00158 }
00159 
00160 template <class T_Element, class T_Alloc>
00161 SparseVector<T_Element,T_Alloc>&
00162 SparseVector<T_Element,T_Alloc>::operator=(
00163   const SparseVector<T_Element,T_Alloc>& sp_vec)
00164 {
00165   if(this == &sp_vec) return *this; // assignment to self
00166 
00167   know_is_sorted_ = sp_vec.know_is_sorted_;
00168   assume_sorted_ = sp_vec.assume_sorted_;
00169 
00170   if( max_nz() < sp_vec.nz() ) {
00171     // need to allocate more storage
00172     resize(0,0);  // free current storage
00173     resize(sp_vec.dim(),sp_vec.nz(),sp_vec.offset());
00174   }
00175   else if( nz() ) {
00176     // Don't allocate new memory, just call distructors on current elements
00177     // and reset to uninitialized.
00178     for(iterator ele_itr = begin(); ele_itr != end();) {
00179 #ifdef _PG_CXX
00180       (ele_itr++)->~element_type();
00181 #else
00182       alloc_.destroy(ele_itr++);
00183 #endif
00184     }
00185     // Set the other data
00186     size_ = sp_vec.size_;
00187   }
00188   
00189   // set nz and offset
00190   index_lookup_.set_sp_vec(index_lookup_.ele(),sp_vec.nz(),sp_vec.offset()); 
00191 
00192   if( sp_vec.nz() ) {
00193     // Perform an uninitialized copy of the elements
00194     iterator    ele_to_itr    = index_lookup_.ele();
00195     const_iterator  ele_from_itr  = sp_vec.begin();
00196     while(ele_from_itr != sp_vec.end()) {
00197 #ifdef _PG_CXX
00198       new (ele_to_itr++) element_type(*ele_from_itr++);
00199 #else
00200       alloc_.construct(ele_to_itr++,*ele_from_itr++);
00201 #endif
00202     }
00203   }
00204 
00205   return *this;
00206 }
00207 
00208 template <class T_Element, class T_Alloc>
00209 SparseVector<T_Element,T_Alloc>&
00210 SparseVector<T_Element,T_Alloc>::operator=(
00211   const SparseVectorSlice<T_Element>& sp_vec_slc )
00212 {
00213   know_is_sorted_ = false;
00214   assume_sorted_ = sp_vec_slc.is_sorted();
00215 
00216   if(max_nz() < sp_vec_slc.nz()) {
00217     // need to allocate more storage
00218     resize(0,0);  // free current storage
00219     resize(sp_vec_slc.dim(),sp_vec_slc.nz(),sp_vec_slc.offset());
00220   }
00221   else if( nz() ) {
00222     // Don't allocate new memory, just call distructors on current elements
00223     // and reset to uninitialized.
00224     for(iterator ele_itr = begin(); ele_itr != end();) {
00225 #ifdef _PG_CXX
00226       (ele_itr++)->~element_type();
00227 #else
00228       alloc_.destroy(ele_itr++);
00229 #endif
00230     }
00231     // Set the other data
00232     size_ = sp_vec_slc.dim();
00233   }
00234   
00235   // set nz and offset
00236   index_lookup_.set_sp_vec(index_lookup_.ele(),sp_vec_slc.nz()
00237     ,sp_vec_slc.offset()); 
00238 
00239   if( sp_vec_slc.nz() ) {
00240     // Perform an uninitialized copy of the elements
00241     iterator    ele_to_itr    = index_lookup_.ele();
00242     const_iterator  ele_from_itr  = sp_vec_slc.begin();
00243     while(ele_from_itr != sp_vec_slc.end()) {
00244 #ifdef _PG_CXX
00245       new (ele_to_itr++) element_type(*ele_from_itr++);
00246 #else
00247       alloc_.construct(ele_to_itr++,*ele_from_itr++);
00248 #endif
00249     }
00250   }
00251 
00252   return *this;
00253 }
00254 
00255 template <class T_Element, class T_Alloc>
00256 EOverLap SparseVector<T_Element,T_Alloc>::overlap(const SparseVectorSlice<T_Element>& sv) const
00257 {
00258   if(!sv.dim()) return DenseLinAlgPack::NO_OVERLAP;
00259 
00260   const_iterator                    this_begin  = begin();
00261   typename SparseVectorSlice<T_Element>::const_iterator   sv_begin  = sv.begin();
00262 
00263   if( this_begin == sv_begin && end() == sv.end() )
00264   {
00265     return DenseLinAlgPack::SAME_MEM;
00266   }
00267 
00268   if(   ( this_begin < sv_begin && end() < sv_begin )
00269     ||  ( sv_begin < this_begin && sv.end() < this_begin )  )
00270   {
00271     return DenseLinAlgPack::NO_OVERLAP;
00272   }
00273 
00274   return DenseLinAlgPack::SOME_OVERLAP;
00275 }
00276 
00277 template <class T_Element, class T_Alloc>
00278 void SparseVector<T_Element,T_Alloc>::resize(size_type size, size_type max_nz
00279   , difference_type offset)
00280 {
00281   // free existing storage
00282   if(index_lookup_.ele()) {
00283     for(element_type* p = index_lookup_.ele(); p < index_lookup_.ele() + index_lookup_.nz(); ++p) {
00284 #ifdef _PG_CXX
00285       p->~element_type();
00286 #else
00287       alloc_.destroy(p);
00288 #endif
00289     }
00290 #ifdef _PG_CXX
00291     delete [] index_lookup_.ele();
00292 #else
00293     alloc_.deallocate(index_lookup_.ele(), max_nz_);
00294 #endif
00295   }
00296   
00297   // reinitialize
00298   max_nz_ = 0;
00299   know_is_sorted_ = false;
00300 
00301   if(max_nz) {
00302     // reallocate storage
00303 #ifdef _PG_CXX
00304     index_lookup_.set_sp_vec( new element_type[max_nz_ = max_nz], 0, offset );
00305 #else
00306     index_lookup_.set_sp_vec( alloc_.allocate( max_nz_ = max_nz, 0 ), 0, offset );
00307 #endif
00308     size_ = size;
00309   }
00310   else {
00311     // reinitialize to no storage
00312     index_lookup_.set_sp_vec( 0, 0, offset );
00313     size_ = size; // allow size to be nonzero with nz = 0
00314   }
00315 }
00316 
00317 template <class T_Element, class T_Alloc>
00318 void SparseVector<T_Element,T_Alloc>::uninitialized_resize(size_type size, size_type nz, size_type max_nz
00319   , difference_type offset)
00320 {
00321 #ifdef TEUCHOS_DEBUG
00322   TEUCHOS_TEST_FOR_EXCEPTION(
00323     nz > max_nz, std::length_error
00324     ,"SparseVector<...>::uninitialized_resize(...) : nz can not be greater"
00325     " than max_nz" );
00326 #endif
00327   resize(size,max_nz,offset);
00328   index_lookup_.set_sp_vec(index_lookup_.ele(), nz, index_lookup_.offset());
00329 }
00330 
00331 template <class T_Element, class T_Alloc>
00332 void SparseVector<T_Element,T_Alloc>::insert_element(element_type ele)
00333 {
00334   assert_space(1);
00335   assert_is_sorted();
00336   // Find the insertion point
00337   if( nz() ) {
00338     typedef SparseVectorUtilityPack::SpVecIndexLookup<T_Element> SpVecIndexLookup;
00339     typedef typename SpVecIndexLookup::poss_type poss_type;
00340     index_lookup_.validate_state();
00341     poss_type poss
00342       = ( nz()
00343           ? index_lookup_.find_poss(ele.index(), SpVecIndexLookup::LOWER_ELE)
00344           : poss_type(0,SpVecIndexLookup::BEFORE_ELE)
00345         );
00346     // Make sure this element does not already exist!
00347 #ifdef TEUCHOS_DEBUG
00348     TEUCHOS_TEST_FOR_EXCEPTION(
00349       nz() && poss.rel == SpVecIndexLookup::EQUAL_TO_ELE, std::length_error
00350       ,"SparseVector<...>::insert_element(...) : Error, this index"
00351       " all ready exists!" );
00352 #endif
00353     const size_type
00354       insert_poss = (poss.rel == SpVecIndexLookup::BEFORE_ELE ? poss.poss : poss.poss+1);
00355     // Copy elements out of the way to make room for inserted element
00356     std::copy_backward( // This assumes element_type supports assignment!
00357       index_lookup_.ele() + insert_poss, index_lookup_.ele() + index_lookup_.nz()
00358       , index_lookup_.ele() + index_lookup_.nz() + 1 );
00359     index_lookup_.ele()[insert_poss] = ele;
00360     index_lookup_.incr_nz();
00361   }
00362   else { // The first element we are adding!
00363     index_lookup_.ele()[0] = ele;
00364     index_lookup_.incr_nz();
00365   }
00366 }
00367 
00368 template <class T_Element, class T_Alloc>
00369 void SparseVector<T_Element,T_Alloc>::sort() {
00370   if( index_lookup_.nz() > 0 )
00371     std::stable_sort(begin(),end(),compare_element_indexes_less<element_type>());
00372   know_is_sorted_ = true;
00373 }
00374 
00375 template <class T_Element, class T_Alloc>
00376 void SparseVector<T_Element,T_Alloc>::assert_valid_and_sorted() const
00377 {
00378 
00379   if(!index_lookup_.nz()) return; // An empty sparse vector is certainly valid
00380 
00381   // Loop through the elements.  If they are sorted then duplicate
00382   // elements will be adjacent to each other so they will be easy to 
00383   // find.
00384   typename T_Element::index_type last_index;
00385   for(T_Element* p = index_lookup_.ele();
00386     p < index_lookup_.ele() + index_lookup_.nz(); ++p)
00387   {
00388     typename T_Element::index_type curr_index = p->index() + offset();
00389 #ifdef TEUCHOS_DEBUG
00390     TEUCHOS_TEST_FOR_EXCEPTION(
00391       (1 > curr_index) || (curr_index > dim()), std::out_of_range
00392       ,"SparseVector<...>::assert_valid_and_sorted():"
00393       << " Error, not in range:  element (0-based) " << p - index_lookup_.ele() - 1
00394       << " with index + offset = " << curr_index
00395       << " is not in range" );
00396 #endif
00397     if(p == index_lookup_.ele()) { // skip these tests for the first element
00398       last_index = curr_index;
00399       continue;
00400     }
00401 #ifdef TEUCHOS_DEBUG
00402     TEUCHOS_TEST_FOR_EXCEPTION(
00403       curr_index < last_index, NotSortedException
00404       ,"SparseVector<...>::assert_valid_and_sorted():"
00405       << " Error, not sorted:  element (0-based) " << p - index_lookup_.ele() - 1
00406       << " and " << p - index_lookup_.ele() << " are not in assending order" );
00407     TEUCHOS_TEST_FOR_EXCEPTION(
00408       curr_index == last_index, DuplicateIndexesException
00409       ,"SparseVector<...>::assert_valid_and_sorted():"
00410       << " Error, duplicate indexes:  element (0-based) " << p - index_lookup_.ele() - 1
00411       << " and " << p - index_lookup_.ele() << " have duplicate indexes" );
00412 #endif
00413     last_index = curr_index;
00414   }
00415 }
00416 
00417 
00418 // /////////////////////////////////////////////////////////////////////////////////////
00419 // Definitions of members for SparseVectorSlice<>
00420 
00421 template <class T_Element>
00422 EOverLap SparseVectorSlice<T_Element>::overlap(const SparseVectorSlice<T_Element>& sv) const
00423 {
00424   if(!sv.dim()) return DenseLinAlgPack::NO_OVERLAP;
00425 
00426   const_iterator          this_begin  = begin(),
00427                   sv_begin  = sv.begin();
00428 
00429   if( this_begin == sv_begin && end() == sv.end() )
00430   {
00431     return DenseLinAlgPack::SAME_MEM;
00432   }
00433 
00434   if(   ( this_begin < sv_begin && end() < sv_begin )
00435     ||  ( sv_begin < this_begin && sv.end() < this_begin )  )
00436   {
00437     return DenseLinAlgPack::NO_OVERLAP;
00438   }
00439 
00440   return DenseLinAlgPack::SOME_OVERLAP;
00441 }
00442 
00443 } // end namespace AbstractLinAlgPack 
00444 
00445 #endif // SPARSE_VECTOR_CLASS_DEF_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends