more explicits, transitioning min_quad_with_fixed (sparse), lu lagrange is broken for constraints, is_symmetric with threshold, full from sparse

Former-commit-id: 6c8fe8299de8064639f19f5749813b35c7e5bd15

include/igl/full.cpp

@@ -1,16 +1,17 @@
 #include "full.h"
 
-template <typename T>
+template <typename T,typename DerivedB>
 IGL_INLINE void igl::full(
   const Eigen::SparseMatrix<T> & A,
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B)
+  Eigen::PlainObjectBase<DerivedB>& B)
 {
-  B = Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>::Zero(A.rows(),A.cols());
+  using namespace Eigen;
+  B = PlainObjectBase<DerivedB >::Zero(A.rows(),A.cols());
   // Iterate over outside
   for(int k=0; k<A.outerSize(); ++k)
   {
     // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
+    for(typename SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
     {
       B(it.row(),it.col()) = it.value();
     }
@@ -30,5 +31,5 @@ IGL_INLINE void igl::full(
 // generated by autoexplicit.sh
 template void igl::full<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
-template void igl::full<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::full<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/full.h

@@ -13,10 +13,10 @@ namespace igl
   //   A  m by n sparse matrix
   // Output:
   //   B  m by n dense/full matrix
-  template <typename T>
+  template <typename T,typename DerivedB>
   IGL_INLINE void full(
     const Eigen::SparseMatrix<T> & A,
-    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & B);
+    Eigen::PlainObjectBase<DerivedB> & B);
   // If already full then this will just be a copy by assignment
   template <typename DerivedA,typename DerivedB>
   IGL_INLINE void full(

include/igl/is_symmetric.cpp

@@ -1,4 +1,5 @@
 #include "is_symmetric.h"
+#include "find.h"
 
 template <typename T>
 IGL_INLINE bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
@@ -33,10 +34,32 @@ IGL_INLINE bool igl::is_symmetric(
   return AmAT.nonZeros() == 0;
 }
 
+template <typename AType, typename epsilonT>
+IGL_INLINE bool igl::is_symmetric(
+  const Eigen::SparseMatrix<AType>& A, 
+  const epsilonT epsilon)
+{
+  using namespace Eigen;
+  using namespace igl;
+  if(A.rows() != A.cols())
+  {
+    return false;
+  }
+  SparseMatrix<AType> AT = A.transpose();
+  SparseMatrix<AType> AmAT = A-AT;
+  VectorXi AmATI,AmATJ;
+  Matrix<AType,Dynamic,1> AmATV;
+  find(AmAT,AmATI,AmATJ,AmATV);
+  
+  return AmATV.maxCoeff() < epsilon && AmATV.minCoeff() > -epsilon;
+}
+
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
 template bool igl::is_symmetric<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
 // generated by autoexplicit.sh
 template bool igl::is_symmetric<double>(Eigen::SparseMatrix<double, 0, int> const&);
+template bool igl::is_symmetric<double, double>(Eigen::SparseMatrix<double, 0, int> const&, double);
+template bool igl::is_symmetric<double, int>(Eigen::SparseMatrix<double, 0, int> const&, int);
 #endif

include/igl/is_symmetric.h

@@ -9,8 +9,12 @@ namespace igl
   // Inputs:
   //   A  m by m matrix
   // Returns true if the matrix is square and symmetric
-  template <typename T>
-  IGL_INLINE bool is_symmetric(const Eigen::SparseMatrix<T>& A);
+  template <typename AT>
+  IGL_INLINE bool is_symmetric(const Eigen::SparseMatrix<AT>& A);
+  // Inputs:
+  //   epsilon threshold on L1 difference between A and A'
+  template <typename AT, typename epsilonT>
+  IGL_INLINE bool is_symmetric(const Eigen::SparseMatrix<AT>& A, const epsilonT epsilon);
   template <typename DerivedA>
   IGL_INLINE bool is_symmetric(
     const Eigen::PlainObjectBase<DerivedA>& A);

include/igl/lu_lagrange.cpp

@@ -143,4 +143,5 @@ IGL_INLINE bool igl::lu_lagrange(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template bool igl::lu_lagrange<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/lu_lagrange.h

@@ -7,6 +7,8 @@
 #include <Eigen/Sparse>
 namespace igl
 {
+  // KNOWN BUGS: This does not seem to be correct for non-empty C
+  //
   // LU_LAGRANGE Compute a LU decomposition for a special type of
   // matrix Q that is symmetric but not positive-definite:
   // Q = [A'*A C

include/igl/matlab/MatlabWorkspace.cpp

@@ -23,6 +23,8 @@ IGL_INLINE igl::MatlabWorkspace::~MatlabWorkspace()
 IGL_INLINE void igl::MatlabWorkspace::clear()
 {
   for_each(data.begin(),data.end(),&mxDestroyArray);
+  data.clear();
+  names.clear();
 }
 
 IGL_INLINE bool igl::MatlabWorkspace::write(const std::string & path) const
@@ -50,6 +52,74 @@ IGL_INLINE bool igl::MatlabWorkspace::write(const std::string & path) const
   return true;
 }
 
+IGL_INLINE bool igl::MatlabWorkspace::read(const std::string & path)
+{
+  using namespace std;
+
+  MATFile * mat_file;
+
+  mat_file = matOpen(path.c_str(), "r");
+  if (mat_file == NULL) 
+  {
+    cerr<<"Error: failed to open "<<path<<endl;
+    return false;
+  }
+
+  int ndir;
+  const char ** dir = (const char **)matGetDir(mat_file, &ndir);
+  if (dir == NULL) {
+    cerr<<"Error reading directory of file "<< path<<endl;
+    return false;
+  }
+  mxFree(dir);
+
+  // Must close and reopen
+  if(matClose(mat_file) != 0)
+  {
+    cerr<<"Error: failed to close file "<<path<<endl;
+    return false;
+  }
+  mat_file = matOpen(path.c_str(), "r");
+  if (mat_file == NULL) 
+  {
+    cerr<<"Error: failed to open "<<path<<endl;
+    return false;
+  }
+  
+
+  /* Read in each array. */
+  for (int i=0; i<ndir; i++) 
+  {
+    const char * name;
+    mxArray * mx_data = matGetNextVariable(mat_file, &name);
+    if (mx_data == NULL) 
+    {
+      cerr<<"Error: matGetNextVariable failed in "<<path<<endl;
+      return false;
+    } 
+    const int dims = mxGetNumberOfDimensions(mx_data);
+    assert(dims == 2);
+    if(dims != 2)
+    {
+      fprintf(stderr,"Variable '%s' has %d ≠ 2 dimensions. Skipping\n",
+          name,dims);
+      mxDestroyArray(mx_data);
+      continue;
+    }
+    // don't destroy
+    names.push_back(name);
+    data.push_back(mx_data);
+  }
+
+  if(matClose(mat_file) != 0)
+  {
+    cerr<<"Error: failed to close file "<<path<<endl;
+    return false;
+  }
+
+  return true;
+}
+
 // Treat everything as a double
 template <typename DerivedM>
 IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
@@ -157,6 +227,102 @@ IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save_index(
   return this->save_index(V,name);
 }
 
+template <typename DerivedM>
+IGL_INLINE bool igl::MatlabWorkspace::find( 
+  const std::string & name,
+  Eigen::PlainObjectBase<DerivedM>& M)
+{
+  using namespace std;
+  const int i = std::find(names.begin(), names.end(), name)-names.begin();
+  if(i>=(int)names.size())
+  {
+    return false;
+  }
+  assert(i<=(int)data.size());
+  mxArray * mx_data = data[i];
+  assert(!mxIsSparse(mx_data));
+  assert(mxGetNumberOfDimensions(mx_data) == 2);
+  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
+  const int m = mxGetM(mx_data);
+  const int n = mxGetN(mx_data);
+
+  // Handle vectors
+  if(DerivedM::IsVectorAtCompileTime)
+  {
+    assert(m==1 || n==1);
+    M.resize(m*n);
+  }else
+  {
+    M.resize(m,n);
+  }
+  copy(
+    mxGetPr(mx_data),
+    mxGetPr(mx_data)+mxGetNumberOfElements(mx_data),
+    M.data());
+  return true;
+}
+
+template <typename MT>
+IGL_INLINE bool igl::MatlabWorkspace::find( 
+  const std::string & name,
+  Eigen::SparseMatrix<MT>& M)
+{
+  using namespace std;
+  using namespace Eigen;
+  const int i = std::find(names.begin(), names.end(), name)-names.begin();
+  if(i>=(int)names.size())
+  {
+    return false;
+  }
+  assert(i<=(int)data.size());
+  mxArray * mx_data = data[i];
+  assert(mxIsSparse(mx_data));
+  assert(mxGetNumberOfDimensions(mx_data) == 2);
+  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
+  const int m = mxGetM(mx_data);
+  const int n = mxGetN(mx_data);
+
+  // Copy data immediately
+  double * pr = mxGetPr(mx_data);
+  mwIndex * ir = mxGetIr(mx_data);
+  mwIndex * jc = mxGetJc(mx_data);
+
+  vector<Triplet<MT> > MIJV;
+  MIJV.reserve(mxGetNumberOfElements(mx_data));
+
+  // Iterate over outside
+  int k = 0;
+  for(int j=0; j<n;j++)
+  {
+    // Iterate over inside
+    while(k<(int)jc[j+1])
+    {
+      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
+      assert((int)ir[k]<m);
+      assert((int)j<n);
+      MIJV.push_back(Triplet<MT >(ir[k],j,pr[k]));
+      k++;
+    }
+  }
+  M.resize(m,n);
+  M.setFromTriplets(MIJV.begin(),MIJV.end());
+  return true;
+
+}
+
+template <typename DerivedM>
+IGL_INLINE bool igl::MatlabWorkspace::find_index( 
+  const std::string & name,
+  Eigen::PlainObjectBase<DerivedM>& M)
+{
+  if(!find(name,M))
+  {
+    return false;
+  }
+  M.array() -= 1;
+  return true;
+}
+
 
 //template <typename Data>
 //bool igl::MatlabWorkspace::save(const Data & M, const std::string & name)

include/igl/matlab/MatlabWorkspace.h

@@ -24,6 +24,8 @@ namespace igl
   class MatlabWorkspace
   {
     private:
+      // KNOWN BUG: Why not use a map? Any reason to allow duplicate names?
+      //
       // List of names
       std::vector<std::string> names;
       // List of data pointers
@@ -39,6 +41,12 @@ namespace igl
       //   path  path to .mat file
       // Returns true on success, false on failure
       IGL_INLINE bool write(const std::string & path) const;
+      // Load list of variables from .mat file
+      //
+      // Inputs:
+      //   path  path to .mat file
+      // Returns true on success, false on failure
+      IGL_INLINE bool read(const std::string & path);
       // Assign data to a variable name in the workspace
       //
       // Template: 
@@ -85,6 +93,30 @@ namespace igl
       IGL_INLINE MatlabWorkspace& save_index(
         const std::vector<ScalarV> & vV,
         const std::string & name);
+      // Find a certain matrix by name.
+      //
+      // KNOWN BUG: Outputs the first found (not necessarily unique lists).
+      //
+      // Template: 
+      //   DerivedM  eigen matrix (e.g. MatrixXd)
+      // Inputs:
+      //   name  exact name of matrix as string
+      // Outputs:
+      //   M  matrix
+      // Returns true only if found.
+      template <typename DerivedM>
+      IGL_INLINE bool find( 
+        const std::string & name,
+        Eigen::PlainObjectBase<DerivedM>& M);
+      template <typename MT>
+      IGL_INLINE bool find( 
+        const std::string & name,
+        Eigen::SparseMatrix<MT>& M);
+      // Subtracts 1 from all entries
+      template <typename DerivedM>
+      IGL_INLINE bool find_index( 
+        const std::string & name,
+        Eigen::PlainObjectBase<DerivedM>& M);
   };
 }
 

include/igl/matlab_format.cpp

@@ -93,4 +93,5 @@ template Eigen::WithFormat<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const igl:
 template Eigen::WithFormat<Eigen::Array<int, -1, -1, 0, -1, -1> > const igl::matlab_format<Eigen::Array<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Array<int, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const igl::matlab_format<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 #endif

include/igl/min_quad_with_fixed.cpp

@@ -1,13 +1,5 @@
 #include "min_quad_with_fixed.h"
 
-//#include <Eigen/SparseExtra>
-// Bug in unsupported/Eigen/SparseExtra needs iostream first
-#include <iostream>
-#include <unsupported/Eigen/SparseExtra>
-#include <cassert>
-#include <cstdio>
-#include <iostream>
-
 #include "slice.h"
 #include "is_symmetric.h"
 #include "find.h"
@@ -15,6 +7,15 @@
 #include "repmat.h"
 #include "lu_lagrange.h"
 #include "full.h"
+#include "EPS.h"
+
+//#include <Eigen/SparseExtra>
+// Bug in unsupported/Eigen/SparseExtra needs iostream first
+#include <iostream>
+#include <unsupported/Eigen/SparseExtra>
+#include <cassert>
+#include <cstdio>
+#include <iostream>
 
 template <typename T>
 IGL_INLINE bool igl::min_quad_with_fixed_precompute(
@@ -25,6 +26,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   igl::min_quad_with_fixed_data<T> & data
   )
 {
+  using namespace Eigen;
+  using namespace std;
   // number of rows
   int n = A.rows();
   // cache problem size
@@ -76,31 +79,31 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   data.unknown_lagrange.resize(data.unknown.size()+data.lagrange.size());
   data.unknown_lagrange << data.unknown, data.lagrange;
 
-  Eigen::SparseMatrix<T> Auu;
+  SparseMatrix<T> Auu;
   igl::slice(A,data.unknown,data.unknown,Auu);
 
   // determine if A(unknown,unknown) is symmetric and/or positive definite
-  data.Auu_sym = igl::is_symmetric(Auu);
+  data.Auu_sym = igl::is_symmetric(Auu,igl::EPS<double>());
   // Positive definiteness is *not* determined, rather it is given as a
   // parameter
   data.Auu_pd = pd;
 
   // Append lagrange multiplier quadratic terms
-  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;
+  SparseMatrix<T> new_A;
+  Matrix<int,Dynamic,1> AI;
+  Matrix<int,Dynamic,1> AJ;
+  Matrix<T,Dynamic,1> AV;
   igl::find(A,AI,AJ,AV);
-  Eigen::Matrix<int,Eigen::Dynamic,1> AeqI(0,1);
-  Eigen::Matrix<int,Eigen::Dynamic,1> AeqJ(0,1);
-  Eigen::Matrix<T,Eigen::Dynamic,1> AeqV(0,1);
+  Matrix<int,Dynamic,1> AeqI(0,1);
+  Matrix<int,Dynamic,1> AeqJ(0,1);
+  Matrix<T,Dynamic,1> AeqV(0,1);
   if(neq > 0)
   {
     igl::find(Aeq,AeqI,AeqJ,AeqV);
   }
-  Eigen::Matrix<int,Eigen::Dynamic,1> new_AI(AV.size()+AeqV.size()*2);
-  Eigen::Matrix<int,Eigen::Dynamic,1> new_AJ(AV.size()+AeqV.size()*2);
-  Eigen::Matrix<T,Eigen::Dynamic,1>   new_AV(AV.size()+AeqV.size()*2);
+  Matrix<int,Dynamic,1> new_AI(AV.size()+AeqV.size()*2);
+  Matrix<int,Dynamic,1> new_AJ(AV.size()+AeqV.size()*2);
+  Matrix<T,Dynamic,1>   new_AV(AV.size()+AeqV.size()*2);
   new_AI << AI, (AeqI.array()+n).matrix(), AeqJ;
   new_AJ << AJ, AeqJ, (AeqI.array()+n).matrix();
   new_AV << AV, AeqV, AeqV;
@@ -109,14 +112,14 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   // precompute RHS builders
   if(kr > 0)
   {
-    Eigen::SparseMatrix<T> Aulk;
+    SparseMatrix<T> Aulk;
     igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
-    Eigen::SparseMatrix<T> Akul;
+    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();
+    SparseMatrix<T> AkulT = Akul.transpose();
     data.preY = Aulk + AkulT;
   }else
   {
@@ -124,46 +127,60 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   }
 
   // Create factorization
-  if(data.Auu_sym && data.Auu_pd)
+  //if(data.Auu_sym && data.Auu_pd)
+  //{
+  //  data.sparse = true;
+  //  SparseMatrix<T> Aequ(0,0);
+  //  if(neq>0)
+  //  {
+  //    Matrix<int,Dynamic,1> Aeq_rows;
+  //    Aeq_rows.resize(neq);
+  //    for(int i = 0;i<neq;i++)
+  //    {
+  //      Aeq_rows(i) = i;
+  //    }
+  //    igl::slice(Aeq,Aeq_rows,data.unknown,Aequ);
+  //  }
+  //  // Get transpose of Aequ
+  //  SparseMatrix<T> AequT = Aequ.transpose();
+  //  // Compute LU decomposition
+  //  // TODO: This should be using SimplicialLDLT
+  //  // Q: Does SimplicialLDLT work on "Hermitian indefinite matrices" the way
+  //  // that matlabs ldl does?
+  //  // A: Maybe not. The eigen documentation writes, "This class provides a
+  //  // LDL^T Cholesky factorizations without square root of sparse matrices
+  //  // that are selfadjoint and positive definite"
+  //  bool lu_success = igl::lu_lagrange(Auu,AequT,data.L,data.U);
+  //  if(!lu_success)
+  //  {
+  //    return false;
+  //  }
+  //}else
+  //{
+  SparseMatrix<T> NA;
+  igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
+  if(data.Auu_sym)
   {
-    data.sparse = true;
-    Eigen::SparseMatrix<T> Aequ(0,0);
-    if(neq>0)
-    {
-      Eigen::Matrix<int,Eigen::Dynamic,1> Aeq_rows;
-      Aeq_rows.resize(neq);
-      for(int i = 0;i<neq;i++)
-      {
-        Aeq_rows(i) = i;
-      }
-      igl::slice(Aeq,Aeq_rows,data.unknown,Aequ);
-    }
-    // Get transpose of Aequ
-    Eigen::SparseMatrix<T> AequT = Aequ.transpose();
-    // Compute LU decomposition
-    // TODO: This should be using SimplicialLDLT
-    // Q: Does SimplicialLDLT work on "Hermitian indefinite matrices" the way
-    // that matlabs ldl does?
-    // A: Maybe not. The eigen documentation writes, "This class provides a
-    // LDL^T Cholesky factorizations without square root of sparse matrices
-    // that are selfadjoint and positive definite"
-    bool lu_success = igl::lu_lagrange(Auu,AequT,data.L,data.U);
-    if(!lu_success)
+    SimplicialLDLT<SparseMatrix<T> > ldlt_solver;
+    ldlt_solver.compute(NA);
+    if(ldlt_solver.info() != Eigen::Success)
     {
       return false;
     }
+    // Implementation incomplete
+    assert(false);
+    cerr<<"ERROR min_quad_with_fixed_precompute implementation incomplete"<<endl;
+    return false;
   }else
   {
-    Eigen::SparseMatrix<T> NA;
-    igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
     // 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;
+    Matrix<T,Dynamic,Dynamic> NAfull;
     igl::full(NA,NAfull);
     data.lu = 
-      Eigen::FullPivLU<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >(NAfull);
+      FullPivLU<Matrix<T,Dynamic,Dynamic> >(NAfull);
     if(!data.lu.isInvertible())
     {
       fprintf(stderr,
@@ -171,17 +188,19 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
       return false;
     }
   }
+  //}
+
   return true;
 }
 
 
-template <typename T>
+template <typename T, typename DerivedY, typename DerivedZ>
 IGL_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,
+  const Eigen::PlainObjectBase<DerivedY> & Y,
   const Eigen::Matrix<T,Eigen::Dynamic,1> & Beq,
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Z)
+  Eigen::PlainObjectBase<DerivedZ> & Z)
 {
   // number of known rows
   int kr = data.known.size();
@@ -257,4 +276,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::min_quad_with_fixed_precompute<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int> const&, bool, igl::min_quad_with_fixed_data<double>&);
 #endif
+

include/igl/min_quad_with_fixed.h

@@ -44,18 +44,23 @@ namespace igl
     );
 
   // Solves a system previously factored using min_quad_with_fixed_precompute
+  // 
+  // Template:
+  //   T  type of sparse matrix (e.g. double) 
+  //   DerivedY  type of Y (e.g. derived from VectorXd or MatrixXd)
+  //   DerivedZ  type of Z (e.g. derived from VectorXd or MatrixXd)
   // Inputs:
   //   data  factorization struct with all necessary precomputation to solve
   // Outputs:
   //   Z  n by cols solution
   // Returns true on success, false on error
-  template <typename T>
+  template <typename T,typename DerivedY,typename DerivedZ>
   IGL_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,
+    const Eigen::PlainObjectBase<DerivedY> & Y,
     const Eigen::Matrix<T,Eigen::Dynamic,1> & Beq,
-    Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> & Z);
+    Eigen::PlainObjectBase<DerivedZ> & Z);
 }
 
 template <typename T>

include/igl/repmat.cpp

@@ -1,11 +1,11 @@
 #include "repmat.h"
 
-template <typename DerivedA>
+template <typename DerivedA, typename DerivedB>
 IGL_INLINE void igl::repmat(
   const Eigen::PlainObjectBase<DerivedA> & A,
   const int r,
   const int c,
-  Eigen::PlainObjectBase<DerivedA> & B)
+  Eigen::PlainObjectBase<DerivedB> & B)
 {
   assert(r>0);
   assert(c>0);
@@ -54,6 +54,5 @@ IGL_INLINE void igl::repmat(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::repmat<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::repmat<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::repmat<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/repmat.h

@@ -22,12 +22,12 @@ namespace igl
   // Outputs:
   //   B  r*m by c*n output matrix
   //
-  template <typename DerivedA>
+  template <typename DerivedA,typename DerivedB>
   IGL_INLINE void repmat(
     const Eigen::PlainObjectBase<DerivedA> & A,
     const int r,
     const int c,
-    Eigen::PlainObjectBase<DerivedA> & B);
+    Eigen::PlainObjectBase<DerivedB> & B);
   template <typename T>
   IGL_INLINE void repmat(
     const Eigen::SparseMatrix<T> & A,

include/igl/sparse.cpp

@@ -2,6 +2,7 @@
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
+#include <vector>
 #include <unsupported/Eigen/SparseExtra>
 
 template <class IndexVector, class ValueVector, typename T>
@@ -26,6 +27,8 @@ IGL_INLINE void igl::sparse(
   const size_t n,
   Eigen::SparseMatrix<T>& X)
 {
+  using namespace std;
+  using namespace Eigen;
   assert((int)I.maxCoeff() < (int)m);
   assert((int)I.minCoeff() >= 0);
   assert((int)J.maxCoeff() < (int)n);
@@ -37,21 +40,64 @@ IGL_INLINE void igl::sparse(
   assert(J.rows() == V.rows());
   assert(I.cols() == J.cols());
   assert(J.cols() == V.cols());
-  // number of values
-  int nv = V.size();
+  //// number of values
+  //int nv = V.size();
 
-  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++)
+  //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);
+  //}
+  //X = Eigen::SparseMatrix<T>(dyn_X);
+  vector<Triplet<T> > IJV;
+  IJV.reserve(I.size());
+  for(int x = 0;x<I.size();x++)
   {
-    dyn_X.coeffRef((int)I(i),(int)J(i)) += (T)V(i);
+    IJV.push_back(Triplet<T >(I(x),J(x),V(x)));
   }
-  X = Eigen::SparseMatrix<T>(dyn_X);
+  X.resize(m,n);
+  X.setFromTriplets(IJV.begin(),IJV.end());
+}
 
+template <typename DerivedD, typename T>
+IGL_INLINE void igl::sparse(
+  const Eigen::PlainObjectBase<DerivedD>& D,
+  Eigen::SparseMatrix<T>& X)
+{
+  using namespace std;
+  using namespace Eigen;
+  vector<Triplet<T> > DIJV;
+  const int m = D.rows();
+  const int n = D.cols();
+  for(int i = 0;i<m;i++)
+  {
+    for(int j = 0;j<n;j++)
+    {
+      if(D(i,j)!=0)
+      {
+        DIJV.push_back(Triplet<T>(i,j,D(i,j)));
+      }
+    }
+  }
+  X.resize(m,n);
+  X.setFromTriplets(DIJV.begin(),DIJV.end());
+}
+
+template <typename DerivedD>
+IGL_INLINE Eigen::SparseMatrix<typename DerivedD::Scalar > igl::sparse(
+  const Eigen::PlainObjectBase<DerivedD>& D)
+{
+  Eigen::SparseMatrix<typename DerivedD::Scalar > X;
+  igl::sparse(D,X);
+  return X;
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 template void igl::sparse<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double>(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, size_t, size_t, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::sparse<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int> igl::sparse<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
+template Eigen::SparseMatrix<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, 0, int> igl::sparse<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&);
 #endif

include/igl/sparse.h

@@ -40,6 +40,23 @@ namespace igl
     const size_t m,
     const size_t n,
     Eigen::SparseMatrix<T>& X);
+  // Convert a full, dense matrix to a sparse one
+  //
+  // Templates:
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   D  m by n full, dense matrix
+  // Output:
+  //   X  m by n sparse matrix
+  template <typename DerivedD, typename T>
+  IGL_INLINE void sparse(
+    const Eigen::PlainObjectBase<DerivedD>& D,
+    Eigen::SparseMatrix<T>& X);
+  // Wrapper with return
+  template <typename DerivedD>
+  IGL_INLINE Eigen::SparseMatrix<typename DerivedD::Scalar > sparse(
+    const Eigen::PlainObjectBase<DerivedD>& D);
+
 }
 
 #ifdef IGL_HEADER_ONLY

todos.txt

@@ -25,3 +25,4 @@
 - svd.* depends on lapack, should use eigen...
 - use preprocessor macros to automatically include .cpp files at end of .h
 - unify include opengl with convenience includes
+- MatrixBase --> PlainObjectBase