add use_lu as parameter to min_quad_dense

Former-commit-id: b0147681f3363b226914ab220c6c957c66cef045

Makefile

@@ -4,6 +4,11 @@ all: lib examples extras
 # Shared flags etc.
 include Makefile.conf
 
+# optimized default settings
+all: LFLAGS +=
+all: CFLAGS += -O3 -DNDEBUG -j 
+debug: CFLAGS += -g -Wall -Werror
+
 EXTRA_DIRS=
 ifeq ($(IGL_WITH_TETGEN),1)
 	# append tetgen extra dir

Makefile.conf

@@ -5,8 +5,3 @@ IGL_WITH_TETGEN=1
 #############################################################################
 AFLAGS += -arch x86_64 -m64
 CFLAGS += -Wall
-# optimized default settings
-all: LFLAGS +=
-all: CFLAGS += -O3 -DNDEBUG -j 
-debug: CFLAGS += -g -Wall -Werror
-

examples/marching_cubes/example.REMOVED.git-id

@@ -1 +1 @@
-5fb18cdaec3ff1d9fee970dc9d3e8852d5d1cf2e
+b3f3382514052ba833b3b5e749a8b16ee2281c1e

include/igl/min_quad_dense.cpp

@@ -8,6 +8,7 @@ template <typename T>
 IGL_INLINE void igl::min_quad_dense_precompute(
   const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
   const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
+  const bool use_lu_decomposition,
   Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S)
 {
   typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Mat;
@@ -28,39 +29,42 @@ IGL_INLINE void igl::min_quad_dense_precompute(
   LM.block(n, 0, m, n) = Aeq;
   LM.block(n, n, m, m).setZero();
 
-//#define USE_LU_DECOMPOSITION
-#ifdef USE_LU_DECOMPOSITION    
-  // if LM is close to singular, use at your own risk :)
-  Mat LMpinv = LM.inverse();
-#else // use SVD
-  typedef Eigen::Matrix<T, Eigen::Dynamic, 1> Vec; 
-  Vec singValues;
-  Eigen::JacobiSVD<Mat> svd;
-  svd.compute(LM, Eigen::ComputeFullU | Eigen::ComputeFullV );
-  const Mat& u = svd.matrixU();
-  const Mat& v = svd.matrixV();
-  const Vec& singVals = svd.singularValues();
-
-  Vec pi_singVals(n + m);
-  int zeroed = 0;
-  for (int i=0; i<n + m; i++)
+  Mat LMpinv;
+  if(use_lu_decomposition)
+  {
+    // if LM is close to singular, use at your own risk :)
+    LMpinv = LM.inverse();
+  }else
   {
-    T sv = singVals(i, 0);
-    assert(sv >= 0);      
-               // printf("sv: %lg ? %lg\n",(double) sv,(double)treshold);
-    if (sv > treshold) pi_singVals(i, 0) = T(1) / sv;
-    else 
+    // use SVD
+    typedef Eigen::Matrix<T, Eigen::Dynamic, 1> Vec; 
+    Vec singValues;
+    Eigen::JacobiSVD<Mat> svd;
+    svd.compute(LM, Eigen::ComputeFullU | Eigen::ComputeFullV );
+    const Mat& u = svd.matrixU();
+    const Mat& v = svd.matrixV();
+    const Vec& singVals = svd.singularValues();
+
+    Vec pi_singVals(n + m);
+    int zeroed = 0;
+    for (int i=0; i<n + m; i++)
     {
-      pi_singVals(i, 0) = T(0);
-      zeroed++;
+      T sv = singVals(i, 0);
+      assert(sv >= 0);      
+                 // printf("sv: %lg ? %lg\n",(double) sv,(double)treshold);
+      if (sv > treshold) pi_singVals(i, 0) = T(1) / sv;
+      else 
+      {
+        pi_singVals(i, 0) = T(0);
+        zeroed++;
+      }
     }
-  }
 
-  printf("min_quad_dense_precompute: %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
-  Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);
+    printf("min_quad_dense_precompute: %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
+    Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);
 
-  Mat LMpinv = v * pi_diag * u.transpose();
-#endif
+    LMpinv = v * pi_diag * u.transpose();
+  }
   S = LMpinv.block(0, 0, n, n + m);
 
   //// debug:
@@ -81,5 +85,5 @@ IGL_INLINE void igl::min_quad_dense_precompute(
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-  template void igl::min_quad_dense_precompute<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::min_quad_dense_precompute<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, bool, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/min_quad_dense.h

@@ -21,6 +21,7 @@ namespace igl
   //   B  n by 1 column of linear coefficients
   //   Aeq  m by n list of linear equality constraint coefficients
   //   Beq  m by 1 list of linear equality constraint constant values
+  //   use_lu_decomposition  use lu rather than SVD
   // Outputs:
   //   S  n by (n + m) "solve" matrix, such that S*[B', Beq'] is a solution
   // Returns true on success, false on error
@@ -28,6 +29,7 @@ namespace igl
   IGL_INLINE void min_quad_dense_precompute(
     const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A,
     const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& Aeq,    
+    const bool use_lu_decomposition,
     Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S);
 }
 

include/igl/verbose.h

@@ -20,7 +20,7 @@ namespace igl
 
 #include <string>
 // http://channel9.msdn.com/forums/techoff/254707-wrapping-printf-in-c/
-inline  int igl::verbose(const char * msg,...)
+inline int igl::verbose(const char * msg,...)
 {
 #ifdef VERBOSE
   va_list argList;