10 #ifndef EIGEN_SPARSEMATRIX_H
11 #define EIGEN_SPARSEMATRIX_H
46 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
47 struct traits<SparseMatrix<_Scalar, _Options, _StorageIndex> >
49 typedef _Scalar Scalar;
50 typedef _StorageIndex StorageIndex;
51 typedef Sparse StorageKind;
52 typedef MatrixXpr XprKind;
59 SupportedAccessPatterns = InnerRandomAccessPattern
63 template<
typename _Scalar,
int _Options,
typename _StorageIndex,
int DiagIndex>
64 struct traits<Diagonal<SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex> >
66 typedef SparseMatrix<_Scalar, _Options, _StorageIndex> MatrixType;
67 typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
68 typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
70 typedef _Scalar Scalar;
71 typedef Dense StorageKind;
72 typedef _StorageIndex StorageIndex;
73 typedef MatrixXpr XprKind;
77 ColsAtCompileTime = 1,
79 MaxColsAtCompileTime = 1,
84 template<
typename _Scalar,
int _Options,
typename _StorageIndex,
int DiagIndex>
85 struct traits<Diagonal<const SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex> >
86 :
public traits<Diagonal<SparseMatrix<_Scalar, _Options, _StorageIndex>, DiagIndex> >
95 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
100 using Base::convert_index;
106 using Base::operator+=;
107 using Base::operator-=;
112 typedef typename Base::InnerIterator InnerIterator;
113 typedef typename Base::ReverseInnerIterator ReverseInnerIterator;
116 using Base::IsRowMajor;
117 typedef internal::CompressedStorage<Scalar,StorageIndex> Storage;
123 typedef typename Base::ScalarVector ScalarVector;
129 StorageIndex* m_outerIndex;
130 StorageIndex* m_innerNonZeros;
136 inline Index rows()
const {
return IsRowMajor ? m_outerSize : m_innerSize; }
138 inline Index cols()
const {
return IsRowMajor ? m_innerSize : m_outerSize; }
148 inline const Scalar*
valuePtr()
const {
return m_data.valuePtr(); }
152 inline Scalar*
valuePtr() {
return m_data.valuePtr(); }
157 inline const StorageIndex*
innerIndexPtr()
const {
return m_data.indexPtr(); }
182 inline Storage& data() {
return m_data; }
184 inline const Storage& data()
const {
return m_data; }
190 eigen_assert(row>=0 && row<
rows() && col>=0 && col<
cols());
192 const Index outer = IsRowMajor ? row : col;
193 const Index inner = IsRowMajor ? col : row;
194 Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer+1];
195 return m_data.atInRange(m_outerIndex[outer], end, StorageIndex(inner));
208 eigen_assert(row>=0 && row<
rows() && col>=0 && col<
cols());
210 const Index outer = IsRowMajor ? row : col;
211 const Index inner = IsRowMajor ? col : row;
213 Index start = m_outerIndex[outer];
214 Index end = m_innerNonZeros ? m_outerIndex[outer] + m_innerNonZeros[outer] : m_outerIndex[outer+1];
215 eigen_assert(end>=start &&
"you probably called coeffRef on a non finalized matrix");
218 const Index p = m_data.searchLowerIndex(start,end-1,StorageIndex(inner));
219 if((p<end) && (m_data.index(p)==inner))
220 return m_data.value(p);
254 memset(m_outerIndex, 0, (m_outerSize+1)*
sizeof(StorageIndex));
256 memset(m_innerNonZeros, 0, (m_outerSize)*
sizeof(StorageIndex));
264 eigen_assert(
isCompressed() &&
"This function does not make sense in non compressed mode.");
265 m_data.reserve(reserveSize);
268 #ifdef EIGEN_PARSED_BY_DOXYGEN
281 template<
class SizesType>
282 inline void reserve(
const SizesType& reserveSizes);
284 template<
class SizesType>
285 inline void reserve(
const SizesType& reserveSizes,
const typename SizesType::value_type& enableif =
286 #
if (!EIGEN_COMP_MSVC) || (EIGEN_COMP_MSVC>=1500)
289 SizesType::value_type())
291 EIGEN_UNUSED_VARIABLE(enableif);
292 reserveInnerVectors(reserveSizes);
294 #endif // EIGEN_PARSED_BY_DOXYGEN
296 template<
class SizesType>
297 inline void reserveInnerVectors(
const SizesType& reserveSizes)
301 Index totalReserveSize = 0;
303 m_innerNonZeros =
static_cast<StorageIndex*
>(std::malloc(m_outerSize *
sizeof(StorageIndex)));
304 if (!m_innerNonZeros) internal::throw_std_bad_alloc();
307 StorageIndex* newOuterIndex = m_innerNonZeros;
309 StorageIndex count = 0;
310 for(
Index j=0; j<m_outerSize; ++j)
312 newOuterIndex[j] = count;
313 count += reserveSizes[j] + (m_outerIndex[j+1]-m_outerIndex[j]);
314 totalReserveSize += reserveSizes[j];
316 m_data.reserve(totalReserveSize);
317 StorageIndex previousOuterIndex = m_outerIndex[m_outerSize];
318 for(
Index j=m_outerSize-1; j>=0; --j)
320 StorageIndex innerNNZ = previousOuterIndex - m_outerIndex[j];
321 for(
Index i=innerNNZ-1; i>=0; --i)
323 m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i);
324 m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i);
326 previousOuterIndex = m_outerIndex[j];
327 m_outerIndex[j] = newOuterIndex[j];
328 m_innerNonZeros[j] = innerNNZ;
331 m_outerIndex[m_outerSize] = m_outerIndex[m_outerSize-1] + m_innerNonZeros[m_outerSize-1] + reserveSizes[m_outerSize-1];
333 m_data.resize(m_outerIndex[m_outerSize]);
337 StorageIndex* newOuterIndex =
static_cast<StorageIndex*
>(std::malloc((m_outerSize+1)*
sizeof(StorageIndex)));
338 if (!newOuterIndex) internal::throw_std_bad_alloc();
340 StorageIndex count = 0;
341 for(
Index j=0; j<m_outerSize; ++j)
343 newOuterIndex[j] = count;
344 StorageIndex alreadyReserved = (m_outerIndex[j+1]-m_outerIndex[j]) - m_innerNonZeros[j];
345 StorageIndex toReserve = std::max<StorageIndex>(reserveSizes[j], alreadyReserved);
346 count += toReserve + m_innerNonZeros[j];
348 newOuterIndex[m_outerSize] = count;
350 m_data.resize(count);
351 for(
Index j=m_outerSize-1; j>=0; --j)
353 Index offset = newOuterIndex[j] - m_outerIndex[j];
356 StorageIndex innerNNZ = m_innerNonZeros[j];
357 for(
Index i=innerNNZ-1; i>=0; --i)
359 m_data.index(newOuterIndex[j]+i) = m_data.index(m_outerIndex[j]+i);
360 m_data.value(newOuterIndex[j]+i) = m_data.value(m_outerIndex[j]+i);
365 std::swap(m_outerIndex, newOuterIndex);
366 std::free(newOuterIndex);
384 inline Scalar& insertBack(
Index row,
Index col)
386 return insertBackByOuterInner(IsRowMajor?row:col, IsRowMajor?col:row);
391 inline Scalar& insertBackByOuterInner(
Index outer,
Index inner)
393 eigen_assert(
Index(m_outerIndex[outer+1]) == m_data.size() &&
"Invalid ordered insertion (invalid outer index)");
394 eigen_assert( (m_outerIndex[outer+1]-m_outerIndex[outer]==0 || m_data.index(m_data.size()-1)<inner) &&
"Invalid ordered insertion (invalid inner index)");
395 Index p = m_outerIndex[outer+1];
396 ++m_outerIndex[outer+1];
397 m_data.append(Scalar(0), inner);
398 return m_data.value(p);
403 inline Scalar& insertBackByOuterInnerUnordered(
Index outer,
Index inner)
405 Index p = m_outerIndex[outer+1];
406 ++m_outerIndex[outer+1];
407 m_data.append(Scalar(0), inner);
408 return m_data.value(p);
413 inline void startVec(
Index outer)
415 eigen_assert(m_outerIndex[outer]==
Index(m_data.size()) &&
"You must call startVec for each inner vector sequentially");
416 eigen_assert(m_outerIndex[outer+1]==0 &&
"You must call startVec for each inner vector sequentially");
417 m_outerIndex[outer+1] = m_outerIndex[outer];
423 inline void finalize()
427 StorageIndex size = internal::convert_index<StorageIndex>(m_data.size());
428 Index i = m_outerSize;
430 while (i>=0 && m_outerIndex[i]==0)
433 while (i<=m_outerSize)
435 m_outerIndex[i] = size;
443 template<
typename InputIterators>
444 void setFromTriplets(
const InputIterators& begin,
const InputIterators& end);
446 template<
typename InputIterators,
typename DupFunctor>
447 void setFromTriplets(
const InputIterators& begin,
const InputIterators& end, DupFunctor dup_func);
449 void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op<Scalar,Scalar>()); }
451 template<
typename DupFunctor>
452 void collapseDuplicates(DupFunctor dup_func = DupFunctor());
460 return insert(IsRowMajor ? j : i, IsRowMajor ? i : j);
470 eigen_internal_assert(m_outerIndex!=0 && m_outerSize>0);
472 Index oldStart = m_outerIndex[1];
473 m_outerIndex[1] = m_innerNonZeros[0];
474 for(
Index j=1; j<m_outerSize; ++j)
476 Index nextOldStart = m_outerIndex[j+1];
477 Index offset = oldStart - m_outerIndex[j];
480 for(
Index k=0; k<m_innerNonZeros[j]; ++k)
482 m_data.index(m_outerIndex[j]+k) = m_data.index(oldStart+k);
483 m_data.value(m_outerIndex[j]+k) = m_data.value(oldStart+k);
486 m_outerIndex[j+1] = m_outerIndex[j] + m_innerNonZeros[j];
487 oldStart = nextOldStart;
489 std::free(m_innerNonZeros);
491 m_data.resize(m_outerIndex[m_outerSize]);
498 if(m_innerNonZeros != 0)
500 m_innerNonZeros =
static_cast<StorageIndex*
>(std::malloc(m_outerSize *
sizeof(StorageIndex)));
501 for (
Index i = 0; i < m_outerSize; i++)
503 m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i];
510 prune(default_prunning_func(reference,epsilon));
520 template<
typename KeepFunc>
521 void prune(
const KeepFunc& keep = KeepFunc())
527 for(
Index j=0; j<m_outerSize; ++j)
529 Index previousStart = m_outerIndex[j];
531 Index end = m_outerIndex[j+1];
532 for(
Index i=previousStart; i<end; ++i)
534 if(keep(IsRowMajor?j:m_data.index(i), IsRowMajor?m_data.index(i):j, m_data.value(i)))
536 m_data.value(k) = m_data.value(i);
537 m_data.index(k) = m_data.index(i);
542 m_outerIndex[m_outerSize] = k;
557 if (this->
rows() == rows && this->
cols() == cols)
return;
564 StorageIndex newInnerSize = convert_index(IsRowMajor ?
cols :
rows);
570 StorageIndex *newInnerNonZeros =
static_cast<StorageIndex*
>(std::realloc(m_innerNonZeros, (m_outerSize + outerChange) *
sizeof(StorageIndex)));
571 if (!newInnerNonZeros) internal::throw_std_bad_alloc();
572 m_innerNonZeros = newInnerNonZeros;
574 for(
Index i=m_outerSize; i<m_outerSize+outerChange; i++)
575 m_innerNonZeros[i] = 0;
577 else if (innerChange < 0)
580 m_innerNonZeros =
static_cast<StorageIndex*
>(std::malloc((m_outerSize+outerChange+1) *
sizeof(StorageIndex)));
581 if (!m_innerNonZeros) internal::throw_std_bad_alloc();
582 for(
Index i = 0; i < m_outerSize; i++)
583 m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i];
587 if (m_innerNonZeros && innerChange < 0)
589 for(
Index i = 0; i < m_outerSize + (std::min)(outerChange,
Index(0)); i++)
591 StorageIndex &n = m_innerNonZeros[i];
592 StorageIndex start = m_outerIndex[i];
593 while (n > 0 && m_data.index(start+n-1) >= newInnerSize) --n;
597 m_innerSize = newInnerSize;
600 if (outerChange == 0)
603 StorageIndex *newOuterIndex =
static_cast<StorageIndex*
>(std::realloc(m_outerIndex, (m_outerSize + outerChange + 1) *
sizeof(StorageIndex)));
604 if (!newOuterIndex) internal::throw_std_bad_alloc();
605 m_outerIndex = newOuterIndex;
608 StorageIndex last = m_outerSize == 0 ? 0 : m_outerIndex[m_outerSize];
609 for(
Index i=m_outerSize; i<m_outerSize+outerChange+1; i++)
610 m_outerIndex[i] = last;
612 m_outerSize += outerChange;
625 m_innerSize = IsRowMajor ?
cols :
rows;
627 if (m_outerSize !=
outerSize || m_outerSize==0)
629 std::free(m_outerIndex);
630 m_outerIndex =
static_cast<StorageIndex*
>(std::malloc((
outerSize + 1) *
sizeof(StorageIndex)));
631 if (!m_outerIndex) internal::throw_std_bad_alloc();
637 std::free(m_innerNonZeros);
640 memset(m_outerIndex, 0, (m_outerSize+1)*
sizeof(StorageIndex));
645 void resizeNonZeros(
Index size)
661 : m_outerSize(-1), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
663 check_template_parameters();
669 : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
671 check_template_parameters();
676 template<
typename OtherDerived>
678 : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
680 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
681 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
682 check_template_parameters();
683 const bool needToTranspose = (Flags &
RowMajorBit) != (internal::evaluator<OtherDerived>::Flags &
RowMajorBit);
688 #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
689 EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
691 internal::call_assignment_no_alias(*
this, other.
derived());
696 template<
typename OtherDerived,
unsigned int UpLo>
698 : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
700 check_template_parameters();
701 Base::operator=(other);
706 :
Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
708 check_template_parameters();
709 *
this = other.derived();
713 template<
typename OtherDerived>
715 :
Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
717 check_template_parameters();
718 initAssignment(other);
723 template<
typename OtherDerived>
725 :
Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
727 check_template_parameters();
728 *
this = other.derived();
736 std::swap(m_outerIndex, other.m_outerIndex);
737 std::swap(m_innerSize, other.m_innerSize);
738 std::swap(m_outerSize, other.m_outerSize);
739 std::swap(m_innerNonZeros, other.m_innerNonZeros);
740 m_data.swap(other.m_data);
747 eigen_assert(
rows() ==
cols() &&
"ONLY FOR SQUARED MATRICES");
748 this->m_data.resize(
rows());
752 std::free(m_innerNonZeros);
757 if (other.isRValue())
759 swap(other.const_cast_derived());
761 else if(
this!=&other)
763 #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
764 EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
766 initAssignment(other);
769 internal::smart_copy(other.m_outerIndex, other.m_outerIndex + m_outerSize + 1, m_outerIndex);
770 m_data = other.m_data;
774 Base::operator=(other);
780 #ifndef EIGEN_PARSED_BY_DOXYGEN
781 template<
typename OtherDerived>
782 inline SparseMatrix& operator=(
const EigenBase<OtherDerived>& other)
783 {
return Base::operator=(other.derived()); }
784 #endif // EIGEN_PARSED_BY_DOXYGEN
786 template<
typename OtherDerived>
787 EIGEN_DONT_INLINE
SparseMatrix& operator=(
const SparseMatrixBase<OtherDerived>& other);
789 friend std::ostream & operator << (std::ostream & s,
const SparseMatrix& m)
792 s <<
"Nonzero entries:\n";
795 for (Index i=0; i<m.nonZeros(); ++i)
796 s <<
"(" << m.m_data.value(i) <<
"," << m.m_data.index(i) <<
") ";
800 for (Index i=0; i<m.outerSize(); ++i)
802 Index p = m.m_outerIndex[i];
803 Index pe = m.m_outerIndex[i]+m.m_innerNonZeros[i];
806 s <<
"(" << m.m_data.value(k) <<
"," << m.m_data.index(k) <<
") ";
808 for (; k<m.m_outerIndex[i+1]; ++k) {
815 s <<
"Outer pointers:\n";
816 for (
Index i=0; i<m.outerSize(); ++i) {
817 s << m.m_outerIndex[i] <<
" ";
819 s <<
" $" << std::endl;
820 if(!m.isCompressed())
822 s <<
"Inner non zeros:\n";
823 for (Index i=0; i<m.outerSize(); ++i) {
824 s << m.m_innerNonZeros[i] <<
" ";
826 s <<
" $" << std::endl;
830 s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
837 std::free(m_outerIndex);
838 std::free(m_innerNonZeros);
844 # ifdef EIGEN_SPARSEMATRIX_PLUGIN
845 # include EIGEN_SPARSEMATRIX_PLUGIN
850 template<
typename Other>
851 void initAssignment(
const Other& other)
853 resize(other.rows(), other.cols());
856 std::free(m_innerNonZeros);
863 EIGEN_DONT_INLINE Scalar& insertCompressed(Index row, Index col);
867 class SingletonVector
869 StorageIndex m_index;
870 StorageIndex m_value;
872 typedef StorageIndex value_type;
873 SingletonVector(Index i, Index v)
874 : m_index(convert_index(i)), m_value(convert_index(v))
877 StorageIndex operator[](Index i)
const {
return i==m_index ? m_value : 0; }
882 EIGEN_DONT_INLINE Scalar& insertUncompressed(Index row, Index col);
887 EIGEN_STRONG_INLINE Scalar& insertBackUncompressed(Index row, Index col)
889 const Index outer = IsRowMajor ? row : col;
890 const Index inner = IsRowMajor ? col : row;
892 eigen_assert(!isCompressed());
893 eigen_assert(m_innerNonZeros[outer]<=(m_outerIndex[outer+1] - m_outerIndex[outer]));
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));
901 static void check_template_parameters()
903 EIGEN_STATIC_ASSERT(NumTraits<StorageIndex>::IsSigned,THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE);
904 EIGEN_STATIC_ASSERT((Options&(
ColMajor|
RowMajor))==Options,INVALID_MATRIX_TEMPLATE_PARAMETERS);
907 struct default_prunning_func {
908 default_prunning_func(
const Scalar& ref,
const RealScalar& eps) : reference(ref), epsilon(eps) {}
909 inline bool operator() (
const Index&,
const Index&,
const Scalar& value)
const
911 return !internal::isMuchSmallerThan(value, reference, epsilon);
920 template<
typename InputIterator,
typename SparseMatrixType,
typename DupFunctor>
921 void set_from_triplets(
const InputIterator& begin,
const InputIterator& end, SparseMatrixType& mat, DupFunctor dup_func)
923 enum { IsRowMajor = SparseMatrixType::IsRowMajor };
924 typedef typename SparseMatrixType::Scalar Scalar;
925 typedef typename SparseMatrixType::StorageIndex StorageIndex;
926 SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor,StorageIndex> trMat(mat.rows(),mat.cols());
931 typename SparseMatrixType::IndexVector wi(trMat.outerSize());
933 for(InputIterator it(begin); it!=end; ++it)
935 eigen_assert(it->row()>=0 && it->row()<mat.rows() && it->col()>=0 && it->col()<mat.cols());
936 wi(IsRowMajor ? it->col() : it->row())++;
941 for(InputIterator it(begin); it!=end; ++it)
942 trMat.insertBackUncompressed(it->row(),it->col()) = it->value();
945 trMat.collapseDuplicates(dup_func);
992 template<
typename Scalar,
int _Options,
typename _StorageIndex>
993 template<
typename InputIterators>
996 internal::set_from_triplets<InputIterators, SparseMatrix<Scalar,_Options,_StorageIndex> >(begin, end, *
this, internal::scalar_sum_op<Scalar,Scalar>());
1008 template<
typename Scalar,
int _Options,
typename _StorageIndex>
1009 template<
typename InputIterators,
typename DupFunctor>
1012 internal::set_from_triplets<InputIterators, SparseMatrix<Scalar,_Options,_StorageIndex>, DupFunctor>(begin, end, *
this, dup_func);
1016 template<
typename Scalar,
int _Options,
typename _StorageIndex>
1017 template<
typename DupFunctor>
1020 eigen_assert(!isCompressed());
1022 IndexVector wi(innerSize());
1024 StorageIndex count = 0;
1026 for(
Index j=0; j<outerSize(); ++j)
1028 StorageIndex start = count;
1029 Index oldEnd = m_outerIndex[j]+m_innerNonZeros[j];
1030 for(
Index k=m_outerIndex[j]; k<oldEnd; ++k)
1032 Index i = m_data.index(k);
1036 m_data.value(wi(i)) = dup_func(m_data.value(wi(i)), m_data.value(k));
1040 m_data.value(count) = m_data.value(k);
1041 m_data.index(count) = m_data.index(k);
1046 m_outerIndex[j] = start;
1048 m_outerIndex[m_outerSize] = count;
1051 std::free(m_innerNonZeros);
1052 m_innerNonZeros = 0;
1053 m_data.resize(m_outerIndex[m_outerSize]);
1056 template<
typename Scalar,
int _Options,
typename _StorageIndex>
1057 template<
typename OtherDerived>
1058 EIGEN_DONT_INLINE SparseMatrix<Scalar,_Options,_StorageIndex>& SparseMatrix<Scalar,_Options,_StorageIndex>::operator=(
const SparseMatrixBase<OtherDerived>& other)
1060 EIGEN_STATIC_ASSERT((internal::is_same<Scalar, typename OtherDerived::Scalar>::value),
1061 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
1063 #ifdef EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
1064 EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN
1067 const bool needToTranspose = (Flags &
RowMajorBit) != (internal::evaluator<OtherDerived>::Flags &
RowMajorBit);
1068 if (needToTranspose)
1070 #ifdef EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
1071 EIGEN_SPARSE_TRANSPOSED_COPY_PLUGIN
1077 typedef typename internal::nested_eval<OtherDerived,2,typename internal::plain_matrix_type<OtherDerived>::type >::type OtherCopy;
1078 typedef typename internal::remove_all<OtherCopy>::type _OtherCopy;
1079 typedef internal::evaluator<_OtherCopy> OtherCopyEval;
1080 OtherCopy otherCopy(other.derived());
1081 OtherCopyEval otherCopyEval(otherCopy);
1083 SparseMatrix dest(other.rows(),other.cols());
1088 for (Index j=0; j<otherCopy.outerSize(); ++j)
1089 for (
typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it)
1090 ++dest.m_outerIndex[it.index()];
1093 StorageIndex count = 0;
1094 IndexVector positions(dest.outerSize());
1095 for (Index j=0; j<dest.outerSize(); ++j)
1097 StorageIndex tmp = dest.m_outerIndex[j];
1098 dest.m_outerIndex[j] = count;
1099 positions[j] = count;
1102 dest.m_outerIndex[dest.outerSize()] = count;
1104 dest.m_data.resize(count);
1106 for (StorageIndex j=0; j<otherCopy.outerSize(); ++j)
1108 for (
typename OtherCopyEval::InnerIterator it(otherCopyEval, j); it; ++it)
1110 Index pos = positions[it.index()]++;
1111 dest.m_data.index(pos) = j;
1112 dest.m_data.value(pos) = it.value();
1120 if(other.isRValue())
1122 initAssignment(other.derived());
1125 return Base::operator=(other.derived());
1129 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
1132 eigen_assert(row>=0 && row<rows() && col>=0 && col<cols());
1134 const Index outer = IsRowMajor ? row : col;
1135 const Index inner = IsRowMajor ? col : row;
1142 if(m_data.allocatedSize()==0)
1143 m_data.reserve(2*m_innerSize);
1146 m_innerNonZeros =
static_cast<StorageIndex*
>(std::malloc(m_outerSize *
sizeof(StorageIndex)));
1147 if(!m_innerNonZeros) internal::throw_std_bad_alloc();
1149 memset(m_innerNonZeros, 0, (m_outerSize)*
sizeof(StorageIndex));
1153 StorageIndex end = convert_index(m_data.allocatedSize());
1154 for(
Index j=1; j<=m_outerSize; ++j)
1155 m_outerIndex[j] = end;
1160 m_innerNonZeros =
static_cast<StorageIndex*
>(std::malloc(m_outerSize *
sizeof(StorageIndex)));
1161 if(!m_innerNonZeros) internal::throw_std_bad_alloc();
1162 for(
Index j=0; j<m_outerSize; ++j)
1163 m_innerNonZeros[j] = m_outerIndex[j+1]-m_outerIndex[j];
1168 Index data_end = m_data.allocatedSize();
1172 if(m_outerIndex[outer]==data_end)
1174 eigen_internal_assert(m_innerNonZeros[outer]==0);
1178 StorageIndex p = convert_index(m_data.size());
1180 while(j>=0 && m_innerNonZeros[j]==0)
1181 m_outerIndex[j--] = p;
1184 ++m_innerNonZeros[outer];
1185 m_data.append(Scalar(0), inner);
1188 if(data_end != m_data.allocatedSize())
1193 eigen_internal_assert(data_end < m_data.allocatedSize());
1194 StorageIndex new_end = convert_index(m_data.allocatedSize());
1195 for(
Index k=outer+1; k<=m_outerSize; ++k)
1196 if(m_outerIndex[k]==data_end)
1197 m_outerIndex[k] = new_end;
1199 return m_data.value(p);
1204 if(m_outerIndex[outer+1]==data_end && m_outerIndex[outer]+m_innerNonZeros[outer]==m_data.size())
1206 eigen_internal_assert(outer+1==m_outerSize || m_innerNonZeros[outer+1]==0);
1209 ++m_innerNonZeros[outer];
1210 m_data.resize(m_data.size()+1);
1213 if(data_end != m_data.allocatedSize())
1218 eigen_internal_assert(data_end < m_data.allocatedSize());
1219 StorageIndex new_end = convert_index(m_data.allocatedSize());
1220 for(
Index k=outer+1; k<=m_outerSize; ++k)
1221 if(m_outerIndex[k]==data_end)
1222 m_outerIndex[k] = new_end;
1226 Index startId = m_outerIndex[outer];
1227 Index p = m_outerIndex[outer]+m_innerNonZeros[outer]-1;
1228 while ( (p > startId) && (m_data.index(p-1) > inner) )
1230 m_data.index(p) = m_data.index(p-1);
1231 m_data.value(p) = m_data.value(p-1);
1235 m_data.index(p) = convert_index(inner);
1236 return (m_data.value(p) = 0);
1239 if(m_data.size() != m_data.allocatedSize())
1242 m_data.resize(m_data.allocatedSize());
1246 return insertUncompressed(row,col);
1249 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
1252 eigen_assert(!isCompressed());
1254 const Index outer = IsRowMajor ? row : col;
1255 const StorageIndex inner = convert_index(IsRowMajor ? col : row);
1257 Index room = m_outerIndex[outer+1] - m_outerIndex[outer];
1258 StorageIndex innerNNZ = m_innerNonZeros[outer];
1262 reserve(SingletonVector(outer,std::max<StorageIndex>(2,innerNNZ)));
1265 Index startId = m_outerIndex[outer];
1266 Index p = startId + m_innerNonZeros[outer];
1267 while ( (p > startId) && (m_data.index(p-1) > inner) )
1269 m_data.index(p) = m_data.index(p-1);
1270 m_data.value(p) = m_data.value(p-1);
1273 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");
1275 m_innerNonZeros[outer]++;
1277 m_data.index(p) = inner;
1278 return (m_data.value(p) = Scalar(0));
1281 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
1282 EIGEN_DONT_INLINE
typename SparseMatrix<_Scalar,_Options,_StorageIndex>::Scalar& SparseMatrix<_Scalar,_Options,_StorageIndex>::insertCompressed(Index row, Index col)
1284 eigen_assert(isCompressed());
1286 const Index outer = IsRowMajor ? row : col;
1287 const Index inner = IsRowMajor ? col : row;
1289 Index previousOuter = outer;
1290 if (m_outerIndex[outer+1]==0)
1293 while (previousOuter>=0 && m_outerIndex[previousOuter]==0)
1295 m_outerIndex[previousOuter] = convert_index(m_data.size());
1298 m_outerIndex[outer+1] = m_outerIndex[outer];
1304 bool isLastVec = (!(previousOuter==-1 && m_data.size()!=0))
1305 && (std::size_t(m_outerIndex[outer+1]) == m_data.size());
1307 std::size_t startId = m_outerIndex[outer];
1309 std::size_t p = m_outerIndex[outer+1];
1310 ++m_outerIndex[outer+1];
1312 double reallocRatio = 1;
1313 if (m_data.allocatedSize()<=m_data.size())
1316 if (m_data.size()==0)
1325 double nnzEstimate = double(m_outerIndex[outer])*double(m_outerSize)/double(outer+1);
1326 reallocRatio = (nnzEstimate-double(m_data.size()))/
double(m_data.size());
1330 reallocRatio = (std::min)((std::max)(reallocRatio,1.5),8.);
1333 m_data.resize(m_data.size()+1,reallocRatio);
1337 if (previousOuter==-1)
1341 for (Index k=0; k<=(outer+1); ++k)
1342 m_outerIndex[k] = 0;
1344 while(m_outerIndex[k]==0)
1345 m_outerIndex[k++] = 1;
1346 while (k<=m_outerSize && m_outerIndex[k]!=0)
1347 m_outerIndex[k++]++;
1350 k = m_outerIndex[k]-1;
1353 m_data.index(k) = m_data.index(k-1);
1354 m_data.value(k) = m_data.value(k-1);
1363 while (j<=m_outerSize && m_outerIndex[j]!=0)
1364 m_outerIndex[j++]++;
1367 Index k = m_outerIndex[j]-1;
1370 m_data.index(k) = m_data.index(k-1);
1371 m_data.value(k) = m_data.value(k-1);
1377 while ( (p > startId) && (m_data.index(p-1) > inner) )
1379 m_data.index(p) = m_data.index(p-1);
1380 m_data.value(p) = m_data.value(p-1);
1384 m_data.index(p) = inner;
1385 return (m_data.value(p) = Scalar(0));
1388 namespace internal {
1390 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
1391 struct evaluator<SparseMatrix<_Scalar,_Options,_StorageIndex> >
1392 : evaluator<SparseCompressedBase<SparseMatrix<_Scalar,_Options,_StorageIndex> > >
1394 typedef evaluator<SparseCompressedBase<SparseMatrix<_Scalar,_Options,_StorageIndex> > > Base;
1395 typedef SparseMatrix<_Scalar,_Options,_StorageIndex> SparseMatrixType;
1396 evaluator() : Base() {}
1397 explicit evaluator(
const SparseMatrixType &mat) : Base(mat) {}
1404 #endif // EIGEN_SPARSEMATRIX_H