fixed min_quad_* to handle non positive definite systems (though falls back on dense decomposition)

FIXED INLINES ON ~70 FILES. Please remember to put inline before each function


Former-commit-id: f304aef7877028f12dedb792703421c1ad283bb7

adjacency_matrix.h

@@ -31,7 +31,7 @@ namespace igl
   //
   // See also: edges, cotmatrix, diag
   template <typename T>
-  void adjacency_matrix(
+  inline void adjacency_matrix(
     const Eigen::MatrixXi & F, 
     Eigen::SparseMatrix<T>& A);
 }
@@ -40,7 +40,7 @@ namespace igl
 #include "verbose.h"
 
 template <typename T>
-void igl::adjacency_matrix(
+inline void igl::adjacency_matrix(
   const Eigen::MatrixXi & F, 
   Eigen::SparseMatrix<T>& A)
 {

basename.h

@@ -11,13 +11,13 @@ namespace igl
   // Returns string containing basename (see php's basename)
   //
   // See also: dirname, pathinfo
-  std::string basename(const std::string & path);
+  inline std::string basename(const std::string & path);
 }
 
 // Implementation
 #include <algorithm>
 
-std::string igl::basename(const std::string & path)
+inline std::string igl::basename(const std::string & path)
 {
   if(path == "")
   {

cocoa_key_to_anttweakbar_key.h

@@ -8,11 +8,11 @@ namespace igl
   // Inputs:
   //   key  unsigned char key from keyboard
   // Returns int of new key code 
-  int cocoa_key_to_anttweakbar_key(int key);
+  inline int cocoa_key_to_anttweakbar_key(int key);
 }
 
 // Implementation
-int igl::cocoa_key_to_anttweakbar_key(int key)
+inline int igl::cocoa_key_to_anttweakbar_key(int key)
 {
   // I've left commented the AntTweakBar key codes that correspond to keys I
   // don't have on my keyboard. Please fill this in if you have those keys

cotmatrix.h

@@ -23,7 +23,7 @@ namespace igl
   //   L  #V by #V cotangent matrix, each row i corresponding to V(i,:)
   //
   // See also: adjacency_matrix
-  void cotmatrix(
+  inline void cotmatrix(
     const Eigen::MatrixXd & V, 
     const Eigen::MatrixXi & F,
     Eigen::SparseMatrix<double>& L);
@@ -38,7 +38,7 @@ namespace igl
   //   cot2 cotangent of angle at corner #2
   //   cot3  cotangent of angle at corner #3
   //   
-  void computeCotWeights(
+  inline void computeCotWeights(
     const Eigen::Vector3d& v1, 
     const Eigen::Vector3d& v2, 
     const Eigen::Vector3d& v3, 
@@ -52,7 +52,7 @@ namespace igl
 // For error printing
 #include <cstdio>
 
-void igl::computeCotWeights(
+inline void igl::computeCotWeights(
   const Eigen::Vector3d& v1, 
   const Eigen::Vector3d& v2, 
   const Eigen::Vector3d& v3, 
@@ -79,13 +79,14 @@ void igl::computeCotWeights(
   cot3 = (-v23.dot(-v13)) / halfArea /2;
 }
 
-void igl::cotmatrix(
+inline void igl::cotmatrix(
   const Eigen::MatrixXd & V, 
   const Eigen::MatrixXi & F, 
   Eigen::SparseMatrix<double>& L)
 {
-  // Assumes vertices are 3D
-  assert(V.cols() == 3);
+  // 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);
 
@@ -99,9 +100,9 @@ void igl::cotmatrix(
     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));
+    Eigen::Vector3d v1(V(vi1,0), V(vi1,1), (dim==2?0:V(vi1,2)));
+    Eigen::Vector3d v2(V(vi2,0), V(vi2,1), (dim==2?0:V(vi2,2)));
+    Eigen::Vector3d v3(V(vi3,0), V(vi3,1), (dim==2?0:V(vi3,2)));
     // Compute cotangent of angles at each corner
     double cot1, cot2, cot3;
     computeCotWeights(v1, v2, v3, cot1, cot2, cot3);

create_index_vbo.h

@@ -20,7 +20,7 @@ namespace igl
   // Outputs:
   //   F_vbo_id  buffer id for face indices
   //
-  void create_index_vbo(
+  inline void create_index_vbo(
     const Eigen::MatrixXi & F,
     GLuint & F_vbo_id);
 }
@@ -28,7 +28,7 @@ namespace igl
 // Implementation
 
 // http://www.songho.ca/opengl/gl_vbo.html#create
-void igl::create_index_vbo(
+inline void igl::create_index_vbo(
   const Eigen::MatrixXi & F,
   GLuint & F_vbo_id)
 {

create_mesh_vbo.h

@@ -26,7 +26,7 @@ namespace igl
   // NOTE: when using glDrawElements VBOs for V and F using MatrixXd and
   // MatrixXi will have types GL_DOUBLE and GL_UNSIGNED_INT respectively
   //
-  void create_mesh_vbo(
+  inline void create_mesh_vbo(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     GLuint & V_vbo_id,
@@ -40,7 +40,7 @@ namespace igl
   //   V_vbo_id  buffer id for vertex positions
   //   F_vbo_id  buffer id for face indices
   //   N_vbo_id  buffer id for vertex positions
-  void create_mesh_vbo(
+  inline void create_mesh_vbo(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const Eigen::MatrixXd & N,
@@ -55,7 +55,7 @@ namespace igl
 #include "create_index_vbo.h"
 
 // http://www.songho.ca/opengl/gl_vbo.html#create
-void igl::create_mesh_vbo(
+inline void igl::create_mesh_vbo(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   GLuint & V_vbo_id,
@@ -68,7 +68,7 @@ void igl::create_mesh_vbo(
 }
 
 // http://www.songho.ca/opengl/gl_vbo.html#create
-void igl::create_mesh_vbo(
+inline void igl::create_mesh_vbo(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   const Eigen::MatrixXd & N,

create_vector_vbo.h

@@ -23,7 +23,7 @@ namespace igl
   //   V_vbo_id  buffer id for vectors
   //
   template <typename T>
-  void create_vector_vbo(
+  inline void create_vector_vbo(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & V,
     GLuint & V_vbo_id);
 }
@@ -33,7 +33,7 @@ namespace igl
 
 // http://www.songho.ca/opengl/gl_vbo.html#create
 template <typename T>
-void igl::create_vector_vbo(
+inline void igl::create_vector_vbo(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & V,
   GLuint & V_vbo_id)
 {

diag.h

@@ -15,7 +15,7 @@ namespace igl
   // Outputs:
   //   V  a min(m,n) sparse vector
   template <typename T>
-  void diag(
+  inline void diag(
     const Eigen::SparseMatrix<T>& X, 
     Eigen::SparseVector<T>& V);
   // Templates:
@@ -25,7 +25,7 @@ namespace igl
   // Outputs:
   //  X  a m by m sparse matrix
   template <typename T>
-  void diag(
+  inline void diag(
     const Eigen::SparseVector<T>& V,
     Eigen::SparseMatrix<T>& X);
 }
@@ -34,7 +34,7 @@ namespace igl
 #include "verbose.h"
 
 template <typename T>
-void igl::diag(
+inline void igl::diag(
   const Eigen::SparseMatrix<T>& X, 
   Eigen::SparseVector<T>& V)
 {
@@ -58,7 +58,7 @@ void igl::diag(
 }
 
 template <typename T>
-void igl::diag(
+inline void igl::diag(
   const Eigen::SparseVector<T>& V,
   Eigen::SparseMatrix<T>& X)
 {
@@ -66,7 +66,7 @@ void igl::diag(
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
 
   // loop over non-zeros
-  for(typename SparseVector<T>::InnerIterator it(V); it; ++it)
+  for(typename Eigen::SparseVector<T>::InnerIterator it(V); it; ++it)
   {
     dyn_X.coeffRef(it.index(),it.index()) += it.value();
   }

edgetopology.h

@@ -21,7 +21,7 @@ namespace igl
     // FE : #Fx3, Stores the Triangle-Edge relation
     // EF : #Ex2: Stores the Edge-Triangle relation (unsorted)
 
-    void edgetopology(Eigen::MatrixXd& V, Eigen::MatrixXi& F, 
+    inline void edgetopology(Eigen::MatrixXd& V, Eigen::MatrixXi& F, 
                       Eigen::MatrixXi& EV, Eigen::MatrixXi& FE, Eigen::MatrixXi& EF)
     {
         assert(igl::isManifold(V,F));
@@ -81,4 +81,4 @@ namespace igl
     }
 }
 
-#endif
+#endif

find.h

@@ -84,7 +84,7 @@ inline void igl::find(
 
   int i = 0;
   // loop over non-zeros
-  for(typename SparseVector<T>::InnerIterator it(X); it; ++it)
+  for(typename Eigen::SparseVector<T>::InnerIterator it(X); it; ++it)
   {
     I(i) = it.index();
     V(i) = it.value();

is_symmetric.h

@@ -19,8 +19,8 @@ inline bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
   {
     return false;
   }
-  SparseMatrix<T> AT = A.transpose();
-  SparseMatrix<T> AmAT = A-AT;
+  Eigen::SparseMatrix<T> AT = A.transpose();
+  Eigen::SparseMatrix<T> AmAT = A-AT;
   //// Eigen screws up something with LLT if you try to do
   //SparseMatrix<T> AmAT = A-A.transpose();
   //// Eigen crashes at runtime if you try to do

ismanifold.h

@@ -14,7 +14,7 @@
 namespace igl 
 {
     // check if the mesh is edge-manifold
-    bool isManifold(Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+    inline bool isManifold(Eigen::MatrixXd& V, Eigen::MatrixXi& F)
     {
         std::vector<std::vector<int> > TTT;
         for(int f=0;f<F.rows();++f)
@@ -47,4 +47,4 @@ namespace igl
     }
 }
 
-#endif
+#endif

load_shader.h

@@ -17,14 +17,14 @@ namespace igl
   //     GL_FRAGMENT_SHADER
   //     GL_GEOMETRY_SHADER
   // Returns  index id of the newly created shader, 0 on error
-  GLuint load_shader(const char *src,const GLenum type);
+  inline GLuint load_shader(const char *src,const GLenum type);
 }
 
 // Implmentation
 // Copyright Denis Kovacs 4/10/08
 #include "print_shader_info_log.h"
 #include <cstdio>
-GLuint igl::load_shader(const char *src,const GLenum type)
+inline GLuint igl::load_shader(const char *src,const GLenum type)
 {
   GLuint s = glCreateShader(type);
   if(s == 0)

lu_lagrange.h

@@ -29,12 +29,24 @@ namespace igl
   //   U  upper triangular matrix such that Q = L*U
   // Returns true on success, false on error
   //
+  // Note: C should *not* have any empty columns. Typically C is the slice of
+  // the linear constraints matrix Aeq concerning the unknown variables of a
+  // quadratic optimization. Generally constraints may deal with unknowns as
+  // well as knowns. Each linear constraint corresponds to a column of Aeq. As
+  // long as each constraint concerns at least one unknown then the
+  // corresponding column in C will have at least one non zero entry. If a
+  // constraint concerns *no* unknowns, you should double check that this is a
+  // valid constraint. How can you constrain known values to each other? This
+  // is either a contradiction to the knowns' values or redundent. In either
+  // case, it's not this functions responsiblilty to handle empty constraints
+  // so you will get an error.
+  //
   template <typename T>
-  bool lu_lagrange(
-    const SparseMatrix<T> & ATA,
-    const SparseMatrix<T> & C,
-    SparseMatrix<T> & L,
-    SparseMatrix<T> & U);
+  inline bool lu_lagrange(
+    const Eigen::SparseMatrix<T> & ATA,
+    const Eigen::SparseMatrix<T> & C,
+    Eigen::SparseMatrix<T> & L,
+    Eigen::SparseMatrix<T> & U);
 
 }
 
@@ -46,11 +58,11 @@ namespace igl
 #include "sparse.h"
 
 template <typename T>
-bool igl::lu_lagrange(
-  const SparseMatrix<T> & ATA,
-  const SparseMatrix<T> & C,
-  SparseMatrix<T> & L,
-  SparseMatrix<T> & U)
+inline bool igl::lu_lagrange(
+  const Eigen::SparseMatrix<T> & ATA,
+  const Eigen::SparseMatrix<T> & C,
+  Eigen::SparseMatrix<T> & L,
+  Eigen::SparseMatrix<T> & U)
 {
   // number of unknowns
   int n = ATA.rows();
@@ -61,14 +73,42 @@ bool igl::lu_lagrange(
   if(m != 0)
   {
     assert(C.rows() == n);
+    if(C.nonZeros() == 0)
+    {
+      // See note above about empty columns in C
+      fprintf(stderr,"Error: lu_lagrange() C has columns but no entries\n");
+      return false;
+    }
   }
 
+  // Check that each column of C has at least one entry
+  std::vector<bool> has_entry; has_entry.resize(C.cols(),false);
+  // Iterate over outside
+  for(int k=0; k<C.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<T>::InnerIterator it (C,k); it; ++it)
+    {
+      has_entry[it.col()] = true;
+    }
+  }
+  for(int i=0;i<(int)has_entry.size();i++)
+  {
+    if(!has_entry[i])
+    {
+      // See note above about empty columns in C
+      fprintf(stderr,"Error: lu_lagrange() C(:,%d) has no entries\n",i);
+      return false;
+    }
+  }
+
+
 
   // Cholesky factorization of ATA
   //// Eigen fails if you give a full view of the matrix like this:
   //Eigen::SparseLLT<SparseMatrix<T> > ATA_LLT(ATA);
-  SparseMatrix<T> ATA_LT = ATA.template triangularView<Eigen::Lower>();
-  Eigen::SparseLLT<SparseMatrix<T> > ATA_LLT(ATA_LT);
+  Eigen::SparseMatrix<T> ATA_LT = ATA.template triangularView<Eigen::Lower>();
+  Eigen::SparseLLT<Eigen::SparseMatrix<T> > ATA_LLT(ATA_LT);
 
   Eigen::SparseMatrix<T> J = ATA_LLT.matrixL();
 
@@ -81,6 +121,8 @@ bool igl::lu_lagrange(
 
   if(m == 0)
   {
+    // If there are no constraints (C is empty) then LU decomposition is just L
+    // and L' from cholesky decomposition
     L = J;
     U = J.transpose();
   }else
@@ -89,16 +131,13 @@ bool igl::lu_lagrange(
     Eigen::SparseMatrix<T> M = C;
     J.template triangularView<Eigen::Lower>().solveInPlace(M);
 
-
     // Compute cholesky factorizaiton of M'*M
     Eigen::SparseMatrix<T> MTM = M.transpose() * M;
 
-
-    Eigen::SparseLLT<SparseMatrix<T> > MTM_LLT(MTM.template triangularView<Eigen::Lower>());
+    Eigen::SparseLLT<Eigen::SparseMatrix<T> > MTM_LLT(MTM.template triangularView<Eigen::Lower>());
 
     Eigen::SparseMatrix<T> K = MTM_LLT.matrixL();
 
-
     //if(!MTM_LLT.succeeded())
     if(!((K*0).eval().nonZeros() == 0))
     {
@@ -139,10 +178,6 @@ bool igl::lu_lagrange(
     LJ << JJ,                        MI, (KJ.array() + n).matrix(); 
     LV << JV,                        MV,                        KV;
     igl::sparse(LI,LJ,LV,L);
-
-    //// Only keep lower and upper parts
-    //L = L.template triangularView<Lower>();
-    //U = U.template triangularView<Upper>();
   }
 
   return true;

mat_max.h

@@ -22,7 +22,7 @@ namespace igl
   //   I  vector the same size as Y containing the indices along dim of maximum
   //     entries
   template <typename T>
-  void mat_max(
+  inline void mat_max(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
     const int dim,
     Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
@@ -30,7 +30,7 @@ namespace igl
 }
 
 template <typename T>
-void igl::mat_max(
+inline void igl::mat_max(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
   const int dim,
   Eigen::Matrix<T,Eigen::Dynamic,1> & Y,

mat_min.h

@@ -24,7 +24,7 @@ namespace igl
   //
   // See also: mat_max
   template <typename T>
-  void mat_min(
+  inline void mat_min(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
     const int dim,
     Eigen::Matrix<T,Eigen::Dynamic,1> & Y,
@@ -34,7 +34,7 @@ namespace igl
 #include "verbose.h"
 
 template <typename T>
-void igl::mat_min(
+inline void igl::mat_min(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & X,
   const int dim,
   Eigen::Matrix<T,Eigen::Dynamic,1> & Y,

matlab-to-eigen.html

@@ -195,6 +195,15 @@ tr.gotcha2 td
         </td>
       </tr>
 
+      <tr class=d0>
+        <td><pre><code>I = low:step:hi</code></pre></td>
+        <td><pre><code>igl::colon(low,step,hi,I)</code></pre></td>
+        <td>IGL version should be templated enough to handle same situations as
+        matlab's colon. The matlab keyword <b>end</b> does not correspond in
+        the C++ version. You'll have to use M.size(),M.rows() or whatever.
+        </td>
+      </tr>
+
       <!-- Insert rows for each command pair -->
 
       <!-- Leave this here for copy and pasting

min_quad_with_fixed.h

@@ -5,6 +5,7 @@
 #include <Eigen/Core>
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
+#include <Eigen/SparseExtra>
 
 namespace igl
 {
@@ -30,7 +31,7 @@ namespace igl
   //     using min_quad_with_fixed_solve
   // Returns true on success, false on error
   template <typename T>
-  bool min_quad_with_fixed_precompute(
+  inline bool min_quad_with_fixed_precompute(
     const Eigen::SparseMatrix<T>& A,
     const Eigen::MatrixXi & known,
     const Eigen::SparseMatrix<T>& Aeq,
@@ -45,7 +46,7 @@ namespace igl
   //   Z  n by cols solution
   // Returns true on success, false on error
   template <typename T>
-  bool min_quad_with_fixed_solve(
+  inline bool min_quad_with_fixed_solve(
     const min_quad_with_fixed_data<T> & data,
     const Eigen::Matrix<T,Eigen::Dynamic,1> & B,
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Y,
@@ -57,30 +58,46 @@ namespace igl
 #include <Eigen/SparseExtra>
 #include <cassert>
 #include <cstdio>
+
 #include "slice.h"
 #include "is_symmetric.h"
-
 #include "find.h"
 #include "sparse.h"
+#include "repmat.h"
 #include "lu_lagrange.h"
+#include "full.h"
 
 template <typename T>
 struct igl::min_quad_with_fixed_data
 {
+  // Size of original system: number of unknowns + number of knowns
   int n;
+  // Whether A(unknown,unknown) is positive definite
   bool Auu_pd;
+  // Whether A(unknown,unknown) is symmetric
   bool Auu_sym;
+  // Indices of known variables
   Eigen::Matrix<int,Eigen::Dynamic,1> known;
+  // Indices of unknown variables
   Eigen::Matrix<int,Eigen::Dynamic,1> unknown;
+  // Indices of lagrange variables
   Eigen::Matrix<int,Eigen::Dynamic,1> lagrange;
+  // Indices of unknown variable followed by Indices of lagrange variables
   Eigen::Matrix<int,Eigen::Dynamic,1> unknown_lagrange;
-  SparseMatrix<T> preY;
-  SparseMatrix<T> L;
-  SparseMatrix<T> U;
+  // Matrix multiplied against Y when constructing right hand side
+  Eigen::SparseMatrix<T> preY;
+  // Tells whether system is sparse
+  bool sparse;
+  // Lower triangle of LU decomposition of final system matrix
+  Eigen::SparseMatrix<T> L;
+  // Upper triangle of LU decomposition of final system matrix
+  Eigen::SparseMatrix<T> U;
+  // Dense LU factorization
+  Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > lu;
 };
 
 template <typename T>
-bool igl::min_quad_with_fixed_precompute(
+inline bool igl::min_quad_with_fixed_precompute(
   const Eigen::SparseMatrix<T>& A,
   const Eigen::Matrix<int,Eigen::Dynamic,1> & known,
   const Eigen::SparseMatrix<T>& Aeq,
@@ -94,10 +111,9 @@ bool igl::min_quad_with_fixed_precompute(
   data.n = n;
 
   int neq = Aeq.rows();
-  // defulat is to have 0 linear equality constraints
+  // default is to have 0 linear equality constraints
   if(Aeq.size() != 0)
   {
-    //Aeq = Eigen::SparseMatrix<T>(0,n);
     assert(n == Aeq.cols());
   }
 
@@ -150,7 +166,7 @@ bool igl::min_quad_with_fixed_precompute(
   data.Auu_pd = pd;
 
   // Append lagrange multiplier quadratic terms
-  SparseMatrix<T> new_A;
+  Eigen::SparseMatrix<T> new_A;
   Eigen::Matrix<int,Eigen::Dynamic,1> AI;
   Eigen::Matrix<int,Eigen::Dynamic,1> AJ;
   Eigen::Matrix<T,Eigen::Dynamic,1> AV;
@@ -168,33 +184,29 @@ bool igl::min_quad_with_fixed_precompute(
   new_AI << AI, (AeqI.array()+n).matrix(), AeqJ;
   new_AJ << AJ, AeqJ, (AeqI.array()+n).matrix();
   new_AV << AV, AeqV, AeqV;
-  //new_AI.block(0,0,n,1) = AI;
-  //new_AJ.block(0,0,n,1) = AJ;
-  //new_AV.block(0,0,n,1) = AV;
-  //new_AI.block(n,0,neq,1) = AeqI+n;
-  //new_AJ.block(n,0,neq,1) = AeqJ;
-  //new_AV.block(n,0,neq,1) = AeqV;
-  //new_AI.block(n+neq,0,neq,1) = AeqJ;
-  //new_AJ.block(n+neq,0,neq,1) = AeqI+n;
-  //new_AV.block(n+neq,0,neq,1) = AeqV;
   igl::sparse(new_AI,new_AJ,new_AV,n+neq,n+neq,new_A);
 
   // precompute RHS builders
-  Eigen::SparseMatrix<T> Aulk;
-  igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
-  Eigen::SparseMatrix<T> Akul;
-  igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
+  if(kr > 0)
+  {
+    Eigen::SparseMatrix<T> Aulk;
+    igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
+    Eigen::SparseMatrix<T> Akul;
+    igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
 
-  //// This doesn't work!!!
-  //data.preY = Aulk + Akul.transpose();
-  Eigen::SparseMatrix<T> AkulT = Akul.transpose();
-  //// Resize preY before assigning
-  //data.preY.resize(data.unknown_lagrange.size(),data.known.size());
-  data.preY = Aulk + AkulT;
+    //// This doesn't work!!!
+    //data.preY = Aulk + Akul.transpose();
+    Eigen::SparseMatrix<T> AkulT = Akul.transpose();
+    data.preY = Aulk + AkulT;
+  }else
+  {
+    data.preY.resize(data.unknown_lagrange.size(),0);
+  }
 
   // Create factorization
   if(data.Auu_sym && data.Auu_pd)
   {
+    data.sparse = true;
     Eigen::SparseMatrix<T> Aequ(0,0);
     if(neq>0)
     {
@@ -218,15 +230,27 @@ bool igl::min_quad_with_fixed_precompute(
   {
     Eigen::SparseMatrix<T> NA;
     igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
-    assert(false);
-    return false;
+    // Build LU decomposition of NA
+    data.sparse = false;
+    fprintf(stderr,
+      "Warning: min_quad_with_fixed_precompute() resorting to dense LU\n");
+    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NAfull;
+    igl::full(NA,NAfull);
+    data.lu = 
+      Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >(NAfull);
+    if(!data.lu.isInvertible())
+    {
+      fprintf(stderr,
+        "Error: min_quad_with_fixed_precompute() LU not invertible\n");
+      return false;
+    }
   }
   return true;
 }
 
 
 template <typename T>
-bool igl::min_quad_with_fixed_solve(
+inline bool igl::min_quad_with_fixed_solve(
   const igl::min_quad_with_fixed_data<T> & data,
   const Eigen::Matrix<T,Eigen::Dynamic,1> & B,
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Y,
@@ -245,10 +269,6 @@ bool igl::min_quad_with_fixed_solve(
   // number of lagrange multipliers aka linear equality constraints
   int neq = data.lagrange.size();
 
-  if(neq != 0)
-  {
-  }
-
   // append lagrange multiplier rhs's
   Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
   BBeq << B, (Beq*-2.0);
@@ -259,7 +279,13 @@ bool igl::min_quad_with_fixed_solve(
   Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> BBequlcols;
   igl::repmat(BBequl,1,cols,BBequlcols);
   Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
-  NB = data.preY * Y + BBequlcols;
+  if(kr == 0)
+  {
+    NB = BBequlcols;
+  }else
+  {
+    NB = data.preY * Y + BBequlcols;
+  }
 
   // resize output
   Z.resize(data.n,cols);
@@ -272,8 +298,23 @@ bool igl::min_quad_with_fixed_solve(
       Z(data.known(i),j) = Y(i,j);
     }
   }
-  data.L.template triangularView<Lower>().solveInPlace(NB);
-  data.U.template triangularView<Upper>().solveInPlace(NB);
+
+  //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
+
+  if(data.sparse)
+  {
+    //std::cout<<"data.LIJV=["<<std::endl;print_ijv(data.L,1);std::cout<<std::endl<<"];"<<
+    //  std::endl<<"data.L=sparse(data.LIJV(:,1),data.LIJV(:,2),data.LIJV(:,3),"<<
+    //  data.L.rows()<<","<<data.L.cols()<<");"<<std::endl;
+    //std::cout<<"data.UIJV=["<<std::endl;print_ijv(data.U,1);std::cout<<std::endl<<"];"<<
+    //  std::endl<<"data.U=sparse(data.UIJV(:,1),data.UIJV(:,2),data.UIJV(:,3),"<<
+    //  data.U.rows()<<","<<data.U.cols()<<");"<<std::endl;
+    data.L.template triangularView<Eigen::Lower>().solveInPlace(NB);
+    data.U.template triangularView<Eigen::Upper>().solveInPlace(NB);
+  }else
+  {
+    NB = data.lu.solve(NB);
+  }
   // Now NB contains sol/-0.5
   NB *= -0.5;
   // Now NB contains solution

normalize_rows.h

@@ -10,10 +10,10 @@ namespace igl
   //   A  #rows by k input matrix
   // Outputs:
   //   B  #rows by k input matrix, can be the same as A
-  void normalize_rows(const Eigen::MatrixXd & A, Eigen::MatrixXd & B);
+  inline void normalize_rows(const Eigen::MatrixXd & A, Eigen::MatrixXd & B);
 }
 
-void igl::normalize_rows(const Eigen::MatrixXd & A, Eigen::MatrixXd & B)
+inline void igl::normalize_rows(const Eigen::MatrixXd & A, Eigen::MatrixXd & B)
 {
   // Resize output
   B.resize(A.rows(),A.cols());

pathinfo.h

@@ -23,7 +23,7 @@ namespace igl
   //    '.')
   //
   //  See also: basename, dirname
-  void pathinfo(
+  inline void pathinfo(
     const std::string & path,
     std::string & dirname,
     std::string & basename,
@@ -38,7 +38,7 @@ namespace igl
 // Verbose should be removed once everythings working correctly
 #include "verbose.h"
 
-void igl::pathinfo(
+inline void igl::pathinfo(
   const std::string & path,
   std::string & dirname,
   std::string & basename,

per_face_normals.h

@@ -9,14 +9,14 @@ namespace igl
   //   F  #F by 3 eigne Matrix of face (triangle) indices
   // Output:
   //   N  #F by 3 eigen Matrix of mesh face (triangle) 3D normals
-  void per_face_normals(
+  inline void per_face_normals(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     Eigen::MatrixXd & N);
 }
 // Implementation
 
-void igl::per_face_normals(
+inline void igl::per_face_normals(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & N)

per_vertex_normals.h

@@ -13,7 +13,7 @@ namespace igl
   //   F  #F by 3 eigne Matrix of face (triangle) indices
   // Output:
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
-  void per_vertex_normals(
+  inline void per_vertex_normals(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     Eigen::MatrixXd & N);
@@ -23,7 +23,7 @@ namespace igl
 #include "per_face_normals.h"
 #include "normalize_rows.h"
 
-void igl::per_vertex_normals(
+inline void igl::per_vertex_normals(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
   Eigen::MatrixXd & N)

print_gl_get.h

@@ -12,13 +12,13 @@ namespace igl
   // Prints the value of pname found by issuing glGet*(pname,*)
   // Inputs:
   //   pname  enum key to gl parameter
-  void print_gl_get(GLenum pname);
+  inline void print_gl_get(GLenum pname);
 }
 
 
 // implementation
 #include <cstdio>
-void igl::print_gl_get(GLenum pname)
+inline void igl::print_gl_get(GLenum pname)
 {
   double dM[16];
 

quat_to_mat.h

@@ -10,13 +10,13 @@ namespace igl
   // Output:
   //   mat  pointer to 16 elements of matrix
   template <typename Q_type>
-  void quat_to_mat(const Q_type * quat, Q_type * mat);
+  inline void quat_to_mat(const Q_type * quat, Q_type * mat);
 }
 
 // Implementation
 
 template <typename Q_type>
-void igl::quat_to_mat(const Q_type * quat, Q_type * mat)
+inline void igl::quat_to_mat(const Q_type * quat, Q_type * mat)
 {
   Q_type yy2 = 2.0f * quat[1] * quat[1];
   Q_type xy2 = 2.0f * quat[0] * quat[1];

read.h

@@ -21,13 +21,13 @@ namespace igl
   // Outputs:
   //   V  eigen double matrix #V by 3
   //   F  eigen int matrix #F by 3
-  bool read(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
+  inline bool read(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
 }
 
 // Implementation
 #include <readOBJ.h>
 #include <readOFF.h>
-bool igl::read(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+inline bool igl::read(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
 {
     const char* p;
     for (p = str.c_str(); *p != '\0'; p++)

readOBJ.h

@@ -29,11 +29,11 @@ namespace igl
   //
   // KNOWN BUG: This only knows how to face lines without normal or texture
   // indices. It will probably crash or give garbage on anything else.
-  bool readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
+  inline bool readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
 }
 
 // Implementation
-bool igl::readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+inline bool igl::readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
 {
     std::ifstream s(str.c_str());
     if (s.is_open() == false)

readOFF.h

@@ -20,12 +20,12 @@ namespace igl
   // Outputs:
   //   V  eigen double matrix #V by 3
   //   F  eigen int matrix #F by 3
-  bool readOFF (const std::string meshfile, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
+  inline bool readOFF (const std::string meshfile, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
 }
 
 // Implementation
 
-bool igl::readOFF (const std::string meshfile, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+inline bool igl::readOFF (const std::string meshfile, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
 {
     int vnum, fnum;
     FILE *fp = fopen (meshfile.c_str(), "r");

repdiag.h

@@ -24,13 +24,13 @@ namespace igl
 
   // Sparse version
   template <typename T>
-  void repdiag(
+  inline void repdiag(
     const Eigen::SparseMatrix<T>& A,
     const int d,
     Eigen::SparseMatrix<T>& B);
   // Dense version
   template <typename T>
-  void repdiag(
+  inline void repdiag(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & A,
     const int d,
     Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
@@ -39,7 +39,7 @@ namespace igl
 // Implementation
 
 template <typename T>
-void igl::repdiag(
+inline void igl::repdiag(
   const Eigen::SparseMatrix<T>& A,
   const int d,
   Eigen::SparseMatrix<T>& B)
@@ -57,7 +57,7 @@ void igl::repdiag(
     for (int k=0; k<A.outerSize(); ++k)
     {
       // loop inner level
-      for (typename SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
+      for (typename Eigen::SparseMatrix<T>::InnerIterator it(A,k); it; ++it)
       {
         dyn_B.coeffRef(i*m+it.row(),i*n+it.col()) += it.value();
       }
@@ -68,7 +68,7 @@ void igl::repdiag(
 }
 
 template <typename T>
-void igl::repdiag(
+inline void igl::repdiag(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & A,
   const int d,
   Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B)

repmat.h

@@ -27,6 +27,12 @@ namespace igl
     const int r,
     const int c,
     Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
+  template <typename T>
+  inline void repmat(
+    const Eigen::SparseMatrix<T> & A,
+    const int r,
+    const int c,
+    Eigen::SparseMatrix<T> & B);
 }
 
 // Implementation
@@ -52,4 +58,33 @@ inline void igl::repmat(
     }
   }
 }
+
+template <typename T>
+inline void igl::repmat(
+  const Eigen::SparseMatrix<T> & A,
+  const int r,
+  const int c,
+  Eigen::SparseMatrix<T> & 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++)
+  {
+    for(int j = 0;j<c;j++)
+    {
+      // 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)
+        {
+          B.insert(i*A.rows()+it.row(),j*A.cols() + it.col()) = it.value();
+        }
+      }
+    }
+  }
+  B.finalize();
+}
 #endif

sparse.h

@@ -1,5 +1,6 @@
 #ifndef IGL_SPARSE_H
 #define IGL_SPARSE_H
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 namespace igl

stdin_to_temp.h

@@ -18,13 +18,13 @@ namespace igl
   //
   // Note: Caller is responsible for closing the file (tmpfile() automatically
   // unlinks the file so there is no need to remove/delete/unlink the file)
-  bool stdin_to_temp(FILE ** temp_file);
+  inline bool stdin_to_temp(FILE ** temp_file);
 }
 
 // IMPLEMENTATION
 #include <iostream>
 
-bool igl::stdin_to_temp(FILE ** temp_file)
+inline bool igl::stdin_to_temp(FILE ** temp_file)
 {
   // get a temporary file
   *temp_file = tmpfile();

sum.h

@@ -18,14 +18,14 @@ namespace igl
   //   or
   //   S  m-long sparse vector (if dim == 2)
   template <typename T>
-  void sum(
+  inline void sum(
     const Eigen::SparseMatrix<T>& X, 
     const int dim,
     Eigen::SparseVector<T>& S);
 }
 
 template <typename T>
-void igl::sum(
+inline void igl::sum(
   const Eigen::SparseMatrix<T>& X, 
   const int dim,
   Eigen::SparseVector<T>& S)

trackball.h

@@ -15,7 +15,7 @@ namespace igl
   // Outputs:
   //   quat  the resulting rotation (as quaternion)
   template <typename Q_type>
-  void trackball(
+  inline void trackball(
     const int w,
     const int h,
     const Q_type speed_factor,
@@ -58,7 +58,7 @@ static inline Q_type _QuatIY(int y, int w, int h)
 // code is straight from its source in TwMgr.cpp
 // http://www.antisphere.com/Wiki/tools:anttweakbar
 template <typename Q_type>
-void igl::trackball(
+inline void igl::trackball(
   const int w,
   const int h,
   const Q_type speed_factor,

transpose_blocks.h

@@ -39,7 +39,7 @@ namespace igl
   // 204 208];
   //   
   template <typename T>
-  void transpose_blocks(
+  inline void transpose_blocks(
     const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & A,
     const size_t k,
     const size_t dim,
@@ -50,7 +50,7 @@ namespace igl
 #include <cassert>
 
 template <typename T>
-void igl::transpose_blocks(
+inline void igl::transpose_blocks(
   const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & A,
   const size_t k,
   const size_t dim,

tt.h

@@ -15,7 +15,7 @@
 namespace igl 
 {
     // Compute triangle-triangle adjacency
-    void tt(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXi& TT)
+    inline void tt(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXi& TT)
     {
         assert(igl::isManifold(V,F));
         std::vector<std::vector<int> > TTT;
@@ -47,4 +47,4 @@ namespace igl
     }
 }
 
-#endif
+#endif

write.h

@@ -15,7 +15,7 @@
 namespace igl 
 {
     // write mesh to an ascii file with automatic detection of file format. supported: obj, off)
-    void write(std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+    inline void write(std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
     {
         const char* p;
         for (p = str.c_str(); *p != '\0'; p++)
@@ -31,4 +31,4 @@ namespace igl
     }
 }
 
-#endif
+#endif

writeOBJ.h

@@ -20,7 +20,7 @@ namespace igl
   //   V  eigen double matrix #V by 3 (mesh vertices)
   //   F  eigen int matrix #F by 3 (mesh indices)
   // Returns true on success, false on error
-  bool writeOBJ(const std::string str, const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
+  inline bool writeOBJ(const std::string str, const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
 }
 
 // Implementation
@@ -28,7 +28,7 @@ namespace igl
 #include <fstream>
 #include <cstdio>
 
-bool igl::writeOBJ(const std::string str, const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
+inline bool igl::writeOBJ(const std::string str, const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
 {
   std::ofstream s(str.c_str());
 

writeOFF.h

@@ -12,7 +12,7 @@
 namespace igl 
 {
     // write mesh to an ascii off file
-    void writeOFF (std::string fname, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+    inline void writeOFF (std::string fname, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
     {
         FILE *fp = fopen (fname.c_str(), "w");
         
@@ -31,4 +31,4 @@ namespace igl
     }
 }
 
-#endif
+#endif