# HG changeset patch
# User David Bateman <dbateman@free.fr>
# Date 1206476322 14400
# Node ID 017e13f0d05629a397dda3f661c0969741452650
# Parent  86bbba911de8d129e9a9145ed7fbc36db2796e1a
Fix fall back from sparse cholesky factorization to LU when matrix detected as not being positive definite

diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- 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 +
diff --git a/liboctave/sparse-base-chol.cc b/liboctave/sparse-base-chol.cc
--- 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");
diff --git a/liboctave/sparse-base-chol.h b/liboctave/sparse-base-chol.h
--- 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; }