1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_COMPRESSED_STORAGE_H 11 #define EIGEN_COMPRESSED_STORAGE_H 12 13 namespace Eigen { 14 15 namespace internal { 16 17 /** \internal 18 * Stores a sparse set of values as a list of values and a list of indices. 19 * 20 */ 21 template<typename _Scalar,typename _Index> 22 class CompressedStorage 23 { 24 public: 25 26 typedef _Scalar Scalar; 27 typedef _Index Index; 28 29 protected: 30 31 typedef typename NumTraits<Scalar>::Real RealScalar; 32 33 public: 34 CompressedStorage()35 CompressedStorage() 36 : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) 37 {} 38 CompressedStorage(size_t size)39 CompressedStorage(size_t size) 40 : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) 41 { 42 resize(size); 43 } 44 CompressedStorage(const CompressedStorage & other)45 CompressedStorage(const CompressedStorage& other) 46 : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0) 47 { 48 *this = other; 49 } 50 51 CompressedStorage& operator=(const CompressedStorage& other) 52 { 53 resize(other.size()); 54 internal::smart_copy(other.m_values, other.m_values + m_size, m_values); 55 internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices); 56 return *this; 57 } 58 swap(CompressedStorage & other)59 void swap(CompressedStorage& other) 60 { 61 std::swap(m_values, other.m_values); 62 std::swap(m_indices, other.m_indices); 63 std::swap(m_size, other.m_size); 64 std::swap(m_allocatedSize, other.m_allocatedSize); 65 } 66 ~CompressedStorage()67 ~CompressedStorage() 68 { 69 delete[] m_values; 70 delete[] m_indices; 71 } 72 reserve(size_t size)73 void reserve(size_t size) 74 { 75 size_t newAllocatedSize = m_size + size; 76 if (newAllocatedSize > m_allocatedSize) 77 reallocate(newAllocatedSize); 78 } 79 squeeze()80 void squeeze() 81 { 82 if (m_allocatedSize>m_size) 83 reallocate(m_size); 84 } 85 86 void resize(size_t size, double reserveSizeFactor = 0) 87 { 88 if (m_allocatedSize<size) 89 reallocate(size + size_t(reserveSizeFactor*double(size))); 90 m_size = size; 91 } 92 append(const Scalar & v,Index i)93 void append(const Scalar& v, Index i) 94 { 95 Index id = static_cast<Index>(m_size); 96 resize(m_size+1, 1); 97 m_values[id] = v; 98 m_indices[id] = i; 99 } 100 size()101 inline size_t size() const { return m_size; } allocatedSize()102 inline size_t allocatedSize() const { return m_allocatedSize; } clear()103 inline void clear() { m_size = 0; } 104 value(size_t i)105 inline Scalar& value(size_t i) { return m_values[i]; } value(size_t i)106 inline const Scalar& value(size_t i) const { return m_values[i]; } 107 index(size_t i)108 inline Index& index(size_t i) { return m_indices[i]; } index(size_t i)109 inline const Index& index(size_t i) const { return m_indices[i]; } 110 Map(Index * indices,Scalar * values,size_t size)111 static CompressedStorage Map(Index* indices, Scalar* values, size_t size) 112 { 113 CompressedStorage res; 114 res.m_indices = indices; 115 res.m_values = values; 116 res.m_allocatedSize = res.m_size = size; 117 return res; 118 } 119 120 /** \returns the largest \c k such that for all \c j in [0,k) index[\c j]\<\a key */ searchLowerIndex(Index key)121 inline Index searchLowerIndex(Index key) const 122 { 123 return searchLowerIndex(0, m_size, key); 124 } 125 126 /** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */ searchLowerIndex(size_t start,size_t end,Index key)127 inline Index searchLowerIndex(size_t start, size_t end, Index key) const 128 { 129 while(end>start) 130 { 131 size_t mid = (end+start)>>1; 132 if (m_indices[mid]<key) 133 start = mid+1; 134 else 135 end = mid; 136 } 137 return static_cast<Index>(start); 138 } 139 140 /** \returns the stored value at index \a key 141 * If the value does not exist, then the value \a defaultValue is returned without any insertion. */ 142 inline Scalar at(Index key, const Scalar& defaultValue = Scalar(0)) const 143 { 144 if (m_size==0) 145 return defaultValue; 146 else if (key==m_indices[m_size-1]) 147 return m_values[m_size-1]; 148 // ^^ optimization: let's first check if it is the last coefficient 149 // (very common in high level algorithms) 150 const size_t id = searchLowerIndex(0,m_size-1,key); 151 return ((id<m_size) && (m_indices[id]==key)) ? m_values[id] : defaultValue; 152 } 153 154 /** Like at(), but the search is performed in the range [start,end) */ 155 inline Scalar atInRange(size_t start, size_t end, Index key, const Scalar& defaultValue = Scalar(0)) const 156 { 157 if (start>=end) 158 return Scalar(0); 159 else if (end>start && key==m_indices[end-1]) 160 return m_values[end-1]; 161 // ^^ optimization: let's first check if it is the last coefficient 162 // (very common in high level algorithms) 163 const size_t id = searchLowerIndex(start,end-1,key); 164 return ((id<end) && (m_indices[id]==key)) ? m_values[id] : defaultValue; 165 } 166 167 /** \returns a reference to the value at index \a key 168 * If the value does not exist, then the value \a defaultValue is inserted 169 * such that the keys are sorted. */ 170 inline Scalar& atWithInsertion(Index key, const Scalar& defaultValue = Scalar(0)) 171 { 172 size_t id = searchLowerIndex(0,m_size,key); 173 if (id>=m_size || m_indices[id]!=key) 174 { 175 resize(m_size+1,1); 176 for (size_t j=m_size-1; j>id; --j) 177 { 178 m_indices[j] = m_indices[j-1]; 179 m_values[j] = m_values[j-1]; 180 } 181 m_indices[id] = key; 182 m_values[id] = defaultValue; 183 } 184 return m_values[id]; 185 } 186 187 void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits<RealScalar>::dummy_precision()) 188 { 189 size_t k = 0; 190 size_t n = size(); 191 for (size_t i=0; i<n; ++i) 192 { 193 if (!internal::isMuchSmallerThan(value(i), reference, epsilon)) 194 { 195 value(k) = value(i); 196 index(k) = index(i); 197 ++k; 198 } 199 } 200 resize(k,0); 201 } 202 203 protected: 204 reallocate(size_t size)205 inline void reallocate(size_t size) 206 { 207 Scalar* newValues = new Scalar[size]; 208 Index* newIndices = new Index[size]; 209 size_t copySize = (std::min)(size, m_size); 210 // copy 211 internal::smart_copy(m_values, m_values+copySize, newValues); 212 internal::smart_copy(m_indices, m_indices+copySize, newIndices); 213 // delete old stuff 214 delete[] m_values; 215 delete[] m_indices; 216 m_values = newValues; 217 m_indices = newIndices; 218 m_allocatedSize = size; 219 } 220 221 protected: 222 Scalar* m_values; 223 Index* m_indices; 224 size_t m_size; 225 size_t m_allocatedSize; 226 227 }; 228 229 } // end namespace internal 230 231 } // end namespace Eigen 232 233 #endif // EIGEN_COMPRESSED_STORAGE_H 234