Mercurial > octave
diff liboctave/sparse-base-chol.cc @ 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 | daff886a8e2a |
children | 2176f2b4599e 72830070a17b |
line wrap: on
line diff
--- a/liboctave/sparse-base-chol.cc Tue Mar 25 15:56:41 2008 -0400 +++ b/liboctave/sparse-base-chol.cc Tue Mar 25 16:17:16 2008 -0400 @@ -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");