patch::: 1.2.9.jcn2

[lilypond.git] / lily / qlpsolve.cc

diff --git a/lily/qlpsolve.cc b/lily/qlpsolve.cc
index 78e07042773efb7ed8b91e7e5025b92d99832a0f..3d64b3779289b4ccd58d559e1caa7375784483e9 100644 (file)
--- a/lily/qlpsolve.cc
+++ b/lily/qlpsolve.cc
@@ -3,7 +3,7 @@
 
   source file of the GNU LilyPond music typesetter
 
-  (c) 1996,  1997--1998 Han-Wen Nienhuys <hanwen@stack.nl>
+  (c) 1996,  1997--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
 
   TODO:
   try fixed point arithmetic, to speed up lily.
@@ -22,13 +22,13 @@ Active_constraints::status() const
   String s ("Active|Inactive [");
   for (int i=0; i< active.size(); i++)
     {
-      s += String (active[i]) + " ";
+      s += to_str (active[i]) + " ";
     }
 
   s+="| ";
   for (int i=0; i< inactive.size(); i++)
     {
-      s += String (inactive[i]) + " ";
+      s += to_str (inactive[i]) + " ";
     }
   s+="]";
 
@@ -70,78 +70,81 @@ Active_constraints::get_lagrange (Vector gradient)
 }
 
 void
-Active_constraints::add (int k)
+Active_constraints::add_constraint (int k)
 {
   // add indices
   int cidx=inactive[k];
-  active.push (cidx);
-
-  inactive.swap (k,inactive.size()-1);
-  inactive.pop();
 
   Vector a (opt->cons_[cidx]);
   // update of matrices
   Vector Ha = H*a;
   Real aHa = a*Ha;
   Vector addrow (Ha.dim());
-  if (abs (aHa) > EPS)
+  bool degenerate =  (abs (aHa) < EPS);
+
+  if (degenerate)
     {
+      warning (String ("Active_constraints::add ():")
+       + _("degenerate constraints"));
+      DOUT << "Ha = "  << Ha.str () << '\n';
       /*
        a != 0, so if Ha = O(EPS), then
        Ha * aH / aHa = O(EPS^2/EPS)
 
        if H*a == 0, the constraints are dependent.
-         */
-      H -= Matrix (Ha/aHa , Ha);
-
+      */
+      degenerate_count_i_ ++;
+    }
+  if (!degenerate)
+    {
+      active.push (cidx);
+      inactive.swap (k,inactive.size()-1);
+      inactive.pop();
 
-      /*
-         sorry, don't know how to justify this. ..
-         */
+      H -= Matrix (Ha/aHa , Ha);
+      
       addrow=Ha;
-      addrow/= aHa;
+      addrow /= aHa;
       A -= Matrix (A*a, addrow);
       A.insert_row (addrow,A.rows());
-    }else
-      WARN << _("degenerate constraints");
+    }
 }
 
 void
-Active_constraints::drop (int k)
+Active_constraints::drop_constraint (int k)
 {
   int q=active.size()-1;
 
-  // drop indices
-  inactive.push (active[k]);
-  active.swap (k,q);
-  A.swap_rows (k,q);
-  active.pop();
 
   Vector a (A.row (q));
   if (a.norm() > EPS)
     {
+      // drop indices
+      inactive.push (active[k]);
+      active.swap (k,q);
+      A.swap_rows (k,q);
+      active.pop();
       /*
 
        */
-      Real q = a*opt->quad_*a;
-      Matrix aaq (a,a/q);
+      Real aqa = a*opt->quad_*a;
+      Matrix aaq (a,a/aqa);
       H += aaq;
       A -= A*opt->quad_*aaq;
-    }else
-      WARN << _("degenerate constraints");
-#ifndef NDEBUG
-  Vector rem_row (A.row (q));
-  assert (rem_row.norm() < EPS);
-#endif
-
-  A.delete_row (q);
+ 
+      A.delete_row (q);
+    }else {
+      degenerate_count_i_ ++;
+      warning (String ("Active_constraints::drop ():")
+       + _("degenerate constraints"));
+    }
 }
 
 
 Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
-  :       A(0,op->dim()),
-         H(op->dim()),
-         opt (op)
+  : A(0,op->dim()),
+    H(op->dim()),
+    opt (op)
 {
   for (int i=0; i < op->cons_.size(); i++)
     inactive.push (i);
@@ -152,6 +155,8 @@ Active_constraints::Active_constraints (Ineq_constrained_qp const *op)
     */
   H=chol.inverse();
   OK();
+
+  degenerate_count_i_ = 0;
 }
 
 /** Find the optimum which is in the planes generated by the active