Mercurial > hg > octave-terminal
changeset 8402:2176f2b4599e
Fix sparse cholesky inversion
author | David Bateman <dbateman@free.fr> |
---|---|
date | Fri, 12 Dec 2008 23:18:20 +0100 |
parents | 712cfdc2e417 |
children | 87cca636a6c6 |
files | liboctave/ChangeLog liboctave/sparse-base-chol.cc src/ChangeLog src/DLD-FUNCTIONS/chol.cc |
diffstat | 4 files changed, 30 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -24,6 +24,11 @@ * CmplxAEPBAL.h, CmplxAEPBAL.cc: Rebase ComplexAEPBAL on base_aepbal. * fCmplxAEPBAL.h, fCmplxAEPBAL.cc: Rebase FloatComplexAEPBAL on base_aepbal. +2008-12-12 David Bateman <dbateman@free.fr> + + * sparse-base-chol.cc (inverse): Fix inversion based on cholesky + factorization. + 2008-12-08 Jaroslav Hajek <highegg@gmail.com> * idx-vector.cc (idx_vector::idx_vector_rep::idx_vector_rep (const
--- a/liboctave/sparse-base-chol.cc +++ b/liboctave/sparse-base-chol.cc @@ -277,15 +277,15 @@ double rcond2; octave_idx_type info; MatrixType mattype (MatrixType::Upper); - chol_type linv = L().transpose().inverse(mattype, info, rcond2, 1, 0); + chol_type linv = L().hermitian().inverse(mattype, info, rcond2, 1, 0); if (perms.length() == n) { p_type Qc = Q(); - retval = Qc * linv.transpose() * linv * Qc.transpose(); + retval = Qc * linv * linv.hermitian() * Qc.transpose(); } else - retval = linv.transpose() * linv; + retval = linv * linv.hermitian (); #endif return retval; }
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -46,6 +46,10 @@ * DLD-FUNCTIONS/balance.cc (Fbalance): Exploit the new AEPBAL functionality. +2008-12-12 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/chol.cc (Fcholinv): Add test. + 2008-12-08 Jaroslav Hajek <highegg@gmail.com> * xpow.cc ( xpow (const DiagMatrix& a, double b),
--- a/src/DLD-FUNCTIONS/chol.cc +++ b/src/DLD-FUNCTIONS/chol.cc @@ -469,6 +469,24 @@ return retval; } +/* + +%!test +%! A = [2,0.2;0.2,1]; +%! issymmetric(A) +%! min(eig(A)) +%! Ainv = inv(A); +%! Ainv1 = cholinv(A); +%! Ainv2 = inv(sparse(A)); +%! Ainv3 = cholinv(sparse(A)); +%! Ainv4 = spcholinv(sparse(A)); +%! assert (norm(Ainv-Ainv1),1e-10) +%! assert (norm(Ainv-Ainv2),1e-10) +%! assert (norm(Ainv-Ainv3),1e-10) +%! assert (norm(Ainv-Ainv4),1e-10) + +*/ + DEFUN_DLD (chol2inv, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} chol2inv (@var{u})\n\