Mercurial > hg > octave-max

changeset 11716:017e13f0d056 release-3-0-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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:18:42 -0400
parents 86bbba911de8
children e5510f2d482a
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-23  David Bateman  <dbateman@free.fr>
 
 	* mx-ops: Definite binary operators for mixed integer array +
--- 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; }