lilypond-0.0.78

diff --git a/flower/choleski.cc b/flower/choleski.cc
index 202e184440cb62f68153b552a8c1038337a41516..d619b97460f724c428e00758a32e94fbb4aa2602 100644 (file)
--- a/flower/choleski.cc
+++ b/flower/choleski.cc
@@ -9,6 +9,9 @@
 #include "choleski.hh"
 const Real EPS = 1e-7;         // so sue me. Hard coded
 
+// for testing new Matrix_storage.
+//#define PARANOID
+
 Vector
 Choleski_decomposition::solve(Vector rhs)const
 {
@@ -37,20 +40,11 @@ Choleski_decomposition::solve(Vector rhs)const
     return x;
 }
 
-/*
-  Standard matrix algorithm.
-  Should add support for banded matrices
-  */
-
-Choleski_decomposition::Choleski_decomposition(Matrix P)
-    : L(P.dim()), D(P.dim())
+void
+Choleski_decomposition::full_matrix_decompose(Matrix const & P)
 {
-    int n = P.dim();
-
-#ifdef PARANOID
-    assert((P-P.transposed()).norm()/P.norm() < EPS);
-#endif
-    
+ 
+   int n = P.dim();
     L.unit();
     for (int k= 0; k < n; k++) {
        for (int j = 0; j < k; j++){
@@ -67,6 +61,52 @@ Choleski_decomposition::Choleski_decomposition(Matrix P)
        D(k) = d;
     }
 
+}
+
+void
+Choleski_decomposition::band_matrix_decompose(Matrix const &P)
+{
+    int n = P.dim();
+    int b = P.band_i();
+    L.unit();
+    
+    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 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());
+}
+
+
+
+
+/*
+  Standard matrix algorithm.
+   */
+
+Choleski_decomposition::Choleski_decomposition(Matrix const & P)
+   : 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);
+    else
+       full_matrix_decompose(P);
+ 
+
 #ifdef PARANOID
     assert((original()-P).norm() / P.norm() < EPS);
 #endif

diff --git a/flower/matrix-storage.cc b/flower/matrix-storage.cc
new file mode 100644 (file)
index 0000000..76b89f7
--- /dev/null
+++ b/flower/matrix-storage.cc
@@ -0,0 +1,139 @@
+/*
+  matrix-storage.cc -- implement Matrix_storage
+
+  source file of the Flower Library
+
+  (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
+*/
+
+#include "full-storage.hh"
+#include "diagonal-storage.hh"
+
+void
+Matrix_storage::set_addition_result(Matrix_storage *&dat, Matrix_storage *right)
+{
+    if (dat && dat->name() == Diagonal_storage::static_name() 
+        && right->name() == Diagonal_storage::static_name()) {
+       Diagonal_storage *L = (Diagonal_storage*)dat;
+       Diagonal_storage* R = (Diagonal_storage*) right;
+
+       if ( R->band_size_i() > L->band_size_i()) { 
+           L->set_band_size(R->band_size_i());
+           return ;
+       }
+    }
+    if  (!dat || !dat->is_type_b(Full_storage::static_name() )) {
+
+       Matrix_storage *new_stor = (dat)? new Full_storage(dat) : 
+           new Full_storage( right->rows(), right->cols());
+       delete dat;
+       dat = new_stor;
+    }
+}
+
+Matrix_storage*
+Matrix_storage::get_product_result(Matrix_storage*left, Matrix_storage*right)
+{
+    Matrix_storage* dest =0;
+    set_product_result(dest, left,right);
+    return dest;
+}
+    
+
+/*
+  hairy
+ */
+void
+Matrix_storage::set_product_result(Matrix_storage*&dest, 
+                          Matrix_storage*left, Matrix_storage*right)
+{
+    if ( left->name() == Diagonal_storage::static_name() 
+        && right->name() == Diagonal_storage::static_name()) {
+       Diagonal_storage *L = (Diagonal_storage*)left;
+       Diagonal_storage* R = (Diagonal_storage*) right;
+
+       if  (L->band_size_i() + R->band_size_i() < L->dim()/2 ) {
+           if (dest ->name() != Diagonal_storage::static_name()){
+               delete dest;
+               dest = new Diagonal_storage;
+           }
+        
+           dest->set_size(L->dim());
+           return;
+       }
+    }
+
+    if ( dest && dest->name() == Full_storage::static_name()) {
+       dest->set_size(left->rows(), right->cols());
+    } else {
+       delete dest;
+       dest = new Full_storage( left->rows(), right->cols());
+    }
+}
+
+IMPLEMENT_IS_TYPE_B(Matrix_storage);
+
+Matrix_storage *
+Matrix_storage::get_full(int n, int m)
+{
+    return new Full_storage(n,m);
+}
+
+
+
+ bool
+Matrix_storage::try_right_multiply(Matrix_storage *, 
+                                  const Matrix_storage *)const
+{
+    return false;
+}
+
+Array<Real>
+Matrix_storage::row(int n) const
+{
+    Array<Real> r;
+    for (int j = 0; j < cols(); j++)
+       r.push(elem(n,j));
+    return r;
+}
+
+Array<Real>
+Matrix_storage::column(int n) const
+{
+    Array<Real> r;
+    for (int i = 0; i < rows(); i++)
+       r.push(elem(i,n));
+    return r;
+}
+
+void
+Matrix_storage::set_size(int rows, int cols)
+{
+        
+    resize(rows,cols);
+}
+       
+void
+Matrix_storage::set_size(int rows)
+{
+        
+       resize(rows);
+}
+       
+
+void
+Matrix_storage::set_band(Matrix_storage *&mat, int b)
+{
+    Matrix_storage* ns = new Diagonal_storage(mat, b);
+    delete mat;
+    mat=ns;
+}
+       
+
+void
+Matrix_storage::set_full(Matrix_storage *&mat)
+{
+    Matrix_storage* ns = new Full_storage(mat);
+    delete mat;
+    mat=ns;
+}