bugfix in Map by Keir Mierle

Eigen/src/Core/Map.h

@@ -85,7 +85,7 @@ template<typename MatrixType, int PacketAccess> class Map
       EIGEN_ONLY_USED_FOR_DEBUG(rows);
       EIGEN_ONLY_USED_FOR_DEBUG(cols);
       ei_assert(rows == this->rows());
-      ei_assert(rows == this->cols());
+      ei_assert(cols == this->cols());
     }
 
     inline void resize(int size)

Eigen/src/Sparse/RandomSetter.h

@@ -25,6 +25,10 @@
 #ifndef EIGEN_RANDOMSETTER_H
 #define EIGEN_RANDOMSETTER_H
 
+/** Represents a std::map
+  *
+  * \see RandomSetter
+  */
 template<typename Scalar> struct StdMapTraits
 {
   typedef int KeyType;
@@ -37,6 +41,10 @@ template<typename Scalar> struct StdMapTraits
 };
 
 #ifdef _HASH_MAP
+/** Represents a __gnu_cxx::hash_map
+  *
+  * \see RandomSetter
+  */
 template<typename Scalar> struct GnuHashMapTraits
 {
   typedef int KeyType;
@@ -50,6 +58,10 @@ template<typename Scalar> struct GnuHashMapTraits
 #endif
 
 #ifdef _DENSE_HASH_MAP_H_
+/** Represents a google::dense_hash_map
+  *
+  * \see RandomSetter
+  */
 template<typename Scalar> struct GoogleDenseHashMapTraits
 {
   typedef int KeyType;
@@ -64,6 +76,10 @@ template<typename Scalar> struct GoogleDenseHashMapTraits
 #endif
 
 #ifdef _SPARSE_HASH_MAP_H_
+/** Represents a google::sparse_hash_map
+  *
+  * \see RandomSetter
+  */
 template<typename Scalar> struct GoogleSparseHashMapTraits
 {
   typedef int KeyType;
@@ -78,7 +94,19 @@ template<typename Scalar> struct GoogleSparseHashMapTraits
 
 /** \class RandomSetter
   *
-  * Typical usage:
+  * \brief The RandomSetter is a wrapper object allowing to set/update a sparse matrix with random access
+  *
+  * \param SparseMatrixType the type of the sparse matrix we are updating
+  * \param MapTraits a traits class representing the map implementation used for the temporary sparse storage.
+  *                  Its default value depends on the system.
+  * \param OuterPacketBits defines the number of rows (or columns) manage by a single map object
+  *                        as a power of two exponent.
+  *
+  * This class temporarily represents a sparse matrix object using a generic map implementation allowing for
+  * efficient random access. The conversion from the compressed representation to a hash_map object is performed
+  * in the RandomSetter constructor, while the sparse matrix is updated back at destruction time. This strategy
+  * suggest the use of nested blocks as in this example:
+  *
   * \code
   * SparseMatrix<double> m(rows,cols);
   * {
@@ -91,11 +119,28 @@ template<typename Scalar> struct GoogleSparseHashMapTraits
   * // and m is ready to use.
   * \endcode
   *
-  * \note for performance and memory consumption reasons it is highly recommended to use
+  * Since hash_map objects are not fully sorted, representing a full matrix as a single hash_map would
-  * Google's hash library. To do so you have two options:
+  * involve a big and costly sort to update the compressed matrix back. To overcome this issue, a RandomSetter
-  *  - include <google/dense_hash_map> yourself \b before Eigen/Sparse header
+  * use multiple hash_map, each representing 2^OuterPacketBits columns or rows according to the storage order.
+  * To reach optimal performance, this value should be adjusted according to the average number of nonzeros
+  * per rows/columns.
+  *
+  * The possible values for the template parameter MapTraits are:
+  *  - \b StdMapTraits: corresponds to std::map. (does not perform very well)
+  *  - \b GnuHashMapTraits: corresponds to __gnu_cxx::hash_map (available only with GCC)
+  *  - \b GoogleDenseHashMapTraits: corresponds to google::dense_hash_map (best efficiency, reasonable memory consumption)
+  *  - \b GoogleSparseHashMapTraits: corresponds to google::sparse_hash_map (best memory consumption, relatively good performance)
+  *
+  * The default map implementation depends on the availability, and the preferred order is:
+  * GoogleSparseHashMapTraits, GnuHashMapTraits, and finally StdMapTraits.
+  *
+  * For performance and memory consumption reasons it is highly recommended to use one of
+  * the Google's hash_map implementation. To enable the support for them, 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.
+  * 
+  * \see http://code.google.com/p/google-sparsehash/
   */
 template<typename SparseMatrixType,
          template <typename T> class MapTraits =
@@ -121,11 +166,19 @@ class RandomSetter
     enum {
       SwapStorage = 1 - MapTraits<ScalarWrapper>::IsSorted,
       TargetRowMajor = (SparseMatrixType::Flags & RowMajorBit) ? 1 : 0,
-      SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor
+      SetterRowMajor = SwapStorage ? 1-TargetRowMajor : TargetRowMajor,
+      IsUpperTriangular = SparseMatrixType::Flags & UpperTriangularBit,
+      IsLowerTriangular = SparseMatrixType::Flags & LowerTriangularBit
     };
 
   public:
 
+    /** Constructs a random setter object from the sparse matrix \a target
+      *
+      * Note that the initial value of \a target are imported. If you want to re-set
+      * a sparse matrix from scratch, then you must set it to zero first using the
+      * setZero() function.
+      */
     inline RandomSetter(SparseMatrixType& target)
       : mp_target(&target)
     {
@@ -153,6 +206,7 @@ class RandomSetter
           (*this)(TargetRowMajor?j:it.index(), TargetRowMajor?it.index():j) = it.value();
     }
 
+    /** Destructor updating back the sparse matrix target */
     ~RandomSetter()
     {
       KeyType keyBitsMask = (1<<m_keyBitsOffset)-1;
@@ -226,8 +280,11 @@ class RandomSetter
       delete[] m_hashmaps;
     }
 
+    /** \returns a reference to the coefficient at given coordinates \a row, \a col */
     Scalar& operator() (int row, int col)
     {
+      ei_assert(((!IsUpperTriangular) || (row<=col)) && "Invalid access to an upper triangular matrix");
+      ei_assert(((!IsLowerTriangular) || (col<=row)) && "Invalid access to an upper triangular matrix");
       const int outer = SetterRowMajor ? row : col;
       const int inner = SetterRowMajor ? col : row;
       const int outerMajor = outer >> OuterPacketBits; // index of the packet/map
@@ -236,7 +293,11 @@ class RandomSetter
       return m_hashmaps[outerMajor][key].value;
     }
 
-    // might be slow
+    /** \returns the number of non zero coefficients 
+      *
+      * \note According to the underlying map/hash_map implementation,
+      * this function might be quite expensive.
+      */
     int nonZeros() const
     {
       int nz = 0;

bench/sparse_setter.cpp

@@ -4,7 +4,7 @@
 // -DNOGMM -DNOMTL -DCSPARSE
 // -I /home/gael/Coding/LinearAlgebra/CSparse/Include/ /home/gael/Coding/LinearAlgebra/CSparse/Lib/libcsparse.a
 #ifndef SIZE
-#define SIZE 1000000
+#define SIZE 100000
 #endif
 
 #ifndef NBPERROW
@@ -22,6 +22,8 @@
 
 #include "BenchSparseUtil.h"
 
+#define CHECK_MEM 
+// #define CHECK_MEM  std/**/::cout << "check mem\n"; getchar();
 
 #define BENCH(X) \
   timer.reset(); \
@@ -34,6 +36,7 @@
 typedef std::vector<Vector2i> Coordinates;
 typedef std::vector<float> Values;
 
+EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals);
 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals);
 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords, const Values& vals);
 EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords, const Values& vals);
@@ -47,17 +50,31 @@ int main(int argc, char *argv[])
 {
   int rows = SIZE;
   int cols = SIZE;
+  bool fullyrand = false;
   //float density = float(NBPERROW)/float(SIZE);
   
   BenchTimer timer;
   Coordinates coords;
   Values values;
+  if(fullyrand)
+  {
     for (int i=0; i<cols*NBPERROW; ++i)
     {
       coords.push_back(Vector2i(ei_random<int>(0,rows-1),ei_random<int>(0,cols-1)));
       values.push_back(ei_random<Scalar>());
     }
+  }
+  else
+  {
+    for (int j=0; j<cols; ++j)
+    for (int i=0; i<NBPERROW; ++i)
+    {
+      coords.push_back(Vector2i(ei_random<int>(0,rows-1),j));
+      values.push_back(ei_random<Scalar>());
+    }
+  }
   std::cout << "nnz = " << coords.size()  << "\n";
+  CHECK_MEM
 
     // dense matrices
     #ifdef DENSEMATRIX
@@ -72,6 +89,15 @@ int main(int argc, char *argv[])
     #endif
 
     // eigen sparse matrices
+    if (!fullyrand)
+    {
+      timer.reset();
+      timer.start();
+      for (int k=0; k<REPEAT; ++k)
+        setinnerrand_eigen(coords,values);
+      timer.stop();
+      std::cout << "Eigen fillrand\t" << timer.value() << "\n";
+    }
     {
       timer.reset();
       timer.start();
@@ -150,6 +176,20 @@ int main(int argc, char *argv[])
   return 0;
 }
 
+EIGEN_DONT_INLINE Scalar* setinnerrand_eigen(const Coordinates& coords, const Values& vals)
+{
+  using namespace Eigen;
+  SparseMatrix<Scalar> mat(SIZE,SIZE);
+  mat.startFill(2000000/*coords.size()*/);
+  for (int i=0; i<coords.size(); ++i)
+  {
+    mat.fillrand(coords[i].x(), coords[i].y()) = vals[i];
+  }
+  mat.endFill();
+  CHECK_MEM;
+  return 0;
+}
+
 EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, const Values& vals)
 {
   using namespace Eigen;
@@ -160,7 +200,7 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_gnu_hash(const Coordinates& coords, cons
     {
       setter(coords[i].x(), coords[i].y()) = vals[i];
     }
-//     std::cout << "check mem\n"; getchar();
+    CHECK_MEM;
   }
   return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
 }
@@ -174,7 +214,7 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_google_dense(const Coordinates& coords,
     RandomSetter<SparseMatrix<Scalar>, GoogleDenseHashMapTraits> setter(mat);
     for (int i=0; i<coords.size(); ++i)
       setter(coords[i].x(), coords[i].y()) = vals[i];
-//     std::cout << "check mem\n"; getchar();
+    CHECK_MEM;
   }
   return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
 }
@@ -187,7 +227,7 @@ EIGEN_DONT_INLINE Scalar* setrand_eigen_google_sparse(const Coordinates& coords,
     RandomSetter<SparseMatrix<Scalar>, GoogleSparseHashMapTraits> setter(mat);
     for (int i=0; i<coords.size(); ++i)
       setter(coords[i].x(), coords[i].y()) = vals[i];
-//     std::cout << "check mem\n"; getchar();
+    CHECK_MEM;
   }
   return 0;//&mat.coeffRef(coords[0].x(), coords[0].y());
 }
@@ -204,7 +244,7 @@ EIGEN_DONT_INLINE Scalar* setrand_ublas_mapped(const Coordinates& coords, const
   {
     aux(coords[i].x(), coords[i].y()) = vals[i];
   }
-//   std::cout << "check mem\n"; getchar();
+  CHECK_MEM;
   compressed_matrix<Scalar> mat(aux);
   return 0;// &mat(coords[0].x(), coords[0].y());
 }
@@ -245,6 +285,7 @@ EIGEN_DONT_INLINE Scalar* setrand_ublas_genvec(const Coordinates& coords, const
   {
     aux(coords[i].x(), coords[i].y()) = vals[i];
   }
+  CHECK_MEM;
   compressed_matrix<Scalar,row_major> mat(aux);
   return 0;//&mat(coords[0].x(), coords[0].y());
 }