# HG changeset patch
# User Daniel Kraft <address@hidden>
# Date 1410265923 -7200
#      Tue Sep 09 14:32:03 2014 +0200
# Node ID be29bd5d63b2316988d1280e314b0443fe956f60
# Parent  afc1a8965664c891630714c00e11bc705a842f12
Rework internal structure of sparse_base_chol_rep

* liboctave/numeric/sparse-base-chol.h: Make sparse_base_chol_rep hold also
  cholmod_factor as a class member.  Add routines to handle it.
* liboctave/numeric/sparse-base-chol.cc: Split the code from
  sparse_base_chol_rep::init into multiple parts for initialising the
  cholmod_common, cleaning up, and converting the cholmod_factor
  to the final cholmod_sparse.

diff -r afc1a8965664 -r be29bd5d63b2 .hgsubstate
--- a/.hgsubstate	Mon Sep 08 08:09:45 2014 +0200
+++ b/.hgsubstate	Tue Sep 09 14:32:03 2014 +0200
@@ -1,1 +1,1 @@
-6d4e36653a40da6604507406f2a97e3e64bf9dbf gnulib-hg
+2539dbbdf52a445c223b4366edac5cde5530be64 gnulib-hg
diff -r afc1a8965664 -r be29bd5d63b2 liboctave/numeric/sparse-base-chol.cc
--- a/liboctave/numeric/sparse-base-chol.cc	Mon Sep 08 08:09:45 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.cc	Tue Sep 09 14:32:03 2014 +0200
@@ -34,6 +34,82 @@
 #include "MatrixType.h"
 
 #ifdef HAVE_CHOLMOD
+
+template <class chol_type, class chol_elt, class p_type>
+void
+sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::init_common
+  ()
+{
+  CHOLMOD_NAME(start) (&Common);
+  Common.prefer_zomplex = false;
+
+  double spu = octave_sparse_params::get_key ("spumoni");
+  if (spu == 0.)
+    {
+      Common.print = -1;
+      Common.print_function = 0;
+    }
+  else
+    {
+      Common.print = static_cast<int> (spu) + 2;
+      Common.print_function =&SparseCholPrint;
+    }
+
+  Common.error_handler = &SparseCholError;
+  Common.complex_divide = CHOLMOD_NAME(divcomplex);
+  Common.hypotenuse = CHOLMOD_NAME(hypot);
+
+  Common.final_asis = false;
+  Common.final_super = false;
+  Common.final_ll = true;
+  Common.final_pack = true;
+  Common.final_monotonic = true;
+  Common.final_resymbol = false;
+}
+
+template <class chol_type, class chol_elt, class p_type>
+void
+sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep
+  ::factor_to_sparse (bool natural)
+{
+  assert (Lfactor);
+
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  cond = CHOLMOD_NAME(rcond) (Lfactor, &Common);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  minor_p = Lfactor->minor;
+
+  clear_sparse ();
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, &Common);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  if (minor_p > 0 && minor_p < Lfactor->n)
+    {
+      size_t n1 = Lfactor->n + 1;
+      Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1,
+                                          sizeof(octave_idx_type),
+                                          Lsparse->p, &n1, &Common);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CHOLMOD_NAME(reallocate_sparse)
+        (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], Lsparse, &Common);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      Lsparse->ncol = minor_p;
+    }
+
+  drop_zeros (Lsparse);
+
+  if (! natural)
+    {
+      perms.resize (Lfactor->n);
+      for (octave_idx_type i = 0; i < Lfactor->n; i++)
+        perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i];
+    }
+
+  clear_factor ();
+}
+
 // Can't use CHOLMOD_NAME(drop)(0.0, S, cm). It doesn't treat complex matrices
 template <class chol_type, class chol_elt, class p_type>
 void
@@ -95,35 +171,6 @@
       return -1;
     }
 
-  cholmod_common *cm = &Common;
-
-  // Setup initial parameters
-  CHOLMOD_NAME(start) (cm);
-  cm->prefer_zomplex = false;
-
-  double spu = octave_sparse_params::get_key ("spumoni");
-  if (spu == 0.)
-    {
-      cm->print = -1;
-      cm->print_function = 0;
-    }
-  else
-    {
-      cm->print = static_cast<int> (spu) + 2;
-      cm->print_function =&SparseCholPrint;
-    }
-
-  cm->error_handler = &SparseCholError;
-  cm->complex_divide = CHOLMOD_NAME(divcomplex);
-  cm->hypotenuse = CHOLMOD_NAME(hypot);
-
-  cm->final_asis = false;
-  cm->final_super = false;
-  cm->final_ll = true;
-  cm->final_pack = true;
-  cm->final_monotonic = true;
-  cm->final_resymbol = false;
-
   cholmod_sparse A;
   cholmod_sparse *ac = &A;
   double dummy;
@@ -157,62 +204,22 @@
   // use natural ordering if no q output parameter
   if (natural)
     {
-      cm->nmethods = 1 ;
-      cm->method[0].ordering = CHOLMOD_NATURAL ;
-      cm->postorder = false ;
+      Common.nmethods = 1 ;
+      Common.method[0].ordering = CHOLMOD_NATURAL ;
+      Common.postorder = false ;
     }
 
-  cholmod_factor *Lfactor;
+  clear_factor ();
   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
-  CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
+  Lfactor = CHOLMOD_NAME(analyze) (ac, &Common);
+  CHOLMOD_NAME(factorize) (ac, Lfactor, &Common);
   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
-  is_pd = cm->status == CHOLMOD_OK;
-  info = (is_pd ? 0 : cm->status);
+  is_pd = Common.status == CHOLMOD_OK;
+  info = (is_pd ? 0 : Common.status);
 
   if (is_pd || force)
-    {
-      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-      cond = CHOLMOD_NAME(rcond) (Lfactor, cm);
-      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-      minor_p = Lfactor->minor;
-
-      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)
-        {
-          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];
-        }
-
-      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;
-    }
+    factor_to_sparse (natural);
 #else
   (*current_liboctave_error_handler)
     ("Missing CHOLMOD. Sparse cholesky factorization disabled");
diff -r afc1a8965664 -r be29bd5d63b2 liboctave/numeric/sparse-base-chol.h
--- a/liboctave/numeric/sparse-base-chol.h	Mon Sep 08 08:09:45 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.h	Tue Sep 09 14:32:03 2014 +0200
@@ -37,31 +37,44 @@
   {
   public:
     sparse_base_chol_rep (void)
-      : count (1), Lsparse (0), Common (), is_pd (false), minor_p (0),
-        perms (), cond (0)
-    { }
+      : count (1), Lsparse (0), Lfactor (0), Common (),
+        is_pd (false), minor_p (0), perms (), cond (0)
+    {
+      init_common ();
+    }
 
     sparse_base_chol_rep (const chol_type& a, bool natural, bool force)
-      : count (1), Lsparse (0), Common (), is_pd (false), minor_p (0),
-        perms (), cond (0)
+      : count (1), Lsparse (0), Lfactor (0), Common (),
+        is_pd (false), minor_p (0), perms (), cond (0)
     {
+      init_common ();
       init (a, natural, force);
     }
 
     sparse_base_chol_rep (const chol_type& a, octave_idx_type& info,
                           bool natural, bool force)
-      : count (1), Lsparse (0), Common (), is_pd (false), minor_p (0),
-        perms (), cond (0)
+      : count (1), Lsparse (0), Lfactor (0), Common (),
+        is_pd (false), minor_p (0), perms (), cond (0)
     {
+      init_common ();
       info = init (a, natural, force);
     }
 
     ~sparse_base_chol_rep (void)
     {
-      if (is_pd)
-        CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
+      clear_factor ();
+      clear_sparse ();
+
+      static char tmp[] = " ";
+      CHOLMOD_NAME (finish) (&Common);
+      CHOLMOD_NAME (print_common) (tmp, &Common);
     }
 
+    /* Convert Lfactor to Lsparse.  This is used after manipulating
+       the factor (using up/downdates, row inserts and so on) is done
+       and before extracting the result.  */
+    void factor_to_sparse (bool natural);
+
     cholmod_sparse * L (void) const { return Lsparse; }
 
     octave_idx_type P (void) const
@@ -82,6 +95,7 @@
 
   private:
     cholmod_sparse *Lsparse;
+    cholmod_factor *Lfactor;
 
     cholmod_common Common;
 
@@ -93,10 +107,30 @@
 
     double cond;
 
+    void init_common ();
+
     octave_idx_type init (const chol_type& a, bool natural, bool force);
 
     void drop_zeros (const cholmod_sparse* S);
 
+    void clear_factor ()
+    {
+      if (Lfactor)
+        {
+          CHOLMOD_NAME (free_factor) (&Lfactor, &Common);
+          Lfactor = 0;
+        }
+    }
+
+    void clear_sparse ()
+    {
+      if (Lsparse)
+        {
+          CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
+          Lsparse = 0;
+        }
+    }
+
     // No copying!
 
     sparse_base_chol_rep (const sparse_base_chol_rep&);
@@ -125,6 +159,8 @@
 
     ~sparse_base_chol_rep (void) { }
 
+    void factor_to_sparse (bool) { }
+
     octave_idx_type P (void) const { return 0; }
 
     ColumnVector perm (void) const { return perms + 1; }
@@ -204,6 +240,8 @@
     return *this;
   }
 
+  void factor_to_sparse (bool natural) { rep->factor_to_sparse (natural); }
+
   chol_type L (void) const;
 
   chol_type R (void) const { return L ().hermitian (); }