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

Generated on Wed May 12 21:50:44 2010 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.4.7