2 qlpsolve.cc -- implement Active_constraints, Inactive_iter
4 source file of the GNU LilyPond music typesetter
6 (c) 1996, 1997 Han-Wen Nienhuys <hanwen@stack.nl>
9 try fixed point arithmetic, to speed up lily.
12 #include "ineq-constrained-qp.hh"
13 #include "qlpsolve.hh"
15 #include "choleski.hh"
17 const Real TOL=1e-1; // roughly 1/30 mm
20 Active_constraints::status() const
22 String s ("Active|Inactive [");
23 for (int i=0; i< active.size(); i++)
25 s += String (active[i]) + " ";
29 for (int i=0; i< inactive.size(); i++)
31 s += String (inactive[i]) + " ";
39 Active_constraints::OK()
44 assert (active.size() +inactive.size () == opt->cons.size ());
45 assert (H.dim() == opt->dim ());
46 assert (active.size() == A.rows ());
49 for (int i=0; i < opt->cons.size(); i++)
51 for (int i=0; i < active.size(); i++)
56 for (int i=0; i < inactive.size(); i++)
61 for (int i=0; i < allcons.size(); i++)
62 assert (allcons[i] == 1);
67 Active_constraints::get_lagrange (Vector gradient)
73 Active_constraints::add (int k)
79 inactive.swap (k,inactive.size()-1);
82 Vector a (opt->cons[cidx]);
86 Vector addrow (Ha.dim());
90 a != 0, so if Ha = O(EPS), then
91 Ha * aH / aHa = O(EPS^2/EPS)
93 if H*a == 0, the constraints are dependent.
95 H -= Matrix (Ha/aHa , Ha);
99 sorry, don't know how to justify this. ..
103 A -= Matrix (A*a, addrow);
104 A.insert_row (addrow,A.rows());
106 WARN << "degenerate constraints";
110 Active_constraints::drop (int k)
112 int q=active.size()-1;
115 inactive.push (active[k]);
120 Vector a (A.row (q));
126 Real q = a*opt->quad*a;
129 A -= A*opt->quad*aaq;
131 WARN << "degenerate constraints";
133 Vector rem_row (A.row (q));
134 assert (rem_row.norm() < EPS);
141 Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
146 for (int i=0; i < op->cons.size(); i++)
148 Choleski_decomposition chol (op->quad);
157 /** Find the optimum which is in the planes generated by the active
161 Active_constraints::find_active_optimum (Vector g)