00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. Eigen itself is part of the KDE project. 00003 // 00004 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> 00005 // 00006 // Eigen is free software; you can redistribute it and/or 00007 // modify it under the terms of the GNU Lesser General Public 00008 // License as published by the Free Software Foundation; either 00009 // version 3 of the License, or (at your option) any later version. 00010 // 00011 // Alternatively, you can redistribute it and/or 00012 // modify it under the terms of the GNU General Public License as 00013 // published by the Free Software Foundation; either version 2 of 00014 // the License, or (at your option) any later version. 00015 // 00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 00019 // GNU General Public License for more details. 00020 // 00021 // You should have received a copy of the GNU Lesser General Public 00022 // License and a copy of the GNU General Public License along with 00023 // Eigen. If not, see <http://www.gnu.org/licenses/>. 00024 00025 #ifndef EIGEN_AMBIVECTOR_H 00026 #define EIGEN_AMBIVECTOR_H 00027 00033 template<typename _Scalar> class AmbiVector 00034 { 00035 public: 00036 typedef _Scalar Scalar; 00037 typedef typename NumTraits<Scalar>::Real RealScalar; 00038 AmbiVector(int size) 00039 : m_buffer(0), m_size(0), m_allocatedSize(0), m_allocatedElements(0), m_mode(-1) 00040 { 00041 resize(size); 00042 } 00043 00044 void init(RealScalar estimatedDensity); 00045 void init(int mode); 00046 00047 void nonZeros() const; 00048 00050 void setBounds(int start, int end) { m_start = start; m_end = end; } 00051 00052 void setZero(); 00053 00054 void restart(); 00055 Scalar& coeffRef(int i); 00056 Scalar coeff(int i); 00057 00058 class Iterator; 00059 00060 ~AmbiVector() { delete[] m_buffer; } 00061 00062 void resize(int size) 00063 { 00064 if (m_allocatedSize < size) 00065 reallocate(size); 00066 m_size = size; 00067 } 00068 00069 int size() const { return m_size; } 00070 00071 protected: 00072 00073 void reallocate(int size) 00074 { 00075 // if the size of the matrix is not too large, let's allocate a bit more than needed such 00076 // that we can handle dense vector even in sparse mode. 00077 delete[] m_buffer; 00078 if (size<1000) 00079 { 00080 int allocSize = (size * sizeof(ListEl))/sizeof(Scalar); 00081 m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl); 00082 m_buffer = new Scalar[allocSize]; 00083 } 00084 else 00085 { 00086 m_allocatedElements = (size*sizeof(Scalar))/sizeof(ListEl); 00087 m_buffer = new Scalar[size]; 00088 } 00089 m_size = size; 00090 m_start = 0; 00091 m_end = m_size; 00092 } 00093 00094 void reallocateSparse() 00095 { 00096 int copyElements = m_allocatedElements; 00097 m_allocatedElements = std::min(int(m_allocatedElements*1.5),m_size); 00098 int allocSize = m_allocatedElements * sizeof(ListEl); 00099 allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0); 00100 Scalar* newBuffer = new Scalar[allocSize]; 00101 std::memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); 00102 delete[] m_buffer; 00103 m_buffer = newBuffer; 00104 } 00105 00106 protected: 00107 // element type of the linked list 00108 struct ListEl 00109 { 00110 int next; 00111 int index; 00112 Scalar value; 00113 }; 00114 00115 // used to store data in both mode 00116 Scalar* m_buffer; 00117 int m_size; 00118 int m_start; 00119 int m_end; 00120 int m_allocatedSize; 00121 int m_allocatedElements; 00122 int m_mode; 00123 00124 // linked list mode 00125 int m_llStart; 00126 int m_llCurrent; 00127 int m_llSize; 00128 00129 private: 00130 AmbiVector(const AmbiVector&); 00131 00132 }; 00133 00135 template<typename Scalar> 00136 void AmbiVector<Scalar>::nonZeros() const 00137 { 00138 if (m_mode==IsSparse) 00139 return m_llSize; 00140 else 00141 return m_end - m_start; 00142 } 00143 00144 template<typename Scalar> 00145 void AmbiVector<Scalar>::init(RealScalar estimatedDensity) 00146 { 00147 if (estimatedDensity>0.1) 00148 init(IsDense); 00149 else 00150 init(IsSparse); 00151 } 00152 00153 template<typename Scalar> 00154 void AmbiVector<Scalar>::init(int mode) 00155 { 00156 m_mode = mode; 00157 if (m_mode==IsSparse) 00158 { 00159 m_llSize = 0; 00160 m_llStart = -1; 00161 } 00162 } 00163 00169 template<typename Scalar> 00170 void AmbiVector<Scalar>::restart() 00171 { 00172 m_llCurrent = m_llStart; 00173 } 00174 00176 template<typename Scalar> 00177 void AmbiVector<Scalar>::setZero() 00178 { 00179 if (m_mode==IsDense) 00180 { 00181 for (int i=m_start; i<m_end; ++i) 00182 m_buffer[i] = Scalar(0); 00183 } 00184 else 00185 { 00186 ei_assert(m_mode==IsSparse); 00187 m_llSize = 0; 00188 m_llStart = -1; 00189 } 00190 } 00191 00192 template<typename Scalar> 00193 Scalar& AmbiVector<Scalar>::coeffRef(int i) 00194 { 00195 if (m_mode==IsDense) 00196 return m_buffer[i]; 00197 else 00198 { 00199 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer); 00200 // TODO factorize the following code to reduce code generation 00201 ei_assert(m_mode==IsSparse); 00202 if (m_llSize==0) 00203 { 00204 // this is the first element 00205 m_llStart = 0; 00206 m_llCurrent = 0; 00207 ++m_llSize; 00208 llElements[0].value = Scalar(0); 00209 llElements[0].index = i; 00210 llElements[0].next = -1; 00211 return llElements[0].value; 00212 } 00213 else if (i<llElements[m_llStart].index) 00214 { 00215 // this is going to be the new first element of the list 00216 ListEl& el = llElements[m_llSize]; 00217 el.value = Scalar(0); 00218 el.index = i; 00219 el.next = m_llStart; 00220 m_llStart = m_llSize; 00221 ++m_llSize; 00222 m_llCurrent = m_llStart; 00223 return el.value; 00224 } 00225 else 00226 { 00227 int nextel = llElements[m_llCurrent].next; 00228 ei_assert(i>=llElements[m_llCurrent].index && "you must call restart() before inserting an element with lower or equal index"); 00229 while (nextel >= 0 && llElements[nextel].index<=i) 00230 { 00231 m_llCurrent = nextel; 00232 nextel = llElements[nextel].next; 00233 } 00234 00235 if (llElements[m_llCurrent].index==i) 00236 { 00237 // the coefficient already exists and we found it ! 00238 return llElements[m_llCurrent].value; 00239 } 00240 else 00241 { 00242 if (m_llSize>=m_allocatedElements) 00243 { 00244 reallocateSparse(); 00245 llElements = reinterpret_cast<ListEl*>(m_buffer); 00246 } 00247 ei_internal_assert(m_llSize<m_allocatedElements && "internal error: overflow in sparse mode"); 00248 // let's insert a new coefficient 00249 ListEl& el = llElements[m_llSize]; 00250 el.value = Scalar(0); 00251 el.index = i; 00252 el.next = llElements[m_llCurrent].next; 00253 llElements[m_llCurrent].next = m_llSize; 00254 ++m_llSize; 00255 return el.value; 00256 } 00257 } 00258 } 00259 } 00260 00261 template<typename Scalar> 00262 Scalar AmbiVector<Scalar>::coeff(int i) 00263 { 00264 if (m_mode==IsDense) 00265 return m_buffer[i]; 00266 else 00267 { 00268 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_buffer); 00269 ei_assert(m_mode==IsSparse); 00270 if ((m_llSize==0) || (i<llElements[m_llStart].index)) 00271 { 00272 return Scalar(0); 00273 } 00274 else 00275 { 00276 int elid = m_llStart; 00277 while (elid >= 0 && llElements[elid].index<i) 00278 elid = llElements[elid].next; 00279 00280 if (llElements[elid].index==i) 00281 return llElements[m_llCurrent].value; 00282 else 00283 return Scalar(0); 00284 } 00285 } 00286 } 00287 00289 template<typename _Scalar> 00290 class AmbiVector<_Scalar>::Iterator 00291 { 00292 public: 00293 typedef _Scalar Scalar; 00294 typedef typename NumTraits<Scalar>::Real RealScalar; 00295 00302 Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision<RealScalar>()) 00303 : m_vector(vec) 00304 { 00305 m_epsilon = epsilon; 00306 m_isDense = m_vector.m_mode==IsDense; 00307 if (m_isDense) 00308 { 00309 m_cachedIndex = m_vector.m_start-1; 00310 ++(*this); 00311 } 00312 else 00313 { 00314 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer); 00315 m_currentEl = m_vector.m_llStart; 00316 while (m_currentEl>=0 && ei_abs(llElements[m_currentEl].value)<m_epsilon) 00317 m_currentEl = llElements[m_currentEl].next; 00318 if (m_currentEl<0) 00319 { 00320 m_cachedIndex = -1; 00321 } 00322 else 00323 { 00324 m_cachedIndex = llElements[m_currentEl].index; 00325 m_cachedValue = llElements[m_currentEl].value; 00326 } 00327 } 00328 } 00329 00330 int index() const { return m_cachedIndex; } 00331 Scalar value() const { return m_cachedValue; } 00332 00333 operator bool() const { return m_cachedIndex>=0; } 00334 00335 Iterator& operator++() 00336 { 00337 if (m_isDense) 00338 { 00339 do { 00340 ++m_cachedIndex; 00341 } while (m_cachedIndex<m_vector.m_end && ei_abs(m_vector.m_buffer[m_cachedIndex])<m_epsilon); 00342 if (m_cachedIndex<m_vector.m_end) 00343 m_cachedValue = m_vector.m_buffer[m_cachedIndex]; 00344 else 00345 m_cachedIndex=-1; 00346 } 00347 else 00348 { 00349 ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer); 00350 do { 00351 m_currentEl = llElements[m_currentEl].next; 00352 } while (m_currentEl>=0 && ei_abs(llElements[m_currentEl].value)<m_epsilon); 00353 if (m_currentEl<0) 00354 { 00355 m_cachedIndex = -1; 00356 } 00357 else 00358 { 00359 m_cachedIndex = llElements[m_currentEl].index; 00360 m_cachedValue = llElements[m_currentEl].value; 00361 } 00362 } 00363 return *this; 00364 } 00365 00366 protected: 00367 const AmbiVector& m_vector; // the target vector 00368 int m_currentEl; // the current element in sparse/linked-list mode 00369 RealScalar m_epsilon; // epsilon used to prune zero coefficients 00370 int m_cachedIndex; // current coordinate 00371 Scalar m_cachedValue; // current value 00372 bool m_isDense; // mode of the vector 00373 00374 private: 00375 Iterator& operator=(const Iterator&); 00376 }; 00377 00378 00379 #endif // EIGEN_AMBIVECTOR_H