From ad26f79eed1c16bee540d9ec52a411912346fc74 Mon Sep 17 00:00:00 2001
From: fred <fred>
Date: Sun, 24 Mar 2002 19:50:19 +0000
Subject: [PATCH] lilypond-0.1.1

---
 flower/include/choleski.hh |  4 +++-
 flower/include/vector.hh   |  5 +++++
 flower/matrix.cc           | 10 ++++++++++
 3 files changed, 18 insertions(+), 1 deletion(-)

diff --git a/flower/include/choleski.hh b/flower/include/choleski.hh
index a3a02a11e0..abd19ee676 100644
--- 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 fdeab97f13..7f75735344 100644
--- 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 78e3d777b0..7e1854827b 100644
--- 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;
-- 
2.39.5