Rework sparse repmat (#1819)

* Add test repmat function

* Rework repmat function

* Add repmat test for colMajor sparse matrices

* derived majorType from params

* missing templates

Co-authored-by: Paul Rötzer <[email protected]>
Co-authored-by: Alec Jacobson <[email protected]>
Co-authored-by: Alec Jacobson <[email protected]>

include/igl/repmat.cpp

@@ -29,39 +29,51 @@ IGL_INLINE void igl::repmat(
   }
 }
 
-template <typename T>
+template <typename T, int majorType>
 IGL_INLINE void igl::repmat(
-  const Eigen::SparseMatrix<T> & A,
+  const Eigen::SparseMatrix<T, majorType> & A,
   const int r,
   const int c,
-  Eigen::SparseMatrix<T> & B)
+  Eigen::SparseMatrix<T, majorType> & B)
 {
   assert(r>0);
   assert(c>0);
-  B.resize(r*A.rows(),c*A.cols());
-  B.reserve(r*c*A.nonZeros());
-  for(int i = 0;i<r;i++)
+  B.resize(r*A.rows(), c*A.cols());
+  std::vector<Eigen::Triplet<T>> b;
+  b.reserve(r*c*A.nonZeros());
+
+  for(int i = 0; i < r; i++)
   {
-    for(int j = 0;j<c;j++)
+    for(int j = 0; j < c; j++)
     {
-      // Loop outer level
-      for (int k=0; k<A.outerSize(); ++k)
+      // loop outer level
+      for (int k = 0; k < A.outerSize(); ++k)
       {
         // loop inner level
-        for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
+        for (typename Eigen::SparseMatrix<T, majorType>::InnerIterator 
+          it(A,k); it; ++it)
         {
-          B.insert(i*A.rows()+it.row(),j*A.cols() + it.col()) = it.value();
+          Eigen::Triplet<T> triplet(i * A.rows() + it.row(), j * A.cols() 
+            + it.col(), it.value());
+          b.push_back(triplet);
         }
       }
     }
   }
-  B.finalize();
+    B.setFromTriplets(b.begin(), b.end());
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::repmat<double, 0>(Eigen::SparseMatrix<double, 0, int> const&, int, int, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
 template void igl::repmat<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::repmat<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::repmat<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::repmat<double, 1>(Eigen::SparseMatrix<double, 1, int> const&, int, int, Eigen::SparseMatrix<double, 1, int>&);
+template void igl::repmat<int, 1>(Eigen::SparseMatrix<int, 1, int> const&, int, int, Eigen::SparseMatrix<int, 1, int>&);
+template void igl::repmat<float, 1>(Eigen::SparseMatrix<float, 1, int> const&, int, int, Eigen::SparseMatrix<float, 1, int>&);
+template void igl::repmat<int, 0>(Eigen::SparseMatrix<int, 0, int> const&, int, int, Eigen::SparseMatrix<int, 0, int>&);
+template void igl::repmat<float, 0>(Eigen::SparseMatrix<float, 0, int> const&, int, int, Eigen::SparseMatrix<float, 0, int>&);
 #endif

include/igl/repmat.h

@@ -38,12 +38,12 @@ namespace igl
     const int r,
     const int c,
     Eigen::PlainObjectBase<DerivedB> & B);
-  template <typename T>
+  template <typename T, int majorType>
   IGL_INLINE void repmat(
-    const Eigen::SparseMatrix<T> & A,
+    const Eigen::SparseMatrix<T, majorType> & A,
     const int r,
     const int c,
-    Eigen::SparseMatrix<T> & B);
+    Eigen::SparseMatrix<T, majorType> & B);
 }
 
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/repmat.cpp

@@ -0,0 +1,61 @@
+#include <test_common.h>
+#include <Eigen/Sparse>
+#include <igl/repmat.h>
+
+
+template <typename T, int majorType>
+void checkRepmat(
+  Eigen::SparseMatrix<T, majorType> &A,
+  Eigen::SparseMatrix<T, majorType> &B,
+  const int rows,
+  const int cols,
+  const int r,
+  const int c)
+{
+  for (int i = 0; i < rows; i++)
+  {
+    for (int j = 0; j < cols; j++)
+    {
+      for (int ii = 0; ii < r; ii++)
+      {
+        for (int jj = 0; jj < c; jj++)
+        {
+          REQUIRE (A.coeff(i, j) == B.coeff(i + ii * rows, j + jj * cols));
+        }
+      }
+    }
+  }
+}
+
+template <typename T, int majorType>
+void testRepmat(
+  const int rows,
+  const int cols,
+  const int r,
+  const int c)
+{
+  Eigen::SparseMatrix<T, majorType> A;
+  A = (Eigen::Matrix<T, Eigen::Dynamic,
+    Eigen::Dynamic>().setRandom(rows, cols).sparseView());
+  Eigen::SparseMatrix<T, majorType> B;
+  igl::repmat(A, r, c, B);
+
+  REQUIRE (B.rows() == r * rows);
+  REQUIRE (B.cols() == c * cols);
+  checkRepmat<T, majorType>(A, B, rows, cols, r, c);
+}
+
+TEST_CASE("repmat: sparse rowMajor", "[igl]")
+{
+  testRepmat<double, Eigen::RowMajor>(4,  5, 2, 4);
+  testRepmat<   int, Eigen::RowMajor>(2,  8, 3, 4);
+  testRepmat< float, Eigen::RowMajor>(6, 10, 2, 2);
+}
+
+TEST_CASE("repmat: sparse colMajor", "[igl]")
+{
+  testRepmat<double, Eigen::ColMajor>(4,  5, 3, 5);
+  testRepmat<   int, Eigen::ColMajor>(2,  8, 3, 5);
+  testRepmat< float, Eigen::ColMajor>(6, 10, 3, 5);
+}
+