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