harmonic templates for python

include/igl/harmonic.cpp

@@ -77,8 +77,8 @@ template <
   typename Derivedbc,
   typename DerivedW>
 IGL_INLINE bool igl::harmonic(
-  const Eigen::SparseMatrix<DerivedL> & L,
-  const Eigen::SparseMatrix<DerivedM> & M,
+  const Eigen::SparseCompressedBase<DerivedL> & L,
+  const Eigen::SparseCompressedBase<DerivedM> & M,
   const Eigen::MatrixBase<Derivedb> & b,
   const Eigen::MatrixBase<Derivedbc> & bc,
   const int k,
@@ -89,15 +89,16 @@ IGL_INLINE bool igl::harmonic(
   assert((k==1 || n == M.cols() ) && "M must be same size as L");
   assert((k==1 || n == M.rows() ) && "M must be square");
   assert((k==1 || igl::isdiag(M))  && "Mass matrix should be diagonal");
+  typedef typename DerivedL::Scalar Scalar;
 
-  Eigen::SparseMatrix<DerivedL> Q;
+  Eigen::SparseMatrix<Scalar> Q;
   igl::harmonic(L,M,k,Q);
 
-  typedef DerivedL Scalar;
+
   min_quad_with_fixed_data<Scalar> data;
   min_quad_with_fixed_precompute(Q,b,Eigen::SparseMatrix<Scalar>(),true,data);
   W.resize(n,bc.cols());
-  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  typedef Eigen::Matrix<typename Derivedbc::Scalar,Eigen::Dynamic,1> VectorXS;
   const VectorXS B = VectorXS::Zero(n,1);
   for(int w = 0;w<bc.cols();w++)
   {
@@ -117,17 +118,17 @@ template <
   typename DerivedM,
   typename DerivedQ>
 IGL_INLINE void igl::harmonic(
-  const Eigen::SparseMatrix<DerivedL> & L,
-  const Eigen::SparseMatrix<DerivedM> & M,
+  const Eigen::SparseCompressedBase<DerivedL> & L,
+  const Eigen::SparseCompressedBase<DerivedM> & M,
   const int k,
-  Eigen::SparseMatrix<DerivedQ> & Q)
+  DerivedQ & Q)
 {
   assert(L.rows() == L.cols()&&"L should be square");
   Q = -L;
   if(k == 1) return;
   assert(L.rows() == M.rows()&&"L should match M's dimensions");
   assert(M.rows() == M.cols()&&"M should be square");
-  Eigen::SparseMatrix<DerivedM> Mi;
+  Eigen::SparseMatrix<typename DerivedM::Scalar> Mi;
   invert_diag(M,Mi);
   // This is **not** robust for k>2. See KKT system in [Jacobson et al. 2010]
   // of the kharmonic function in gptoolbox
@@ -146,9 +147,9 @@ IGL_INLINE void igl::harmonic(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   const int k,
-  Eigen::SparseMatrix<DerivedQ> & Q)
+  DerivedQ & Q)
 {
-  Eigen::SparseMatrix<DerivedQ> L,M;
+  DerivedQ L,M;
   cotmatrix(V,F,L);
   if(k>1)
   {
@@ -166,7 +167,7 @@ template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Mat
 // generated by autoexplicit.sh
 template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
-template void igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 // generated by autoexplicit.sh
 template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh

include/igl/harmonic.h

@@ -75,8 +75,8 @@ namespace igl
     typename Derivedbc,
     typename DerivedW>
   IGL_INLINE bool harmonic(
-    const Eigen::SparseMatrix<DerivedL> & L,
-    const Eigen::SparseMatrix<DerivedM> & M,
+    const Eigen::SparseCompressedBase<DerivedL> & L,
+    const Eigen::SparseCompressedBase<DerivedM> & M,
     const Eigen::MatrixBase<Derivedb> & b,
     const Eigen::MatrixBase<Derivedbc> & bc,
     const int k,
@@ -95,10 +95,10 @@ namespace igl
     typename DerivedM,
     typename DerivedQ>
   IGL_INLINE void harmonic(
-    const Eigen::SparseMatrix<DerivedL> & L,
-    const Eigen::SparseMatrix<DerivedM> & M,
+    const Eigen::SparseCompressedBase<DerivedL> & L,
+    const Eigen::SparseCompressedBase<DerivedM> & M,
     const int k,
-    Eigen::SparseMatrix<DerivedQ> & Q);
+    DerivedQ & Q);
   // Inputs:
   //   V  #V by dim vertex positions
   //   F  #F by simplex-size list of element indices
@@ -113,7 +113,7 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F,
     const int k,
-    Eigen::SparseMatrix<DerivedQ> & Q);
+    DerivedQ & Q);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/invert_diag.cpp

@@ -7,18 +7,22 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "invert_diag.h"
 
-template <typename T>
+template <typename DerivedX, typename MatY>
 IGL_INLINE void igl::invert_diag(
-  const Eigen::SparseMatrix<T>& X, 
-  Eigen::SparseMatrix<T>& Y)
+  const Eigen::SparseCompressedBase<DerivedX>& X,
+  MatY& Y)
 {
+  typedef typename DerivedX::Scalar Scalar;
 #ifndef NDEBUG
-  typename Eigen::SparseVector<T> dX = X.diagonal().sparseView();
+  Eigen::SparseMatrix<Scalar> tmp = X;
+  Eigen::SparseVector<Scalar> dX = tmp.diagonal().sparseView();
   // Check that there are no zeros along the diagonal
   assert(dX.nonZeros() == dX.size());
 #endif
   // http://www.alecjacobson.com/weblog/?p=2552
-  if(&Y != &X)
+
+
+  if((void *)&Y != (void *)&X)
   {
     Y = X;
   }
@@ -26,13 +30,13 @@ IGL_INLINE void igl::invert_diag(
   for(int k=0; k<Y.outerSize(); ++k)
   {
     // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (Y,k); it; ++it)
+    for(typename MatY::InnerIterator it (Y,k); it; ++it)
     {
       if(it.col() == it.row())
       {
-        T v = it.value();
+        Scalar v = it.value();
         assert(v != 0);
-        v = ((T)1.0)/v;
+        v = ((Scalar)1.0)/v;
         Y.coeffRef(it.row(),it.col()) = v;
       }
     }
@@ -41,6 +45,6 @@ IGL_INLINE void igl::invert_diag(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::invert_diag<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&);
-template void igl::invert_diag<float>(Eigen::SparseMatrix<float, 0, int> const&, Eigen::SparseMatrix<float, 0, int>&);
+template void igl::invert_diag<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseCompressedBase<Eigen::SparseMatrix<double, 0, int>> const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::invert_diag<Eigen::SparseMatrix<float, 0, int>, Eigen::SparseMatrix<float, 0, int> >(Eigen::SparseCompressedBase<Eigen::SparseMatrix<float, 0, int>> const&, Eigen::SparseMatrix<float, 0, int>&);
 #endif

include/igl/invert_diag.h

@@ -22,10 +22,10 @@ namespace igl
   //   X  an m by n sparse matrix
   // Outputs:
   //   Y  an m by n sparse matrix
-  template <typename T>
+  template <typename DerivedX, typename MatY>
   IGL_INLINE void invert_diag(
-    const Eigen::SparseMatrix<T>& X, 
-    Eigen::SparseMatrix<T>& Y);
+    const Eigen::SparseCompressedBase<DerivedX>& X,
+    MatY& Y);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/isdiag.cpp

@@ -7,14 +7,14 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "isdiag.h"
 
-template <typename Scalar>
-IGL_INLINE bool igl::isdiag(const Eigen::SparseMatrix<Scalar> & A)
+template <typename Derived>
+IGL_INLINE bool igl::isdiag(const Eigen::SparseCompressedBase<Derived> & A)
 {
   // Iterate over outside of A
   for(int k=0; k<A.outerSize(); ++k)
   {
     // Iterate over inside
-    for(typename Eigen::SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+    for(typename Derived::InnerIterator it (A,k); it; ++it)
     {
       if(it.row() != it.col() && it.value()!=0)
       {
@@ -28,5 +28,5 @@ IGL_INLINE bool igl::isdiag(const Eigen::SparseMatrix<Scalar> & A)
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template bool igl::isdiag<double>(class Eigen::SparseMatrix<double,0,int> const &);
+template bool igl::isdiag<Eigen::SparseMatrix<double,0,int>>(Eigen::SparseCompressedBase<Eigen::SparseMatrix<double,0,int>> const &);
 #endif

include/igl/isdiag.h

@@ -17,8 +17,8 @@ namespace igl
   // Inputs:
   //   A  m by n sparse matrix
   // Returns true iff and only if the matrix is diagonal.
-  template <typename Scalar>
-  IGL_INLINE bool isdiag(const Eigen::SparseMatrix<Scalar> & A);
+  template <typename Derived>
+  IGL_INLINE bool isdiag(const Eigen::SparseCompressedBase<Derived> & A);
 };
 #ifndef IGL_STATIC_LIBRARY
 #  include "isdiag.cpp"

include/igl/min_quad_with_fixed.cpp

@@ -51,6 +51,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     assert(n == Aeq.cols() && "#Aeq.cols() should match A.rows()");
   }
 
+  assert(known.cols() == 1 && "known should be a vector");
   assert(A.rows() == n && "A should be square");
   assert(A.cols() == n && "A should be square");
 
@@ -70,7 +71,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   unknown_mask.resize(n,true);
   for(int i = 0;i<kr;i++)
   {
-    unknown_mask[known(i)] = false;
+    unknown_mask[known(i, 0)] = false;
   }
   int u = 0;
   for(int i = 0;i<n;i++)