# HG changeset patch
# User Daniel Kraft <address@hidden>
# Date 1410437803 -7200
#      Thu Sep 11 14:16:43 2014 +0200
# Node ID 27362c8204a8b87320e5bbea6f22a12157beb7c2
# Parent  be29bd5d63b2316988d1280e314b0443fe956f60
Support sparse matrices in cholinsert

* libinterp/dldfcn/chol.cc: Add support for sparse matrices to cholinsert.
  Update the docstring and add tests.
* liboctave/numeric/sparse-base-chol.h: Extend to support row additions.
* liboctave/numeric/sparse-base-chol.cc: Ditto.

diff -r be29bd5d63b2 -r 27362c8204a8 libinterp/dldfcn/chol.cc
--- a/libinterp/dldfcn/chol.cc	Tue Sep 09 14:32:03 2014 +0200
+++ b/libinterp/dldfcn/chol.cc	Thu Sep 11 14:16:43 2014 +0200
@@ -855,6 +855,8 @@
 @end itemize\n\
 \n\
 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
+\n\
+Sparse matrices are supported, but only if @var{R} and @var{u} are real.\n\
 @seealso{chol, cholupdate, choldelete, cholshift}\n\
 @end deftypefn")
 {
@@ -883,7 +885,28 @@
           if (j > 0 && j <= n+1)
             {
               int err = 0;
-              if (argr.is_single_type () || argu.is_single_type ())
+              if (argr.is_sparse_type () || argu.is_sparse_type ())
+                {
+                  if (argr.is_real_type () && argu.is_real_type ())
+                    {
+                      // real case
+                      const SparseMatrix R = argr.sparse_matrix_value ();
+                      const SparseMatrix u = argu.sparse_matrix_value ();
+
+                      SparseCHOL fact;
+                      fact.set (R, j - 1);
+                      fact.insert_sym (u, j - 1);
+
+                      fact.factor_to_sparse (true);
+                      retval(0) = fact.R ();
+                    }
+                  else
+                    {
+                      // complex case: not supported by CHOLMOD
+                      error ("cholinsert: sparse complex matrix not supported");
+                    }
+                }
+              else if (argr.is_single_type () || argu.is_single_type ())
                 {
                   if (argr.is_real_type () && argu.is_real_type ())
                     {
@@ -1096,6 +1119,23 @@
 %! assert (cca, ccau,   16*eps);
 %! assert (cca, ccal2', 16*eps);
 %! assert (cca, ccau2,  16*eps);
+
+%!testif HAVE_CHOLMOD
+%! a = sparse ([4, 2, 0; 2, 3, 0; 0, 0, 1]);
+%! ap = a([1, 3], [1, 3]);
+%! u = a(:, 2);
+%!
+%! r = chol (a);
+%! rp = chol (ap);
+%! rp = cholinsert (rp, 2, u);
+%! assert (rp, r, 16 * eps);
+
+%!error <sparse complex matrix not supported>
+%! cholinsert (sparse (1i), 2, [1; 1]);
+%!error <sparse complex matrix not supported>
+%! cholinsert (sparse (1), 2, [1i; 1]);
+%!error <Invalid Cholesky factor>
+%! cholinsert (sparse ([1, 0; 1, 1]), 3, [1; 1; 1]);
 */
 
 DEFUN_DLD (choldelete, args, ,
diff -r be29bd5d63b2 -r 27362c8204a8 liboctave/numeric/sparse-base-chol.cc
--- a/liboctave/numeric/sparse-base-chol.cc	Tue Sep 09 14:32:03 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.cc	Thu Sep 11 14:16:43 2014 +0200
@@ -35,6 +35,41 @@
 
 #ifdef HAVE_CHOLMOD
 
+/* Initialise a cholmod_sparse object from an Octave sparse matrix.
+   No memory is copied, everything is shared.  */
+template <class chol_type>
+static void
+octave_to_cholmod_sparse (cholmod_sparse& sp, const chol_type& mat)
+{
+  sp.nrow = mat.rows ();
+  sp.ncol = mat.cols ();
+
+  sp.p = mat.cidx ();
+  sp.i = mat.ridx ();
+  sp.nzmax = mat.nnz ();
+  sp.packed = true;
+  sp.sorted = true;
+  sp.nz = 0;
+#ifdef USE_64_BIT_IDX_T
+  sp.itype = CHOLMOD_LONG;
+#else
+  sp.itype = CHOLMOD_INT;
+#endif
+  sp.dtype = CHOLMOD_DOUBLE;
+  sp.stype = 1;
+#ifdef OCTAVE_CHOLMOD_TYPE
+  sp.xtype = OCTAVE_CHOLMOD_TYPE;
+#else
+  sp.xtype = CHOLMOD_REAL;
+#endif
+
+  static double dummy;
+  if (mat.rows () < 1)
+    sp.x = &dummy;
+  else
+    sp.x = mat.data ();
+}
+
 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
@@ -69,11 +104,136 @@
 
 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)
+sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::set
+  (const chol_type& matL, octave_idx_type k)
+{
+  const octave_idx_type nr = matL.rows ();
+  const octave_idx_type nc = matL.cols ();
+  if (nr != nc)
+    {
+      (*current_liboctave_error_handler) ("Cholesky factor should be square");
+      return;
+    }
+  assert (k == -1 || (k >= 0 && k <= nr));
+
+  /* Size of the result may be larger by one if we add an identity row.  */
+  const octave_idx_type n = (k == -1 ? nr : nr + 1);
+  
+  clear_factor ();
+
+  Lfactor = CHOLMOD_NAME (allocate_factor) (n, &Common);
+  Lfactor->is_ll = true;
+  Lfactor->is_super = false;
+  Lfactor->is_monotonic = true;
+#ifdef USE_64_BIT_IDX_T
+  Lfactor->itype = CHOLMOD_LONG;
+#else
+  Lfactor->itype = CHOLMOD_INT;
+#endif
+  Lfactor->dtype = CHOLMOD_DOUBLE;
+#ifdef OCTAVE_CHOLMOD_TYPE
+  Lfactor->xtype = OCTAVE_CHOLMOD_TYPE;
+#else
+  Lfactor->xtype = CHOLMOD_REAL;
+#endif
+
+  const octave_idx_type nnzIn = matL.nnz ();
+  const octave_idx_type nnz = (k == -1 ? nnzIn : nnzIn + 1);
+
+  Lfactor->nzmax = nnz;
+  octave_idx_type *p = new octave_idx_type[n + 1];
+  octave_idx_type *i = new octave_idx_type[nnz];
+  chol_elt *x = new chol_elt[nnz];
+  octave_idx_type *nz = new octave_idx_type[n];
+  octave_idx_type *next = new octave_idx_type[n + 2];
+  octave_idx_type *prev = new octave_idx_type[n + 2];
+
+  for (octave_idx_type c = 0; c <= n; ++c)
+    if (k == -1 || c <= k)
+      p[c] = matL.cidx (c);
+    else
+      p[c] = matL.cidx (c - 1) + 1;
+  for (octave_idx_type c = 0; c < n; ++c)
+    nz[c] = p[c + 1] - p[c];
+
+  /* Copy the data and take care of inserting the identity row
+     (and column) if applicable.  */
+  octave_idx_type in = 0;
+  octave_idx_type out = 0;
+  for (octave_idx_type c = 0; c < n; ++c)
+    if (c == k)
+      {
+        i[out] = c;
+        x[out] = 1.;
+        ++out;
+      }
+    else
+      for (octave_idx_type el = 0; el < nz[c]; ++el)
+        {
+          i[out] = matL.ridx (in);
+          if (k != -1 && i[out] >= k)
+            ++i[out];
+          x[out] = matL.data (in);
+          ++in;
+          ++out;
+        }
+  assert (in == nnzIn && out == nnz);
+
+  /* Columns are monotonic.  Reflect this in the next/prev lists.  */
+  for (octave_idx_type j = 0; j < n; ++j)
+    next[j] = j + 1;
+  next[n] = -1;
+  next[n + 1] = 0;
+  for (octave_idx_type j = 1; j <= n; ++j)
+    prev[j] = j - 1;
+  prev[0] = n + 1;
+  prev[n + 1] = -1;
+
+  Lfactor->p = p;
+  Lfactor->i = i;
+  Lfactor->x = x;
+  Lfactor->nz = nz;
+  Lfactor->next = next;
+  Lfactor->prev = prev;
+
+  if (! CHOLMOD_NAME (check_factor) (Lfactor, &Common))
+    (*current_liboctave_error_handler) ("Invalid Cholesky factor");
+}
+
+template <class chol_type, class chol_elt, class p_type>
+octave_idx_type
+sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::insert_sym
+  (const chol_type& u, octave_idx_type k)
 {
   assert (Lfactor);
 
+  cholmod_sparse uc;
+  octave_to_cholmod_sparse (uc, u);
+
+  volatile octave_idx_type res = 0;
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  res = CHOLMOD_NAME (rowadd) (k, &uc, Lfactor, &Common);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  return res;
+}
+
+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);
+
+  /* If we have an LDL' factor due to row-addition, make sure to convert
+     it back to LL'.  Otherwise, the sparse matrix constructed will represent
+     also an LDL' factor -- this is not what we want.  */
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  CHOLMOD_NAME (change_factor) (Lfactor->xtype, true, false, false, false,
+                                Lfactor, &Common);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  assert (Lfactor->is_ll);
+
   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
   cond = CHOLMOD_NAME(rcond) (Lfactor, &Common);
   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
@@ -171,35 +331,8 @@
       return -1;
     }
 
-  cholmod_sparse A;
-  cholmod_sparse *ac = &A;
-  double dummy;
-  ac->nrow = a_nr;
-  ac->ncol = a_nc;
-
-  ac->p = a.cidx ();
-  ac->i = a.ridx ();
-  ac->nzmax = a.nnz ();
-  ac->packed = true;
-  ac->sorted = true;
-  ac->nz = 0;
-#ifdef USE_64_BIT_IDX_T
-  ac->itype = CHOLMOD_LONG;
-#else
-  ac->itype = CHOLMOD_INT;
-#endif
-  ac->dtype = CHOLMOD_DOUBLE;
-  ac->stype = 1;
-#ifdef OCTAVE_CHOLMOD_TYPE
-  ac->xtype = OCTAVE_CHOLMOD_TYPE;
-#else
-  ac->xtype = CHOLMOD_REAL;
-#endif
-
-  if (a_nr < 1)
-    ac->x = &dummy;
-  else
-    ac->x = a.data ();
+  cholmod_sparse ac;
+  octave_to_cholmod_sparse (ac, a);
 
   // use natural ordering if no q output parameter
   if (natural)
@@ -211,8 +344,8 @@
 
   clear_factor ();
   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  Lfactor = CHOLMOD_NAME(analyze) (ac, &Common);
-  CHOLMOD_NAME(factorize) (ac, Lfactor, &Common);
+  Lfactor = CHOLMOD_NAME(analyze) (&ac, &Common);
+  CHOLMOD_NAME(factorize) (&ac, Lfactor, &Common);
   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 
   is_pd = Common.status == CHOLMOD_OK;
@@ -228,6 +361,19 @@
 }
 
 template <class chol_type, class chol_elt, class p_type>
+void
+sparse_base_chol<chol_type, chol_elt, p_type>::set
+  (const chol_type& matR, octave_idx_type k)
+{
+#ifdef HAVE_CHOLMOD
+  rep->set (matR.hermitian (), k);
+#else
+  (*current_liboctave_error_handler)
+    ("Missing CHOLMOD. Sparse cholesky factorization disabled");
+#endif
+}
+
+template <class chol_type, class chol_elt, class p_type>
 chol_type
 sparse_base_chol<chol_type, chol_elt, p_type>::L (void) const
 {
diff -r be29bd5d63b2 -r 27362c8204a8 liboctave/numeric/sparse-base-chol.h
--- a/liboctave/numeric/sparse-base-chol.h	Tue Sep 09 14:32:03 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.h	Thu Sep 11 14:16:43 2014 +0200
@@ -70,6 +70,14 @@
       CHOLMOD_NAME (print_common) (tmp, &Common);
     }
 
+    /* Set the internal factor from the given Octave matrix.  The matrix
+       should be lower-triangular, as per the structure of cholmod_factor
+       objects.  It should be in the L L' factorisation type.  If k is
+       given, insert an identity-row at that index in addition.  */
+    void set (const chol_type& matL, octave_idx_type k = -1);
+
+    octave_idx_type insert_sym (const chol_type& u, octave_idx_type k);
+
     /* 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.  */
@@ -159,6 +167,10 @@
 
     ~sparse_base_chol_rep (void) { }
 
+    void set (const chol_type& mat, octave_idx_type k = -1) { }
+
+    octave_idx_type insert_sym (const chol_type& u, octave_idx_type k) { }
+
     void factor_to_sparse (bool) { }
 
     octave_idx_type P (void) const { return 0; }
@@ -240,6 +252,17 @@
     return *this;
   }
 
+  /* Set the internal factor from the given Octave matrix.  The matrix
+     should be upper-triangular (as necessary for cholinsert and friends)
+     and in the L L' factorisation type.  */
+  void set (const chol_type& matR, octave_idx_type k = -1);
+
+  octave_idx_type
+  insert_sym (const chol_type& u, octave_idx_type k)
+  {
+    return rep->insert_sym (u, k);
+  }
+
   void factor_to_sparse (bool natural) { rep->factor_to_sparse (natural); }
 
   chol_type L (void) const;