faster cat

Former-commit-id: 56bdf2c2a4a6aab2cfe7e3c3c2b55f7a6d8b4516

include/igl/cat.cpp

@@ -6,12 +6,14 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cat.h"
+
 #include <cstdio>
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
 #include <unsupported/Eigen/SparseExtra>
 
+
 // Sparse matrices need to be handled carefully. Because C++ does not 
 // Template:
 //   Scalar  sparse matrix scalar type, e.g. double
@@ -22,6 +24,7 @@ IGL_INLINE void igl::cat(
     const Eigen::SparseMatrix<Scalar> & B, 
     Eigen::SparseMatrix<Scalar> & C)
 {
+
   assert(dim == 1 || dim == 2);
   using namespace Eigen;
   // Special case if B or A is empty
@@ -36,6 +39,7 @@ IGL_INLINE void igl::cat(
     return;
   }
 
+#if false
   // This **must** be DynamicSparseMatrix, otherwise this implementation is
   // insanely slow
   DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
@@ -77,6 +81,109 @@ IGL_INLINE void igl::cat(
   }
 
   C = SparseMatrix<Scalar>(dyn_C);
+#elif false
+  std::vector<Triplet<Scalar> > CIJV;
+  CIJV.reserve(A.nonZeros() + B.nonZeros());
+  {
+    // Iterate over outside of A
+    for(int k=0; k<A.outerSize(); ++k)
+    {
+      // Iterate over inside
+      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+      {
+        CIJV.emplace_back(it.row(),it.col(),it.value());
+      }
+    }
+    // Iterate over outside of B
+    for(int k=0; k<B.outerSize(); ++k)
+    {
+      // Iterate over inside
+      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+      {
+        int r = (dim == 1 ? A.rows()+it.row() : it.row());
+        int c = (dim == 2 ? A.cols()+it.col() : it.col());
+        CIJV.emplace_back(r,c,it.value());
+      }
+    }
+
+  }
+
+  C = SparseMatrix<Scalar>( 
+      dim == 1 ? A.rows()+B.rows() : A.rows(),
+      dim == 1 ? A.cols()          : A.cols()+B.cols());
+  C.reserve(A.nonZeros() + B.nonZeros());
+  C.setFromTriplets(CIJV.begin(),CIJV.end());
+#else
+  C = SparseMatrix<Scalar>( 
+      dim == 1 ? A.rows()+B.rows() : A.rows(),
+      dim == 1 ? A.cols()          : A.cols()+B.cols());
+  Eigen::VectorXi per_col = Eigen::VectorXi::Zero(C.cols());
+  if(dim == 1)
+  {
+    assert(A.outerSize() == B.outerSize());
+    for(int k = 0;k<A.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+      {
+        per_col(k)++;
+      }
+      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+      {
+        per_col(k)++;
+      }
+    }
+  }else
+  {
+    for(int k = 0;k<A.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+      {
+        per_col(k)++;
+      }
+    }
+    for(int k = 0;k<B.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+      {
+        per_col(A.cols() + k)++;
+      }
+    }
+  }
+  C.reserve(per_col);
+  if(dim == 1)
+  {
+    for(int k = 0;k<A.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+      {
+        C.insert(it.row(),k) = it.value();
+      }
+      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+      {
+        C.insert(A.rows()+it.row(),k) = it.value();
+      }
+    }
+  }else
+  {
+    for(int k = 0;k<A.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+      {
+        C.insert(it.row(),k) = it.value();
+      }
+    }
+    for(int k = 0;k<B.outerSize();++k)
+    {
+      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
+      {
+        C.insert(it.row(),A.cols()+k) = it.value();
+      }
+    }
+  }
+  C.makeCompressed();
+
+#endif
+
 }
 
 template <typename Derived, class MatC>