Overhaul Sparse Core

Eigen/SparseCore

@@ -17,6 +17,7 @@
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
+#include <numeric>
 
 /** 
   * \defgroup SparseCore_Module SparseCore module

Eigen/src/SparseCore/CompressedStorage.h

@@ -132,15 +132,7 @@ class CompressedStorage
     /** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */
     inline Index searchLowerIndex(Index start, Index end, Index key) const
     {
-      while(end>start)
+      return static_cast<Index>(std::distance(m_indices, std::lower_bound(m_indices + start, m_indices + end, key)));
-      {
-        Index mid = (end+start)>>1;
-        if (m_indices[mid]<key)
-          start = mid+1;
-        else
-          end = mid;
-      }
-      return start;
     }
 
     /** \returns the stored value at index \a key
@@ -198,19 +190,12 @@ class CompressedStorage
       return m_values[id];
     }
 
-    void moveChunk(Index from, Index to, Index chunkSize)
+    inline void moveChunk(Index from, Index to, Index chunkSize)
     {
-      eigen_internal_assert(to+chunkSize <= m_size);
+      eigen_internal_assert(chunkSize >= 0 && to+chunkSize <= m_size);
-      if(to>from && from+chunkSize>to)
+      if (chunkSize > 0) {
-      {
+        internal::smart_memmove(m_values + from, m_values + from + chunkSize, m_values + to);
-        // move backward
+        internal::smart_memmove(m_indices + from, m_indices + from + chunkSize, m_indices + to);
-        internal::smart_memmove(m_values+from,  m_values+from+chunkSize,  m_values+to);
-        internal::smart_memmove(m_indices+from, m_indices+from+chunkSize, m_indices+to);
-      }
-      else
-      {
-        internal::smart_copy(m_values+from,  m_values+from+chunkSize,  m_values+to);
-        internal::smart_copy(m_indices+from, m_indices+from+chunkSize, m_indices+to);
       }
     }
 

test/sparse_basic.cpp

@@ -477,14 +477,45 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
         refMat_prod(r,c) *= v;
       refMat_last(r,c) = v;
     }
+
     SparseMatrixType m(rows,cols);
     m.setFromTriplets(triplets.begin(), triplets.end());
     VERIFY_IS_APPROX(m, refMat_sum);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
 
     m.setFromTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>());
     VERIFY_IS_APPROX(m, refMat_prod);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
     m.setFromTriplets(triplets.begin(), triplets.end(), [] (Scalar,Scalar b) { return b; });
     VERIFY_IS_APPROX(m, refMat_last);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
+
+    // test setFromSortedTriplets
+
+    struct triplet_comp {
+      inline bool operator()(const TripletType& a, const TripletType& b) {
+        return SparseMatrixType::IsRowMajor ? ((a.row() != b.row()) ? (a.row() < b.row()) : (a.col() < b.col()))
+                                            : ((a.col() != b.col()) ? (a.col() < b.col()) : (a.row() < b.row()));
+      }
+    };
+
+    // stable_sort is only necessary when the reduction functor is dependent on the order of the triplets
+    // this is the case with refMat_last
+    // for most cases, std::sort is sufficient and preferred
+    std::stable_sort(triplets.begin(), triplets.end(), triplet_comp());
+
+    m.setZero();
+    m.setFromSortedTriplets(triplets.begin(), triplets.end());
+    VERIFY_IS_APPROX(m, refMat_sum);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
+
+    m.setFromSortedTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>());
+    VERIFY_IS_APPROX(m, refMat_prod);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
+
+    m.setFromSortedTriplets(triplets.begin(), triplets.end(), [](Scalar, Scalar b) { return b; });
+    VERIFY_IS_APPROX(m, refMat_last);
+    VERIFY_IS_EQUAL(m.innerIndicesAreSorted(), m.outerSize());
   }
   
   // test Map