From: fred Date: Fri, 11 Oct 1996 20:42:20 +0000 (+0000) Subject: flower-1.0.2 X-Git-Tag: release/1.5.59~7093 X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=714d8266107a76f6e44ba6485b9c1cfd2806bcfd;p=lilypond.git flower-1.0.2 --- diff --git a/flower/Sources.make b/flower/Sources.make new file mode 100644 index 0000000000..18fb18cb24 --- /dev/null +++ b/flower/Sources.make @@ -0,0 +1,12 @@ + +cc=lgetopt.cc string.cc dataf.cc textdb.cc unionfind.cc \ + smat.cc matrix.cc choleski.cc vector.cc dstream.cc\ + matdebug.cc + +templatecc=cursor.cc list.cc tsmat.cc +inl=findcurs.inl link.inl list.inl +hh=cursor.hh cursor.inl lgetopt.hh link.hh list.hh dstream.hh \ + string.hh stringutil.hh vray.hh textdb.hh textstr.hh assoc.hh\ + findcurs.hh unionfind.hh compare.hh handle.hh matrix.hh\ + smat.hh vsmat.hh vector.hh real.hh choleski.hh\ + tsmat.hh tvsmat.hh diff --git a/flower/choleski.cc b/flower/choleski.cc new file mode 100644 index 0000000000..82740a27c1 --- /dev/null +++ b/flower/choleski.cc @@ -0,0 +1,97 @@ +#include "choleski.hh" +const Real EPS = 1e-7; // so sue me. Hard coded + +Vector +Choleski_decomposition::solve(Vector rhs)const +{ + int n= rhs.dim(); + assert(n == L.dim()); + Vector y(n); + + // 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); + } + for (int i=0; i < n; i++) { + assert(D(i)); + y(i) /= D(i); + } + + // backward subst + Vector x(n); + 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); + } + return x; +} + +/* + Standard matrix algorithm. + */ + +Choleski_decomposition::Choleski_decomposition(Matrix P) + : L(P.dim()), D(P.dim()) +{ + int n = P.dim(); + assert((P-P.transposed()).norm()/P.norm() < EPS); + + 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 l=0; l < k; l++) + sum += sqr(L(k,l))*D(l); + Real d = P(k,k) - sum; + D(k) = d; + } + +#ifdef NDEBUG + assert((original()-P).norm() / P.norm() < EPS); +#endif +} + +Matrix +Choleski_decomposition::original() const +{ + Matrix T(L.dim()); + T.set_diag(D); + return L*T*L.transposed(); +} + +Matrix +Choleski_decomposition::inverse() const +{ + int n=L.dim(); + Matrix invm(n); + Vector e_i(n); + for (int i = 0; i < n; i++) { + e_i.set_unit(i); + Vector inv(solve(e_i)); + for (int j = 0 ; jmult_ok(i,j); dat->mult_next(i,j)) + r += sqr(dat->elem(i,j)); + return sqrt(r); +} + +//inline +Real +Matrix::operator()(int i,int j) const +{ + assert(i >= 0 && j >= 0); + assert(i < rows() && j < cols()); + return dat->elem(i,j); +} + +//inline +Real & +Matrix::operator()(int i, int j) +{ + assert(i >= 0 && j >= 0); + assert(i < rows() && j < cols()); + return dat->elem(i,j); +} + +void +Matrix::fill(Real r) +{ + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j)=r; +} + +void +Matrix::set_diag(Real r) +{ + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j)=(i==j) ? r: 0.0; +} + +void +Matrix::set_diag(Vector d) +{ + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j)=(i==j) ? d(i): 0.0; +} + +void +Matrix::operator+=(const Matrix&m) +{ + assert(m.cols() == cols()); + assert(m.rows() == rows()); + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j) += m(i,j); +} + +void +Matrix::operator-=(const Matrix&m) +{ + assert(m.cols() == cols()); + assert(m.rows() == rows()); + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j) -= m(i,j); +} + + +void +Matrix::operator*=(Real a) +{ + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j) *= a; +} + +void +Matrix::operator=(const Matrix&m) +{ + if (&m == this) + return ; + delete dat; + dat = m.dat->clone(); +} + +Matrix::Matrix(const Matrix &m) +{ + m.OK(); + + dat = m.dat->clone(); +} + + +Matrix::Matrix(int n, int m) +{ + dat = virtual_smat::get_full(n,m); + fill(0); +} + +Matrix::Matrix(int n) +{ + dat = virtual_smat::get_full(n,n); + fill(0); +} + +Matrix::Matrix(Vector v, Vector w) +{ + dat = virtual_smat::get_full(v.dim(), w.dim()); + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dat->elem(i,j)=v(i)*w(j); +} + + +Vector +Matrix::row(int k) const +{ + int n=cols(); + + + Vector v(n); + for(int i=0; i < n; i++) + v(i)=dat->elem(k,i); + + return v; +} + +Vector +Matrix::col(int k) const +{ + int n=rows(); + Vector v(n); + for(int i=0; i < n; i++) + v(i)=dat->elem(i,k); + return v; +} + +Vector +Matrix::left_multiply(const Vector& v) const +{ + Vector dest(v.dim()); + assert(dat->cols()==v.dim()); + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dest(i)+= dat->elem(j,i)*v(j); + return dest; +} + +Vector +Matrix::operator *(const Vector& v) const +{ + Vector dest(rows()); + assert(dat->cols()==v.dim()); + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) + dest(i)+= dat->elem(i,j)*v(j); + return dest; +} + +Matrix +operator /(Matrix const& m1,Real a) +{ + Matrix m(m1); + m /= a; + return m; +} + +void +Matrix::transpose() // delegate to storage? +{ + for (int i=0, j=0; dat->mult_ok(i,j); dat->mult_next(i,j)) { + if (i >= j) + continue; + Real r=dat->elem(i,j); + dat->elem(i,j) = dat->elem(j,i); + dat->elem(j,i)=r; + } +} + +Matrix +Matrix::operator-() const +{ + OK(); + Matrix m(*this); + m*=-1.0; + return m; +} + +Matrix +Matrix::transposed() const +{ + Matrix m(*this); + m.transpose(); + return m; +} + + +/* should do something smarter: bandmatrix * bandmatrix is also banded matrix. */ +Matrix +operator *(const Matrix &m1, const Matrix &m2) +{ + Matrix result(m1.rows(), m2.cols()); + result.set_product(m1,m2); + return result; +} + +void +Matrix::set_product(const Matrix &m1, const Matrix &m2) +{ + assert(m1.cols()==m2.rows()); + assert(cols()==m2.cols() && rows()==m1.rows()); + + for (int i=0, j=0; dat->mult_ok(i,j); + dat->mult_next(i,j)) { + Real r=0.0; + for (int k = 0; k < m1.cols(); k++) + r += m1(i,k)*m2(k,j); + dat->elem(i,j)=r; + } +} + +void +Matrix::insert_row(Vector v, int k) +{ + assert(v.dim()==cols()); + dat->insert_row(k); + for (int j=0; j < cols(); j++) + dat->elem(k,j)=v(j); +} + + +void +Matrix::swap_columns(int c1, int c2) +{ + assert(c1>=0&& c1 < cols()&&c2 < cols() && c2 >=0); + for (int i=0; i< rows(); i++) { + Real r=dat->elem(i,c1); + dat->elem(i,c1) = dat->elem(i,c2); + dat->elem(i,c2)=r; + } +} + +void +Matrix::swap_rows(int c1, int c2) +{ + assert(c1>=0&& c1 < rows()&&c2 < rows() && c2 >=0); + for (int i=0; i< cols(); i++) { + Real r=dat->elem(c1,i); + dat->elem(c1,i) = dat->elem(c2,i); + dat->elem(c2,i)=r; + } +} + + +int +Matrix::dim() const +{ + assert(cols() == rows()); + return rows(); +} + diff --git a/flower/vector.cc b/flower/vector.cc new file mode 100644 index 0000000000..2cbfdfc899 --- /dev/null +++ b/flower/vector.cc @@ -0,0 +1,20 @@ +#include "vector.hh" +#include "string.hh" + +Vector::Vector(const Vector&n) + :dat(n.dat) +{ +} + +Vector +Vector::operator-() const +{ + Vector v(*this); v*=-1; return v; +} + +void +Vector::set_unit(int j) +{ + fill(0.0); + dat[j] = 1.0; +} diff --git a/flower/vector.hh b/flower/vector.hh new file mode 100644 index 0000000000..1929c674a9 --- /dev/null +++ b/flower/vector.hh @@ -0,0 +1,111 @@ +#ifndef VECTOR_HH +#define VECTOR_HH + +#include +#include "real.hh" +#include "vray.hh" + +class Dstream; +class String; +void set_matrix_debug(Dstream&ds); + +/// a row of numbers +class Vector { + svec dat; +public: + void OK() const { dat.OK();} + int dim() const { return dat.sz(); } + Vector() { } + Vector(const Vector&n); + Vector(int n) { + dat.set_size(n); + fill(0); + } + void insert(Real v, int i) { + dat.insert(v,i); + } + void del(int i) { dat.del(i); } + operator String() const; + void fill(Real r) { + for (int i=0; i < dim(); i++) + dat[i] =r; + } + + void operator +=(Vector v) { + assert(v.dim() == dim()); + for (int i=0; i < dim(); i++) + dat[i] += v.dat[i]; + } + + void operator /=(Real a) { + (*this) *= 1/a; + } + + void operator *=(Real a) { + for (int i=0; i < dim(); i++) + dat[i] *= a; + } + + void operator -=(Vector v) { + assert(v.dim() == dim()); + for (int i=0; i < dim(); i++) + dat[i] -= v(i); + } + + Real &operator()(int i) { return dat[i]; } + Real operator()(int i) const { return dat[i]; } + Real elem(int i) { return dat[i]; } + Real operator *(Vector v) const { + Real ip=0; + assert(v.dim() == dim()); + for (int i=0; i < dim(); i++) + ip += dat[i] *v(i); + return ip; + } + Vector operator-() const; + Real norm() { + return sqrt(norm_sq() ); + } + Real norm_sq() { + return ((*this) * (*this)); + } + operator svec () { return dat; } + void print() const; + /// set to j-th element of unit-base + void set_unit(int j) ; +}; +/** + a vector. Storage is handled in svec, Vector only does the mathematics. + */ + +inline Vector +operator+(Vector a, Vector const &b) { + a += b; + return a; +} + +inline Vector +operator-(Vector a, Vector const &b) { + a -= b; + return a; +} + +inline Vector +operator*(Vector v, Real a) { + v *= a; + return v; +} + +inline Vector +operator*( Real a,Vector v) { + v *= a; + return v; +} + +inline Vector +operator/(Vector v,Real a) { + v *= 1/a; + return v; +} + +#endif