Mercurial > hg > octave-lojdl
changeset 7637:2be056f03720
Fix fall back from sparse cholesky factorization to LU when matrix detected as not being positive definite
author | David Bateman <dbateman@free.fr> |
---|---|
date | Tue, 25 Mar 2008 16:17:16 -0400 |
parents | 99c410f7f0b0 |
children | 2df457529cfa |
files | liboctave/ChangeLog liboctave/sparse-base-chol.cc liboctave/sparse-base-chol.h |
diffstat | 3 files changed, 46 insertions(+), 30 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,11 @@ +2008-03-25 David Bateman <dbateman@free.fr> + + * sparse-base-chol.h (sparse_base_chol_rep::~sparse_base_chol_rep + (void)): Only free the factorization if it was created + * spase-base-chol.cc (sparse_base_chol_rep::init): Don't attempt + to factorize a matrix that has been flagged as not being positive + definite. + 2008-03-25 John W. Eaton <jwe@octave.org> * lo-mappers.cc (xtrunc): New function.
--- a/liboctave/sparse-base-chol.cc +++ b/liboctave/sparse-base-chol.cc @@ -165,45 +165,53 @@ BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; Lfactor = CHOLMOD_NAME(analyze) (ac, cm); CHOLMOD_NAME(factorize) (ac, Lfactor, cm); - cond = CHOLMOD_NAME(rcond) (Lfactor, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - minor_p = Lfactor->minor; is_pd = cm->status == CHOLMOD_OK; info = (is_pd ? 0 : cm->status); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - if (minor_p > 0 && minor_p < a_nr) + if (is_pd) { - size_t n1 = a_nr + 1; - Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1, sizeof(octave_idx_type), - Lsparse->p, &n1, cm); BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CHOLMOD_NAME(reallocate_sparse) - (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm); + cond = CHOLMOD_NAME(rcond) (Lfactor, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - Lsparse->ncol = minor_p; - } + + minor_p = Lfactor->minor; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - drop_zeros (Lsparse); + if (minor_p > 0 && minor_p < a_nr) + { + size_t n1 = a_nr + 1; + Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1, + sizeof(octave_idx_type), + Lsparse->p, &n1, cm); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(reallocate_sparse) + (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + Lsparse->ncol = minor_p; + } + + drop_zeros (Lsparse); - if (! natural) - { - perms.resize (a_nr); - for (octave_idx_type i = 0; i < a_nr; i++) - perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i]; + if (! natural) + { + perms.resize (a_nr); + for (octave_idx_type i = 0; i < a_nr; i++) + perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i]; + } + + static char tmp[] = " "; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CHOLMOD_NAME(free_factor) (&Lfactor, cm); + CHOLMOD_NAME(finish) (cm); + CHOLMOD_NAME(print_common) (tmp, cm); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; } - - static char tmp[] = " "; - - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CHOLMOD_NAME(free_factor) (&Lfactor, cm); - CHOLMOD_NAME(finish) (cm); - CHOLMOD_NAME(print_common) (tmp, cm); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; #else (*current_liboctave_error_handler) ("Missing CHOLMOD. Sparse cholesky factorization disabled");
--- a/liboctave/sparse-base-chol.h +++ b/liboctave/sparse-base-chol.h @@ -51,8 +51,8 @@ const bool natural) : count (1) { info = init (a, natural); } - ~sparse_base_chol_rep (void) - { CHOLMOD_NAME(free_sparse) (&Lsparse, &Common); } + ~sparse_base_chol_rep (void) + { if (is_pd) CHOLMOD_NAME(free_sparse) (&Lsparse, &Common); } cholmod_sparse * L (void) const { return Lsparse; }