lilypond-0.0.1

diff --git a/qlpsolve.hh b/qlpsolve.hh
new file mode 100644 (file)
index 0000000..4d29bcc
--- /dev/null
+++ b/qlpsolve.hh
@@ -0,0 +1,78 @@
+#include "qlp.hh"
+#include "matrix.hh"
+
+
+class Active_constraints {
+    friend class Inactive_iter;
+    
+
+    Matrix A,H;
+    svec<int> active;
+    svec<int> inactive;                // actually this is a set, not an array.
+    const Ineq_constrained_qp *opt;
+
+public:
+    String status()const;
+    
+    Vector vec(int k) const { return opt->cons[k]; }
+    Real rhs(int k) const { return opt->consrhs[k]; }
+    
+    /// drop constraint
+    void drop (int k);
+    /** drop constraint k from the active set. k is the index of the
+    constraint in #active#
+    
+    */
+    
+    ///     add constraint j
+    void add(int j);
+    /**
+    add constraint j to the active set j is the index of the
+    constraint in #inactive#   
+    */
+
+    /// exchange in and out.
+    void exchange(int in, int out) { add(in); drop (out); }
+    
+    /// 
+    Vector find_active_optimum(Vector g);
+
+    /// get lagrange multipliers.
+    Vector get_lagrange(Vector v);
+
+    Active_constraints(Ineq_constrained_qp const *op);
+    /** construct: no constraints active, n vars. Put the equalities
+     into the constraints.  */
+
+    /// check invariants
+    void OK();
+};
+
+/**
+    This class represents the set of active (binding) constraints
+    which can be active while the QLP algorithm is in a feasible
+    point. The active constraints are numbered.
+    If the constraints are of the form
+
+      A^T*x >= b
+
+    then the binding constraints are those where the >= is equality.
+    
+  */
+
+///
+class Inactive_iter {
+    int j;
+    Active_constraints const* ac;
+public:
+    Inactive_iter(Active_constraints const &c) { ac=&c; j=0; }
+    int idx() const { return j; }
+    void operator ++(int) { j++; }
+    int constraint_id() const { return ac->inactive[j]; }
+    Vector vec() const { return ac->vec(constraint_id()); }
+    Real rhs() const { return ac->rhs(constraint_id()); }
+    bool ok() const { return j < ac->inactive.sz(); }
+};
+/**
+    loop through the inactive constraints.
+  */