mv min_quad_dense --> kkt_inverse

include/igl/arap_dof.cpp

@@ -27,7 +27,7 @@
 
 #include "get_seconds_hires.h"
 //#include "MKLEigenInterface.h"
-#include "min_quad_dense.h"
+#include "kkt_inverse.h"
 #include "get_seconds.h"
 #include "columnize.h"
 
@@ -543,7 +543,7 @@ IGL_INLINE bool igl::arap_dof_recomputation(
 #endif
 
   // Compute dense solve matrix (alternative of matrix factorization)
-  //printf("min_quad_dense_precompute()\n");
+  //printf("kkt_inverse()\n");
   MatrixXd Qfull(*Q);
   MatrixXd A_eqfull(A_eq);
   MatrixXd M_Solve;
@@ -552,7 +552,7 @@ IGL_INLINE bool igl::arap_dof_recomputation(
   bool use_lu = data.effective_dim != 2;
   //use_lu = false;
   //printf("use_lu: %s\n",(use_lu?"TRUE":"FALSE"));
-  min_quad_dense_precompute(Qfull, A_eqfull, use_lu,M_Solve);
+  kkt_inverse(Qfull, A_eqfull, use_lu,M_Solve);
   double timer0_end = get_seconds_hires();
   verbose("Bob timing: %.20f\n", (timer0_end - timer0_start)*1000.0);
 

include/igl/kkt_inverse.cpp

@@ -0,0 +1,97 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <[email protected]>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "kkt_inverse.h"
+
+#include <Eigen/Core>
+#include <Eigen/LU>
+#include "EPS.h"
+#include <cstdio>
+
+template <typename T>
+IGL_INLINE void igl::kkt_inverse(
+  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;
+        // This threshold seems to matter a lot but I'm not sure how to
+        // set it
+  const T treshold = igl::FLOAT_EPS;
+  //const T treshold = igl::DOUBLE_EPS;
+
+  const int n = A.rows();
+  assert(A.cols() == n);
+  const int m = Aeq.rows();
+  assert(Aeq.cols() == n);
+
+  // Lagrange multipliers method:
+  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> LM(n + m, n + m);
+  LM.block(0, 0, n, n) = A;
+  LM.block(0, n, n, m) = Aeq.transpose();
+  LM.block(n, 0, m, n) = Aeq;
+  LM.block(n, n, m, m).setZero();
+
+  Mat LMpinv;
+  if(use_lu_decomposition)
+  {
+    // if LM is close to singular, use at your own risk :)
+    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++)
+    {
+      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("kkt_inverse : %i singular values zeroed (threshold = %e)\n", zeroed, treshold);
+    Eigen::DiagonalMatrix<T, Eigen::Dynamic> pi_diag(pi_singVals);
+
+    LMpinv = v * pi_diag * u.transpose();
+  }
+  S = LMpinv.block(0, 0, n, n + m);
+
+  //// debug:
+  //mlinit(&g_pEngine);
+  //
+  //mlsetmatrix(&g_pEngine, "A", A);
+  //mlsetmatrix(&g_pEngine, "Aeq", Aeq);
+  //mlsetmatrix(&g_pEngine, "LM", LM);
+  //mlsetmatrix(&g_pEngine, "u", u);
+  //mlsetmatrix(&g_pEngine, "v", v);
+  //MatrixXd svMat = singVals;
+  //mlsetmatrix(&g_pEngine, "singVals", svMat);
+  //mlsetmatrix(&g_pEngine, "LMpinv", LMpinv);
+  //mlsetmatrix(&g_pEngine, "S", S);
+
+  //int hu = 1;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::kkt_inverse<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/kkt_inverse.h

@@ -0,0 +1,65 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <[email protected]>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_KKT_INVERSE_H
+#define IGL_KKT_INVERSE_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+//// debug
+//#include <matlabinterface.h>
+//Engine *g_pEngine;
+
+
+namespace igl
+{
+  // Systems of the form:
+  //  
+  //  / A   Aeqᵀ \  / x \ = / b   \
+  //  \ Aeq    0 /  \ λ /   \ beq /
+  // \_____.______/\__.__/ \___.___/
+  //       M          z        c
+  //
+  // Arise, for example, when solve convex, linear equality constrained
+  // quadratic minimization problems:
+  //
+  // min ½ xᵀ A x - xᵀb  subject to Aeq x = beq
+  //
+  // This function constructs a matrix S such that x = S c solves the system
+  // above. That is:
+  //
+  //  S = [In 0] M⁻¹
+  //
+  //  so that 
+  //
+  //  x = S c
+  //
+  // Templates:
+  //   T  should be a eigen matrix primitive type like float or double
+  // Inputs:
+  //   A  n by n matrix of quadratic coefficients
+  //   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
+  template <typename T>
+  IGL_INLINE void kkt_inverse(
+    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);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "kkt_inverse.cpp"
+#endif
+
+#endif