Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace Eigen {
16
47namespace internal {
48template<typename Scalar_, int Options_, typename StorageIndex_>
49struct 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
65template<typename Scalar_, int Options_, typename StorageIndex_, int DiagIndex>
66struct 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 typename remove_reference<MatrixTypeNested>::type _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
86template<typename Scalar_, int Options_, typename StorageIndex_, int DiagIndex>
87struct 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
97template<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:
129 typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(IsRowMajor?RowMajorBit:0),StorageIndex> TransposedSparseMatrix;
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
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
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
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
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
848protected:
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
884public:
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 }
899protected:
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
1007private:
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
1022namespace internal {
1023
1024template<typename InputIterator, typename SparseMatrixType, typename DupFunctor>
1025void 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
1096template<typename Scalar, int Options_, typename StorageIndex_>
1097template<typename InputIterators>
1098void 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
1112template<typename Scalar, int Options_, typename StorageIndex_>
1113template<typename InputIterators,typename DupFunctor>
1114void 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
1120template<typename Scalar, int Options_, typename StorageIndex_>
1121template<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
1160template<typename Scalar, int Options_, typename StorageIndex_>
1161template<typename OtherDerived>
1162EIGEN_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 typename internal::remove_all<OtherCopy>::type _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
1233template<typename Scalar_, int Options_, typename StorageIndex_>
1234typename 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
1353template<typename Scalar_, int Options_, typename StorageIndex_>
1354EIGEN_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
1385template<typename Scalar_, int Options_, typename StorageIndex_>
1386EIGEN_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
1492namespace internal {
1493
1494template<typename Scalar_, int Options_, typename StorageIndex_>
1495struct 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
const ConstDiagonalReturnType diagonal() const
Definition: SparseMatrix.h:655
void swap(SparseMatrix &other)
Definition: SparseMatrix.h:730
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:161
void setZero()
Definition: SparseMatrix.h:255
Index cols() const
Definition: SparseMatrix.h:142
StorageIndex * innerIndexPtr()
Definition: SparseMatrix.h:165
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
const Scalar * valuePtr() const
Definition: SparseMatrix.h:152
void makeCompressed()
Definition: SparseMatrix.h:467
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
StorageIndex * innerNonZeroPtr()
Definition: SparseMatrix.h:183
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
StorageIndex * outerIndexPtr()
Definition: SparseMatrix.h:174
const StorageIndex * innerNonZeroPtr() const
Definition: SparseMatrix.h:179
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
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:170
SparseMatrix(const SparseMatrixBase< OtherDerived > &other)
Definition: SparseMatrix.h:679
~SparseMatrix()
Definition: SparseMatrix.h:835
Scalar * valuePtr()
Definition: SparseMatrix.h:156
void prune(const KeepFunc &keep=KeepFunc())
Definition: SparseMatrix.h:523
void reserve(Index reserveSize)
Definition: SparseMatrix.h:267
Scalar & coeffRef(Index row, Index col)
Definition: SparseMatrix.h:210
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: B01_Experimental.dox:1
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
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Derived & derived()
Definition: EigenBase.h:48
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235