quadprog for small dense programs

include/igl/min_quad_with_fixed.cpp

@@ -582,9 +582,325 @@ IGL_INLINE bool igl::min_quad_with_fixed(
   return min_quad_with_fixed_solve(data,B,Y,Beq,Z);
 }
 
+//#include <igl/matlab_format.h>
+#include <iostream>
+#include <type_traits>
+
+template <typename Scalar, int n, int m, bool Hpd>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc,
+  const Eigen::Matrix<Scalar,m,n> & A,
+  const Eigen::Matrix<Scalar,m,1> & b)
+{
+  //printf("igl::min_quad_with_fixed<S,n,m,Hpd>\n");
+  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
+  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
+  //std::cout<<igl::matlab_format(k,"k")<<std::endl;
+  //std::cout<<igl::matlab_format(bc,"bc")<<std::endl;
+  //std::cout<<igl::matlab_format(A,"A")<<std::endl;
+  //std::cout<<igl::matlab_format(b,"b")<<std::endl;
+  const auto dyn_n = n == Eigen::Dynamic ? H.rows() : n;
+  const auto dyn_m = m == Eigen::Dynamic ? A.rows() : m;
+  const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+m;
+  const auto dyn_nn = nn == Eigen::Dynamic ? dyn_n+dyn_m : nn;
+  if(dyn_m == 0)
+  {
+    return igl::min_quad_with_fixed<Scalar,n,Hpd>(H,f,k,bc);
+  }
+  // min_x ½ xᵀ H x + xᵀ f   subject to A x = b and x(k) = bc(k)
+  // let zᵀ = [xᵀ λᵀ]
+  // min_z ½ zᵀ [H Aᵀ;A 0] z + zᵀ [f;-b]   z(k) = bc(k)
+  const auto make_HH = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,nn> HH =
+      Eigen::Matrix<Scalar,nn,nn>::Zero(dyn_nn,dyn_nn);
+    HH.topLeftCorner(dyn_n,dyn_n) = H;
+    HH.bottomLeftCorner(dyn_m,dyn_n) = A;
+    HH.topRightCorner(dyn_n,dyn_m) = A.transpose();
+    return HH;
+  };
+  const Eigen::Matrix<Scalar,nn,nn> HH = make_HH();
+  const auto make_ff  = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,1> ff(dyn_nn);
+    ff.head(dyn_n) =  f;
+    ff.tail(dyn_m) = -b;
+    return ff;
+  };
+  const Eigen::Matrix<Scalar,nn,1> ff = make_ff();
+  const auto make_kk  = [&]()
+  {
+    Eigen::Array<bool,nn,1> kk =
+      Eigen::Array<bool,nn,1>::Constant(dyn_nn,1,false);
+    kk.head(dyn_n) =  k;
+    return kk;
+  };
+  const Eigen::Array<bool,nn,1> kk = make_kk();
+  const auto make_bcbc= [&]()
+  {
+    Eigen::Matrix<Scalar,nn,1> bcbc(dyn_nn);
+    bcbc.head(dyn_n) =  bc;
+    return bcbc;
+  };
+  const Eigen::Matrix<double,nn,1> bcbc = make_bcbc();
+  const Eigen::Matrix<Scalar,nn,1> xx = 
+    min_quad_with_fixed<Scalar,nn,false>(HH,ff,kk,bcbc);
+  //std::cout<<igl::matlab_format(xx.head(dyn_n).eval(),"x")<<std::endl;
+  return xx.head(dyn_n);
+}
+
+template <typename Scalar, int n, bool Hpd>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc)
+{
+  //printf("igl::min_quad_with_fixed<S,n,Hpd>\n");
+  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
+  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
+  //std::cout<<igl::matlab_format(k,"k")<<std::endl;
+  //std::cout<<igl::matlab_format(bc,"bc")<<std::endl;
+  assert(H.isApprox(H.transpose(),1e-7)); 
+  assert(H.rows() == H.cols());
+  assert(H.rows() == f.size());
+  assert(H.rows() == k.size());
+  assert(H.rows() == bc.size());
+  const auto kcount = k.count();
+  // Everything fixed
+  if(kcount == (Eigen::Dynamic?H.rows():n))
+  {
+    //std::cout<<igl::matlab_format(bc,"x")<<std::endl;
+    return bc;
+  }
+  // Nothing fixed
+  if(kcount == 0)
+  {
+    // avoid function call
+    typedef Eigen::Matrix<Scalar,n,n> MatrixSn;
+    typedef typename 
+      std::conditional<Hpd,Eigen::LLT<MatrixSn>,Eigen::CompleteOrthogonalDecomposition<MatrixSn>>::type
+      Solver;
+    //std::cout<<igl::matlab_format( Solver(H).solve(-f).eval() ,"x")<<std::endl;
+    return Solver(H).solve(-f);
+  }
+  // All-but-one fixed
+  if( (Eigen::Dynamic?H.rows():n)-kcount == 1)
+  {
+    // which one is not fixed?
+    int u = -1;
+    for(int i=0;i<k.size();i++){ if(!k(i)){ u=i; break; } } 
+    assert(u>=0);
+    // min ½ x(u) Huu x(u) + x(u)(fu + H(u,k)bc(k))
+    // Huu x(u) = -(fu + H(u,k) bc(k))
+    // x(u) = (-fu + ∑ -Huj bcj)/Huu
+    Eigen::Matrix<Scalar,n,1> x = bc;
+    x(u) = -f(u);
+    for(int i=0;i<k.size();i++){ if(i!=u){ x(u)-=bc(i)*H(i,u); } } 
+    x(u) /= H(u,u);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    return x;
+  }
+  // Is there a smart template way to do this?
+  switch(kcount)
+  {
+    case 0: assert(false && "Handled above."); return Eigen::Matrix<Scalar,n,1>();
+    // maxi=16;for i=1:maxi;fprintf('    case %d:\n    {\n     const bool D = (n-%d<=0)||(%d>=n)||(n>%d);\n     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:%d,Hpd>(H,f,k,bc);\n    }\n',[i i i maxi i]);end
+    case 1:
+    {
+     const bool D = (n-1<=0)||(1>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:1,Hpd>(H,f,k,bc);
+    }
+    case 2:
+    {
+     const bool D = (n-2<=0)||(2>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:2,Hpd>(H,f,k,bc);
+    }
+    case 3:
+    {
+     const bool D = (n-3<=0)||(3>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:3,Hpd>(H,f,k,bc);
+    }
+    case 4:
+    {
+     const bool D = (n-4<=0)||(4>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:4,Hpd>(H,f,k,bc);
+    }
+    case 5:
+    {
+     const bool D = (n-5<=0)||(5>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:5,Hpd>(H,f,k,bc);
+    }
+    case 6:
+    {
+     const bool D = (n-6<=0)||(6>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:6,Hpd>(H,f,k,bc);
+    }
+    case 7:
+    {
+     const bool D = (n-7<=0)||(7>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:7,Hpd>(H,f,k,bc);
+    }
+    case 8:
+    {
+     const bool D = (n-8<=0)||(8>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:8,Hpd>(H,f,k,bc);
+    }
+    case 9:
+    {
+     const bool D = (n-9<=0)||(9>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:9,Hpd>(H,f,k,bc);
+    }
+    case 10:
+    {
+     const bool D = (n-10<=0)||(10>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:10,Hpd>(H,f,k,bc);
+    }
+    case 11:
+    {
+     const bool D = (n-11<=0)||(11>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:11,Hpd>(H,f,k,bc);
+    }
+    case 12:
+    {
+     const bool D = (n-12<=0)||(12>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:12,Hpd>(H,f,k,bc);
+    }
+    case 13:
+    {
+     const bool D = (n-13<=0)||(13>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:13,Hpd>(H,f,k,bc);
+    }
+    case 14:
+    {
+     const bool D = (n-14<=0)||(14>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:14,Hpd>(H,f,k,bc);
+    }
+    case 15:
+    {
+     const bool D = (n-15<=0)||(15>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:15,Hpd>(H,f,k,bc);
+    }
+    case 16:
+    {
+     const bool D = (n-16<=0)||(16>=n)||(n>16);
+     return min_quad_with_fixed<Scalar,D?Eigen::Dynamic:n,D?Eigen::Dynamic:16,Hpd>(H,f,k,bc);
+    }
+    default:
+      return min_quad_with_fixed<Scalar,Eigen::Dynamic,Eigen::Dynamic,Hpd>(H,f,k,bc);
+  }
+}
+
+IGL_INLINE template <typename Scalar, int n, int kcount, bool Hpd>
+Eigen::Matrix<Scalar,n,1> igl::min_quad_with_fixed(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Array<bool,n,1> & k,
+  const Eigen::Matrix<Scalar,n,1> & bc)
+{
+  //printf("igl::min_quad_with_fixed<S,n,kcount,Hpd>\n");
+  // 0 and n should be handle outside this function
+  static_assert(kcount==Eigen::Dynamic || kcount>0                  ,"");
+  static_assert(kcount==Eigen::Dynamic || kcount<n                  ,"");
+  const int ucount = n==Eigen::Dynamic ? Eigen::Dynamic : n-kcount;
+  static_assert(kcount==Eigen::Dynamic || ucount+kcount == n        ,"");
+  static_assert((n==Eigen::Dynamic) == (ucount==Eigen::Dynamic),"");
+  static_assert((kcount==Eigen::Dynamic) == (ucount==Eigen::Dynamic),"");
+  assert((n==Eigen::Dynamic) || n == H.rows());
+  assert((kcount==Eigen::Dynamic) || kcount == k.count());
+  typedef Eigen::Matrix<Scalar,ucount,ucount> MatrixSuu;
+  typedef Eigen::Matrix<Scalar,ucount,kcount> MatrixSuk;
+  typedef Eigen::Matrix<Scalar,n,1>      VectorSn;
+  typedef Eigen::Matrix<Scalar,ucount,1> VectorSu;
+  typedef Eigen::Matrix<Scalar,kcount,1> VectorSk;
+  const auto dyn_n = n==Eigen::Dynamic ? H.rows() : n;
+  const auto dyn_kcount = kcount==Eigen::Dynamic ? k.count() : kcount;
+  const auto dyn_ucount = ucount==Eigen::Dynamic ? dyn_n- dyn_kcount : ucount;
+  // For ucount==2 or kcount==2 this calls the coefficient initiliazer rather
+  // than the size initilizer, but I guess that's ok.
+  MatrixSuu Huu(dyn_ucount,dyn_ucount);
+  MatrixSuk Huk(dyn_ucount,dyn_kcount);
+  VectorSu mrhs(dyn_ucount);
+  VectorSk  bck(dyn_kcount);
+  {
+    int ui = 0;
+    int ki = 0;
+    for(int i = 0;i<dyn_n;i++)
+    {
+      if(k(i))
+      {
+        bck(ki) = bc(i);
+        ki++;
+      }else
+      {
+        mrhs(ui) = f(i);
+        int uj = 0;
+        int kj = 0;
+        for(int j = 0;j<dyn_n;j++)
+        {
+          if(k(j))
+          {
+            Huk(ui,kj) = H(i,j);
+            kj++;
+          }else
+          {
+            Huu(ui,uj) = H(i,j);
+            uj++;
+          }
+        }
+        ui++;
+      }
+    }
+  }
+  mrhs += Huk * bck;
+  typedef typename 
+    std::conditional<Hpd,
+      Eigen::LLT<MatrixSuu>,
+      // LDLT should be faster for indefinite problems but already found some
+      // cases where it was too inaccurate when called via quadprog_primal.
+      // Ideally this function takes LLT,LDLT, or
+      // CompleteOrthogonalDecomposition as a template parameter.  "template
+      // template" parameters did work because LLT,LDLT have different number of
+      // template parameters from CompleteOrthogonalDecomposition.  Perhaps
+      // there's a way to take advantage of LLT and LDLT's default template
+      // parameters (I couldn't figure out how).
+      Eigen::CompleteOrthogonalDecomposition<MatrixSuu>>::type
+    Solver;
+  VectorSu xu = Solver(Huu).solve(-mrhs);
+  VectorSn x(dyn_n);
+  {
+    int ui = 0;
+    int ki = 0;
+    for(int i = 0;i<dyn_n;i++)
+    {
+      if(k(i))
+      {
+        x(i) = bck(ki);
+        ki++;
+      }else
+      {
+        x(i) = xu(ui);
+        ui++;
+      }
+    }
+  }
+  //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+  return x;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> igl::min_quad_with_fixed<double, -1, -1, true>(Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Array<bool, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> igl::min_quad_with_fixed<double, 8, 2, true>(Eigen::Matrix<double, 8, 8, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)1) : ((((8) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 8> const&, Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Array<bool, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Matrix<double, 8, 1, ((Eigen::StorageOptions)0) | ((((8) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 8, 1> const&, Eigen::Matrix<double, 2, 8, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((8) != (1))) ? ((Eigen::StorageOptions)1) : ((((8) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 8> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::min_quad_with_fixed<double, 3, 1, true>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Array<bool, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 1, 3, ((Eigen::StorageOptions)0) | ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 1, 3> const&, Eigen::Matrix<double, 1, 1, ((Eigen::StorageOptions)0) | ((((1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 1, 1> const&);
+template Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> igl::min_quad_with_fixed<double, 2, 0, true>(Eigen::Matrix<double, 2, 2, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 2> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Array<bool, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 0, 2, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 2> const&, Eigen::Matrix<double, 0, 1, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 1> const&);
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::min_quad_with_fixed<double, 3, 0, true>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Array<bool, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 0, 3, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 3> const&, Eigen::Matrix<double, 0, 1, ((Eigen::StorageOptions)0) | ((((0) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((0) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 0, 1> const&);
 template bool igl::min_quad_with_fixed<double, 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::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, 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<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #if EIGEN_VERSION_AT_LEAST(3,3,0)
 #else

include/igl/min_quad_with_fixed.h

@@ -117,6 +117,51 @@ namespace igl
     const Eigen::MatrixBase<DerivedBeq> & Beq,
     const bool pd,
     Eigen::PlainObjectBase<DerivedZ> & Z);
+
+  // Dense version optimized for very small, known at compile time sizes. Still
+  // works for Eigen::Dynamic (and then everything needs to be Dynamic).
+  //
+  // min_x ½ xᵀ H x + xᵀ f  
+  // subject to 
+  //   A x = b
+  //   x(i) = bc(i) iff k(i)==true
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   m  #A or Eigen::Dynamic if not known at compile time
+  //   Hpd  whether H is positive definite (LLT used) or not (QR used)
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   k  #H list of flags whether to fix value
+  //   bc  #H value to fix to (if !k(i) then bc(i) is ignored)
+  //   A  #A by #H list of linear equality constraint coefficients
+  //   b  #A list of lienar equality right-hand sides
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int m, bool Hpd=true>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc,
+    const Eigen::Matrix<Scalar,m,n> & A,
+    const Eigen::Matrix<Scalar,m,1> & b);
+  template <typename Scalar, int n, bool Hpd=true>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc);
+  // Special wrapper where the number of constrained values (i.e., true values
+  // in k) is exposed as a template parameter. Not intended to be called
+  // directly. The overhead of calling the overloads above is already minimal.
+  template <typename Scalar, int n, int kcount, bool Hpd/*no default*/>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> min_quad_with_fixed(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Array<bool,n,1> & k,
+    const Eigen::Matrix<Scalar,n,1> & bc);
 }
 
 template <typename T>

include/igl/quadprog.cpp

@@ -0,0 +1,202 @@
+#include "quadprog.h"
+#include "min_quad_with_fixed.h"
+//#include <igl/matlab_format.h>
+//#include <igl/eigen_format.h>
+#include <iostream>
+
+template <typename Scalar, int n, int ni>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Matrix<Scalar,ni,n> & Ai,
+  const Eigen::Matrix<Scalar,ni,1> & lbi,
+  const Eigen::Matrix<Scalar,ni,1> & ubi,
+  const Eigen::Matrix<Scalar,n,1> & lb,
+  const Eigen::Matrix<Scalar,n,1> & ub)
+{
+  const auto dyn_n = n == Eigen::Dynamic ? H.rows() : n;
+  const auto dyn_ni = ni == Eigen::Dynamic ? Ai.rows() : ni;
+  // min_x ½ xᵀ H x + xᵀ f   
+  // subject to lbi ≤ Ai x ≤ ubi, lb≤x≤ub
+  //
+  // min_x,xi ½ xᵀ H x + xᵀ f
+  // subject to Ai x - xi = 0, lbi≤xi≤ubi, lb≤x≤ub 
+  //
+  // min_z ½ zᵀ [H 0;0 0] z + zᵀ [f;0]
+  // subject to [Ai -I] z = 0, [lb;lbi]≤z≤[ub;ubi]
+  const int nn = n == Eigen::Dynamic ? Eigen::Dynamic : n+ni;
+  const auto dyn_nn = nn == Eigen::Dynamic ? dyn_n+dyn_ni : nn;
+  const auto make_HH = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,nn> HH =
+      Eigen::Matrix<Scalar,nn,nn>::Zero(dyn_nn,dyn_nn);
+    HH.topLeftCorner(dyn_n,dyn_n) = H;
+    return HH;
+  };
+  const Eigen::Matrix<Scalar,nn,nn> HH = make_HH();
+  const auto make_ff = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,1> ff = Eigen::Matrix<Scalar,nn,1>::Zero(dyn_nn,1);
+    ff.head(dyn_n) = f;
+    return ff;
+  };
+  const Eigen::Matrix<Scalar,nn,1> ff = make_ff();
+  const auto make_lblb = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,1> lblb =
+      Eigen::Matrix<Scalar,nn,1>::Zero(dyn_nn,1);
+    lblb.head(dyn_n) = lb;
+    lblb.tail(dyn_ni) = lbi;
+    return lblb;
+  };
+  const Eigen::Matrix<Scalar,nn,1> lblb = make_lblb();
+  const auto make_ubub = [&]()
+  {
+    Eigen::Matrix<Scalar,nn,1> ubub =
+      Eigen::Matrix<Scalar,nn,1>::Zero(dyn_nn,1);
+    ubub.head(dyn_n) = ub;
+    ubub.tail(dyn_ni) = ubi;
+    return ubub;
+  };
+  const Eigen::Matrix<Scalar,nn,1> ubub = make_ubub();
+  const auto make_AA = [&]()
+  {
+    Eigen::Matrix<Scalar,ni,nn> AA =
+      Eigen::Matrix<Scalar,ni,nn>::Zero(dyn_ni,dyn_nn);
+    AA.leftCols(dyn_n) = Ai;
+    AA.rightCols(dyn_ni) =
+      -Eigen::Matrix<Scalar,ni,ni>::Identity(dyn_ni,dyn_ni);
+    return AA;
+  };
+  const Eigen::Matrix<Scalar,ni,nn> AA = make_AA();
+  const Eigen::Matrix<Scalar,ni,1> bb = 
+   Eigen::Matrix<Scalar,ni,1>::Zero(dyn_ni,1);
+  // This ends make a system that looks like [H 0 Aᵀ;0 0 -I;A I 0] which is rank
+  // deficient.
+  const Eigen::Matrix<Scalar,nn,1> xx = 
+    igl::quadprog(HH,ff,AA,bb,lblb,ubub);
+  return xx.head(dyn_n);
+}
+
+
+template <typename Scalar, int n, int m>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Matrix<Scalar,m,n> & A,
+  const Eigen::Matrix<Scalar,m,1> & b,
+  const Eigen::Matrix<Scalar,n,1> & lb,
+  const Eigen::Matrix<Scalar,n,1> & ub)
+{
+  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
+  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
+  //std::cout<<igl::matlab_format(A,"A")<<std::endl;
+  //std::cout<<igl::matlab_format(b,"b")<<std::endl;
+  //std::cout<<igl::matlab_format(lb,"lb")<<std::endl;
+  //std::cout<<igl::matlab_format(ub,"ub")<<std::endl;
+  typedef Eigen::Matrix<Scalar,n,1> VectorSn;
+  typedef Eigen::Array<bool,n,1>    Arraybn;
+  assert( (lb.array() < ub.array() ).all() );
+  const int dyn_n = n == Eigen::Dynamic ? H.rows() : n;
+  VectorSn x(dyn_n);
+  VectorSn bc = VectorSn::Constant(dyn_n,1,-1e26);
+  Arraybn k = Arraybn::Constant(dyn_n,1,false);
+  Eigen::Index iter;
+  // n³ is probably way too conservative. 
+  for(iter = 0;iter<dyn_n*dyn_n*dyn_n;iter++)
+  {
+    // For dual contouring 99% of the time the active set is empty.
+    // Optimize for this common case.
+    x = min_quad_with_fixed(H,f,k,bc,A,b);
+    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
+    // constraint violations 
+    VectorSn vl = lb-x;
+    VectorSn vu = x-ub;
+
+    // try to add/remove constraints
+    Eigen::Index best_add = -1; Scalar worst_offense = 0;
+    bool add_lower;
+    Eigen::Index best_remove = -1; Scalar worst_lambda = 0;
+    for(Eigen::Index i = 0;i<dyn_n;i++)
+    {
+      if(vl(i)>worst_offense)
+      {
+        best_add = i;
+        add_lower = true;
+        worst_offense = vl(i);
+      }
+      if(vu(i)>worst_offense)
+      {
+        best_add = i;
+        add_lower = false;
+        worst_offense = vu(i);
+      }
+      // bias toward adding constraints
+      if(best_add<0 && k(i))
+      {
+        const Scalar sign = bc(i)==ub(i)?1:-1;
+        const Scalar lambda_i = sign * (H.row(i)*x+f(i));
+        //printf("  considering k(%d) (λ = %g)\n",i,lambda_i);
+        if(lambda_i > worst_lambda)
+        {
+          //printf("  removing k(%d) (λ = %g)\n",i,lambda_i);
+          best_remove = i;
+          worst_lambda = lambda_i;
+        }
+      }
+    }
+    // bias toward adding constraints
+    if(best_add >= 0)
+    {
+      const auto i = best_add;
+      assert(!k(i));
+      //add_lower ? printf("  adding lb(%d)\n",i) : printf("  adding lb(%d)\n",i);
+      bc(i) = add_lower ? lb(i) : ub(i);
+      k(i) = true;
+    }else if(best_remove >= 0)
+    {
+      const auto i = best_remove;
+      assert(k(i));
+      //printf("  removing k(%d) (λ = %g)\n",i,worst_lambda);
+      k(i) = false;
+    }else /*if(best_add < 0 && best_remove < 0)*/
+    {
+      return x;
+    }
+  }
+  // Should never happen.
+  assert(false && "quadprog failed after too many iterations");
+  //std::cout<<igl::eigen_format(H,"H")<<std::endl;
+  //std::cout<<igl::eigen_format(f,"f")<<std::endl;
+  //std::cout<<igl::eigen_format(lb,"lb")<<std::endl;
+  //std::cout<<igl::eigen_format(ub,"ub")<<std::endl;
+  return VectorSn::Zero(dyn_n);
+}
+
+template <typename Scalar, int n>
+IGL_INLINE Eigen::Matrix<Scalar,n,1> igl::quadprog(
+  const Eigen::Matrix<Scalar,n,n> & H,
+  const Eigen::Matrix<Scalar,n,1> & f,
+  const Eigen::Matrix<Scalar,n,1> & lb,
+  const Eigen::Matrix<Scalar,n,1> & ub)
+{
+  const int m = n == Eigen::Dynamic ? Eigen::Dynamic : 0;
+  return quadprog(
+    H,f,
+    Eigen::Matrix<Scalar,m,n>(0,H.cols()),
+    Eigen::Matrix<Scalar,m,1>(0,1),
+    lb,ub);
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> igl::quadprog<double, -1, -1>(Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, -1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)1) : ((((-1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, -1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&, Eigen::Matrix<double, -1, 1, ((Eigen::StorageOptions)0) | ((((-1) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((-1) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), -1, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 6, 1, ((Eigen::StorageOptions)0) | ((((6) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 6, 1> igl::quadprog<double, 6, 2>(Eigen::Matrix<double, 6, 6, ((Eigen::StorageOptions)0) | ((((6) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)1) : ((((6) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 6, 6> const&, Eigen::Matrix<double, 6, 1, ((Eigen::StorageOptions)0) | ((((6) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 6, 1> const&, Eigen::Matrix<double, 2, 6, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)1) : ((((6) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 6> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 6, 1, ((Eigen::StorageOptions)0) | ((((6) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 6, 1> const&, Eigen::Matrix<double, 6, 1, ((Eigen::StorageOptions)0) | ((((6) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((6) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 6, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> igl::quadprog<double, 2>(Eigen::Matrix<double, 2, 2, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)1) : ((((2) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 2> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&, Eigen::Matrix<double, 2, 1, ((Eigen::StorageOptions)0) | ((((2) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((2) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 2, 1> const&);
+// generated by autoexplicit.sh
+template Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> igl::quadprog<double, 3>(Eigen::Matrix<double, 3, 3, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)1) : ((((3) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 3> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&, Eigen::Matrix<double, 3, 1, ((Eigen::StorageOptions)0) | ((((3) == (1)) && ((1) != (1))) ? ((Eigen::StorageOptions)1) : ((((1) == (1)) && ((3) != (1))) ? ((Eigen::StorageOptions)0) : ((Eigen::StorageOptions)0))), 3, 1> const&);
+#endif

include/igl/quadprog.h

@@ -0,0 +1,87 @@
+#ifndef IGL_QUADPROG_H
+#define IGL_QUADPROG_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Solve a convex quadratic program. Optimized for small dense problems.
+  // Still works for Eigen::Dynamic (and then everything needs to be Dynamic).
+  //
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   lbi ≤ Ai x ≤ ubi
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   ni  #Ai or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   Ai  #Ai by #H list of linear equality constraint coefficients
+  //   lbi  #Ai list of linear equality lower bounds
+  //   ubi  #Ai list of linear equality upper bounds
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int ni>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,ni,n> & Ai,
+    const Eigen::Matrix<Scalar,ni,1> & lbi,
+    const Eigen::Matrix<Scalar,ni,1> & ubi,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   A x = b
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  //   m  #A or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   A  #A by #H list of linear equality constraint coefficients
+  //   b  #A list of linear equality lower bounds
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n, int m>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,m,n> & A,
+    const Eigen::Matrix<Scalar,m,1> & b,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+  // min_x ½ xᵀ H x + xᵀf
+  // subject to:
+  //   lb ≤ x ≤ u
+  //
+  // Templates:
+  //   Scalar  (e.g., double)
+  //   n  #H or Eigen::Dynamic if not known at compile time
+  // Inputs:
+  //   H  #H by #H quadratic coefficients (only lower triangle used)
+  //   f  #H linear coefficients
+  //   lb  #H list of lower bounds
+  //   ub  #H list of lower bounds
+  // Returns #H-long solution x
+  template <typename Scalar, int n>
+  IGL_INLINE Eigen::Matrix<Scalar,n,1> quadprog(
+    const Eigen::Matrix<Scalar,n,n> & H,
+    const Eigen::Matrix<Scalar,n,1> & f,
+    const Eigen::Matrix<Scalar,n,1> & lb,
+    const Eigen::Matrix<Scalar,n,1> & ub);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "quadprog.cpp"
+#endif 
+
+#endif 

tests/include/igl/min_quad_with_fixed.cpp

@@ -0,0 +1,17 @@
+#include <test_common.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/EPS.h>
+
+TEST_CASE("min_quad_with_fixed: dense", "[igl]" )
+{
+  const Eigen::Matrix<double,3,3> H = (Eigen::Matrix<double,3,3>(3,3)<<62,43,76,43,69,62,76,62,101).finished();
+  const Eigen::Matrix<double,3,1> f = (Eigen::Matrix<double,3,1>(3,1)<<9,8,5).finished();
+  const Eigen::Matrix<double,1,3> A = (Eigen::Matrix<double,1,3>(1,3)<<1,1,1).finished();
+  const Eigen::Matrix<double,1,1> b = (Eigen::Matrix<double,1,1>(1,1)<<2).finished();
+  const Eigen::Array<bool,3,1> k = (Eigen::Array<bool,3,1>()<<true,false,false).finished();
+  const Eigen::Matrix<double,3,1> bc = (Eigen::Matrix<double,3,1>(3,1)<<1,0,0).finished();
+  const Eigen::Matrix<double,3,1> x = igl::min_quad_with_fixed(H,f,k,bc,A,b);
+  REQUIRE(abs(x(0)- 1.0)<igl::EPS<double>());
+  REQUIRE(abs(x(1)- 1.5)<igl::EPS<double>());
+  REQUIRE(abs(x(2)- -.5)<igl::EPS<double>());
+}