better assertions, diag is obsolete, min_quad_with_fixed supports linearly dependent Aeq constraints, but Eigen's sparseQR is not giving a sparse Q so this is a huge bottleneck (dense)

Former-commit-id: dd3d02fb6aad5122f92c8de7e72d7fc1de182fda

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+0.3.1  Linearly dependent constraints in min_quad_with_fixed, SparseQR buggy
 0.3.0  Better active set method support
 0.2.3  More explicits, active set method, opengl/anttweakbar guards
 0.2.2  More explicit instanciations, faster sorts and uniques

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.3.0
+0.3.1

include/igl/active_set.cpp

@@ -41,15 +41,15 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   using namespace std;
   SolverStatus ret = SOLVER_STATUS_ERROR;
   const int n = A.rows();
-  assert(n == A.cols());
+  assert(n == A.cols() && "A must be square");
   // Discard const qualifiers
   //if(B.size() == 0)
   //{
   //  B = Eigen::PlainObjectBase<DerivedB>::Zero(n,1);
   //}
-  assert(n == B.rows());
-  assert(B.cols() == 1);
-  assert(Y.cols() == 1);
+  assert(n == B.rows() && "B.rows() must match A.rows()");
+  assert(B.cols() == 1 && "B must be a column vector");
+  assert(Y.cols() == 1 && "Y must be a column vector");
   assert((Aeq.size() == 0 && Beq.size() == 0) || Aeq.cols() == n);
   assert((Aeq.size() == 0 && Beq.size() == 0) || Aeq.rows() == Beq.rows());
   assert((Aeq.size() == 0 && Beq.size() == 0) || Beq.cols() == 1);
@@ -74,18 +74,18 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   {
     ux = p_ux;
   }
-  assert(lx.rows() == n);
-  assert(ux.rows() == n);
-  assert(ux.cols() == 1);
-  assert(lx.cols() == 1);
-  assert((ux.array()-lx.array()).minCoeff() > 0);
+  assert(lx.rows() == n && "lx must have n rows");
+  assert(ux.rows() == n && "ux must have n rows");
+  assert(ux.cols() == 1 && "lx must be a column vector");
+  assert(lx.cols() == 1 && "ux must be a column vector");
+  assert((ux.array()-lx.array()).minCoeff() > 0 && "ux(i) must be > lx(i)");
   if(Z.size() != 0)
   {
     // Initial guess should have correct size
-    assert(Z.rows() == n);
-    assert(Z.cols() == 1);
+    assert(Z.rows() == n && "Z must have n rows");
+    assert(Z.cols() == 1 && "Z must be a column vector");
   }
-  assert(known.cols() == 1);
+  assert(known.cols() == 1 && "known must be a column vector");
   // Number of knowns
   const int nk = known.size();
 
@@ -105,6 +105,8 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   int iter = 0;
   while(true)
   {
+    //cout<<iter<<":"<<endl;
+    //cout<<"  pre"<<endl;
     // FIND BREACHES OF CONSTRAINTS
     int new_as_lx = 0;
     int new_as_ux = 0;
@@ -239,6 +241,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 #endif
     
+    //cout<<"  min_quad_with_fixed_precompute"<<endl;
     if(!min_quad_with_fixed_precompute(A,known_i,Aeq_i,params.Auu_pd,data))
     {
       cerr<<"Error: min_quad_with_fixed precomputation failed."<<endl;
@@ -250,6 +253,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
+    //cout<<"  min_quad_with_fixed_solve"<<endl;
     Eigen::PlainObjectBase<DerivedZ> sol;
     if(!min_quad_with_fixed_solve(data,B,Y_i,Beq_i,Z,sol))
     {
@@ -257,6 +261,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
+    //cout<<"  post"<<endl;
 
     // Compute Lagrange multiplier values for known_i
     // This needs to be adjusted slightly if A is not symmetric

include/igl/diag.cpp

@@ -11,24 +11,26 @@ IGL_INLINE void igl::diag(
   const Eigen::SparseMatrix<T>& X, 
   Eigen::SparseVector<T>& V)
 {
-  // Get size of input
-  int m = X.rows();
-  int n = X.cols();
-  V = Eigen::SparseVector<T>((m>n?n:m));
-  V.reserve(V.size());
+  assert(false && "Just call X.diagonal().sparseView() directly");
+  V = X.diagonal().sparseView();
+  //// Get size of input
+  //int m = X.rows();
+  //int n = X.cols();
+  //V = Eigen::SparseVector<T>((m>n?n:m));
+  //V.reserve(V.size());
 
-  // Iterate over outside
-  for(int k=0; k<X.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-    {
-      if(it.col() == it.row())
-      {
-        V.coeffRef(it.col()) += it.value();
-      }
-    }
-  }
+  //// Iterate over outside
+  //for(int k=0; k<X.outerSize(); ++k)
+  //{
+  //  // Iterate over inside
+  //  for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+  //  {
+  //    if(it.col() == it.row())
+  //    {
+  //      V.coeffRef(it.col()) += it.value();
+  //    }
+  //  }
+  //}
 }
 
 template <typename T,typename DerivedV>
@@ -36,23 +38,25 @@ IGL_INLINE void igl::diag(
   const Eigen::SparseMatrix<T>& X, 
   Eigen::MatrixBase<DerivedV> & V)
 {
-  // Get size of input
-  int m = X.rows();
-  int n = X.cols();
-  V.derived().resize((m>n?n:m),1);
+  assert(false && "Just call X.diagonal() directly");
+  V = X.diagonal();
+  //// Get size of input
+  //int m = X.rows();
+  //int n = X.cols();
+  //V.derived().resize((m>n?n:m),1);
 
-  // Iterate over outside
-  for(int k=0; k<X.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-    {
-      if(it.col() == it.row())
-      {
-        V(it.col()) = it.value();
-      }
-    }
-  }
+  //// Iterate over outside
+  //for(int k=0; k<X.outerSize(); ++k)
+  //{
+  //  // Iterate over inside
+  //  for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+  //  {
+  //    if(it.col() == it.row())
+  //    {
+  //      V(it.col()) = it.value();
+  //    }
+  //  }
+  //}
 }
 
 template <typename T>

include/igl/diag.h

@@ -6,6 +6,13 @@
 
 namespace igl
 {
+  // http://forum.kde.org/viewtopic.php?f=74&t=117476&p=292388#p292388
+  //
+  // This is superceded by 
+  //   VectorXd d = X.diagonal() and 
+  //   SparseVector<double> d = X.diagonal().sparseView()
+  //
+  //
   // Either extracts the main diagonal of a matrix as a vector OR converts a
   // vector into a matrix with vector along the main diagonal. Like matlab's
   // diag function

include/igl/full.cpp

@@ -5,8 +5,7 @@ IGL_INLINE void igl::full(
   const Eigen::SparseMatrix<T> & A,
   Eigen::PlainObjectBase<DerivedB>& B)
 {
-#warning "Obsolete. Just call B = Matrix(A)"
-  assert(false);
+  assert(false && "Obsolete. Just call B = Matrix(A)");
   using namespace Eigen;
   B = PlainObjectBase<DerivedB >::Zero(A.rows(),A.cols());
   // Iterate over outside
@@ -25,8 +24,7 @@ IGL_INLINE void igl::full(
   const Eigen::PlainObjectBase<DerivedA>& A,
   Eigen::PlainObjectBase<DerivedB>& B)
 {
-#warning "Obsolete. Just call B = Matrix(A)"
-  assert(false);
+  assert(false && "Obsolete. Just call B = Matrix(A)");
   B = A;
 }
 

include/igl/gl_type_size.cpp

@@ -17,7 +17,7 @@ IGL_INLINE int igl::gl_type_size(const GLenum type)
       break;
     default:
       // should handle all other GL_[types]
-      assert(false);
+      assert(false && "Implementation incomplete.");
       break;
   }
   return -1;

include/igl/massmatrix.cpp

@@ -105,10 +105,10 @@ IGL_INLINE void igl::massmatrix(
           break;
         }
       case MASSMATRIX_FULL:
-        assert(false);
+        assert(false && "Implementation incomplete");
         break;
       default:
-        assert(false);
+        assert(false && "Unknown Mass matrix type");
     }
 
   }else if(simplex_size == 4)
@@ -137,7 +137,7 @@ IGL_INLINE void igl::massmatrix(
   }else
   {
     // Unsupported simplex size
-    assert(false);
+    assert(false && "Unsupported simplex size");
   }
   sparse(MI,MJ,MV,n,n,M);
 }

include/igl/min_quad_with_fixed.cpp

@@ -31,6 +31,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   using namespace Eigen;
   using namespace std;
   using namespace igl;
+  //cout<<"    pre"<<endl;
   // number of rows
   int n = A.rows();
   // cache problem size
@@ -83,7 +84,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   data.unknown_lagrange << data.unknown, data.lagrange;
 
   SparseMatrix<T> Auu;
-  igl::slice(A,data.unknown,data.unknown,Auu);
+  slice(A,data.unknown,data.unknown,Auu);
 
   // Positive definiteness is *not* determined, rather it is given as a
   // parameter
@@ -92,100 +93,86 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   {
     // PD implies symmetric
     data.Auu_sym = true;
-    assert(igl::is_symmetric(Auu,igl::EPS<double>()));
+    assert(is_symmetric(Auu,EPS<double>()));
   }else
   {
     // determine if A(unknown,unknown) is symmetric and/or positive definite
-    data.Auu_sym = igl::is_symmetric(Auu,igl::EPS<double>());
+    data.Auu_sym = is_symmetric(Auu,EPS<double>());
   }
 
-  // Append lagrange multiplier quadratic terms
-  SparseMatrix<T> 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,Akul;
-    // Slow
-    igl::slice(new_A,data.unknown_lagrange,data.known,Aulk);
-
-    //// This doesn't work!!!
-    //data.preY = Aulk + Akul.transpose();
-    // 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);
-  }
-
-  // Positive definite and no equality constraints (Postive definiteness
-  // implies symmetric)
-  if(data.Auu_pd && neq == 0)
+  //cout<<"    qr"<<endl;
+  // Determine number of linearly independent constraints
+  int nc = 0;
+  if(neq>0)
   {
-    data.llt.compute(Auu);
-    switch(data.llt.info())
+    // QR decomposition to determine row rank in Aequ
+    slice(Aeq,data.unknown,2,data.Aequ);
+    assert(data.Aequ.rows() == neq);
+    assert(data.Aequ.cols() == data.unknown.size());
+    data.AeqTQR.compute(data.Aequ.transpose().eval());
+    switch(data.AeqTQR.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>::LLT; 
+    nc = data.AeqTQR.rank();
+    assert(nc<=neq);
+    data.Aeq_li = nc == neq;
   }else
   {
-    // 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)
+    data.Aeq_li = true;
+  }
+
+  if(data.Aeq_li)
+  {
+    //cout<<"    Aeq_li=true"<<endl;
+    // Append lagrange multiplier quadratic terms
+    SparseMatrix<T> 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 ));
+
+    // precompute RHS builders
+    if(kr > 0)
+    {
+      SparseMatrix<T> Aulk,Akul;
+      // Slow
+      slice(new_A,data.unknown_lagrange,data.known,Aulk);
+      //// This doesn't work!!!
+      //data.preY = Aulk + Akul.transpose();
+      // Slow
+      if(data.Auu_sym)
+      {
+        data.preY = Aulk*2;
+      }else
+      {
+        slice(new_A,data.known,data.unknown_lagrange,Akul);
+        SparseMatrix<T> AkulT = Akul.transpose();
+        data.preY = Aulk + AkulT;
+      }
+    }else
     {
-      data.ldlt.compute(NA);
-      switch(data.ldlt.info())
+      data.preY.resize(data.unknown_lagrange.size(),0);
+    }
+
+    // Positive definite and no equality constraints (Postive definiteness
+    // implies symmetric)
+    //cout<<"    factorize"<<endl;
+    if(data.Auu_pd && neq == 0)
+    {
+      data.llt.compute(Auu);
+      switch(data.llt.info())
       {
         case Eigen::Success:
           break;
@@ -196,40 +183,137 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
           cerr<<"Error: Other."<<endl;
           return false;
       }
-      data.solver_type = igl::min_quad_with_fixed_data<T>::LDLT; 
+      data.solver_type = min_quad_with_fixed_data<T>::LLT;
     }else
     {
-      // Resort to LU
-      // Bottleneck >1/2
-      data.lu.compute(NA); 
-      //std::cout<<"NA=["<<std::endl<<NA<<std::endl<<"];"<<std::endl;
-      switch(data.lu.info())
+      // Either not PD or there are equality constraints
+      SparseMatrix<T> NA;
+      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)
+      {
+        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 = min_quad_with_fixed_data<T>::LDLT;
+      }else
+      {
+        // Resort to LU
+        // Bottleneck >1/2
+        data.lu.compute(NA);
+        //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 = min_quad_with_fixed_data<T>::LU;
+      }
+    }
+  }else
+  {
+    //cout<<"    Aeq_li=false"<<endl;
+    data.neq = neq;
+    const int nu = data.unknown.size();
+    //cout<<"nu: "<<nu<<endl;
+    //cout<<"neq: "<<neq<<endl;
+    //cout<<"nc: "<<nc<<endl;
+    //cout<<"    matrixR"<<endl;
+    SparseMatrix<T> AeqTR,AeqTQ;
+    AeqTR = data.AeqTQR.matrixR();
+    // This shouldn't be necessary
+    AeqTR.prune(0.0);
+    //cout<<"    matrixQ"<<endl;
+    // THIS IS ESSENTIALLY DENSE AND THIS IS BY FAR THE BOTTLENECK
+    // http://forum.kde.org/viewtopic.php?f=74&t=117500
+    AeqTQ = data.AeqTQR.matrixQ();
+    // This shouldn't be necessary
+    AeqTQ.prune(0.0);
+    //cout<<"AeqTQ: "<<AeqTQ.rows()<<" "<<AeqTQ.cols()<<endl;
+    //cout<<matlab_format(AeqTQ,"AeqTQ")<<endl;
+    //cout<<"    perms"<<endl;
+    SparseMatrix<double> I(neq,neq);
+    I.setIdentity();
+    data.AeqTE = data.AeqTQR.colsPermutation() * I;
+    data.AeqTET = data.AeqTQR.colsPermutation().transpose() * I;
+    assert(AeqTR.rows() == neq);
+    assert(AeqTQ.rows() == nu);
+    assert(AeqTQ.cols() == nu);
+    assert(AeqTR.cols() == neq);
+    //cout<<"    slice"<<endl;
+    data.AeqTQ1 = AeqTQ.topLeftCorner(nu,nc);
+    data.AeqTQ1T = data.AeqTQ1.transpose().eval();
+    // ALREADY TRIM (Not 100% sure about this)
+    data.AeqTR1 = AeqTR.topLeftCorner(nc,nc);
+    data.AeqTR1T = data.AeqTR1.transpose().eval();
+    //cout<<"AeqTR1T.size() "<<data.AeqTR1T.rows()<<" "<<data.AeqTR1T.cols()<<endl;
+    // Null space
+    data.AeqTQ2 = AeqTQ.bottomRightCorner(nu,nu-nc);
+    data.AeqTQ2T = data.AeqTQ2.transpose().eval();
+    //cout<<"    proj"<<endl;
+    // Projected hessian
+    SparseMatrix<T> QRAuu = data.AeqTQ2T * Auu * data.AeqTQ2;
+    {
+      //cout<<"    factorize"<<endl;
+      // QRAuu should always be PD
+      data.llt.compute(QRAuu);
+      switch(data.llt.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; 
+      data.solver_type = min_quad_with_fixed_data<T>::QR_LLT;
     }
+    //cout<<"    smash"<<endl;
+    // Known value multiplier
+    SparseMatrix<T> Auk;
+    slice(A,data.unknown,data.known,Auk);
+    SparseMatrix<T> Aku;
+    slice(A,data.known,data.unknown,Aku);
+    SparseMatrix<T> AkuT = Aku.transpose();
+    data.preY = Auk + AkuT;
+    // Needed during solve
+    data.Auu = Auu;
+    slice(Aeq,data.known,2,data.Aeqk);
+    assert(data.Aeqk.rows() == neq);
+    assert(data.Aeqk.cols() == data.known.size());
   }
-
   return true;
 }
 
 
 template <
   typename T,
-  typename DerivedB, 
+  typename DerivedB,
   typename DerivedY,
-  typename DerivedBeq, 
+  typename DerivedBeq,
   typename DerivedZ,
   typename Derivedsol>
 IGL_INLINE bool igl::min_quad_with_fixed_solve(
@@ -241,6 +325,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   Eigen::PlainObjectBase<Derivedsol> & sol)
 {
   using namespace std;
+  using namespace Eigen;
+  using namespace igl;
   // number of known rows
   int kr = data.known.size();
   if(kr!=0)
@@ -252,13 +338,25 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   assert(B.cols() == 1);
   assert(Beq.size() == 0 || Beq.cols() == 1);
 
+  // resize output
+  Z.resize(data.n,cols);
+  // Set known values
+  for(int i = 0;i < kr;i++)
+  {
+    for(int j = 0;j < cols;j++)
+    {
+      Z(data.known(i),j) = Y(i,j);
+    }
+  }
+
+  if(data.Aeq_li)
+  {
+
   // number of lagrange multipliers aka linear equality constraints
   int neq = data.lagrange.size();
-
   // append lagrange multiplier rhs's
   Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
   BBeq << B, (Beq*-2.0);
-
   // Build right hand side
   Eigen::Matrix<T,Eigen::Dynamic,1> BBequl;
   igl::slice(BBeq,data.unknown_lagrange,BBequl);
@@ -273,20 +371,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
     NB = data.preY * Y + BBequlcols;
   }
 
-  // resize output
-  Z.resize(data.n,cols);
-
-  // Set known values
-  for(int i = 0;i < kr;i++)
-  {
-    for(int j = 0;j < cols;j++)
-    {
-      Z(data.known(i),j) = Y(i,j);
-    }
-  }
-
   //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
-
   //cout<<matlab_format(NB,"NB")<<endl;
   switch(data.solver_type)
   {
@@ -316,14 +401,72 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
       Z(data.unknown_lagrange(i),j) = sol(i,j);
     }
   }
+  }else
+  {
+    assert(data.solver_type == min_quad_with_fixed_data<T>::QR_LLT);
+    PlainObjectBase<DerivedBeq> eff_Beq;
+    // Adjust Aeq rhs to include known parts
+    eff_Beq = 
+      //data.AeqTQR.colsPermutation().transpose() * (-data.Aeqk * Y + Beq);
+      data.AeqTET * (-data.Aeqk * Y + Beq);
+    // Where did this -0.5 come from? Probably the same place as above.
+    PlainObjectBase<DerivedB> Bu;
+    slice(B,data.unknown,Bu);
+    PlainObjectBase<DerivedB> NB;
+    NB = -0.5*(Bu + data.preY * Y);
+    // Trim eff_Beq
+    const int nc = data.AeqTQR.rank();
+    const int neq = Beq.rows();
+    eff_Beq = eff_Beq.topLeftCorner(nc,1).eval();
+    data.AeqTR1T.template triangularView<Lower>().solveInPlace(eff_Beq);
+    // Now eff_Beq = (data.AeqTR1T \ (data.AeqTET * (-data.Aeqk * Y + Beq)))
+    PlainObjectBase<DerivedB> lambda_0;
+    lambda_0 = data.AeqTQ1 * eff_Beq;
+    //cout<<matlab_format(lambda_0,"lambda_0")<<endl;
+    PlainObjectBase<DerivedB> QRB;
+    QRB = -data.AeqTQ2T * (data.Auu * lambda_0) + data.AeqTQ2T * NB;
+    PlainObjectBase<Derivedsol> lambda;
+    lambda = data.llt.solve(QRB);
+    // prepare output
+    PlainObjectBase<Derivedsol> solu;
+    solu = data.AeqTQ2 * lambda + lambda_0;
+    //  http://www.math.uh.edu/~rohop/fall_06/Chapter3.pdf
+    PlainObjectBase<Derivedsol> solLambda;
+    {
+      PlainObjectBase<Derivedsol> temp1,temp2;
+      temp1 = (data.AeqTQ1T * NB - data.AeqTQ1T * data.Auu * solu);
+      data.AeqTR1.template triangularView<Upper>().solveInPlace(temp1);
+      //cout<<matlab_format(temp1,"temp1")<<endl;
+      temp2 = PlainObjectBase<Derivedsol>::Zero(neq,1);
+      temp2.topLeftCorner(nc,1) = temp1;
+      //solLambda = data.AeqTQR.colsPermutation() * temp2;
+      solLambda = data.AeqTE * temp2;
+    }
+    // sol is [Z(unknown);Lambda]
+    assert(data.unknown.size() == solu.rows());
+    assert(cols == solu.cols());
+    assert(data.neq == neq);
+    assert(data.neq == solLambda.rows());
+    assert(cols == solLambda.cols());
+    sol.resize(data.unknown.size()+data.neq,cols);
+    sol.block(0,0,solu.rows(),solu.cols()) = solu;
+    sol.block(solu.rows(),0,solLambda.rows(),solLambda.cols()) = solLambda;
+    for(int u = 0;u<data.unknown.size();u++)
+    {
+      for(int j = 0;j<Z.cols();j++)
+      {
+        Z(data.unknown(u),j) = solu(u,j);
+      }
+    }
+  }
   return true;
 }
 
 template <
   typename T,
-  typename DerivedB, 
+  typename DerivedB,
   typename DerivedY,
-  typename DerivedBeq, 
+  typename DerivedBeq,
   typename DerivedZ>
 IGL_INLINE bool igl::min_quad_with_fixed_solve(
   const min_quad_with_fixed_data<T> & data,

include/igl/min_quad_with_fixed.h

@@ -116,12 +116,30 @@ struct igl::min_quad_with_fixed_data
     LLT = 0,
     LDLT = 1,
     LU = 2,
-    NUM_SOLVER_TYPES = 3
+    QR_LLT = 3,
+    NUM_SOLVER_TYPES = 4
   } 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;
+  // QR factorization
+  // Are rows of Aeq linearly independent?
+  bool Aeq_li;
+  // Columns of Aeq corresponding to unknowns
+  int neq;
+  Eigen::SparseQR<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int> >  AeqTQR;
+  Eigen::SparseMatrix<T> Aeqk;
+  Eigen::SparseMatrix<T> Aequ;
+  Eigen::SparseMatrix<T> Auu;
+  Eigen::SparseMatrix<T> AeqTQ1;
+  Eigen::SparseMatrix<T> AeqTQ1T;
+  Eigen::SparseMatrix<T> AeqTQ2;
+  Eigen::SparseMatrix<T> AeqTQ2T;
+  Eigen::SparseMatrix<T> AeqTR1;
+  Eigen::SparseMatrix<T> AeqTR1T;
+  Eigen::SparseMatrix<T> AeqTE;
+  Eigen::SparseMatrix<T> AeqTET;
   // Debug
   Eigen::SparseMatrix<T> NA;
   Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;

include/igl/slice.cpp

@@ -87,7 +87,7 @@ IGL_INLINE void igl::slice(
       igl::colon(0,X.rows()-1,C);
       return slice(X,C,R,Y);
     default:
-      assert(false);
+      assert(false && "Unsupported dimension");
       return;
   }
 }

include/igl/sparse.cpp

@@ -66,8 +66,7 @@ IGL_INLINE void igl::sparse(
   const Eigen::PlainObjectBase<DerivedD>& D,
   Eigen::SparseMatrix<T>& X)
 {
-#warning "This is obsolete. Call .sparseView() instead"
-  assert(false);
+  assert(false && "Obsolete. Just call D.sparseView() directly");
   using namespace std;
   using namespace Eigen;
   vector<Triplet<T> > DIJV;
@@ -91,8 +90,7 @@ template <typename DerivedD>
 IGL_INLINE Eigen::SparseMatrix<typename DerivedD::Scalar > igl::sparse(
   const Eigen::PlainObjectBase<DerivedD>& D)
 {
-#warning "This is obsolete. Call .sparseView() instead"
-  assert(false);
+  assert(false && "Obsolete. Just call D.sparseView() directly");
   Eigen::SparseMatrix<typename DerivedD::Scalar > X;
   igl::sparse(D,X);
   return X;

style.css

@@ -190,6 +190,7 @@ ul
   padding: 0px 20px;
   list-style-type:none;
 }
+li { padding-left: 30px;text-indent: -30px;}
 li { font-size: 100% }
 li li { font-size: 90% }
 li li li { font-size: 80% }

style_guidelines.html

@@ -3,8 +3,9 @@
   <head>
     <title>libigl - Style Guidelines</title>
     <link href="./style.css" rel="stylesheet" type="text/css">
-  <body class=article_body>
-  <div class=article>
+  <body>
+    <div id=container>
+    <div class=article_inner>
     <a href=.><img src=libigl-logo.jpg alt="igl logo" class=center></a>
     <h1>libigl Style Guidelines</h1>
     <p>
@@ -103,7 +104,8 @@ External dependencies can go in the external/ directory
       </li>
       <li>
 Do not use the using namespace directive anywhere outside a local scope. This
-means never write: <code>using namespace std;</code> or <code>using namespace igl;</code> etc.
+means never write: <code>using namespace std;</code> or <code>using namespace
+igl;</code> etc. at the top of a file.
       </li>
       <li>
 Function names should either be the same as the corresponding MATLAB function
@@ -112,6 +114,8 @@ or should use all lowercase separated by underscores: e.g.
       </li>
       <li> Classes should be in <code>CamelCase</code></li>
       <li> No tabs, only two spaces per indentation level </li>
+      <li> Be generous with assertions, but always identify them with strings:
+      <br> e.g. <code>assert(m&lt;n && "m must be less than n");</code></li>
   </ul>
 
     <h2>Specific rules</h2>
@@ -156,8 +160,8 @@ and outputs.
         <th>script</th>
         <th>Description</th>
       <tr class=d0>
-        <td><code>grep -L inline *</code></td>
-        <td>Find files that aren't using "inline"</td>
+        <td><code>grep -L IGL_INLINE *</code></td>
+        <td>Find files that aren't using "IGL_INLINE"</td>
       </tr>
       <tr class=d1>
         <td><code>grep -L "namespace igl" *</code></td>
@@ -167,27 +171,28 @@ and outputs.
         <td><code>grep -P '\t' *</code></td>
         <td>Find files using [TAB] character</td>
       </tr>
-      <tr class=d1>
-        <td><code>grep -L "^\/\/ Implementation" *</code></td>
-        <td>Find files that don't contain // Implementation</td>
-      </tr>
-      <tr class=d0>
-        <td><code>grep -L "^#ifndef IGL_" *</code></td>
-        <td>Find files that don't contain #ifndef IGL_</td>
-      </tr>
       <tr class=d1>
         <td><code>grep -l "^using namespace" *.cpp</code></td>
         <td>Find .cpp files that contain ^using namespace</td>
       </tr>
       <tr class=d0>
+        <td><code>grep -l 'assert([^"]*);' *</code></td>
+        <td>Find files using asserts without identifier strings</td>
+      </tr>
+      <tr class=d1>
         <td><code>grep -l "ifndef IGL_.*[^H] *$" *</code></td>
         <td>Find .h files that contain ifndef IGL_*[^H]</td>
       </tr>
+      <tr class=d0>
+        <td><code>grep -L "^#ifndef IGL_" *</code></td>
+        <td>Find files that don't contain #ifndef IGL_</td>
+      </tr>
     </table>
 
     <p>See also: <a href=tutorial.html>tutorial</a>, <a
     href=doc.html>auto-documentation</a>, <a href=file-formats/index.html>file
     formats</a></p>
   </div>
+  </div>
   </body>
 </html>