sparse module:

- remove some useless stuff => let's focus on a single sparse matrix format
 - finalize the new RandomSetter

Eigen/src/Sparse/RandomSetter.h

@@ -29,6 +29,9 @@ template<typename Scalar> struct StdMapTraits
 {
   typedef int KeyType;
   typedef std::map<KeyType,Scalar> Type;
+  enum {
+    IsSorted = 1
+  };
 
   static void setInvalidKey(Type&, const KeyType&) {}
 };
@@ -38,6 +41,9 @@ template<typename Scalar> struct GnuHashMapTraits
 {
   typedef int KeyType;
   typedef __gnu_cxx::hash_map<KeyType,Scalar> Type;
+  enum {
+    IsSorted = 0
+  };
 
   static void setInvalidKey(Type&, const KeyType&) {}
 };
@@ -48,6 +54,9 @@ template<typename Scalar> struct GoogleDenseHashMapTraits
 {
   typedef int KeyType;
   typedef google::dense_hash_map<KeyType,Scalar> Type;
+  enum {
+    IsSorted = 0
+  };
 
   static void setInvalidKey(Type& map, const KeyType& k)
   { map.set_empty_key(k); }
@@ -59,6 +68,9 @@ template<typename Scalar> struct GoogleSparseHashMapTraits
 {
   typedef int KeyType;
   typedef google::sparse_hash_map<KeyType,Scalar> Type;
+  enum {
+    IsSorted = 0
+  };
 
   static void setInvalidKey(Type&, const KeyType&) {}
 };
@@ -66,10 +78,35 @@ template<typename Scalar> struct GoogleSparseHashMapTraits
 
 /** \class RandomSetter
   *
+  * Typical usage:
+  * \code
+  * SparseMatrix<double> m(rows,cols);
+  * {
+  *   RandomSetter<SparseMatrix<double> > w(m);
+  *   // don't use m but w instead with read/write random access to the coefficients:
+  *   for(;;)
+  *     w(rand(),rand()) = rand;
+  * }
+  * // when w is deleted, the data are copied back to m
+  * // and m is ready to use.
+  * \endcode
+  *
+  * \note for performance and memory consumption reasons it is highly recommended to use
+  * Google's hash library. To do so you have two options:
+  *  - include <google/dense_hash_map> yourself \b before Eigen/Sparse header
+  *  - define EIGEN_GOOGLEHASH_SUPPORT
+  * In the later case the inclusion of <google/dense_hash_map> is made for you.
   */
 template<typename SparseMatrixType,
-         template <typename T> class HashMapTraits = StdMapTraits,
-         int OuterPacketBits = 6>
+         template <typename T> class MapTraits =
+#if defined _DENSE_HASH_MAP_H_
+          GoogleDenseHashMapTraits
+#elif defined _HASH_MAP
+          GnuHashMapTraits
+#else
+          StdMapTraits
+#endif
+         ,int OuterPacketBits = 6>
 class RandomSetter
 {
     typedef typename ei_traits<SparseMatrixType>::Scalar Scalar;
@@ -78,11 +115,13 @@ class RandomSetter
       ScalarWrapper() : value(0) {}
       Scalar value;
     };
-    typedef typename HashMapTraits<ScalarWrapper>::KeyType KeyType;
-    typedef typename HashMapTraits<ScalarWrapper>::Type HashMapType;
+    typedef typename MapTraits<ScalarWrapper>::KeyType KeyType;
+    typedef typename MapTraits<ScalarWrapper>::Type HashMapType;
     static const int OuterPacketMask = (1 << OuterPacketBits) - 1;
     enum {
-      RowMajor = SparseMatrixType::Flags & RowMajorBit
+      SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
+      TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
+      SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor
     };
 
   public:
@@ -90,31 +129,114 @@ class RandomSetter
     inline RandomSetter(SparseMatrixType& target)
       : mp_target(&target)
     {
-      m_outerPackets = target.outerSize() >> OuterPacketBits;
-      if (target.outerSize()&OuterPacketMask)
+      const int outerSize = SwapStorage ? target.innerSize() : target.outerSize();
+      const int innerSize = SwapStorage ? target.outerSize() : target.innerSize();
+      m_outerPackets = outerSize >> OuterPacketBits;
+      if (outerSize&OuterPacketMask)
         m_outerPackets += 1;
       m_hashmaps = new HashMapType[m_outerPackets];
-      KeyType ik = (1<<OuterPacketBits)*mp_target->innerSize()+1;
+      // compute number of bits needed to store inner indices
+      int aux = innerSize - 1;
+      m_keyBitsOffset = 0;
+      while (aux)
+      {
+        m_keyBitsOffset++;
+        aux = aux >> 1;
+      }
+      KeyType ik = (1<<(OuterPacketBits+m_keyBitsOffset));
       for (int k=0; k<m_outerPackets; ++k)
-        HashMapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
+        MapTraits<ScalarWrapper>::setInvalidKey(m_hashmaps[k],ik);
+
+      // insert current coeffs
+      for (int j=0; j<mp_target->outerSize(); ++j)
+        for (typename SparseMatrixType::InnerIterator it(*mp_target,j); it; ++it)
+          (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
     }
 
     ~RandomSetter()
     {
+      KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
+      if (!SwapStorage) // also means the map is sorted
+      {
+        mp_target->startFill(nonZeros());
+        for (int k=0; k<m_outerPackets; ++k)
+        {
+          const int outerOffset = (1<<OuterPacketBits) * k;
+          typename HashMapType::iterator end = m_hashmaps[k].end();
+          for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
+          {
+            const int outer = (it->first >> m_keyBitsOffset) + outerOffset;
+            const int inner = it->first & keyBitsMask;
+            mp_target->fill(TargetRowMajor ? outer : inner, TargetRowMajor ? inner : outer) = it->second.value;
+          }
+        }
+        mp_target->endFill();
+      }
+      else
+      {
+        VectorXi positions(mp_target->outerSize());
+        positions.setZero();
+        // pass 1
+        for (int k=0; k<m_outerPackets; ++k)
+        {
+          typename HashMapType::iterator end = m_hashmaps[k].end();
+          for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
+          {
+            const int outer = it->first & keyBitsMask;
+            positions[outer]++;
+          }
+        }
+        // prefix sum
+        int count = 0;
+        for (int j=0; j<mp_target->outerSize(); ++j)
+        {
+          int tmp = positions[j];
+          mp_target->_outerIndexPtr()[j] = count;
+          positions[j] = count;
+          count += tmp;
+        }
+        mp_target->_outerIndexPtr()[mp_target->outerSize()] = count;
+        mp_target->resizeNonZeros(count);
+        // pass 2
+        for (int k=0; k<m_outerPackets; ++k)
+        {
+          const int outerOffset = (1<<OuterPacketBits) * k;
+          typename HashMapType::iterator end = m_hashmaps[k].end();
+          for (typename HashMapType::iterator it = m_hashmaps[k].begin(); it!=end; ++it)
+          {
+            const int inner = (it->first >> m_keyBitsOffset) + outerOffset;
+            const int outer = it->first & keyBitsMask;
+            // sorted insertion
+            // Note that we have to deal with at most 2^OuterPacketBits unsorted coefficients,
+            // moreover those 2^OuterPacketBits coeffs are likely to be sparse, an so only a
+            // small fraction of them have to be sorted, whence the following simple procedure:
+            int posStart = mp_target->_outerIndexPtr()[outer];
+            int i = (positions[outer]++) - 1;
+            while ( (i >= posStart) && (mp_target->_innerIndexPtr()[i] > inner) )
+            {
+              mp_target->_valuePtr()[i+1] = mp_target->_valuePtr()[i];
+              mp_target->_innerIndexPtr()[i+1] = mp_target->_innerIndexPtr()[i];
+              --i;
+            }
+            mp_target->_innerIndexPtr()[i+1] = inner;
+            mp_target->_valuePtr()[i+1] = it->second.value;
+          }
+        }
+      }
       delete[] m_hashmaps;
     }
 
     Scalar& operator() (int row, int col)
     {
-      const int outer = RowMajor ? row : col;
-      const int inner = RowMajor ? col : row;
-      const int outerMajor = outer >> OuterPacketBits;
-      const int outerMinor = outer & OuterPacketMask;
-      const KeyType key = inner + outerMinor * mp_target->innerSize();
-
+      const int outer = SetterRowMajor ? row : col;
+      const int inner = SetterRowMajor ? col : row;
+      const int outerMajor = outer >> OuterPacketBits; // index of the packet/map
+      const int outerMinor = outer & OuterPacketMask;  // index of the inner vector in the packet
+      const KeyType key = (KeyType(outerMinor)<<m_keyBitsOffset) | inner;
       return m_hashmaps[outerMajor][key].value;
     }
 
+    // might be slow
     int nonZeros() const
     {
       int nz = 0;
@@ -129,6 +251,7 @@ class RandomSetter
     HashMapType* m_hashmaps;
     SparseMatrixType* mp_target;
     int m_outerPackets;
+    unsigned char m_keyBitsOffset;
 };
 
 #endif // EIGEN_RANDOMSETTER_H