working on speeding up the slim iterations. First version of cached AtA and sparse

Former-commit-id: 94fce5c08be75027508fa3aa7556d440a79f8525

include/igl/slim.cpp

@@ -34,6 +34,14 @@
 #include <Eigen/SparseCholesky>
 #include <Eigen/IterativeLinearSolvers>
 
+#include <igl/Timer.h>
+#include <igl/sparse_fast.h>
+#include <igl/sparse_AtA_fast.h>
+
+#ifdef CHOLMOD
+#include <Eigen/CholmodSupport>
+#endif
+
 namespace igl
 {
   namespace slim
@@ -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,23 +498,90 @@ namespace igl
     IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
     {
       // formula (35) in paper
+
+      igl::Timer t;
+
+      t.start();
       Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
       buildA(s,A);
+      t.stop();
+      std::cerr << "buildA: " << t.getElapsedTime() << std::endl;
 
+      // // // TEST CODE
+      // std::vector<Eigen::Triplet<double> > test;
+      // test.push_back(Eigen::Triplet<double>(4,3,10));
+      // test.push_back(Eigen::Triplet<double>(1,0,20));
+      // test.push_back(Eigen::Triplet<double>(2,2,30));
+      // test.push_back(Eigen::Triplet<double>(4,4,2.5));
+
+      // A = Eigen::SparseMatrix<double>(5,5);
+      // A.setFromTriplets(test.begin(),test.end());
+      A.makeCompressed();
+
+      //std::cin.get();
+
+      
+
+      t.start();
       Eigen::SparseMatrix<double> At = A.transpose();
       At.makeCompressed();
+      t.stop();
+      std::cerr << "At: " << t.getElapsedTime() << std::endl;
 
+      t.start();
       Eigen::SparseMatrix<double> id_m(At.rows(), At.rows());
       id_m.setIdentity();
+      t.stop();
+      std::cerr << "idm: " << t.getElapsedTime() << std::endl;
+
+      // std::cerr << "\n TEST CODE" << std::endl;
+      //   s.WGL_M.resize(5);
+      //   s.WGL_M[0] = 1;
+      //   s.WGL_M[1] = 2;
+      //   s.WGL_M[2] = 3;
+      //   s.WGL_M[3] = 3;
+      //   s.WGL_M[4] = 10;
+      t.start();
+      Eigen::SparseMatrix<double> AtA_slow = At * s.WGL_M.asDiagonal() * A;
+      t.stop();
+      std::cerr << "AtA slow: " << t.getElapsedTime() << std::endl;
+      // std::cerr << "AtA slow: " << AtA_slow << std::endl;
+
+      t.start();
+      Eigen::SparseMatrix<double> AtA_fast;
+      igl::sparse_AtA_fast_data data;
+      data.W = s.WGL_M;
+      igl::sparse_AtA_fast_precompute(A,AtA_fast,data);
+      t.stop();
+      std::cerr << "AtA fast pre: " << t.getElapsedTime() << std::endl;
+      t.start();
+      igl::sparse_AtA_fast(A,AtA_fast,data);
+      t.stop();
+      std::cerr << "AtA fast: " << t.getElapsedTime() << std::endl;
+      // std::cerr << "AtA fast: " << AtA_fast  << std::endl;
+
+      std::cerr << "Difference: " << (AtA_fast - AtA_slow).norm() << std::endl;
+
+        exit(0);
 
       // add proximal penalty
+      t.start();
       L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; //add also a proximal term
       L.makeCompressed();
+      t.stop();
+      std::cerr << "proximal: " << t.getElapsedTime() << std::endl;
 
+      t.start();
       buildRhs(s, At);
+      t.stop();
+      std::cerr << "rhs: " << t.getElapsedTime() << std::endl;
+
+      t.start();
       Eigen::SparseMatrix<double> OldL = L;
       add_soft_constraints(s,L);
       L.makeCompressed();
+      t.stop();
+      std::cerr << "soft: " << t.getElapsedTime() << std::endl;
     }
 
     IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
@@ -780,7 +867,26 @@ namespace igl
           }
         }
       }
+
+      // igl::Timer t;
+      // t.start();
       A.setFromTriplets(IJV.begin(), IJV.end());
+      // t.stop();
+      // std::cerr << "setFromTriplets: " << t.getElapsedTime() << std::endl;
+
+      // // Debug Code
+      // Eigen::VectorXi data;
+      // t.start();
+      // igl::fast_sparse_precompute(IJV,A,data);
+      // t.stop();
+      // std::cerr << "fast_precompute: " << t.getElapsedTime() << std::endl;
+
+      // t.start();
+      // igl::fast_sparse(IJV,A,data);
+      // t.stop();
+      // std::cerr << "fast_precompute: " << t.getElapsedTime() << std::endl;
+
+
     }
 
     IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At)

include/igl/sparse_AtA_fast.cpp

@@ -0,0 +1,115 @@
+// 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_AtA_fast.h"
+
+#include <iostream>
+#include <vector>
+#include <unordered_map>
+#include <map>
+#include <utility>
+
+IGL_INLINE void igl::sparse_AtA_fast_precompute(
+    const Eigen::SparseMatrix<double>& A,
+    Eigen::SparseMatrix<double>& AtA,
+    igl::sparse_AtA_fast_data& data)
+{
+  // 1 Compute At
+  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.at(col).push_back(row);
+      Col_IndexPtr.at(col).push_back(value_index);
+    }
+  }
+
+  Eigen::SparseMatrix<double> 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());
+
+  // 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 corresponding indices
+      for (unsigned i=0;i<Col_RowPtr.at(row).size();++i)
+      {
+        for (unsigned j=0;j<Col_RowPtr.at(col).size();++j)
+        {
+          if (Col_RowPtr.at(row)[i] == Col_RowPtr.at(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(data.W[Col_RowPtr[col][j]]);
+          }
+        }
+      }
+    }
+  }
+  data.I_outer.push_back(data.I_row.size()); // makes it more efficient to iterate later on
+}
+
+IGL_INLINE void igl::sparse_AtA_fast(
+    const Eigen::SparseMatrix<double>& A,
+    Eigen::SparseMatrix<double>& AtA,
+    const igl::sparse_AtA_fast_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.I_w[j] * *(A.valuePtr() + data.I_col[j]);
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/sparse_AtA_fast.h

@@ -0,0 +1,47 @@
+// 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_ATA_FAST_H
+#define IGL_SPARSE_ATA_FAST_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 sparse_AtA_fast_data
+  {
+    // Weights
+    Eigen::VectorXd W;
+
+    // Flatten composition rules
+    std::vector<int> I_row;
+    std::vector<int> I_col;
+    std::vector<double> I_w;
+
+    // For each entry of AtA, points to the beginning
+    // of the composition rules
+    std::vector<int> I_outer;
+  };
+
+  IGL_INLINE void sparse_AtA_fast_precompute(
+    const Eigen::SparseMatrix<double>& A,
+    Eigen::SparseMatrix<double>& AtA,
+    sparse_AtA_fast_data& data);
+
+  IGL_INLINE void sparse_AtA_fast(
+    const Eigen::SparseMatrix<double>& A,
+    Eigen::SparseMatrix<double>& AtA,
+    const sparse_AtA_fast_data& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "sparse_AtA_fast.cpp"
+#endif
+
+#endif

include/igl/sparse_fast.cpp

@@ -0,0 +1,115 @@
+// 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_fast.h"
+
+#include <iostream>
+#include <vector>
+#include <unordered_map>
+#include <map>
+#include <utility>
+
+// Hashing function
+// namespace std {
+//     template<> 
+//     struct hash<pair<int, int> > {
+//         size_t operator()(const pair<int, int>& p) const {
+//             size_t seed = 0;
+//             hash<int> h;
+//             seed ^= h(p.first) + 0x9e3779b9 + (seed << 6) + (seed >> 2); 
+//             seed ^= h(p.second) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
+//             return seed;
+//         }
+//     };
+// }
+
+
+IGL_INLINE void igl::sparse_fast_precompute(
+  const Eigen::VectorXi & I,
+  const Eigen::VectorXi & J,
+  Eigen::SparseMatrix<double>& 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_fast_precompute(t,X,data);
+}
+
+IGL_INLINE void igl::sparse_fast_precompute(
+  const std::vector<Eigen::Triplet<double> >& triplets,
+  Eigen::SparseMatrix<double>& X,
+  Eigen::VectorXi& data)
+{
+    // Construct an empty sparse matrix
+    X.setFromTriplets(triplets.begin(),triplets.end());
+    X.makeCompressed();
+
+    // Build hash table for all nnz entries
+    // TODO: this is slow and could be done in nlogn
+    std::map<std::pair<int,int>,int> id;
+
+    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;
+
+        std::pair<int,int> rc = std::make_pair(row,col);
+        id[rc] = value_index;
+      }
+    }
+
+    // Compute the indices
+    data.resize(triplets.size());
+    for (unsigned i=0; i<triplets.size(); ++i)
+      data[i] = id[std::make_pair(triplets[i].row(),triplets[i].col())];
+}
+  
+IGL_INLINE void igl::sparse_fast(
+  const std::vector<Eigen::Triplet<double> >& triplets,
+  Eigen::SparseMatrix<double>& 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();
+}
+
+IGL_INLINE void igl::sparse_fast(
+  const Eigen::VectorXd & V,
+  Eigen::SparseMatrix<double>& 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
+#endif

include/igl/sparse_fast.h

@@ -0,0 +1,64 @@
+// 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_FAST_H
+#define IGL_SPARSE_FAST_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), functions
+  // like the sparse function in matlab. Divides the construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values.
+  //
+  // 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 
+  //
+  IGL_INLINE void sparse_fast_precompute(
+    const Eigen::VectorXi & I,
+    const Eigen::VectorXi & J,
+    Eigen::SparseMatrix<double>& X,
+    Eigen::VectorXi& data);
+  
+  IGL_INLINE void sparse_fast_precompute(
+    const std::vector<Eigen::Triplet<double> >& triplets,
+    Eigen::SparseMatrix<double>& X,
+    Eigen::VectorXi& data);
+
+  IGL_INLINE void sparse_fast(
+    const std::vector<Eigen::Triplet<double> >& triplets,
+    Eigen::SparseMatrix<double>& X,
+    const Eigen::VectorXi& data);
+
+  IGL_INLINE void sparse_fast(
+    const Eigen::VectorXd & V,
+    Eigen::SparseMatrix<double>& X,
+    const Eigen::VectorXi& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "sparse_fast.cpp"
+#endif
+
+#endif