X-Git-Url: https://git.donarmstrong.com/lilypond.git?a=blobdiff_plain;f=flower%2Fcholeski.cc;h=eae3c49cea5a999f756af105a364160419382ff4;hb=1a66290a98e7de8d6d41485b5b71a9f7e1fe35c7;hp=086759193fb57b4152ddbb3331168bd05d2d5c61;hpb=6cea332dd176e75dd6c2ab2d92505cadfdfcd99f;p=lilypond.git diff --git a/flower/choleski.cc b/flower/choleski.cc index 086759193f..eae3c49cea 100644 --- a/flower/choleski.cc +++ b/flower/choleski.cc @@ -13,93 +13,93 @@ const Real EPS = 1e-7; // so sue me. Hard coded //#define PARANOID void -Choleski_decomposition::full_matrix_solve(Vector &out, Vector const &rhs)const +Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs)const { int n= rhs.dim(); - assert(n == L.dim()); + assert (n == L.dim()); Vector y; - y.set_dim( n); - out.set_dim(n); + y.set_dim (n); + out.set_dim (n); // forward substitution for (int i=0; i < n; i++) { - Real sum(0.0); + 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); + sum += y (j) * L(i,j); + y (i) = (rhs (i) - sum)/L(i,i); } for (int i=0; i < n; i++) - y(i) /= D(i); + y (i) /= D(i); // backward subst - Vector &x(out); // using input as return val. + Vector &x (out); // using input as return val. for (int i=n-1; i >= 0; i--) { - Real sum(0.0); + 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); + sum += L(j,i)*x (j); + x (i) = (y (i) - sum)/L(i,i); } } void -Choleski_decomposition::band_matrix_solve(Vector &out, Vector const &rhs)const +Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs)const { int n= rhs.dim(); int b = L.band_i(); - assert(n == L.dim()); + assert (n == L.dim()); - out.set_dim(n); + out.set_dim (n); Vector y; - y.set_dim(n); + y.set_dim (n); // forward substitution for (int i=0; i < n; i++) { - Real sum(0.0); + Real sum (0.0); for (int j= 0 >? i - b; j < i; j++) - sum += y(j) * L(i,j); - y(i) = (rhs(i) - sum)/L(i,i); + sum += y (j) * L(i,j); + y (i) = (rhs (i) - sum)/L(i,i); } for (int i=0; i < n; i++) - y(i) /= D(i); + y (i) /= D(i); // backward subst - Vector &x(out); // using input as return val. + Vector &x (out); // using input as return val. for (int i=n-1; i >= 0; i--) { - Real sum(0.0); + Real sum (0.0); for (int j=i+1; j <= i + b&&j < n ; j++) - sum += L(j,i)*x(j); - x(i) = (y(i) - sum)/L(i,i); + sum += L(j,i)*x (j); + x (i) = (y (i) - sum)/L(i,i); } } void -Choleski_decomposition::solve(Vector &x, Vector const &rhs)const +Choleski_decomposition::solve (Vector &x, Vector const &rhs)const { if (L.band_b()) { - band_matrix_solve(x,rhs); + band_matrix_solve (x,rhs); } else - full_matrix_solve(x,rhs); + full_matrix_solve (x,rhs); } Vector -Choleski_decomposition::solve(Vector rhs)const +Choleski_decomposition::solve (Vector rhs)const { Vector r; - solve(r, rhs); + solve (r, rhs); return r; } void -Choleski_decomposition::full_matrix_decompose(Matrix const & P) +Choleski_decomposition::full_matrix_decompose (Matrix const & P) { int n = P.dim(); L.unit(); for (int k= 0; k < n; k++) { for (int j = 0; j < k; j++){ - Real sum(0.0); + 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); @@ -107,7 +107,7 @@ Choleski_decomposition::full_matrix_decompose(Matrix const & P) Real sum=0.0; for (int l=0; l < k; l++) - sum += sqr(L(k,l))*D(l); + sum += sqr (L(k,l))*D(l); Real d = P(k,k) - sum; D(k) = d; } @@ -115,7 +115,7 @@ Choleski_decomposition::full_matrix_decompose(Matrix const & P) } void -Choleski_decomposition::band_matrix_decompose(Matrix const &P) +Choleski_decomposition::band_matrix_decompose (Matrix const &P) { int n = P.dim(); int b = P.band_i(); @@ -123,7 +123,7 @@ Choleski_decomposition::band_matrix_decompose(Matrix const &P) for (int i= 0; i < n; i++) { for (int j = 0 >? i - b; j < i; j++){ - Real sum(0.0); + 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); @@ -131,12 +131,12 @@ Choleski_decomposition::band_matrix_decompose(Matrix const &P) Real sum=0.0; for (int l=0 >? i - b; l < i; l++) - sum += sqr(L(i,l))*D(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()); + assert ( L.band_i() == P.band_i ()); } @@ -146,20 +146,20 @@ Choleski_decomposition::band_matrix_decompose(Matrix const &P) Standard matrix algorithm. */ -Choleski_decomposition::Choleski_decomposition(Matrix const & P) - : L(P.dim()), D(P.dim()) +Choleski_decomposition::Choleski_decomposition (Matrix const & P) + : L(P.dim()), D(P.dim ()) { #ifdef PARANOID - assert((P-P.transposed()).norm()/P.norm() < EPS); + assert ((P-P.transposed()).norm ()/P.norm () < EPS); #endif if (P.band_b()) - band_matrix_decompose(P); + band_matrix_decompose (P); else - full_matrix_decompose(P); + full_matrix_decompose (P); #ifdef PARANOID - assert((original()-P).norm() / P.norm() < EPS); + assert ((original()-P).norm () / P.norm () < EPS); #endif } @@ -167,7 +167,7 @@ Matrix Choleski_decomposition::original() const { Matrix T(L.dim()); - T.set_diag(D); + T.set_diag (D); return L*T*L.transposed(); } @@ -175,20 +175,20 @@ Matrix Choleski_decomposition::inverse() const { int n=L.dim(); - Matrix invm(n); - Vector e_i(n); - Vector inv(n); + Matrix invm (n); + Vector e_i (n); + Vector inv (n); for (int i = 0; i < n; i++) { - e_i.set_unit(i); - solve(inv, e_i); + e_i.set_unit (i); + solve (inv, e_i); for (int j = 0 ; j