lilypond-0.1.1

diff --git a/flower/include/choleski.hh b/flower/include/choleski.hh
index a3a02a11e0782d7d7c49ec41794fd1a547bbea9f..abd19ee67635a42476b06dca6944d9691a2d55d3 100644 (file)
--- a/flower/include/choleski.hh
+++ b/flower/include/choleski.hh
@@ -30,7 +30,7 @@ struct Choleski_decomposition {
     solve Px = rhs
     */
     Vector solve(Vector rhs) const;
-
+    void solve (Vector &dest, Vector const &rhs)const;
     Vector operator * (Vector rhs) const { return solve (rhs); }
   /**
     return the inverse of the matrix P.
@@ -41,6 +41,8 @@ struct Choleski_decomposition {
     */
     Matrix original() const;
 private:
+    void full_matrix_solve(Vector &,Vector const&)const;
+    void band_matrix_solve(Vector &, Vector const&)const;
     void full_matrix_decompose(Matrix const & P);
     void band_matrix_decompose(Matrix const &P);
         

diff --git a/flower/include/vector.hh b/flower/include/vector.hh
index fdeab97f1391b84997d7528d9e0334588a5b6aa7..7f7573534498fdd125edccaf5a2dcaa51ece4187 100644 (file)
--- a/flower/include/vector.hh
+++ b/flower/include/vector.hh
@@ -23,6 +23,11 @@ public:
        dat.set_size(n);
        fill(0);
     }
+    void set_dim(int i)
+    {
+       dat.set_size(i);
+    }
+       
     void insert(Real v, int i) {
        dat.insert(v,i);
     }

diff --git a/flower/matrix.cc b/flower/matrix.cc
index 78e3d777b01237593e144a6c10851257c95f55b5..7e1854827b223149729514389bb4d261e13f5833 100644 (file)
--- a/flower/matrix.cc
+++ b/flower/matrix.cc
@@ -108,6 +108,10 @@ Matrix::operator=(Matrix const &m)
 int
 Matrix::band_i()const
 {
+    if ( band_b() ) {
+       Diagonal_storage const * diag = (Diagonal_storage*) dat;
+       return diag->band_size_i();
+    }
     int starty = dim();
     while (starty >= 0 ) {
        for ( int i = starty, j = 0; i < dim(); i++, j++ )
@@ -206,9 +210,13 @@ operator /(Matrix const& m1,Real a)
     return  m;
 }
 
+/*
+  ugh. Only works for square matrices.
+ */
 void
 Matrix::transpose()            // delegate to storage?
 {
+#if 1
     for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) {
        if (i >= j)
            continue;
@@ -216,6 +224,7 @@ Matrix::transpose()         // delegate to storage?
        dat->elem(i,j) = dat->elem(j,i);
        dat->elem(j,i)=r;
     }
+#endif
 }
 
 Matrix
@@ -252,6 +261,7 @@ Matrix::set_product(Matrix const &m1, Matrix const &m2)
     
     if (m1.dat->try_right_multiply(dat, m2.dat))
        return; 
+    
     for (int i=0, j=0; dat->mult_ok(i,j);
         dat->mult_next(i,j)) {
        Real r=0.0;