release: 1.0.1

[lilypond.git] / flower / choleski.cc

diff --git a/flower/choleski.cc b/flower/choleski.cc
index 4a98301af3ec39e63d61ebee2590e5a2dce6b22a..8fa18c67f52433881da099ae157871c8cc633701 100644 (file)
--- a/flower/choleski.cc
+++ b/flower/choleski.cc
@@ -3,7 +3,7 @@
 
   source file of the Flower Library
 
-  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+  (c)  1997--1998 Han-Wen Nienhuys <hanwen@cs.uu.nl>
 */
 
 #include "choleski.hh"
@@ -13,7 +13,7 @@ const Real EPS = 1e-7;                // so sue me. Hard coded
 //#define PARANOID
 
 void
-Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs)const
+Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs) const
 {
   int n= rhs.dim();
   assert (n == L.dim());
@@ -24,34 +24,34 @@ Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs)const
   // forward substitution
   for (int i=0; i < n; i++) 
     {
-       Real sum (0.0);
-       for (int j=0; j < i; j++)
-           sum += y (j) * L(i,j);
-       y (i) = (rhs (i) - sum)/L(i,i);
+      Real sum (0.0);
+      for (int j=0; j < i; j++)
+       sum += y (j) * L(i,j);
+      y (i) = (rhs (i) - sum)/L(i,i);
     }
   
   for (int i=0; i < n; i++)
-       y (i) /= D(i);
+    y (i) /= D(i);
 
   // backward subst
   Vector &x (out);             // using input as return val.
   for (int i=n-1; i >= 0; i--) 
     {
-       Real sum (0.0);
-       for (int j=i+1; j < n; j++)
-           sum += L(j,i)*x (j);
-       x (i) = (y (i) - sum)/L(i,i);
+      Real sum (0.0);
+      for (int j=i+1; j < n; j++)
+       sum += L(j,i)*x (j);
+      x (i) = (y (i) - sum)/L(i,i);
     }
 }
 
 void
-Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs)const
+Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs) const
 {
   int n= rhs.dim();
   int b = L.band_i();
   assert (n == L.dim());
 
-    out.set_dim (n);
+  out.set_dim (n);
   
   Vector y;
   y.set_dim (n);
@@ -59,38 +59,38 @@ Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs)const
   // forward substitution
   for (int i=0; i < n; i++) 
     {
-       Real sum (0.0);
-       for (int j= 0 >? i - b; j < i; j++)
-           sum += y (j) * L(i,j);
-       y (i) = (rhs (i) - sum)/L(i,i);
+      Real sum (0.0);
+      for (int j= 0 >? i - b; j < i; j++)
+       sum += y (j) * L(i,j);
+      y (i) = (rhs (i) - sum)/L(i,i);
     }
   for (int i=0; i < n; i++)
-       y (i) /= D(i);
+    y (i) /= D(i);
 
   // backward subst
   Vector &x (out);             // using input as return val.
   for (int i=n-1; i >= 0; i--) 
     {
-       Real sum (0.0);
-       for (int j=i+1; j <= i + b&&j < n ; j++)
-           sum += L(j,i)*x (j);
-       x (i) = (y (i) - sum)/L(i,i);
+      Real sum (0.0);
+      for (int j=i+1; j <= i + b&&j < n ; j++)
+       sum += L(j,i)*x (j);
+      x (i) = (y (i) - sum)/L(i,i);
     }
 }
 
 void
-Choleski_decomposition::solve (Vector &x, Vector const &rhs)const
+Choleski_decomposition::solve (Vector &x, Vector const &rhs) const
 {
-  if (L.band_b()) 
+  if (band_b_) 
     {
-       band_matrix_solve (x,rhs);
+      band_matrix_solve (x,rhs);
     }
   else
-       full_matrix_solve (x,rhs);
+    full_matrix_solve (x,rhs);
 }
 
 Vector
-Choleski_decomposition::solve (Vector rhs)const
+Choleski_decomposition::solve (Vector rhs) const
 {
   Vector r;
   solve (r, rhs);
@@ -101,23 +101,23 @@ void
 Choleski_decomposition::full_matrix_decompose (Matrix const & P)
 {
  
-   int n = P.dim();
+  int n = P.dim();
   L.unit();
   for (int k= 0; k < n; k++) 
     {
-       for (int j = 0; j < k; j++)
-         {
-           Real sum (0.0);
-           for (int l=0; l < j; l++)
-               sum += L(k,l)*L(j,l)*D(l);
-           L(k,j) = (P(k,j) - sum)/D(j);
-         }
-       Real sum=0.0;
+      for (int j = 0; j < k; j++)
+       {
+         Real sum (0.0);
+         for (int l=0; l < j; l++)
+           sum += L(k,l)*L(j,l)*D(l);
+         L(k,j) = (P(k,j) - sum)/D(j);
+       }
+      Real sum=0.0;
        
-       for (int l=0; l < k; l++)
-           sum += sqr (L(k,l))*D(l);
-       Real d = P(k,k) - sum;
-       D(k) = d;
+      for (int l=0; l < k; l++)
+       sum += sqr (L(k,l))*D(l);
+      Real d = P(k,k) - sum;
+      D(k) = d;
     }
 
 }
@@ -131,22 +131,22 @@ Choleski_decomposition::band_matrix_decompose (Matrix const &P)
   
   for (int i= 0; i < n; i++) 
     {
-       for (int j = 0 >? i - b; j < i; j++)
-         {
-           Real sum (0.0);
-           for (int l=0 >? i - b; l < j; l++)
-               sum += L(i,l)*L(j,l)*D(l);
-           L(i,j) = (P(i,j) - sum)/D(j);
-         }
-       Real sum=0.0;
+      for (int j = 0 >? i - b; j < i; j++)
+       {
+         Real sum (0.0);
+         for (int l=0 >? i - b; l < j; l++)
+           sum += L(i,l)*L(j,l)*D(l);
+         L(i,j) = (P(i,j) - sum)/D(j);
+       }
+      Real sum=0.0;
        
-       for (int l=0 >? i - b; l < i; l++)
-           sum += sqr (L(i,l))*D(l);
-       Real d = P(i,i) - sum;
-       D(i) = d;
+      for (int l=0 >? i - b; l < i; l++)
+       sum += sqr (L(i,l))*D(l);
+      Real d = P(i,i) - sum;
+      D(i) = d;
     }
-  L.try_set_band();
-  assert ( L.band_i() == P.band_i ());
+  L.set_band();
+  band_b_ = true;
 }
 
 
@@ -157,16 +157,19 @@ Choleski_decomposition::band_matrix_decompose (Matrix const &P)
    */
 
 Choleski_decomposition::Choleski_decomposition (Matrix const & P)
-   : L(P.dim()), D(P.dim ())
+  : L(P.dim()), D(P.dim ())
 { 
 #ifdef PARANOID
   assert ((P-P.transposed()).norm ()/P.norm () < EPS);
 #endif
-  if  (P.band_b()) 
-       band_matrix_decompose (P);
+  band_b_ = false;
+
+  int b = P.calc_band_i ();
+
+  if (b <= P.dim ()/2)  
+    band_matrix_decompose (P);
   else
-       full_matrix_decompose (P);
- 
+    full_matrix_decompose (P);
 
 #ifdef PARANOID
   assert ((original()-P).norm () / P.norm () < EPS);
@@ -190,10 +193,10 @@ Choleski_decomposition::inverse() const
   Vector inv (n);
   for (int i = 0; i < n; i++) 
     {
-       e_i.set_unit (i);
-       solve (inv, e_i);
-       for (int j = 0 ; j<n; j++)
-           invm (i,j) = inv (j);
+      e_i.set_unit (i);
+      solve (inv, e_i);
+      for (int j = 0 ; j<n; j++)
+       invm (i,j) = inv (j);
     }
   
 #ifdef PARANOID