1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 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_DYNAMIC_SPARSEMATRIX_H 11 #define EIGEN_DYNAMIC_SPARSEMATRIX_H 12 13 namespace Eigen { 14 15 /** \deprecated use a SparseMatrix in an uncompressed mode 16 * 17 * \class DynamicSparseMatrix 18 * 19 * \brief A sparse matrix class designed for matrix assembly purpose 20 * 21 * \param _Scalar the scalar type, i.e. the type of the coefficients 22 * 23 * Unlike SparseMatrix, this class provides a much higher degree of flexibility. In particular, it allows 24 * random read/write accesses in log(rho*outer_size) where \c rho is the probability that a coefficient is 25 * nonzero and outer_size is the number of columns if the matrix is column-major and the number of rows 26 * otherwise. 27 * 28 * Internally, the data are stored as a std::vector of compressed vector. The performances of random writes might 29 * decrease as the number of nonzeros per inner-vector increase. In practice, we observed very good performance 30 * till about 100 nonzeros/vector, and the performance remains relatively good till 500 nonzeros/vectors. 31 * 32 * \see SparseMatrix 33 */ 34 35 namespace internal { 36 template<typename _Scalar, int _Options, typename _StorageIndex> 37 struct traits<DynamicSparseMatrix<_Scalar, _Options, _StorageIndex> > 38 { 39 typedef _Scalar Scalar; 40 typedef _StorageIndex StorageIndex; 41 typedef Sparse StorageKind; 42 typedef MatrixXpr XprKind; 43 enum { 44 RowsAtCompileTime = Dynamic, 45 ColsAtCompileTime = Dynamic, 46 MaxRowsAtCompileTime = Dynamic, 47 MaxColsAtCompileTime = Dynamic, 48 Flags = _Options | NestByRefBit | LvalueBit, 49 CoeffReadCost = NumTraits<Scalar>::ReadCost, 50 SupportedAccessPatterns = OuterRandomAccessPattern 51 }; 52 }; 53 } 54 55 template<typename _Scalar, int _Options, typename _StorageIndex> 56 class DynamicSparseMatrix 57 : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Options, _StorageIndex> > 58 { 59 typedef SparseMatrixBase<DynamicSparseMatrix> Base; 60 using Base::convert_index; 61 public: 62 EIGEN_SPARSE_PUBLIC_INTERFACE(DynamicSparseMatrix) 63 // FIXME: why are these operator already alvailable ??? 64 // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, +=) 65 // EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(DynamicSparseMatrix, -=) 66 typedef MappedSparseMatrix<Scalar,Flags> Map; 67 using Base::IsRowMajor; 68 using Base::operator=; 69 enum { 70 Options = _Options 71 }; 72 73 protected: 74 75 typedef DynamicSparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0), StorageIndex> TransposedSparseMatrix; 76 77 Index m_innerSize; 78 std::vector<internal::CompressedStorage<Scalar,StorageIndex> > m_data; 79 80 public: 81 82 inline Index rows() const { return IsRowMajor ? outerSize() : m_innerSize; } 83 inline Index cols() const { return IsRowMajor ? m_innerSize : outerSize(); } 84 inline Index innerSize() const { return m_innerSize; } 85 inline Index outerSize() const { return convert_index(m_data.size()); } 86 inline Index innerNonZeros(Index j) const { return m_data[j].size(); } 87 88 std::vector<internal::CompressedStorage<Scalar,StorageIndex> >& _data() { return m_data; } 89 const std::vector<internal::CompressedStorage<Scalar,StorageIndex> >& _data() const { return m_data; } 90 91 /** \returns the coefficient value at given position \a row, \a col 92 * This operation involes a log(rho*outer_size) binary search. 93 */ 94 inline Scalar coeff(Index row, Index col) const 95 { 96 const Index outer = IsRowMajor ? row : col; 97 const Index inner = IsRowMajor ? col : row; 98 return m_data[outer].at(inner); 99 } 100 101 /** \returns a reference to the coefficient value at given position \a row, \a col 102 * This operation involes a log(rho*outer_size) binary search. If the coefficient does not 103 * exist yet, then a sorted insertion into a sequential buffer is performed. 104 */ 105 inline Scalar& coeffRef(Index row, Index col) 106 { 107 const Index outer = IsRowMajor ? row : col; 108 const Index inner = IsRowMajor ? col : row; 109 return m_data[outer].atWithInsertion(inner); 110 } 111 112 class InnerIterator; 113 class ReverseInnerIterator; 114 115 void setZero() 116 { 117 for (Index j=0; j<outerSize(); ++j) 118 m_data[j].clear(); 119 } 120 121 /** \returns the number of non zero coefficients */ 122 Index nonZeros() const 123 { 124 Index res = 0; 125 for (Index j=0; j<outerSize(); ++j) 126 res += m_data[j].size(); 127 return res; 128 } 129 130 131 132 void reserve(Index reserveSize = 1000) 133 { 134 if (outerSize()>0) 135 { 136 Index reserveSizePerVector = (std::max)(reserveSize/outerSize(),Index(4)); 137 for (Index j=0; j<outerSize(); ++j) 138 { 139 m_data[j].reserve(reserveSizePerVector); 140 } 141 } 142 } 143 144 /** Does nothing: provided for compatibility with SparseMatrix */ 145 inline void startVec(Index /*outer*/) {} 146 147 /** \returns a reference to the non zero coefficient at position \a row, \a col assuming that: 148 * - the nonzero does not already exist 149 * - the new coefficient is the last one of the given inner vector. 150 * 151 * \sa insert, insertBackByOuterInner */ 152 inline Scalar& insertBack(Index row, Index col) 153 { 154 return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row); 155 } 156 157 /** \sa insertBack */ 158 inline Scalar& insertBackByOuterInner(Index outer, Index inner) 159 { 160 eigen_assert(outer<Index(m_data.size()) && inner<m_innerSize && "out of range"); 161 eigen_assert(((m_data[outer].size()==0) || (m_data[outer].index(m_data[outer].size()-1)<inner)) 162 && "wrong sorted insertion"); 163 m_data[outer].append(0, inner); 164 return m_data[outer].value(m_data[outer].size()-1); 165 } 166 167 inline Scalar& insert(Index row, Index col) 168 { 169 const Index outer = IsRowMajor ? row : col; 170 const Index inner = IsRowMajor ? col : row; 171 172 Index startId = 0; 173 Index id = static_cast<Index>(m_data[outer].size()) - 1; 174 m_data[outer].resize(id+2,1); 175 176 while ( (id >= startId) && (m_data[outer].index(id) > inner) ) 177 { 178 m_data[outer].index(id+1) = m_data[outer].index(id); 179 m_data[outer].value(id+1) = m_data[outer].value(id); 180 --id; 181 } 182 m_data[outer].index(id+1) = inner; 183 m_data[outer].value(id+1) = 0; 184 return m_data[outer].value(id+1); 185 } 186 187 /** Does nothing: provided for compatibility with SparseMatrix */ 188 inline void finalize() {} 189 190 /** Suppress all nonzeros which are smaller than \a reference under the tolerence \a epsilon */ 191 void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision()) 192 { 193 for (Index j=0; j<outerSize(); ++j) 194 m_data[j].prune(reference,epsilon); 195 } 196 197 /** Resize the matrix without preserving the data (the matrix is set to zero) 198 */ 199 void resize(Index rows, Index cols) 200 { 201 const Index outerSize = IsRowMajor ? rows : cols; 202 m_innerSize = convert_index(IsRowMajor ? cols : rows); 203 setZero(); 204 if (Index(m_data.size()) != outerSize) 205 { 206 m_data.resize(outerSize); 207 } 208 } 209 210 void resizeAndKeepData(Index rows, Index cols) 211 { 212 const Index outerSize = IsRowMajor ? rows : cols; 213 const Index innerSize = IsRowMajor ? cols : rows; 214 if (m_innerSize>innerSize) 215 { 216 // remove all coefficients with innerCoord>=innerSize 217 // TODO 218 //std::cerr << "not implemented yet\n"; 219 exit(2); 220 } 221 if (m_data.size() != outerSize) 222 { 223 m_data.resize(outerSize); 224 } 225 } 226 227 /** The class DynamicSparseMatrix is deprectaed */ 228 EIGEN_DEPRECATED inline DynamicSparseMatrix() 229 : m_innerSize(0), m_data(0) 230 { 231 eigen_assert(innerSize()==0 && outerSize()==0); 232 } 233 234 /** The class DynamicSparseMatrix is deprectaed */ 235 EIGEN_DEPRECATED inline DynamicSparseMatrix(Index rows, Index cols) 236 : m_innerSize(0) 237 { 238 resize(rows, cols); 239 } 240 241 /** The class DynamicSparseMatrix is deprectaed */ 242 template<typename OtherDerived> 243 EIGEN_DEPRECATED explicit inline DynamicSparseMatrix(const SparseMatrixBase<OtherDerived>& other) 244 : m_innerSize(0) 245 { 246 Base::operator=(other.derived()); 247 } 248 249 inline DynamicSparseMatrix(const DynamicSparseMatrix& other) 250 : Base(), m_innerSize(0) 251 { 252 *this = other.derived(); 253 } 254 255 inline void swap(DynamicSparseMatrix& other) 256 { 257 //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n"); 258 std::swap(m_innerSize, other.m_innerSize); 259 //std::swap(m_outerSize, other.m_outerSize); 260 m_data.swap(other.m_data); 261 } 262 263 inline DynamicSparseMatrix& operator=(const DynamicSparseMatrix& other) 264 { 265 if (other.isRValue()) 266 { 267 swap(other.const_cast_derived()); 268 } 269 else 270 { 271 resize(other.rows(), other.cols()); 272 m_data = other.m_data; 273 } 274 return *this; 275 } 276 277 /** Destructor */ 278 inline ~DynamicSparseMatrix() {} 279 280 public: 281 282 /** \deprecated 283 * Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */ 284 EIGEN_DEPRECATED void startFill(Index reserveSize = 1000) 285 { 286 setZero(); 287 reserve(reserveSize); 288 } 289 290 /** \deprecated use insert() 291 * inserts a nonzero coefficient at given coordinates \a row, \a col and returns its reference assuming that: 292 * 1 - the coefficient does not exist yet 293 * 2 - this the coefficient with greater inner coordinate for the given outer coordinate. 294 * In other words, assuming \c *this is column-major, then there must not exists any nonzero coefficient of coordinates 295 * \c i \c x \a col such that \c i >= \a row. Otherwise the matrix is invalid. 296 * 297 * \see fillrand(), coeffRef() 298 */ 299 EIGEN_DEPRECATED Scalar& fill(Index row, Index col) 300 { 301 const Index outer = IsRowMajor ? row : col; 302 const Index inner = IsRowMajor ? col : row; 303 return insertBack(outer,inner); 304 } 305 306 /** \deprecated use insert() 307 * Like fill() but with random inner coordinates. 308 * Compared to the generic coeffRef(), the unique limitation is that we assume 309 * the coefficient does not exist yet. 310 */ 311 EIGEN_DEPRECATED Scalar& fillrand(Index row, Index col) 312 { 313 return insert(row,col); 314 } 315 316 /** \deprecated use finalize() 317 * Does nothing. Provided for compatibility with SparseMatrix. */ 318 EIGEN_DEPRECATED void endFill() {} 319 320 # ifdef EIGEN_DYNAMICSPARSEMATRIX_PLUGIN 321 # include EIGEN_DYNAMICSPARSEMATRIX_PLUGIN 322 # endif 323 }; 324 325 template<typename Scalar, int _Options, typename _StorageIndex> 326 class DynamicSparseMatrix<Scalar,_Options,_StorageIndex>::InnerIterator : public SparseVector<Scalar,_Options,_StorageIndex>::InnerIterator 327 { 328 typedef typename SparseVector<Scalar,_Options,_StorageIndex>::InnerIterator Base; 329 public: 330 InnerIterator(const DynamicSparseMatrix& mat, Index outer) 331 : Base(mat.m_data[outer]), m_outer(outer) 332 {} 333 334 inline Index row() const { return IsRowMajor ? m_outer : Base::index(); } 335 inline Index col() const { return IsRowMajor ? Base::index() : m_outer; } 336 inline Index outer() const { return m_outer; } 337 338 protected: 339 const Index m_outer; 340 }; 341 342 template<typename Scalar, int _Options, typename _StorageIndex> 343 class DynamicSparseMatrix<Scalar,_Options,_StorageIndex>::ReverseInnerIterator : public SparseVector<Scalar,_Options,_StorageIndex>::ReverseInnerIterator 344 { 345 typedef typename SparseVector<Scalar,_Options,_StorageIndex>::ReverseInnerIterator Base; 346 public: 347 ReverseInnerIterator(const DynamicSparseMatrix& mat, Index outer) 348 : Base(mat.m_data[outer]), m_outer(outer) 349 {} 350 351 inline Index row() const { return IsRowMajor ? m_outer : Base::index(); } 352 inline Index col() const { return IsRowMajor ? Base::index() : m_outer; } 353 inline Index outer() const { return m_outer; } 354 355 protected: 356 const Index m_outer; 357 }; 358 359 namespace internal { 360 361 template<typename _Scalar, int _Options, typename _StorageIndex> 362 struct evaluator<DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> > 363 : evaluator_base<DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> > 364 { 365 typedef _Scalar Scalar; 366 typedef DynamicSparseMatrix<_Scalar,_Options,_StorageIndex> SparseMatrixType; 367 typedef typename SparseMatrixType::InnerIterator InnerIterator; 368 typedef typename SparseMatrixType::ReverseInnerIterator ReverseInnerIterator; 369 370 enum { 371 CoeffReadCost = NumTraits<_Scalar>::ReadCost, 372 Flags = SparseMatrixType::Flags 373 }; 374 375 evaluator() : m_matrix(0) {} 376 evaluator(const SparseMatrixType &mat) : m_matrix(&mat) {} 377 378 operator SparseMatrixType&() { return m_matrix->const_cast_derived(); } 379 operator const SparseMatrixType&() const { return *m_matrix; } 380 381 Scalar coeff(Index row, Index col) const { return m_matrix->coeff(row,col); } 382 383 Index nonZerosEstimate() const { return m_matrix->nonZeros(); } 384 385 const SparseMatrixType *m_matrix; 386 }; 387 388 } 389 390 } // end namespace Eigen 391 392 #endif // EIGEN_DYNAMIC_SPARSEMATRIX_H 393