From 71b1d9ad2c8494181b85291d9f9aac064aac58d1 Mon Sep 17 00:00:00 2001 From: fred Date: Sun, 24 Mar 2002 19:49:53 +0000 Subject: [PATCH] lilypond-0.0.78 --- flower/choleski.cc | 66 +++++++++++++++---- flower/matrix-storage.cc | 139 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 192 insertions(+), 13 deletions(-) create mode 100644 flower/matrix-storage.cc diff --git a/flower/choleski.cc b/flower/choleski.cc index 202e184440..d619b97460 100644 --- 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 index 0000000000..76b89f7f36 --- /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 +*/ + +#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 +Matrix_storage::row(int n) const +{ + Array r; + for (int j = 0; j < cols(); j++) + r.push(elem(n,j)); + return r; +} + +Array +Matrix_storage::column(int n) const +{ + Array 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; +} -- 2.39.5