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