Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
SparseMatrix.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2014 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_SPARSEMATRIX_H
11 #define EIGEN_SPARSEMATRIX_H
12 
13 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
47 namespace internal {
48 template<typename Scalar_, int Options_, typename StorageIndex_>
49 struct traits<SparseMatrix<Scalar_, Options_, StorageIndex_> >
50 {
51  typedef Scalar_ Scalar;
52  typedef StorageIndex_ StorageIndex;
53  typedef Sparse StorageKind;
54  typedef MatrixXpr XprKind;
55  enum {
56  RowsAtCompileTime = Dynamic,
57  ColsAtCompileTime = Dynamic,
58  MaxRowsAtCompileTime = Dynamic,
59  MaxColsAtCompileTime = Dynamic,
60  Flags = Options_ | NestByRefBit | LvalueBit | CompressedAccessBit,
61  SupportedAccessPatterns = InnerRandomAccessPattern
62  };
63 };
64 
65 template<typename Scalar_, int Options_, typename StorageIndex_, int DiagIndex>
66 struct traits<Diagonal<SparseMatrix<Scalar_, Options_, StorageIndex_>, DiagIndex> >
67 {
68  typedef SparseMatrix<Scalar_, Options_, StorageIndex_> MatrixType;
69  typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
70  typedef std::remove_reference_t<MatrixTypeNested> MatrixTypeNested_;
71 
72  typedef Scalar_ Scalar;
73  typedef Dense StorageKind;
74  typedef StorageIndex_ StorageIndex;
75  typedef MatrixXpr XprKind;
76 
77  enum {
78  RowsAtCompileTime = Dynamic,
79  ColsAtCompileTime = 1,
80  MaxRowsAtCompileTime = Dynamic,
81  MaxColsAtCompileTime = 1,
82  Flags = LvalueBit
83  };
84 };
85 
86 template<typename Scalar_, int Options_, typename StorageIndex_, int DiagIndex>
87 struct traits<Diagonal<const SparseMatrix<Scalar_, Options_, StorageIndex_>, DiagIndex> >
88  : public traits<Diagonal<SparseMatrix<Scalar_, Options_, StorageIndex_>, DiagIndex> >
89 {
90  enum {
91  Flags = 0
92  };
93 };
94 
95 } // end namespace internal
96 
97 template<typename Scalar_, int Options_, typename StorageIndex_>
99  : public SparseCompressedBase<SparseMatrix<Scalar_, Options_, StorageIndex_> >
100 {
102  using Base::convert_index;
103  friend class SparseVector<Scalar_,0,StorageIndex_>;
104  template<typename, typename, typename, typename, typename>
105  friend struct internal::Assignment;
106  public:
107  using Base::isCompressed;
108  using Base::nonZeros;
109  EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix)
110  using Base::operator+=;
111  using Base::operator-=;
112 
116  typedef typename Base::InnerIterator InnerIterator;
117  typedef typename Base::ReverseInnerIterator ReverseInnerIterator;
118 
119 
120  using Base::IsRowMajor;
121  typedef internal::CompressedStorage<Scalar,StorageIndex> Storage;
122  enum {
123  Options = Options_
124  };
125 
126  typedef typename Base::IndexVector IndexVector;
127  typedef typename Base::ScalarVector ScalarVector;
128  protected:
130 
131  Index m_outerSize;
132  Index m_innerSize;
133  StorageIndex* m_outerIndex;
134  StorageIndex* m_innerNonZeros; // optional, if null then the data is compressed
135  Storage m_data;
136 
137  public:
138 
140  inline Index rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
142  inline Index cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }
143 
145  inline Index innerSize() const { return m_innerSize; }
147  inline Index outerSize() const { return m_outerSize; }
148 
152  inline const Scalar* valuePtr() const { return m_data.valuePtr(); }
156  inline Scalar* valuePtr() { return m_data.valuePtr(); }
157 
161  inline const StorageIndex* innerIndexPtr() const { return m_data.indexPtr(); }
165  inline StorageIndex* innerIndexPtr() { return m_data.indexPtr(); }
166 
170  inline const StorageIndex* outerIndexPtr() const { return m_outerIndex; }
174  inline StorageIndex* outerIndexPtr() { return m_outerIndex; }
175 
179  inline const StorageIndex* innerNonZeroPtr() const { return m_innerNonZeros; }
183  inline StorageIndex* innerNonZeroPtr() { return m_innerNonZeros; }
184 
186  inline Storage& data() { return m_data; }
188  inline const Storage& data() const { return m_data; }
189 
192  inline Scalar coeff(Index row, Index col) const
193  {
194  eigen_assert(row>=0 && row<rows() && col>=0 && col<cols());
195 
196  const Index outer = IsRowMajor ? row : col;
197  const Index inner = IsRowMajor ? col : row;
198  Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer+1];
199  return m_data.atInRange(m_outerIndex[outer], end, StorageIndex(inner));
200  }
201 
210  inline Scalar& coeffRef(Index row, Index col)
211  {
212  eigen_assert(row>=0 && row<rows() && col>=0 && col<cols());
213 
214  const Index outer = IsRowMajor ? row : col;
215  const Index inner = IsRowMajor ? col : row;
216 
217  Index start = m_outerIndex[outer];
218  Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer+1];
219  eigen_assert(end>=start && "you probably called coeffRef on a non finalized matrix");
220  if(end<=start)
221  return insert(row,col);
222  const Index p = m_data.searchLowerIndex(start,end-1,StorageIndex(inner));
223  if((p<end) && (m_data.index(p)==inner))
224  return m_data.value(p);
225  else
226  return insert(row,col);
227  }
228 
244  Scalar& insert(Index row, Index col);
245 
246  public:
247 
255  inline void setZero()
256  {
257  m_data.clear();
258  std::fill_n(m_outerIndex, m_outerSize + 1, StorageIndex(0));
259  if(m_innerNonZeros) {
260  std::fill_n(m_innerNonZeros, m_outerSize, StorageIndex(0));
261  }
262  }
263 
267  inline void reserve(Index reserveSize)
268  {
269  eigen_assert(isCompressed() && "This function does not make sense in non compressed mode.");
270  m_data.reserve(reserveSize);
271  }
272 
273  #ifdef EIGEN_PARSED_BY_DOXYGEN
286  template<class SizesType>
287  inline void reserve(const SizesType& reserveSizes);
288  #else
289  template<class SizesType>
290  inline void reserve(const SizesType& reserveSizes, const typename SizesType::value_type& enableif =
291  typename SizesType::value_type())
292  {
293  EIGEN_UNUSED_VARIABLE(enableif);
294  reserveInnerVectors(reserveSizes);
295  }
296  #endif // EIGEN_PARSED_BY_DOXYGEN
297  protected:
298  template<class SizesType>
299  inline void reserveInnerVectors(const SizesType& reserveSizes)
300  {
301  if(isCompressed())
302  {
303  Index totalReserveSize = 0;
304  // turn the matrix into non-compressed mode
305  m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
306  if (!m_innerNonZeros) internal::throw_std_bad_alloc();
307 
308  // temporarily use m_innerSizes to hold the new starting points.
309  StorageIndex* newOuterIndex = m_innerNonZeros;
310 
311  StorageIndex count = 0;
312  for(Index j=0; j<m_outerSize; ++j)
313  {
314  newOuterIndex[j] = count;
315  count += reserveSizes[j] + (m_outerIndex[j+1]-m_outerIndex[j]);
316  totalReserveSize += reserveSizes[j];
317  }
318  m_data.reserve(totalReserveSize);
319  StorageIndex previousOuterIndex = m_outerIndex[m_outerSize];
320  for(Index j=m_outerSize-1; j>=0; --j)
321  {
322  StorageIndex innerNNZ = previousOuterIndex - m_outerIndex[j];
323  for(Index i=innerNNZ-1; i>=0; --i)
324  {
325  m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i);
326  m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i);
327  }
328  previousOuterIndex = m_outerIndex[j];
329  m_outerIndex[j] = newOuterIndex[j];
330  m_innerNonZeros[j] = innerNNZ;
331  }
332  if(m_outerSize>0)
333  m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize-1] + m_innerNonZeros[m_outerSize-1] + reserveSizes[m_outerSize-1];
334 
335  m_data.resize(m_outerIndex[m_outerSize]);
336  }
337  else
338  {
339  StorageIndex* newOuterIndex = static_cast<StorageIndex*>(std::malloc((m_outerSize+1)*sizeof(StorageIndex)));
340  if (!newOuterIndex) internal::throw_std_bad_alloc();
341 
342  StorageIndex count = 0;
343  for(Index j=0; j<m_outerSize; ++j)
344  {
345  newOuterIndex[j] = count;
346  StorageIndex alreadyReserved = (m_outerIndex[j+1]-m_outerIndex[j]) - m_innerNonZeros[j];
347  StorageIndex toReserve = std::max<StorageIndex>(reserveSizes[j], alreadyReserved);
348  count += toReserve + m_innerNonZeros[j];
349  }
350  newOuterIndex[m_outerSize] = count;
351 
352  m_data.resize(count);
353  for(Index j=m_outerSize-1; j>=0; --j)
354  {
355  Index offset = newOuterIndex[j] - m_outerIndex[j];
356  if(offset>0)
357  {
358  StorageIndex innerNNZ = m_innerNonZeros[j];
359  for(Index i=innerNNZ-1; i>=0; --i)
360  {
361  m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i);
362  m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i);
363  }
364  }
365  }
366 
367  std::swap(m_outerIndex, newOuterIndex);
368  std::free(newOuterIndex);
369  }
370 
371  }
372  public:
373 
374  //--- low level purely coherent filling ---
375 
386  inline Scalar& insertBack(Index row, Index col)
387  {
388  return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row);
389  }
390 
393  inline Scalar& insertBackByOuterInner(Index outer, Index inner)
394  {
395  eigen_assert(Index(m_outerIndex[outer+1]) == m_data.size() && "Invalid ordered insertion (invalid outer index)");
396  eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)<inner) && "Invalid ordered insertion (invalid inner index)");
397  Index p = m_outerIndex[outer+1];
398  ++m_outerIndex[outer+1];
399  m_data.append(Scalar(0), inner);
400  return m_data.value(p);
401  }
402 
405  inline Scalar& insertBackByOuterInnerUnordered(Index outer, Index inner)
406  {
407  Index p = m_outerIndex[outer+1];
408  ++m_outerIndex[outer+1];
409  m_data.append(Scalar(0), inner);
410  return m_data.value(p);
411  }
412 
415  inline void startVec(Index outer)
416  {
417  eigen_assert(m_outerIndex[outer]==Index(m_data.size()) && "You must call startVec for each inner vector sequentially");
418  eigen_assert(m_outerIndex[outer+1]==0 && "You must call startVec for each inner vector sequentially");
419  m_outerIndex[outer+1] = m_outerIndex[outer];
420  }
421 
425  inline void finalize()
426  {
427  if(isCompressed())
428  {
429  StorageIndex size = internal::convert_index<StorageIndex>(m_data.size());
430  Index i = m_outerSize;
431  // find the last filled column
432  while (i>=0 && m_outerIndex[i]==0)
433  --i;
434  ++i;
435  while (i<=m_outerSize)
436  {
437  m_outerIndex[i] = size;
438  ++i;
439  }
440  }
441  }
442 
443  //---
444 
445  template<typename InputIterators>
446  void setFromTriplets(const InputIterators& begin, const InputIterators& end);
447 
448  template<typename InputIterators,typename DupFunctor>
449  void setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func);
450 
451  void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op<Scalar,Scalar>()); }
452 
453  template<typename DupFunctor>
454  void collapseDuplicates(DupFunctor dup_func = DupFunctor());
455 
456  //---
457 
460  Scalar& insertByOuterInner(Index j, Index i)
461  {
462  return insert(IsRowMajor ? j : i, IsRowMajor ? i : j);
463  }
464 
468  {
469  if(isCompressed())
470  return;
471 
472  eigen_internal_assert(m_outerIndex!=0 && m_outerSize>0);
473 
474  Index oldStart = m_outerIndex[1];
475  m_outerIndex[1] = m_innerNonZeros[0];
476  for(Index j=1; j<m_outerSize; ++j)
477  {
478  Index nextOldStart = m_outerIndex[j+1];
479  Index offset = oldStart - m_outerIndex[j];
480  if(offset>0)
481  {
482  for(Index k=0; k<m_innerNonZeros[j]; ++k)
483  {
484  m_data.index(m_outerIndex[j]+k) = m_data.index(oldStart+k);
485  m_data.value(m_outerIndex[j]+k) = m_data.value(oldStart+k);
486  }
487  }
488  m_outerIndex[j+1] = m_outerIndex[j] + m_innerNonZeros[j];
489  oldStart = nextOldStart;
490  }
491  std::free(m_innerNonZeros);
492  m_innerNonZeros = 0;
493  m_data.resize(m_outerIndex[m_outerSize]);
494  m_data.squeeze();
495  }
496 
498  void uncompress()
499  {
500  if(m_innerNonZeros != 0)
501  return;
502  m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
503  for (Index i = 0; i < m_outerSize; i++)
504  {
505  m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i];
506  }
507  }
508 
510  void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits<RealScalar>::dummy_precision())
511  {
512  prune(default_prunning_func(reference,epsilon));
513  }
514 
522  template<typename KeepFunc>
523  void prune(const KeepFunc& keep = KeepFunc())
524  {
525  // TODO optimize the uncompressed mode to avoid moving and allocating the data twice
526  makeCompressed();
527 
528  StorageIndex k = 0;
529  for(Index j=0; j<m_outerSize; ++j)
530  {
531  Index previousStart = m_outerIndex[j];
532  m_outerIndex[j] = k;
533  Index end = m_outerIndex[j+1];
534  for(Index i=previousStart; i<end; ++i)
535  {
536  if(keep(IsRowMajor?j:m_data.index(i), IsRowMajor?m_data.index(i):j, m_data.value(i)))
537  {
538  m_data.value(k) = m_data.value(i);
539  m_data.index(k) = m_data.index(i);
540  ++k;
541  }
542  }
543  }
544  m_outerIndex[m_outerSize] = k;
545  m_data.resize(k,0);
546  }
547 
557  {
558  // No change
559  if (this->rows() == rows && this->cols() == cols) return;
560 
561  // If one dimension is null, then there is nothing to be preserved
562  if(rows==0 || cols==0) return resize(rows,cols);
563 
564  Index innerChange = IsRowMajor ? cols - this->cols() : rows - this->rows();
565  Index outerChange = IsRowMajor ? rows - this->rows() : cols - this->cols();
566  StorageIndex newInnerSize = convert_index(IsRowMajor ? cols : rows);
567 
568  // Deals with inner non zeros
569  if (m_innerNonZeros)
570  {
571  // Resize m_innerNonZeros
572  StorageIndex *newInnerNonZeros = static_cast<StorageIndex*>(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) * sizeof(StorageIndex)));
573  if (!newInnerNonZeros) internal::throw_std_bad_alloc();
574  m_innerNonZeros = newInnerNonZeros;
575 
576  for(Index i=m_outerSize; i<m_outerSize+outerChange; i++)
577  m_innerNonZeros[i] = 0;
578  }
579  else if (innerChange < 0)
580  {
581  // Inner size decreased: allocate a new m_innerNonZeros
582  m_innerNonZeros = static_cast<StorageIndex*>(std::malloc((m_outerSize + outerChange) * sizeof(StorageIndex)));
583  if (!m_innerNonZeros) internal::throw_std_bad_alloc();
584  for(Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++)
585  m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i];
586  for(Index i = m_outerSize; i < m_outerSize + outerChange; i++)
587  m_innerNonZeros[i] = 0;
588  }
589 
590  // Change the m_innerNonZeros in case of a decrease of inner size
591  if (m_innerNonZeros && innerChange < 0)
592  {
593  for(Index i = 0; i < m_outerSize + (std::min)(outerChange, Index(0)); i++)
594  {
595  StorageIndex &n = m_innerNonZeros[i];
596  StorageIndex start = m_outerIndex[i];
597  while (n > 0 && m_data.index(start+n-1) >= newInnerSize) --n;
598  }
599  }
600 
601  m_innerSize = newInnerSize;
602 
603  // Re-allocate outer index structure if necessary
604  if (outerChange == 0)
605  return;
606 
607  StorageIndex *newOuterIndex = static_cast<StorageIndex*>(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) * sizeof(StorageIndex)));
608  if (!newOuterIndex) internal::throw_std_bad_alloc();
609  m_outerIndex = newOuterIndex;
610  if (outerChange > 0)
611  {
612  StorageIndex lastIdx = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize];
613  for(Index i=m_outerSize; i<m_outerSize+outerChange+1; i++)
614  m_outerIndex[i] = lastIdx;
615  }
616  m_outerSize += outerChange;
617  }
618 
627  {
628  const Index outerSize = IsRowMajor ? rows : cols;
629  m_innerSize = IsRowMajor ? cols : rows;
630  m_data.clear();
631  if (m_outerSize != outerSize || m_outerSize==0)
632  {
633  std::free(m_outerIndex);
634  m_outerIndex = static_cast<StorageIndex*>(std::malloc((outerSize + 1) * sizeof(StorageIndex)));
635  if (!m_outerIndex) internal::throw_std_bad_alloc();
636 
637  m_outerSize = outerSize;
638  }
639  if(m_innerNonZeros)
640  {
641  std::free(m_innerNonZeros);
642  m_innerNonZeros = 0;
643  }
644  std::fill_n(m_outerIndex, m_outerSize + 1, StorageIndex(0));
645  }
646 
649  void resizeNonZeros(Index size)
650  {
651  m_data.resize(size);
652  }
653 
656 
662 
664  inline SparseMatrix()
665  : m_outerSize(-1), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
666  {
667  resize(0, 0);
668  }
669 
672  : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
673  {
674  resize(rows, cols);
675  }
676 
678  template<typename OtherDerived>
680  : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
681  {
682  EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
683  YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
684  const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator<OtherDerived>::Flags & RowMajorBit);
685  if (needToTranspose)
686  *this = other.derived();
687  else
688  {
689  #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
690  EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
691  #endif
692  internal::call_assignment_no_alias(*this, other.derived());
693  }
694  }
695 
697  template<typename OtherDerived, unsigned int UpLo>
699  : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
700  {
701  Base::operator=(other);
702  }
703 
705  inline SparseMatrix(const SparseMatrix& other)
706  : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
707  {
708  *this = other.derived();
709  }
710 
712  template<typename OtherDerived>
713  SparseMatrix(const ReturnByValue<OtherDerived>& other)
714  : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
715  {
716  initAssignment(other);
717  other.evalTo(*this);
718  }
719 
721  template<typename OtherDerived>
722  explicit SparseMatrix(const DiagonalBase<OtherDerived>& other)
723  : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
724  {
725  *this = other.derived();
726  }
727 
730  inline void swap(SparseMatrix& other)
731  {
732  //EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
733  std::swap(m_outerIndex, other.m_outerIndex);
734  std::swap(m_innerSize, other.m_innerSize);
735  std::swap(m_outerSize, other.m_outerSize);
736  std::swap(m_innerNonZeros, other.m_innerNonZeros);
737  m_data.swap(other.m_data);
738  }
739 
742  inline void setIdentity()
743  {
744  eigen_assert(rows() == cols() && "ONLY FOR SQUARED MATRICES");
745  this->m_data.resize(rows());
746  Eigen::Map<IndexVector>(this->m_data.indexPtr(), rows()).setLinSpaced(0, StorageIndex(rows()-1));
747  Eigen::Map<ScalarVector>(this->m_data.valuePtr(), rows()).setOnes();
748  Eigen::Map<IndexVector>(this->m_outerIndex, rows()+1).setLinSpaced(0, StorageIndex(rows()));
749  std::free(m_innerNonZeros);
750  m_innerNonZeros = 0;
751  }
752  inline SparseMatrix& operator=(const SparseMatrix& other)
753  {
754  if (other.isRValue())
755  {
756  swap(other.const_cast_derived());
757  }
758  else if(this!=&other)
759  {
760  #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
761  EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
762  #endif
763  initAssignment(other);
764  if(other.isCompressed())
765  {
766  internal::smart_copy(other.m_outerIndex, other.m_outerIndex + m_outerSize + 1, m_outerIndex);
767  m_data = other.m_data;
768  }
769  else
770  {
771  Base::operator=(other);
772  }
773  }
774  return *this;
775  }
776 
777 #ifndef EIGEN_PARSED_BY_DOXYGEN
778  template<typename OtherDerived>
779  inline SparseMatrix& operator=(const EigenBase<OtherDerived>& other)
780  { return Base::operator=(other.derived()); }
781 
782  template<typename Lhs, typename Rhs>
783  inline SparseMatrix& operator=(const Product<Lhs,Rhs,AliasFreeProduct>& other);
784 #endif // EIGEN_PARSED_BY_DOXYGEN
785 
786  template<typename OtherDerived>
787  EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other);
788 
789  friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m)
790  {
791  EIGEN_DBG_SPARSE(
792  s << "Nonzero entries:\n";
793  if(m.isCompressed())
794  {
795  for (Index i=0; i<m.nonZeros(); ++i)
796  s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
797  }
798  else
799  {
800  for (Index i=0; i<m.outerSize(); ++i)
801  {
802  Index p = m.m_outerIndex[i];
803  Index pe = m.m_outerIndex[i]+m.m_innerNonZeros[i];
804  Index k=p;
805  for (; k<pe; ++k) {
806  s << "(" << m.m_data.value(k) << "," << m.m_data.index(k) << ") ";
807  }
808  for (; k<m.m_outerIndex[i+1]; ++k) {
809  s << "(_,_) ";
810  }
811  }
812  }
813  s << std::endl;
814  s << std::endl;
815  s << "Outer pointers:\n";
816  for (Index i=0; i<m.outerSize(); ++i) {
817  s << m.m_outerIndex[i] << " ";
818  }
819  s << " $" << std::endl;
820  if(!m.isCompressed())
821  {
822  s << "Inner non zeros:\n";
823  for (Index i=0; i<m.outerSize(); ++i) {
824  s << m.m_innerNonZeros[i] << " ";
825  }
826  s << " $" << std::endl;
827  }
828  s << std::endl;
829  );
830  s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
831  return s;
832  }
833 
835  inline ~SparseMatrix()
836  {
837  std::free(m_outerIndex);
838  std::free(m_innerNonZeros);
839  }
840 
842  Scalar sum() const;
843 
844 # ifdef EIGEN_SPARSEMATRIX_PLUGIN
845 # include EIGEN_SPARSEMATRIX_PLUGIN
846 # endif
847 
848 protected:
849 
850  template<typename Other>
851  void initAssignment(const Other& other)
852  {
853  resize(other.rows(), other.cols());
854  if(m_innerNonZeros)
855  {
856  std::free(m_innerNonZeros);
857  m_innerNonZeros = 0;
858  }
859  }
860 
863  EIGEN_DONT_INLINE Scalar& insertCompressed(Index row, Index col);
864 
867  class SingletonVector
868  {
869  StorageIndex m_index;
870  StorageIndex m_value;
871  public:
872  typedef StorageIndex value_type;
873  SingletonVector(Index i, Index v)
874  : m_index(convert_index(i)), m_value(convert_index(v))
875  {}
876 
877  StorageIndex operator[](Index i) const { return i==m_index ? m_value : 0; }
878  };
879 
882  EIGEN_DONT_INLINE Scalar& insertUncompressed(Index row, Index col);
883 
884 public:
887  EIGEN_STRONG_INLINE Scalar& insertBackUncompressed(Index row, Index col)
888  {
889  const Index outer = IsRowMajor ? row : col;
890  const Index inner = IsRowMajor ? col : row;
891 
892  eigen_assert(!isCompressed());
893  eigen_assert(m_innerNonZeros[outer]<=(m_outerIndex[outer+1] - m_outerIndex[outer]));
894 
895  Index p = m_outerIndex[outer] + m_innerNonZeros[outer]++;
896  m_data.index(p) = convert_index(inner);
897  return (m_data.value(p) = Scalar(0));
898  }
899 protected:
900  struct IndexPosPair {
901  IndexPosPair(Index a_i, Index a_p) : i(a_i), p(a_p) {}
902  Index i;
903  Index p;
904  };
905 
919  template<typename DiagXpr, typename Func>
920  void assignDiagonal(const DiagXpr diagXpr, const Func& assignFunc)
921  {
922  Index n = diagXpr.size();
923 
924  const bool overwrite = internal::is_same<Func, internal::assign_op<Scalar,Scalar> >::value;
925  if(overwrite)
926  {
927  if((this->rows()!=n) || (this->cols()!=n))
928  this->resize(n, n);
929  }
930 
931  if(m_data.size()==0 || overwrite)
932  {
933  typedef Array<StorageIndex,Dynamic,1> ArrayXI;
934  this->makeCompressed();
935  this->resizeNonZeros(n);
936  Eigen::Map<ArrayXI>(this->innerIndexPtr(), n).setLinSpaced(0,StorageIndex(n)-1);
937  Eigen::Map<ArrayXI>(this->outerIndexPtr(), n+1).setLinSpaced(0,StorageIndex(n));
938  Eigen::Map<Array<Scalar,Dynamic,1> > values = this->coeffs();
939  values.setZero();
940  internal::call_assignment_no_alias(values, diagXpr, assignFunc);
941  }
942  else
943  {
944  bool isComp = isCompressed();
945  internal::evaluator<DiagXpr> diaEval(diagXpr);
946  std::vector<IndexPosPair> newEntries;
947 
948  // 1 - try in-place update and record insertion failures
949  for(Index i = 0; i<n; ++i)
950  {
951  internal::LowerBoundIndex lb = this->lower_bound(i,i);
952  Index p = lb.value;
953  if(lb.found)
954  {
955  // the coeff already exists
956  assignFunc.assignCoeff(m_data.value(p), diaEval.coeff(i));
957  }
958  else if((!isComp) && m_innerNonZeros[i] < (m_outerIndex[i+1]-m_outerIndex[i]))
959  {
960  // non compressed mode with local room for inserting one element
961  m_data.moveChunk(p, p+1, m_outerIndex[i]+m_innerNonZeros[i]-p);
962  m_innerNonZeros[i]++;
963  m_data.value(p) = Scalar(0);
964  m_data.index(p) = StorageIndex(i);
965  assignFunc.assignCoeff(m_data.value(p), diaEval.coeff(i));
966  }
967  else
968  {
969  // defer insertion
970  newEntries.push_back(IndexPosPair(i,p));
971  }
972  }
973  // 2 - insert deferred entries
974  Index n_entries = Index(newEntries.size());
975  if(n_entries>0)
976  {
977  Storage newData(m_data.size()+n_entries);
978  Index prev_p = 0;
979  Index prev_i = 0;
980  for(Index k=0; k<n_entries;++k)
981  {
982  Index i = newEntries[k].i;
983  Index p = newEntries[k].p;
984  internal::smart_copy(m_data.valuePtr()+prev_p, m_data.valuePtr()+p, newData.valuePtr()+prev_p+k);
985  internal::smart_copy(m_data.indexPtr()+prev_p, m_data.indexPtr()+p, newData.indexPtr()+prev_p+k);
986  for(Index j=prev_i;j<i;++j)
987  m_outerIndex[j+1] += k;
988  if(!isComp)
989  m_innerNonZeros[i]++;
990  prev_p = p;
991  prev_i = i;
992  newData.value(p+k) = Scalar(0);
993  newData.index(p+k) = StorageIndex(i);
994  assignFunc.assignCoeff(newData.value(p+k), diaEval.coeff(i));
995  }
996  {
997  internal::smart_copy(m_data.valuePtr()+prev_p, m_data.valuePtr()+m_data.size(), newData.valuePtr()+prev_p+n_entries);
998  internal::smart_copy(m_data.indexPtr()+prev_p, m_data.indexPtr()+m_data.size(), newData.indexPtr()+prev_p+n_entries);
999  for(Index j=prev_i+1;j<=m_outerSize;++j)
1000  m_outerIndex[j] += n_entries;
1001  }
1002  m_data.swap(newData);
1003  }
1004  }
1005  }
1006 
1007 private:
1008  EIGEN_STATIC_ASSERT(NumTraits<StorageIndex>::IsSigned,THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE)
1009  EIGEN_STATIC_ASSERT((Options&(ColMajor|RowMajor))==Options,INVALID_MATRIX_TEMPLATE_PARAMETERS)
1010 
1011  struct default_prunning_func {
1012  default_prunning_func(const Scalar& ref, const RealScalar& eps) : reference(ref), epsilon(eps) {}
1013  inline bool operator() (const Index&, const Index&, const Scalar& value) const
1014  {
1015  return !internal::isMuchSmallerThan(value, reference, epsilon);
1016  }
1017  Scalar reference;
1018  RealScalar epsilon;
1019  };
1020 };
1021 
1022 namespace internal {
1023 
1024 template<typename InputIterator, typename SparseMatrixType, typename DupFunctor>
1025 void set_from_triplets(const InputIterator& begin, const InputIterator& end, SparseMatrixType& mat, DupFunctor dup_func)
1026 {
1027  enum { IsRowMajor = SparseMatrixType::IsRowMajor };
1028  typedef typename SparseMatrixType::Scalar Scalar;
1029  typedef typename SparseMatrixType::StorageIndex StorageIndex;
1030  SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor,StorageIndex> trMat(mat.rows(),mat.cols());
1031 
1032  if(begin!=end)
1033  {
1034  // pass 1: count the nnz per inner-vector
1035  typename SparseMatrixType::IndexVector wi(trMat.outerSize());
1036  wi.setZero();
1037  for(InputIterator it(begin); it!=end; ++it)
1038  {
1039  eigen_assert(it->row()>=0 && it->row()<mat.rows() && it->col()>=0 && it->col()<mat.cols());
1040  wi(IsRowMajor ? it->col() : it->row())++;
1041  }
1042 
1043  // pass 2: insert all the elements into trMat
1044  trMat.reserve(wi);
1045  for(InputIterator it(begin); it!=end; ++it)
1046  trMat.insertBackUncompressed(it->row(),it->col()) = it->value();
1047 
1048  // pass 3:
1049  trMat.collapseDuplicates(dup_func);
1050  }
1051 
1052  // pass 4: transposed copy -> implicit sorting
1053  mat = trMat;
1054 }
1055 
1056 }
1057 
1058 
1096 template<typename Scalar, int Options_, typename StorageIndex_>
1097 template<typename InputIterators>
1098 void SparseMatrix<Scalar,Options_,StorageIndex_>::setFromTriplets(const InputIterators& begin, const InputIterators& end)
1099 {
1100  internal::set_from_triplets<InputIterators, SparseMatrix<Scalar,Options_,StorageIndex_> >(begin, end, *this, internal::scalar_sum_op<Scalar,Scalar>());
1101 }
1102 
1112 template<typename Scalar, int Options_, typename StorageIndex_>
1113 template<typename InputIterators,typename DupFunctor>
1114 void SparseMatrix<Scalar,Options_,StorageIndex_>::setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func)
1115 {
1116  internal::set_from_triplets<InputIterators, SparseMatrix<Scalar,Options_,StorageIndex_>, DupFunctor>(begin, end, *this, dup_func);
1117 }
1118 
1120 template<typename Scalar, int Options_, typename StorageIndex_>
1121 template<typename DupFunctor>
1123 {
1124  eigen_assert(!isCompressed());
1125  // TODO, in practice we should be able to use m_innerNonZeros for that task
1126  IndexVector wi(innerSize());
1127  wi.fill(-1);
1128  StorageIndex count = 0;
1129  // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers
1130  for(Index j=0; j<outerSize(); ++j)
1131  {
1132  StorageIndex start = count;
1133  Index oldEnd = m_outerIndex[j]+m_innerNonZeros[j];
1134  for(Index k=m_outerIndex[j]; k<oldEnd; ++k)
1135  {
1136  Index i = m_data.index(k);
1137  if(wi(i)>=start)
1138  {
1139  // we already meet this entry => accumulate it
1140  m_data.value(wi(i)) = dup_func(m_data.value(wi(i)), m_data.value(k));
1141  }
1142  else
1143  {
1144  m_data.value(count) = m_data.value(k);
1145  m_data.index(count) = m_data.index(k);
1146  wi(i) = count;
1147  ++count;
1148  }
1149  }
1150  m_outerIndex[j] = start;
1151  }
1152  m_outerIndex[m_outerSize] = count;
1153 
1154  // turn the matrix into compressed form
1155  std::free(m_innerNonZeros);
1156  m_innerNonZeros = 0;
1157  m_data.resize(m_outerIndex[m_outerSize]);
1158 }
1159 
1160 template<typename Scalar, int Options_, typename StorageIndex_>
1161 template<typename OtherDerived>
1162 EIGEN_DONT_INLINE SparseMatrix<Scalar,Options_,StorageIndex_>& SparseMatrix<Scalar,Options_,StorageIndex_>::operator=(const SparseMatrixBase<OtherDerived>& other)
1163 {
1164  EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
1165  YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
1166 
1167  #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
1168  EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
1169  #endif
1170 
1171  const bool needToTranspose = (Flags & RowMajorBit) != (internal::evaluator<OtherDerived>::Flags & RowMajorBit);
1172  if (needToTranspose)
1173  {
1174  #ifdef EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
1175  EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
1176  #endif
1177  // two passes algorithm:
1178  // 1 - compute the number of coeffs per dest inner vector
1179  // 2 - do the actual copy/eval
1180  // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed
1181  typedef typename internal::nested_eval<OtherDerived,2,typename internal::plain_matrix_type<OtherDerived>::type >::type OtherCopy;
1182  typedef internal::remove_all_t<OtherCopy> OtherCopy_;
1183  typedef internal::evaluator<OtherCopy_> OtherCopyEval;
1184  OtherCopy otherCopy(other.derived());
1185  OtherCopyEval otherCopyEval(otherCopy);
1186 
1187  SparseMatrix dest(other.rows(),other.cols());
1188  Eigen::Map<IndexVector> (dest.m_outerIndex,dest.outerSize()).setZero();
1189 
1190  // pass 1
1191  // FIXME the above copy could be merged with that pass
1192  for (Index j=0; j<otherCopy.outerSize(); ++j)
1193  for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it)
1194  ++dest.m_outerIndex[it.index()];
1195 
1196  // prefix sum
1197  StorageIndex count = 0;
1198  IndexVector positions(dest.outerSize());
1199  for (Index j=0; j<dest.outerSize(); ++j)
1200  {
1201  StorageIndex tmp = dest.m_outerIndex[j];
1202  dest.m_outerIndex[j] = count;
1203  positions[j] = count;
1204  count += tmp;
1205  }
1206  dest.m_outerIndex[dest.outerSize()] = count;
1207  // alloc
1208  dest.m_data.resize(count);
1209  // pass 2
1210  for (StorageIndex j=0; j<otherCopy.outerSize(); ++j)
1211  {
1212  for (typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it)
1213  {
1214  Index pos = positions[it.index()]++;
1215  dest.m_data.index(pos) = j;
1216  dest.m_data.value(pos) = it.value();
1217  }
1218  }
1219  this->swap(dest);
1220  return *this;
1221  }
1222  else
1223  {
1224  if(other.isRValue())
1225  {
1226  initAssignment(other.derived());
1227  }
1228  // there is no special optimization
1229  return Base::operator=(other.derived());
1230  }
1231 }
1232 
1233 template<typename Scalar_, int Options_, typename StorageIndex_>
1234 typename SparseMatrix<Scalar_,Options_,StorageIndex_>::Scalar& SparseMatrix<Scalar_,Options_,StorageIndex_>::insert(Index row, Index col)
1235 {
1236  eigen_assert(row>=0 && row<rows() && col>=0 && col<cols());
1237 
1238  const Index outer = IsRowMajor ? row : col;
1239  const Index inner = IsRowMajor ? col : row;
1240 
1241  if(isCompressed())
1242  {
1243  if(nonZeros()==0)
1244  {
1245  // reserve space if not already done
1246  if(m_data.allocatedSize()==0)
1247  m_data.reserve(2*m_innerSize);
1248 
1249  // turn the matrix into non-compressed mode
1250  m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
1251  if(!m_innerNonZeros) internal::throw_std_bad_alloc();
1252 
1253  std::fill(m_innerNonZeros, m_innerNonZeros + m_outerSize, StorageIndex(0));
1254 
1255  // pack all inner-vectors to the end of the pre-allocated space
1256  // and allocate the entire free-space to the first inner-vector
1257  StorageIndex end = convert_index(m_data.allocatedSize());
1258  for(Index j=1; j<=m_outerSize; ++j)
1259  m_outerIndex[j] = end;
1260  }
1261  else
1262  {
1263  // turn the matrix into non-compressed mode
1264  m_innerNonZeros = static_cast<StorageIndex*>(std::malloc(m_outerSize * sizeof(StorageIndex)));
1265  if(!m_innerNonZeros) internal::throw_std_bad_alloc();
1266  for(Index j=0; j<m_outerSize; ++j)
1267  m_innerNonZeros[j] = m_outerIndex[j+1]-m_outerIndex[j];
1268  }
1269  }
1270 
1271  // check whether we can do a fast "push back" insertion
1272  Index data_end = m_data.allocatedSize();
1273 
1274  // First case: we are filling a new inner vector which is packed at the end.
1275  // We assume that all remaining inner-vectors are also empty and packed to the end.
1276  if(m_outerIndex[outer]==data_end)
1277  {
1278  eigen_internal_assert(m_innerNonZeros[outer]==0);
1279 
1280  // pack previous empty inner-vectors to end of the used-space
1281  // and allocate the entire free-space to the current inner-vector.
1282  StorageIndex p = convert_index(m_data.size());
1283  Index j = outer;
1284  while(j>=0 && m_innerNonZeros[j]==0)
1285  m_outerIndex[j--] = p;
1286 
1287  // push back the new element
1288  ++m_innerNonZeros[outer];
1289  m_data.append(Scalar(0), inner);
1290 
1291  // check for reallocation
1292  if(data_end != m_data.allocatedSize())
1293  {
1294  // m_data has been reallocated
1295  // -> move remaining inner-vectors back to the end of the free-space
1296  // so that the entire free-space is allocated to the current inner-vector.
1297  eigen_internal_assert(data_end < m_data.allocatedSize());
1298  StorageIndex new_end = convert_index(m_data.allocatedSize());
1299  for(Index k=outer+1; k<=m_outerSize; ++k)
1300  if(m_outerIndex[k]==data_end)
1301  m_outerIndex[k] = new_end;
1302  }
1303  return m_data.value(p);
1304  }
1305 
1306  // Second case: the next inner-vector is packed to the end
1307  // and the current inner-vector end match the used-space.
1308  if(m_outerIndex[outer+1]==data_end && m_outerIndex[outer]+m_innerNonZeros[outer]==m_data.size())
1309  {
1310  eigen_internal_assert(outer+1==m_outerSize || m_innerNonZeros[outer+1]==0);
1311 
1312  // add space for the new element
1313  ++m_innerNonZeros[outer];
1314  m_data.resize(m_data.size()+1);
1315 
1316  // check for reallocation
1317  if(data_end != m_data.allocatedSize())
1318  {
1319  // m_data has been reallocated
1320  // -> move remaining inner-vectors back to the end of the free-space
1321  // so that the entire free-space is allocated to the current inner-vector.
1322  eigen_internal_assert(data_end < m_data.allocatedSize());
1323  StorageIndex new_end = convert_index(m_data.allocatedSize());
1324  for(Index k=outer+1; k<=m_outerSize; ++k)
1325  if(m_outerIndex[k]==data_end)
1326  m_outerIndex[k] = new_end;
1327  }
1328 
1329  // and insert it at the right position (sorted insertion)
1330  Index startId = m_outerIndex[outer];
1331  Index p = m_outerIndex[outer]+m_innerNonZeros[outer]-1;
1332  while ( (p > startId) && (m_data.index(p-1) > inner) )
1333  {
1334  m_data.index(p) = m_data.index(p-1);
1335  m_data.value(p) = m_data.value(p-1);
1336  --p;
1337  }
1338 
1339  m_data.index(p) = convert_index(inner);
1340  return (m_data.value(p) = Scalar(0));
1341  }
1342 
1343  if(m_data.size() != m_data.allocatedSize())
1344  {
1345  // make sure the matrix is compatible to random un-compressed insertion:
1346  m_data.resize(m_data.allocatedSize());
1347  this->reserveInnerVectors(Array<StorageIndex,Dynamic,1>::Constant(m_outerSize, 2));
1348  }
1349 
1350  return insertUncompressed(row,col);
1351 }
1352 
1353 template<typename Scalar_, int Options_, typename StorageIndex_>
1354 EIGEN_DONT_INLINE typename SparseMatrix<Scalar_,Options_,StorageIndex_>::Scalar& SparseMatrix<Scalar_,Options_,StorageIndex_>::insertUncompressed(Index row, Index col)
1355 {
1356  eigen_assert(!isCompressed());
1357 
1358  const Index outer = IsRowMajor ? row : col;
1359  const StorageIndex inner = convert_index(IsRowMajor ? col : row);
1360 
1361  Index room = m_outerIndex[outer+1] - m_outerIndex[outer];
1362  StorageIndex innerNNZ = m_innerNonZeros[outer];
1363  if(innerNNZ>=room)
1364  {
1365  // this inner vector is full, we need to reallocate the whole buffer :(
1366  reserve(SingletonVector(outer,std::max<StorageIndex>(2,innerNNZ)));
1367  }
1368 
1369  Index startId = m_outerIndex[outer];
1370  Index p = startId + m_innerNonZeros[outer];
1371  while ( (p > startId) && (m_data.index(p-1) > inner) )
1372  {
1373  m_data.index(p) = m_data.index(p-1);
1374  m_data.value(p) = m_data.value(p-1);
1375  --p;
1376  }
1377  eigen_assert((p<=startId || m_data.index(p-1)!=inner) && "you cannot insert an element that already exists, you must call coeffRef to this end");
1378 
1379  m_innerNonZeros[outer]++;
1380 
1381  m_data.index(p) = inner;
1382  return (m_data.value(p) = Scalar(0));
1383 }
1384 
1385 template<typename Scalar_, int Options_, typename StorageIndex_>
1386 EIGEN_DONT_INLINE typename SparseMatrix<Scalar_,Options_,StorageIndex_>::Scalar& SparseMatrix<Scalar_,Options_,StorageIndex_>::insertCompressed(Index row, Index col)
1387 {
1388  eigen_assert(isCompressed());
1389 
1390  const Index outer = IsRowMajor ? row : col;
1391  const Index inner = IsRowMajor ? col : row;
1392 
1393  Index previousOuter = outer;
1394  if (m_outerIndex[outer+1]==0)
1395  {
1396  // we start a new inner vector
1397  while (previousOuter>=0 && m_outerIndex[previousOuter]==0)
1398  {
1399  m_outerIndex[previousOuter] = convert_index(m_data.size());
1400  --previousOuter;
1401  }
1402  m_outerIndex[outer+1] = m_outerIndex[outer];
1403  }
1404 
1405  // here we have to handle the tricky case where the outerIndex array
1406  // starts with: [ 0 0 0 0 0 1 ...] and we are inserted in, e.g.,
1407  // the 2nd inner vector...
1408  bool isLastVec = (!(previousOuter==-1 && m_data.size()!=0))
1409  && (std::size_t(m_outerIndex[outer+1]) == m_data.size());
1410 
1411  std::size_t startId = m_outerIndex[outer];
1412  // FIXME let's make sure sizeof(long int) == sizeof(std::size_t)
1413  std::size_t p = m_outerIndex[outer+1];
1414  ++m_outerIndex[outer+1];
1415 
1416  double reallocRatio = 1;
1417  if (m_data.allocatedSize()<=m_data.size())
1418  {
1419  // if there is no preallocated memory, let's reserve a minimum of 32 elements
1420  if (m_data.size()==0)
1421  {
1422  m_data.reserve(32);
1423  }
1424  else
1425  {
1426  // we need to reallocate the data, to reduce multiple reallocations
1427  // we use a smart resize algorithm based on the current filling ratio
1428  // in addition, we use double to avoid integers overflows
1429  double nnzEstimate = double(m_outerIndex[outer])*double(m_outerSize)/double(outer+1);
1430  reallocRatio = (nnzEstimate-double(m_data.size()))/double(m_data.size());
1431  // furthermore we bound the realloc ratio to:
1432  // 1) reduce multiple minor realloc when the matrix is almost filled
1433  // 2) avoid to allocate too much memory when the matrix is almost empty
1434  reallocRatio = (std::min)((std::max)(reallocRatio,1.5),8.);
1435  }
1436  }
1437  m_data.resize(m_data.size()+1,reallocRatio);
1438 
1439  if (!isLastVec)
1440  {
1441  if (previousOuter==-1)
1442  {
1443  // oops wrong guess.
1444  // let's correct the outer offsets
1445  for (Index k=0; k<=(outer+1); ++k)
1446  m_outerIndex[k] = 0;
1447  Index k=outer+1;
1448  while(m_outerIndex[k]==0)
1449  m_outerIndex[k++] = 1;
1450  while (k<=m_outerSize && m_outerIndex[k]!=0)
1451  m_outerIndex[k++]++;
1452  p = 0;
1453  --k;
1454  k = m_outerIndex[k]-1;
1455  while (k>0)
1456  {
1457  m_data.index(k) = m_data.index(k-1);
1458  m_data.value(k) = m_data.value(k-1);
1459  k--;
1460  }
1461  }
1462  else
1463  {
1464  // we are not inserting into the last inner vec
1465  // update outer indices:
1466  Index j = outer+2;
1467  while (j<=m_outerSize && m_outerIndex[j]!=0)
1468  m_outerIndex[j++]++;
1469  --j;
1470  // shift data of last vecs:
1471  Index k = m_outerIndex[j]-1;
1472  while (k>=Index(p))
1473  {
1474  m_data.index(k) = m_data.index(k-1);
1475  m_data.value(k) = m_data.value(k-1);
1476  k--;
1477  }
1478  }
1479  }
1480 
1481  while ( (p > startId) && (m_data.index(p-1) > inner) )
1482  {
1483  m_data.index(p) = m_data.index(p-1);
1484  m_data.value(p) = m_data.value(p-1);
1485  --p;
1486  }
1487 
1488  m_data.index(p) = inner;
1489  return (m_data.value(p) = Scalar(0));
1490 }
1491 
1492 namespace internal {
1493 
1494 template<typename Scalar_, int Options_, typename StorageIndex_>
1495 struct evaluator<SparseMatrix<Scalar_,Options_,StorageIndex_> >
1496  : evaluator<SparseCompressedBase<SparseMatrix<Scalar_,Options_,StorageIndex_> > >
1497 {
1498  typedef evaluator<SparseCompressedBase<SparseMatrix<Scalar_,Options_,StorageIndex_> > > Base;
1499  typedef SparseMatrix<Scalar_,Options_,StorageIndex_> SparseMatrixType;
1500  evaluator() : Base() {}
1501  explicit evaluator(const SparseMatrixType &mat) : Base(mat) {}
1502 };
1503 
1504 }
1505 
1506 } // end namespace Eigen
1507 
1508 #endif // EIGEN_SPARSEMATRIX_H
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:49
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:67
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
Common base class for sparse [compressed]-{row|column}-storage format.
Definition: SparseCompressedBase.h:40
Index nonZeros() const
Definition: SparseCompressedBase.h:58
bool isCompressed() const
Definition: SparseCompressedBase.h:109
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:30
internal::traits< SparseMatrix< Scalar_, Options_, StorageIndex_ > >::StorageIndex StorageIndex
Definition: SparseMatrixBase.h:45
A versatible sparse matrix representation.
Definition: SparseMatrix.h:100
Scalar coeff(Index row, Index col) const
Definition: SparseMatrix.h:192
StorageIndex * innerIndexPtr()
Definition: SparseMatrix.h:165
StorageIndex * innerNonZeroPtr()
Definition: SparseMatrix.h:183
const ConstDiagonalReturnType diagonal() const
Definition: SparseMatrix.h:655
void swap(SparseMatrix &other)
Definition: SparseMatrix.h:730
void setZero()
Definition: SparseMatrix.h:255
Index cols() const
Definition: SparseMatrix.h:142
SparseMatrix(const ReturnByValue< OtherDerived > &other)
Copy constructor with in-place evaluation.
Definition: SparseMatrix.h:713
void setFromTriplets(const InputIterators &begin, const InputIterators &end, DupFunctor dup_func)
Definition: SparseMatrix.h:1114
Index outerSize() const
Definition: SparseMatrix.h:147
SparseMatrix()
Definition: SparseMatrix.h:664
void uncompress()
Definition: SparseMatrix.h:498
Scalar * valuePtr()
Definition: SparseMatrix.h:156
const StorageIndex * innerNonZeroPtr() const
Definition: SparseMatrix.h:179
void makeCompressed()
Definition: SparseMatrix.h:467
StorageIndex * outerIndexPtr()
Definition: SparseMatrix.h:174
Scalar & coeffRef(Index row, Index col)
Definition: SparseMatrix.h:210
const Scalar * valuePtr() const
Definition: SparseMatrix.h:152
SparseMatrix(const SparseSelfAdjointView< OtherDerived, UpLo > &other)
Definition: SparseMatrix.h:698
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:626
bool isCompressed() const
Definition: SparseCompressedBase.h:109
Index rows() const
Definition: SparseMatrix.h:140
SparseMatrix(const SparseMatrix &other)
Definition: SparseMatrix.h:705
void setFromTriplets(const InputIterators &begin, const InputIterators &end)
Definition: SparseMatrix.h:1098
SparseMatrix(const DiagonalBase< OtherDerived > &other)
Copy constructor with in-place evaluation.
Definition: SparseMatrix.h:722
void setIdentity()
Definition: SparseMatrix.h:742
Index innerSize() const
Definition: SparseMatrix.h:145
SparseMatrix(Index rows, Index cols)
Definition: SparseMatrix.h:671
void conservativeResize(Index rows, Index cols)
Definition: SparseMatrix.h:556
void prune(const Scalar &reference, const RealScalar &epsilon=NumTraits< RealScalar >::dummy_precision())
Definition: SparseMatrix.h:510
SparseMatrix(const SparseMatrixBase< OtherDerived > &other)
Definition: SparseMatrix.h:679
~SparseMatrix()
Definition: SparseMatrix.h:835
void prune(const KeepFunc &keep=KeepFunc())
Definition: SparseMatrix.h:523
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:170
void reserve(Index reserveSize)
Definition: SparseMatrix.h:267
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:161
Scalar & insert(Index row, Index col)
Definition: SparseMatrix.h:1234
void reserve(const SizesType &reserveSizes)
DiagonalReturnType diagonal()
Definition: SparseMatrix.h:661
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:47
a sparse vector class
Definition: SparseVector.h:68
static const lastp1_t end
Definition: IndexedViewHelper.h:183
@ ColMajor
Definition: Constants.h:321
@ RowMajor
Definition: Constants.h:323
const unsigned int LvalueBit
Definition: Constants.h:146
const unsigned int RowMajorBit
Definition: Constants.h:68
const unsigned int CompressedAccessBit
Definition: Constants.h:193
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231