Merge pull request #667 from danielepanozzo/fast-slim

Improved slim

Former-commit-id: 7d5a907ff65bac9ba9fd7e8440572328d5eb6986

include/igl/AtA_cached.cpp

@@ -0,0 +1,130 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// 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 "AtA_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    igl::AtA_cached_data& data)
+{
+  // 1 Compute At (this could be avoided, but performance-wise it will not make a difference)
+  std::vector<std::vector<int> > Col_RowPtr;
+  std::vector<std::vector<int> > Col_IndexPtr;
+
+  Col_RowPtr.resize(A.cols());
+  Col_IndexPtr.resize(A.cols());
+
+  for (unsigned k=0; k<A.outerSize(); ++k)
+  {
+    unsigned outer_index = *(A.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == A.outerSize()) ? A.nonZeros() : *(A.outerIndexPtr()+k+1); 
+    
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(A.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < A.cols());
+      assert(col >= 0);
+      assert(row < A.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < A.nonZeros());
+
+      Col_RowPtr[col].push_back(row);
+      Col_IndexPtr[col].push_back(value_index);
+    }
+  }
+
+  Eigen::SparseMatrix<Scalar> At = A.transpose();
+  At.makeCompressed();
+  AtA = At * A;
+  AtA.makeCompressed();
+
+  assert(AtA.isCompressed());
+
+  // If weights are not provided, use 1
+  if (data.W.size() == 0)
+    data.W = Eigen::VectorXd::Ones(A.rows());
+  assert(data.W.size() == A.rows());
+
+  data.I_outer.reserve(AtA.outerSize());
+  data.I_row.reserve(2*AtA.nonZeros());
+  data.I_col.reserve(2*AtA.nonZeros());
+  data.I_w.reserve(2*AtA.nonZeros());
+
+  // 2 Construct the rules
+  for (unsigned k=0; k<AtA.outerSize(); ++k)
+  {
+    unsigned outer_index = *(AtA.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == AtA.outerSize()) ? AtA.nonZeros() : *(AtA.outerIndexPtr()+k+1); 
+    
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(AtA.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < AtA.cols());
+      assert(col >= 0);
+      assert(row < AtA.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < AtA.nonZeros());
+
+      data.I_outer.push_back(data.I_row.size());
+
+      // Find correspondences
+      unsigned i=0;
+      unsigned j=0;
+      while (i<Col_RowPtr[row].size() && j<Col_RowPtr[col].size())
+      {
+          if (Col_RowPtr[row][i] == Col_RowPtr[col][j])
+          {
+            data.I_row.push_back(Col_IndexPtr[row][i]);
+            data.I_col.push_back(Col_IndexPtr[col][j]);
+            data.I_w.push_back(Col_RowPtr[col][j]);
+            ++i;
+            ++j;
+          } else 
+          if (Col_RowPtr[row][i] > Col_RowPtr[col][j])
+            ++j;
+          else
+            ++i;
+
+      }
+    }
+  }
+  data.I_outer.push_back(data.I_row.size()); // makes it more efficient to iterate later on
+
+  igl::AtA_cached(A,AtA,data);
+}
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const igl::AtA_cached_data& data)
+{
+  for (unsigned i=0; i<data.I_outer.size()-1; ++i)
+  {
+    *(AtA.valuePtr() + i) = 0;
+    for (unsigned j=data.I_outer[i]; j<data.I_outer[i+1]; ++j)
+      *(AtA.valuePtr() + i) += *(A.valuePtr() + data.I_row[j]) * data.W[data.I_w[j]] * *(A.valuePtr() + data.I_col[j]);
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::AtA_cached<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data const&);
+template void igl::AtA_cached_precompute<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data&);
+#endif

include/igl/AtA_cached.h

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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/.
+#ifndef IGL_ATA_CACHED_H
+#define IGL_ATA_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+  struct AtA_cached_data
+  {
+    // Weights
+    Eigen::VectorXd W;
+
+    // Flatten composition rules
+    std::vector<int> I_row;
+    std::vector<int> I_col;
+    std::vector<int> I_w;
+
+    // For each entry of AtA, points to the beginning
+    // of the composition rules
+    std::vector<int> I_outer;
+  };
+
+  // Computes At * W * A, where A is sparse and W is diagonal. Divides the 
+  // construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values. Compared to
+  // evaluating it directly, this version is slower for the first time (since it requires a
+  // precomputation), but faster to the subsequent evaluations.
+  //
+  // Input:
+  //   A m x n sparse matrix
+  //   data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
+  // Outputs:
+  //   AtA  m by m matrix computed as AtA * W * A
+  //
+  // Example:
+  // AtA_data = igl::AtA_cached_data();
+  // AtA_data.W = W;
+  // if (s.AtA.rows() == 0)
+  //   igl::AtA_cached_precompute(s.A,s.AtA,s.AtA_data);
+  // else
+  //   igl::AtA_cached(s.A,s.AtA,s.AtA_data);
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    AtA_cached_data& data);
+
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const AtA_cached_data& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "AtA_cached.cpp"
+#endif
+
+#endif

include/igl/slice.cpp

@@ -350,6 +350,7 @@ template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<in
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
 template void igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::slice<unsigned int, unsigned int>(Eigen::SparseMatrix<unsigned int, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<unsigned int, 0, int>&);
 #ifdef WIN32
 template void igl::slice<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>>>(class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> &);
 template void igl::slice<class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);

include/igl/slice_cached.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// 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 "slice_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <igl/slice.h>
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached_precompute(
+  const Eigen::SparseMatrix<TX>& X,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+  Eigen::SparseMatrix<TY>& Y,
+  Eigen::VectorXi& data)
+{
+  // Create a sparse matrix whose entries are the ids
+  Eigen::SparseMatrix<unsigned> TS = X.template cast<unsigned>();
+
+  TS.makeCompressed();
+  for (unsigned i=0;i<TS.nonZeros();++i)
+    *(TS.valuePtr() + i) = i;
+
+  Eigen::SparseMatrix<unsigned> TS_sliced;
+  igl::slice(TS,R,C,TS_sliced);
+  Y = TS_sliced.cast<TY>();
+
+  data.resize(TS_sliced.nonZeros());
+  for (unsigned i=0;i<data.size();++i)
+  {
+    data[i] = *(TS_sliced.valuePtr() + i);
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+  }
+}
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached(
+  const Eigen::SparseMatrix<TX>& X,
+  Eigen::SparseMatrix<TY>& Y,
+  const Eigen::VectorXi& data)
+{
+  for (unsigned i=0; i<data.size(); ++i)
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::slice_cached_precompute<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::slice_cached<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+#endif

include/igl/slice_cached.h

@@ -0,0 +1,66 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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/.
+#ifndef IGL_SLICE_CACHED_H
+#define IGL_SLICE_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+
+  // Act like the matlab X(row_indices,col_indices) operator, where
+  // row_indices, col_indices are non-negative integer indices. This is a fast version
+  // of igl::slice that can analyze and store the sparsity structure. It is slower at the // first evaluation (slice_cached_precompute), but faster on the subsequent ones.
+  // 
+  // Inputs:
+  //   X  m by n matrix
+  //   R  list of row indices
+  //   C  list of column indices
+  //   
+  // Output:
+  //   Y  #R by #C matrix
+  //   data Temporary data used by slice_cached to repeat this operation
+  //
+  // Usage:
+  //
+  // // Construct and slice up Laplacian
+  // SparseMatrix<double> L,L_sliced;
+  // igl::cotmatrix(V,F,L);
+
+  // // Normal igl::slice call
+  // igl::slice(L,in,in,L_in_in);
+
+  // // Fast version
+  // static VectorXi data; // static or saved in a global state
+  // if (data.size() == 0)
+  //   igl::slice_cached_precompute(L,in,in,L_sliced,data);
+  // else
+  //   igl::slice_cached(L,L_sliced,temp);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached_precompute(
+    const Eigen::SparseMatrix<TX>& X,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+    Eigen::SparseMatrix<TY>& Y,
+    Eigen::VectorXi& data);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached(
+    const Eigen::SparseMatrix<TX>& X,
+    Eigen::SparseMatrix<TY>& Y,
+    const Eigen::VectorXi& data);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slice_cached.cpp"
+#endif
+
+#endif

include/igl/slim.cpp

@@ -34,6 +34,14 @@
 #include <Eigen/SparseCholesky>
 #include <Eigen/IterativeLinearSolvers>
 
+#include <igl/Timer.h>
+#include <igl/sparse_cached.h>
+#include <igl/AtA_cached.h>
+
+#ifdef CHOLMOD
+#include <Eigen/CholmodSupport>
+#endif
+
 namespace igl
 {
   namespace slim
@@ -42,8 +50,8 @@ namespace igl
     IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
                                                     const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2,
                                                     Eigen::SparseMatrix<double> &D1, Eigen::SparseMatrix<double> &D2);
-    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A);
-    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At);
+    IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV);
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A);
     IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
     IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
     IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
@@ -395,8 +403,12 @@ namespace igl
       Eigen::SparseMatrix<double> L;
       build_linear_system(s,L);
 
+      igl::Timer t;
+      
+      //t.start();
       // solve
       Eigen::VectorXd Uc;
+#ifndef CHOLMOD
       if (s.dim == 2)
       {
         SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
@@ -406,13 +418,21 @@ namespace igl
       { // seems like CG performs much worse for 2D and way better for 3D
         Eigen::VectorXd guess(uv.rows() * s.dim);
         for (int i = 0; i < s.v_num; i++) for (int j = 0; j < s.dim; j++) guess(uv.rows() * j + i) = uv(i, j); // flatten vector
-        ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower | Upper> solver;
-        solver.setTolerance(1e-8);
-        Uc = solver.compute(L).solveWithGuess(s.rhs, guess);
+        ConjugateGradient<Eigen::SparseMatrix<double>, Lower | Upper> cg;
+        cg.setTolerance(1e-8);
+        cg.compute(L);
+        Uc = cg.solveWithGuess(s.rhs, guess);
       }
-
+#else
+        CholmodSimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+        Uc = solver.compute(L).solve(s.rhs);
+#endif
       for (int i = 0; i < s.dim; i++)
         uv.col(i) = Uc.block(i * s.v_n, 0, s.v_n, 1);
+
+      // t.stop();
+      // std::cerr << "solve: " << t.getElapsedTime() << std::endl;
+
     }
 
 
@@ -478,20 +498,60 @@ namespace igl
     IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
     {
       // formula (35) in paper
+      std::vector<Eigen::Triplet<double> > IJV;
+      
+      #ifdef SLIM_CACHED
+      buildA(s,IJV);
+      if (s.A.rows() == 0)
+      {
+        s.A = Eigen::SparseMatrix<double>(s.dim * s.dim * s.f_n, s.dim * s.v_n);
+        igl::sparse_cached_precompute(IJV,s.A,s.A_data);
+      }
+      else
+        igl::sparse_cached(IJV,s.A,s.A_data);
+      #else
       Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
-      buildA(s,A);
+      buildA(s,IJV);
+      A.setFromTriplets(IJV.begin(),IJV.end());
+      A.makeCompressed();
+      #endif
 
+      #ifdef SLIM_CACHED
+      #else
       Eigen::SparseMatrix<double> At = A.transpose();
       At.makeCompressed();
+      #endif
+
+      #ifdef SLIM_CACHED
+      Eigen::SparseMatrix<double> id_m(s.A.cols(), s.A.cols());
+      #else
+      Eigen::SparseMatrix<double> id_m(A.cols(), A.cols());
+      #endif
 
-      Eigen::SparseMatrix<double> id_m(At.rows(), At.rows());
       id_m.setIdentity();
 
       // add proximal penalty
+      #ifdef SLIM_CACHED
+      s.AtA_data.W = s.WGL_M;
+      if (s.AtA.rows() == 0)
+        igl::AtA_cached_precompute(s.A,s.AtA,s.AtA_data);
+      else
+        igl::AtA_cached(s.A,s.AtA,s.AtA_data);
+
+      L = s.AtA + s.proximal_p * id_m; //add also a proximal 
+      L.makeCompressed();
+
+      #else
       L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; //add also a proximal term
       L.makeCompressed();
+      #endif
+
+      #ifdef SLIM_CACHED
+      buildRhs(s, s.A);
+      #else
+      buildRhs(s, A);
+      #endif
 
-      buildRhs(s, At);
       Eigen::SparseMatrix<double> OldL = L;
       add_soft_constraints(s,L);
       L.makeCompressed();
@@ -657,10 +717,9 @@ namespace igl
       return energy;
     }
 
-    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A)
+    IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV)
     {
       // formula (35) in paper
-      std::vector<Eigen::Triplet<double> > IJV;
       if (s.dim == 2)
       {
         IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize()));
@@ -780,10 +839,9 @@ namespace igl
           }
         }
       }
-      A.setFromTriplets(IJV.begin(), IJV.end());
     }
 
-    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At)
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A)
     {
       Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n);
       f_rhs.setZero();
@@ -830,7 +888,7 @@ namespace igl
         for (int j = 0; j < s.v_n; j++)
           uv_flat(s.v_n * i + j) = s.V_o(j, i);
 
-      s.rhs = (At * s.WGL_M.asDiagonal() * f_rhs + s.proximal_p * uv_flat);
+      s.rhs = (f_rhs.transpose() * s.WGL_M.asDiagonal() * A).transpose() + s.proximal_p * uv_flat;
     }
 
   }

include/igl/slim.h

@@ -12,6 +12,14 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
+// This option makes the iterations faster (all except the first) by caching the 
+// sparsity pattern of the matrix involved in the assembly. It should be on if you plan to do many iterations, off if you have to change the matrix structure at every iteration.
+#define SLIM_CACHED 
+
+#ifdef SLIM_CACHED
+#include <igl/AtA_cached.h>
+#endif
+
 namespace igl
 {
 
@@ -64,6 +72,13 @@ struct SLIMData
   bool first_solve;
   bool has_pre_calc = false;
   int dim;
+
+  #ifdef SLIM_CACHED
+  Eigen::SparseMatrix<double> A;
+  Eigen::VectorXi A_data;
+  Eigen::SparseMatrix<double> AtA;
+  igl::AtA_cached_data AtA_data;
+  #endif
 };
 
 // Compute necessary information to start using SLIM

include/igl/sparse_cached.cpp

@@ -0,0 +1,122 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// 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 "sparse_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <array>
+#include <unordered_map>
+#include <map>
+#include <utility>
+
+template <typename DerivedI, typename Scalar>
+IGL_INLINE void igl::sparse_cached_precompute(
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedI> & J,
+  Eigen::SparseMatrix<Scalar>& X,
+  Eigen::VectorXi& data)
+{
+  // Generates the triplets
+  std::vector<Eigen::Triplet<double> > t(I.size());
+  for (unsigned i = 0; i<I.size(); ++i)
+    t[i] = Eigen::Triplet<double>(I[i],J[i],1);
+
+  // Call the triplets version
+  sparse_cached_precompute(t,X,data);
+}
+
+template <typename Scalar>
+IGL_INLINE void igl::sparse_cached_precompute(
+  const std::vector<Eigen::Triplet<Scalar> >& triplets,
+  Eigen::SparseMatrix<Scalar>& X,
+  Eigen::VectorXi& data)
+{
+    // Construct an empty sparse matrix
+    X.setFromTriplets(triplets.begin(),triplets.end());
+    X.makeCompressed();
+
+    std::vector<std::array<int,3> > T(triplets.size());
+    for (unsigned i=0; i<triplets.size(); ++i)
+    {
+      T[i][0] = triplets[i].col();
+      T[i][1] = triplets[i].row();
+      T[i][2] = i;
+    }
+
+    std::sort(T.begin(), T.end());
+
+    data.resize(triplets.size());
+
+    int t = 0;
+
+    for (unsigned k=0; k<X.outerSize(); ++k)
+    {
+      unsigned outer_index = *(X.outerIndexPtr()+k);
+      unsigned next_outer_index = (k+1 == X.outerSize()) ? X.nonZeros() : *(X.outerIndexPtr()+k+1); 
+      
+      for (unsigned l=outer_index; l<next_outer_index; ++l)
+      {
+        int col = k;
+        int row = *(X.innerIndexPtr()+l);
+        int value_index = l;
+        assert(col < X.cols());
+        assert(col >= 0);
+        assert(row < X.rows());
+        assert(row >= 0);
+        assert(value_index >= 0);
+        assert(value_index < X.nonZeros());
+
+        std::pair<int,int> p_m = std::make_pair(row,col);
+
+        while (t<T.size() && (p_m == std::make_pair(T[t][1],T[t][0])))
+          data[T[t++][2]] = value_index;
+      }
+    }
+    assert(t==T.size());
+
+}
+  
+template <typename Scalar>
+IGL_INLINE void igl::sparse_cached(
+  const std::vector<Eigen::Triplet<Scalar> >& triplets,
+  Eigen::SparseMatrix<Scalar>& X,
+  const Eigen::VectorXi& data)
+{
+  assert(triplets.size() == data.size());
+
+  // Clear it first
+  for (unsigned i = 0; i<data.size(); ++i)
+    *(X.valuePtr() + data[i]) = 0;
+ 
+  // Then sum them up
+  for (unsigned i = 0; i<triplets.size(); ++i)
+    *(X.valuePtr() + data[i]) += triplets[i].value();
+}
+
+template <typename DerivedV, typename Scalar>
+IGL_INLINE void igl::sparse_cached(
+  const Eigen::MatrixBase<DerivedV>& V,
+  Eigen::SparseMatrix<Scalar>& X,
+  const Eigen::VectorXi& data)
+{
+  assert(V.size() == data.size());
+
+  // Clear it first
+  for (unsigned i = 0; i<data.size(); ++i)
+    *(X.valuePtr() + data[i]) = 0;
+ 
+  // Then sum them up
+  for (unsigned i = 0; i<V.size(); ++i)
+    *(X.valuePtr() + data[i]) += V[i];
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::sparse_cached_precompute<double>(std::vector<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index>, std::allocator<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index> > > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::sparse_cached<double>(std::vector<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index>, std::allocator<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index> > > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+#endif

include/igl/sparse_cached.h

@@ -0,0 +1,82 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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/.
+#ifndef IGL_SPARSE_CACHED_H
+#define IGL_SPARSE_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Build a sparse matrix from list of indices and values (I,J,V), similarly to 
+  // the sparse function in matlab. Divides the construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values. Compared to
+  // igl::sparse, this version is slower for the first time (since it requires a
+  // precomputation), but faster to the subsequent evaluations.
+  //
+  // Templates:
+  //   IndexVector  list of indices, value should be non-negative and should
+  //     expect to be cast to an index. Must implement operator(i) to retrieve
+  //     ith element
+  //   ValueVector  list of values, value should be expect to be cast to type
+  //     T. Must implement operator(i) to retrieve ith element
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   I  nnz vector of row indices of non zeros entries in X
+  //   J  nnz vector of column indices of non zeros entries in X
+  //   V  nnz vector of non-zeros entries in X
+  //   Optional:
+  //     m  number of rows
+  //     n  number of cols
+  // Outputs:
+  //   X  m by n matrix of type T whose entries are to be found 
+  //
+  // Example:
+  //   Eigen::SparseMatrix<double> A;
+  //   std::vector<Eigen::Triplet<double> > IJV;
+  //   buildA(IJV);
+  //   if (A.rows() == 0)
+  //   {
+  //     A = Eigen::SparseMatrix<double>(rows,cols);
+  //     igl::sparse_cached_precompute(IJV,A,A_data);
+  //   }
+  //   else
+  //     igl::sparse_cached(IJV,s.A,s.A_data);
+
+  template <typename DerivedI, typename Scalar>
+  IGL_INLINE void sparse_cached_precompute(
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedI> & J,
+    Eigen::SparseMatrix<Scalar>& X,
+    Eigen::VectorXi& data);
+  
+  template <typename Scalar>
+  IGL_INLINE void sparse_cached_precompute(
+    const std::vector<Eigen::Triplet<Scalar> >& triplets,
+    Eigen::SparseMatrix<Scalar>& X,
+    Eigen::VectorXi& data);
+
+  template <typename Scalar>
+  IGL_INLINE void sparse_cached(
+    const std::vector<Eigen::Triplet<Scalar> >& triplets,
+    Eigen::SparseMatrix<Scalar>& X,
+    const Eigen::VectorXi& data);
+
+  template <typename DerivedV, typename Scalar>
+  IGL_INLINE void sparse_cached(
+    const Eigen::MatrixBase<DerivedV>& V,
+    Eigen::SparseMatrix<Scalar>& X,
+    const Eigen::VectorXi& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "sparse_cached.cpp"
+#endif
+
+#endif

tutorial/710_SLIM/main.cpp

@@ -152,7 +152,7 @@ void deform_3d_demo_iter(igl::viewer::Viewer& viewer) {
     double soft_const_p = 1e5;
     sData.exp_factor = 5.0;
     slim_precompute(V,F,V_0,sData,igl::SLIMData::EXP_CONFORMAL,b,bc,soft_const_p);
-    cout << "precomputed" << endl;
+    //cout << "precomputed" << endl;
 
     first_iter = false;
     display_3d_mesh(viewer);