added reserves to functions using DynamicSparseMatrix, where appropriate

Former-commit-id: a6a6b35550a56ee65f0afd2f08f6995b4640be0e

adjacency_matrix.h

@@ -46,6 +46,7 @@ inline void igl::adjacency_matrix(
 {
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
     dyn_A(F.maxCoeff()+1, F.maxCoeff()+1);
+  dyn_A.reserve(6*(F.maxCoeff()+1));
 
   // Loop over faces
   for(int i = 0;i<F.rows();i++)

cotmatrix.h

@@ -91,6 +91,10 @@ inline void igl::cotmatrix(
   assert(F.cols() == 3);
 
   Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
+  // This is important! it could decrease the comptuation time by a factor of 2
+  // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+  // row
+  dyn_L.reserve(7*V.rows());
   
   // Loop over triangles
   for (unsigned i = 0; i < F.rows(); i++)

examples/sparse/3x3.obj

@@ -0,0 +1,17 @@
+v 40.000000 40.000000 0.000000
+v 40.000000 200.000000 0.000000
+v 40.000000 360.000000 0.000000
+v 200.000000 40.000000 0.000000
+v 200.000000 200.000000 0.000000
+v 200.000000 360.000000 0.000000
+v 360.000000 40.000000 0.000000
+v 360.000000 200.000000 0.000000
+v 360.000000 360.000000 0.000000
+f 1 4 2
+f 4 5 2
+f 2 5 3
+f 5 6 3
+f 4 7 5
+f 7 8 5
+f 5 8 6
+f 8 9 6

examples/sparse/IJV.h

@@ -0,0 +1,40 @@
+// Templates: 
+template <typename I,typename V>
+struct IJV
+{
+  // Subscripts
+  I i,j;
+  // Value
+  V v;
+  IJV(int i, int j, double v)
+  {
+    this->i = i;
+    this->j = j;
+    this->v = v;
+  }
+  // Compare based on i then j
+  bool operator<(const IJV & B) const
+  {
+    if(this->i < B.i)
+    {
+      return true;
+    }else if(this->i > B.i)
+    {
+      return false;
+    }else
+    {
+      if(this->j < B.j)
+      {
+        return true;
+      }else if(this->j > B.j)
+      {
+        return false;
+      }else
+      {
+        // i and j are equal
+        return false;
+      }
+    }
+    return false;
+  }
+};

examples/sparse/Makefile

@@ -8,8 +8,8 @@ all: example
 igl_lib=-I../../
 eigen=-I/usr/local/include/eigen3 -I/usr/local/include/eigen3/unsupported
 
-CFLAGS=-g
-#CFLAGS=-Os -DNDEBUG
+#CFLAGS=-g
+CFLAGS=-Os -DNDEBUG
 inc=$(igl_lib) $(eigen)
 lib=
 

examples/sparse/example.cpp

@@ -16,6 +16,7 @@ using namespace std;
 using namespace igl;
 
 #include "SparseSolver.h"
+#include "IJV.h"
 
 // Build cotangent matrix by looping over triangles filling in eigen's dynamic
 // sparse matrix then casting to sparse matrix
@@ -31,6 +32,9 @@ inline void cotmatrix_dyncast(
   assert(F.cols() == 3);
 
   Eigen::DynamicSparseMatrix<double, Eigen::RowMajor> dyn_L (V.rows(), V.rows());
+
+  // THIS MAKES A HUGE DIFFERENCE 2x
+  dyn_L.reserve(7*V.rows());
   
   // Loop over triangles
   for (unsigned i = 0; i < F.rows(); i++)
@@ -60,6 +64,7 @@ inline void cotmatrix_dyncast(
     //dyn_L.coeffRef(vi3, vi3) -= cot2;
     //dyn_L.coeffRef(vi1, vi1) -= cot2;
   }
+
   for (int k=0; k < dyn_L.outerSize(); ++k)
   {
     double tmp = 0.0f;
@@ -123,6 +128,8 @@ inline void cotmatrix_sparsesolver(
   }
 }
 
+// First build adjacency list then use it to build a cotan matrix loop over
+// vertices in order
 inline void cotmatrix_adjacencylist(
   const Eigen::MatrixXd & V, 
   const Eigen::MatrixXi & F, 
@@ -133,6 +140,105 @@ inline void cotmatrix_adjacencylist(
   //adjacency_list(F,A);
 }
 
+// Build a cotmatrix by building a list of values to add in an IJV list, then
+// sort the IJV list and *add* them to the sparsematrix in sorted order
+inline void cotmatrix_sort_ijv(
+  const Eigen::MatrixXd & V, 
+  const Eigen::MatrixXi & F, 
+  Eigen::SparseMatrix<double>& L)
+{
+  // Assumes vertices are 3D or 2D
+  assert((V.cols() == 3) || (V.cols() == 2));
+  int dim = V.cols();
+  // Assumes faces are triangles
+  assert(F.cols() == 3);
+
+  vector<IJV<int,double> > Lijv;
+  // We know that we will have 12 entries per face
+  Lijv.reserve(6*F.rows());
+
+  // Loop over triangles
+  for (unsigned i = 0; i < F.rows(); i++)
+  {
+    // Corner indices of this triangle
+    int vi1 = F(i,0);
+    int vi2 = F(i,1);
+    int vi3 = F(i,2);
+    // Grab corner positions of this triangle
+    Eigen::Vector3d v1(V(vi1,0), V(vi1,1), V(vi1,2));
+    Eigen::Vector3d v2(V(vi2,0), V(vi2,1), V(vi2,2));
+    Eigen::Vector3d v3(V(vi3,0), V(vi3,1), V(vi3,2));
+    // Compute cotangent of angles at each corner
+    double cot1, cot2, cot3;
+    computeCotWeights(v1, v2, v3, cot1, cot2, cot3);
+    // Throw each corner's cotangent at opposite edge (in both directions)
+    Lijv.push_back(IJV<int,double>(vi1, vi2,cot3));
+    //Lijv.push_back(IJV<int,double>(vi2, vi1,cot3));
+    Lijv.push_back(IJV<int,double>(vi2, vi3,cot1));
+    //Lijv.push_back(IJV<int,double>(vi3, vi2,cot1));
+    Lijv.push_back(IJV<int,double>(vi3, vi1,cot2));
+    //Lijv.push_back(IJV<int,double>(vi1, vi3,cot2));
+    //// diagonal entries
+    Lijv.push_back(IJV<int,double>(vi1, vi1,-cot3));
+    //Lijv.push_back(IJV<int,double>(vi2, vi2,-cot3));
+    Lijv.push_back(IJV<int,double>(vi2, vi2,-cot1));
+    //Lijv.push_back(IJV<int,double>(vi3, vi3,-cot1));
+    Lijv.push_back(IJV<int,double>(vi3, vi3,-cot2));
+    //Lijv.push_back(IJV<int,double>(vi1, vi1,-cot2));
+  }
+  sort(Lijv.begin(),Lijv.end());
+
+  L = Eigen::SparseMatrix<double>(V.rows(),V.rows());
+  L.reserve(7*V.rows());
+
+  // Now that all the (i,j,v) triplets are sorted we can add them to the
+  // SparseMatrix directly
+  int i = -1;
+  int j = -1;
+  double v = 0;
+  vector<IJV<int,double> >::iterator last;
+  for(
+    vector<IJV<int,double> >::iterator Lit = Lijv.begin();
+    Lit != Lijv.end();
+    Lit++)
+  {
+    // Sum up non-unique entries
+    if(i == Lit->i && j == Lit->j)
+    {
+      last->v += Lit->v;
+      Lit->v = 0;
+    }else
+    {
+      i = Lit->i;
+      j = Lit->j;
+      last = Lit;
+    }
+  }
+
+  i = -1;
+  j = -1;
+  // Now our ijv list is sorted and we don't have dups
+  for(
+    vector<IJV<int,double> >::iterator Lit = Lijv.begin();
+    Lit != Lijv.end();
+    Lit++)
+  {
+    // only consider non-zeros
+    if(Lit->v != 0)
+    {
+      // If new outer then need to start
+      if(i != Lit->i)
+      {
+        i = Lit->i;
+        L.startVec(i);
+      }
+      L.insertBack(Lit->j,Lit->i) = Lit->v;
+    }
+  }
+  L.finalize();
+  L = L + SparseMatrix<double>(L.transpose());
+}
+
 int main(int argc, char * argv[])
 {
   if(argc < 2)
@@ -155,9 +261,21 @@ int main(int argc, char * argv[])
 
   double before;
   int trials;
-  int max_trials = 5;
+  int max_trials = 10;
   SparseMatrix<double> L;
 
+  printf("cotmatrix:\n  ");
+  before = get_seconds();
+  for(trials = 0;trials<max_trials;trials++)
+  {
+    cotmatrix(V,F,L);
+  }
+  printf("%g\n",(get_seconds()-before)/(double)trials);
+  if(L.rows() < 10)
+  {
+    print_ijv(L);
+  }
+
   printf("cotmatrix_dyncast:\n  ");
   before = get_seconds();
   for(trials = 0;trials<max_trials;trials++)
@@ -165,14 +283,30 @@ int main(int argc, char * argv[])
     cotmatrix_dyncast(V,F,L);
   }
   printf("%g\n",(get_seconds()-before)/(double)trials);
+  if(L.rows() < 10)
+  {
+    print_ijv(L);
+  }
 
-  printf("cotmatrix_sparsesolver:\n  ");
+  //printf("cotmatrix_sparsesolver:\n  ");
+  //before = get_seconds();
+  //for(trials = 0;trials<max_trials;trials++)
+  //{
+  //  cotmatrix_sparsesolver(V,F,L);
+  //}
+  //printf("%g\n",(get_seconds()-before)/(double)trials);
+
+  printf("cotmatrix_sort_ijv:\n  ");
   before = get_seconds();
   for(trials = 0;trials<max_trials;trials++)
   {
-    cotmatrix_sparsesolver(V,F,L);
+    cotmatrix_sort_ijv(V,F,L);
   }
   printf("%g\n",(get_seconds()-before)/(double)trials);
+  if(L.rows() < 10)
+  {
+    print_ijv(L);
+  }
 
   return 0;
 }

examples/sparse/test.cpp

@@ -0,0 +1,37 @@
+#!/bin/bash
+//usr/bin/tail -n +2 $0 | g++ -o main -x c++ - && ./main && rm main && exit
+#include "IJV.h"
+
+#include <cstdio>
+#include <vector>
+#include <algorithm>
+using namespace std;
+
+int print_ijv(const IJV<int,double> & ijv)
+{
+  printf("%d %d %g\n",
+    ijv.i,
+    ijv.j,
+    ijv.v);
+}
+
+int main(int argc, char * argv[])
+{
+  vector<IJV<int,double> > Aijv;
+  Aijv.push_back(IJV<int,double>(1,2,10));
+  Aijv.push_back(IJV<int,double>(4,2,10));
+  Aijv.push_back(IJV<int,double>(4,3,10));
+  Aijv.push_back(IJV<int,double>(9,2,10));
+  Aijv.push_back(IJV<int,double>(1,2,10));
+  Aijv.push_back(IJV<int,double>(1,3,10));
+  Aijv.push_back(IJV<int,double>(1,1,10));
+  Aijv.push_back(IJV<int,double>(3,2,10));
+
+  printf("Original:\n");
+  for_each(Aijv.begin(),Aijv.end(),print_ijv);
+  sort(Aijv.begin(),Aijv.end());
+  printf("Sorted:\n");
+  for_each(Aijv.begin(),Aijv.end(),print_ijv);
+
+  return 0;
+}

readme.txt

@@ -59,7 +59,9 @@ description of each template, input and output in each prototype.
   = Specific style conventions =
 
   Each file should be wrapped in an ifndef compiler directive. If the
-  file's (and function's) name is example.h, then the ifndef should be
+  file's (and function's) name is example.h, then the ifndef should
+  always begin with IGL_, then the function/file name in all caps then
+  _H. As in:
 #ifndef IGL_EXAMPLE_H
 #define IGL_EXAMPLE_H
 ...
@@ -80,3 +82,4 @@ description of each template, input and output in each prototype.
   Be generous by templating your inputs and outputs. If you do use
   templates, you must document the template just as you document inputs
   and outputs.
+

repdiag.h

@@ -47,8 +47,13 @@ inline void igl::repdiag(
   int m = A.rows();
   int n = A.cols();
 
+  // Should be able to *easily* do this in coherent order without
+  // dynamicsparsematrix
+
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
     dyn_B(m*d,n*d);
+  // Reserve enough space for new non zeros
+  dyn_B.reserve(d*A.nonZeros());
 
   // loop over reps
   for(int i=0;i<d;i++)

slice.h

@@ -78,6 +78,9 @@ inline void igl::slice(
   // Resize output
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
     dyn_Y(ym,yn);
+  // Take a guess at the number of nonzeros (this assumes uniform distribution
+  // not banded or heavily diagonal)
+  dyn_Y.reserve((X.nonzeros()/(X.rows()*X.cols())) * (ym*yn));
   // Iterate over outside
   for(int k=0; k<X.outerSize(); ++k)
   {

sparse.h

@@ -79,8 +79,9 @@ inline void igl::sparse(
   // number of values
   int nv = V.size();
 
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
-    dyn_X(m,n);
+  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(m,n);
+  // over estimate the number of entries
+  dyn_X.reserve(I.size());
   for(int i = 0;i < nv;i++)
   {
     dyn_X.coeffRef((int)I(i),(int)J(i)) += (T)V(i);