min_quad_with_fixed is now truly sparse (still at least 4x slower than matlab)

Former-commit-id: 0162bbfbad7e648760dbc16f6752b08566463e88

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-0.2.2
+0.2.3

include/igl/min_quad_with_fixed.cpp

@@ -7,7 +7,9 @@
 #include "repmat.h"
 #include "lu_lagrange.h"
 #include "full.h"
+#include "matlab_format.h"
 #include "EPS.h"
+#include "cat.h"
 
 //#include <Eigen/SparseExtra>
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
@@ -28,6 +30,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 {
   using namespace Eigen;
   using namespace std;
+  using namespace igl;
   // number of rows
   int n = A.rows();
   // cache problem size
@@ -82,113 +85,143 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   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,igl::EPS<double>());
   // Positive definiteness is *not* determined, rather it is given as a
   // parameter
   data.Auu_pd = pd;
+  if(data.Auu_pd)
+  {
+    // PD implies symmetric
+    data.Auu_sym = true;
+    assert(igl::is_symmetric(Auu,igl::EPS<double>()));
+  }else
+  {
+    // determine if A(unknown,unknown) is symmetric and/or positive definite
+    data.Auu_sym = igl::is_symmetric(Auu,igl::EPS<double>());
+  }
 
   // Append lagrange multiplier quadratic terms
   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);
-  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);
-  }
-  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;
-  igl::sparse(new_AI,new_AJ,new_AV,n+neq,n+neq,new_A);
+  SparseMatrix<T> AeqT = Aeq.transpose();
+  SparseMatrix<T> Z(neq,neq);
+  // This is a bit slower. But why isn't cat fast?
+  new_A = 
+    cat(1,
+      cat(2,   A, AeqT ),
+      cat(2, Aeq,    Z ));
+  //cout<<matlab_format(new_A,"new_A")<<endl;
+  //Matrix<int,Dynamic,1> AI;
+  //Matrix<int,Dynamic,1> AJ;
+  //Matrix<T,Dynamic,1> AV;
+  //// Slow
+  //igl::find(A,AI,AJ,AV);
+  //Matrix<int,Dynamic,1> AeqI(0,1);
+  //Matrix<int,Dynamic,1> AeqJ(0,1);
+  //Matrix<T,Dynamic,1> AeqV(0,1);
+  //if(neq > 0)
+  //{
+  //  // Slow
+  //  igl::find(Aeq,AeqI,AeqJ,AeqV);
+  //}
+  //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;
+  //// Slow
+  //igl::sparse(new_AI,new_AJ,new_AV,n+neq,n+neq,new_A);
+  //cout<<matlab_format(new_A,"new_A")<<endl;
 
   // precompute RHS builders
   if(kr > 0)
   {
-    SparseMatrix<T> Aulk;
+    SparseMatrix<T> Aulk,Akul;
+    // Slow
     igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
-    SparseMatrix<T> Akul;
-    igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
 
     //// This doesn't work!!!
     //data.preY = Aulk + Akul.transpose();
-    SparseMatrix<T> AkulT = Akul.transpose();
-    data.preY = Aulk + AkulT;
+    // Slow
+    if(data.Auu_sym)
+    {
+      data.preY = Aulk*2;
+    }else
+    {
+      igl::slice(new_A,data.known,data.unknown_lagrange,Akul);
+      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;
-  //  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)
+  // Positive definite and no equality constraints (Postive definiteness
+  // implies symmetric)
+  if(data.Auu_pd && neq == 0)
   {
-    SimplicialLDLT<SparseMatrix<T> > ldlt_solver;
-    ldlt_solver.compute(NA);
-    if(ldlt_solver.info() != Eigen::Success)
+    data.llt.compute(Auu);
+    switch(data.ldlt.info())
     {
-      return false;
+      case Eigen::Success:
+        break;
+      case Eigen::NumericalIssue:
+        cerr<<"Error: Numerical issue."<<endl;
+        return false;
+      default:
+        cerr<<"Error: Other."<<endl;
+        return false;
     }
-    // Implementation incomplete
-    assert(false);
-    cerr<<"ERROR min_quad_with_fixed_precompute implementation incomplete"<<endl;
-    return false;
+    data.solver_type = igl::min_quad_with_fixed_data<T>::LLT; 
   }else
   {
-    // Build LU decomposition of NA
-    data.sparse = false;
-    fprintf(stderr,
-      "Warning: min_quad_with_fixed_precompute() resorting to dense LU\n");
-    Matrix<T,Dynamic,Dynamic> NAfull;
-    igl::full(NA,NAfull);
-    data.lu = 
-      FullPivLU<Matrix<T,Dynamic,Dynamic> >(NAfull);
-    if(!data.lu.isInvertible())
+    // Either not PD or there are equality constraints
+    SparseMatrix<T> NA;
+    igl::slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
+    data.NA = NA;
+    // Ideally we'd use LDLT but Eigen doesn't support positive semi-definite
+    // matrices:
+    // http://forum.kde.org/viewtopic.php?f=74&t=106962&p=291990#p291990
+    if(data.Auu_sym && false)
     {
-      fprintf(stderr,
-        "Error: min_quad_with_fixed_precompute() LU not invertible\n");
-      return false;
+      data.ldlt.compute(NA);
+      switch(data.ldlt.info())
+      {
+        case Eigen::Success:
+          break;
+        case Eigen::NumericalIssue:
+          cerr<<"Error: Numerical issue."<<endl;
+          return false;
+        default:
+          cerr<<"Error: Other."<<endl;
+          return false;
+      }
+      data.solver_type = igl::min_quad_with_fixed_data<T>::LDLT; 
+    }else
+    {
+      // Resort to LU
+      // Bottleneck >1/2
+      data.lu.compute(NA); 
+      MatrixXd NAf;
+      full(NA,NAf);
+      //std::cout<<"NA=["<<std::endl<<NA<<std::endl<<"];"<<std::endl;
+      switch(data.lu.info())
+      {
+        case Eigen::Success:
+          break;
+        case Eigen::NumericalIssue:
+          cerr<<"Error: Numerical issue."<<endl;
+          return false;
+        case Eigen::InvalidInput:
+          cerr<<"Error: Invalid Input."<<endl;
+          return false;
+        default:
+          cerr<<"Error: Other."<<endl;
+          return false;
+      }
+      data.solver_type = igl::min_quad_with_fixed_data<T>::LU; 
     }
   }
-  //}
 
   return true;
 }
@@ -202,6 +235,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   const Eigen::Matrix<T,Eigen::Dynamic,1> & Beq,
   Eigen::PlainObjectBase<DerivedZ> & Z)
 {
+  using namespace std;
   // number of known rows
   int kr = data.known.size();
   if(kr!=0)
@@ -246,29 +280,34 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 
   //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
 
-  if(data.sparse)
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> sol;
+  //cout<<matlab_format(NB,"NB")<<endl;
+  switch(data.solver_type)
   {
-    //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);
+    case igl::min_quad_with_fixed_data<T>::LLT:
+      sol = data.llt.solve(NB);
+      break;
+    case igl::min_quad_with_fixed_data<T>::LDLT:
+      sol = data.ldlt.solve(NB);
+      break;
+    case igl::min_quad_with_fixed_data<T>::LU:
+      // Not a bottleneck
+      sol = data.lu.solve(NB);
+      break;
+    default:
+      cerr<<"Error: invalid solver type"<<endl;
+      return false;
   }
-  // Now NB contains sol/-0.5
-  NB *= -0.5;
-  // Now NB contains solution
+  //std::cout<<"sol=["<<std::endl<<sol<<std::endl<<"];"<<std::endl;
+  // Now sol contains sol/-0.5
+  sol *= -0.5;
+  // Now sol contains solution
   // Place solution in Z
-  for(int i = 0;i<(NB.rows()-neq);i++)
+  for(int i = 0;i<(sol.rows()-neq);i++)
   {
-    for(int j = 0;j<NB.cols();j++)
+    for(int j = 0;j<sol.cols();j++)
     {
-      Z(data.unknown_lagrange(i),j) = NB(i,j);
+      Z(data.unknown_lagrange(i),j) = sol(i,j);
     }
   }
   return true;

include/igl/min_quad_with_fixed.h

@@ -73,23 +73,29 @@ struct igl::min_quad_with_fixed_data
   // Whether A(unknown,unknown) is symmetric
   bool Auu_sym;
   // Indices of known variables
-  Eigen::Matrix<int,Eigen::Dynamic,1> known;
+  Eigen::VectorXi known;
   // Indices of unknown variables
-  Eigen::Matrix<int,Eigen::Dynamic,1> unknown;
+  Eigen::VectorXi unknown;
   // Indices of lagrange variables
-  Eigen::Matrix<int,Eigen::Dynamic,1> lagrange;
+  Eigen::VectorXi lagrange;
   // Indices of unknown variable followed by Indices of lagrange variables
-  Eigen::Matrix<int,Eigen::Dynamic,1> unknown_lagrange;
+  Eigen::VectorXi unknown_lagrange;
   // 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;
+  enum SolverType
+  {
+    LLT = 0,
+    LDLT = 1,
+    LU = 2,
+    NUM_SOLVER_TYPES = 3
+  } solver_type;
+  // Solvers
+  Eigen::SimplicialLLT <Eigen::SparseMatrix<T > > llt;
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<T > > ldlt;
+  Eigen::SparseLU<Eigen::SparseMatrix<T, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> >   lu;
+  // Debug
+  Eigen::SparseMatrix<T> NA;
+  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
 };
 
 #ifdef IGL_HEADER_ONLY

include/igl/sparse.h

@@ -40,6 +40,7 @@ namespace igl
     const size_t m,
     const size_t n,
     Eigen::SparseMatrix<T>& X);
+  // THIS MAY BE SUPERSEDED BY EIGEN'S .sparseView
   // Convert a full, dense matrix to a sparse one
   //
   // Templates:

todos.txt

@@ -26,3 +26,4 @@
 - use preprocessor macros to automatically include .cpp files at end of .h
 - unify include opengl with convenience includes
 - MatrixBase --> PlainObjectBase
+- min_quad_with_fixed should probably be a class