# HG changeset patch
# User Daniel Kraft <d@domob.eu>
# Date 1410515468 -7200
#      Fri Sep 12 11:51:08 2014 +0200
# Node ID f1cad7943333b959a474af0c03eea40b84ad28e7
# Parent  b92bf24ae3b32ea5b0ddb8bcf0bcdaa1249d50db
Support multiple inserts at once in cholinsert

* libinterp/dldfcn/chol.cc: Support to insert multiple columns at once in
  cholinsert with sparse matrices.  Document and test this.
* liboctave/numeric/sparse-base-chol.cc: Implement multiple added identity
  rows/columns for factor.
* liboctave/numeric/sparse-base-chol.h: Ditto.

diff -r b92bf24ae3b3 -r f1cad7943333 libinterp/dldfcn/chol.cc
--- a/libinterp/dldfcn/chol.cc	Fri Sep 12 08:48:47 2014 +0200
+++ b/libinterp/dldfcn/chol.cc	Fri Sep 12 11:51:08 2014 +0200
@@ -857,6 +857,9 @@
 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\
+In this case, @var{j} can be a row vector of multiple insert positions\n\
+(sorted in ascending order) and @var{u} should contain as many columns.\n\
+Using this is more efficient that calling cholinsert multiple times.\n\
 @seealso{chol, cholupdate, choldelete, cholshift}\n\
 @end deftypefn")
 {
@@ -870,113 +873,158 @@
       return retval;
     }
 
-  octave_value argr = args(0);
-  octave_value argj = args(1);
-  octave_value argu = args(2);
+  const octave_value argr = args(0);
+  const octave_value argj = args(1);
+  const octave_value argu = args(2);
 
-  if (argr.is_numeric_type () && argu.is_numeric_type ()
-      && argj.is_real_scalar ())
+  if (!argr.is_numeric_type () || !argu.is_numeric_type ()
+      || !argj.is_numeric_type ())
     {
-      octave_idx_type n = argr.rows ();
-      octave_idx_type j = argj.scalar_value ();
+      print_usage ();
+      return retval;
+    }
 
-      if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
+  const octave_idx_type n = argr.rows ();
+  const octave_idx_type inserts = argu.columns ();
+  if (argr.columns () != n || argu.rows () != n + inserts)
+    {
+      error ("cholinsert: dimension mismatch between R and U");
+      return retval;
+    }
+
+  /* The sparse case allows multiple inserts to be performed at once
+     (since this is more efficient than calling cholinsert multiple times).
+     Handle this already here.  */
+  if (argr.is_sparse_type () || argu.is_sparse_type ())
+    {
+      const RowVector js = argj.row_vector_value () - 1;
+
+      if (js.length () == inserts)
         {
-          if (j > 0 && j <= n+1)
+          if (argr.is_real_type () && argu.is_real_type ())
             {
-              int err = 0;
-              if (argr.is_sparse_type () || argu.is_sparse_type ())
+              // real case
+              const SparseMatrix R = argr.sparse_matrix_value ();
+              const SparseMatrix u = argu.sparse_matrix_value ();
+
+              SparseCHOL fact;
+              fact.set (R, js);
+
+              if (inserts == 1)
+                fact.insert_sym (u, js(0));
+              else
+                for (octave_idx_type c = 0; c < js.length (); ++c)
+                  {
+                    const octave_idx_type offs = u.cidx (c);
+                    const octave_idx_type nnz = u.cidx (c + 1) - offs;
+                    SparseMatrix col (u.rows (), 1, nnz);
+                    col.xcidx (0) = 0;
+                    col.xcidx (1) = nnz;
+                    for (octave_idx_type i = 0; i < nnz; ++i)
+                      {
+                        col.xridx (i) = u.ridx (i + offs);
+                        col.xdata (i) = u.data (i + offs);
+                      }
+
+                    fact.insert_sym (col, js(c));
+                  }
+
+              fact.factor_to_sparse (true);
+              retval(0) = fact.R ();
+            }
+          else
+            {
+              // complex case: not supported by CHOLMOD
+              error ("cholinsert: sparse complex matrix not supported");
+            }
+        }
+      else
+        error ("cholinsert: dimension mismatch between J and U");
+
+      return retval;
+    }
+
+  if (inserts != 1)
+    {
+      error ("cholinsert: dimension mismatch for U");
+      return retval;
+    }
+
+  if (argj.is_real_scalar ())
+    {
+      const octave_idx_type j = argj.scalar_value ();
+
+      if (j > 0 && j <= n+1)
+        {
+          int err = 0;
+          if (argr.is_single_type () || argu.is_single_type ())
+            {
+              if (argr.is_real_type () && argu.is_real_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 ();
+                  // real case
+                  FloatMatrix R = argr.float_matrix_value ();
+                  FloatColumnVector u = argu.float_column_vector_value ();
 
-                      SparseCHOL fact;
-                      fact.set (R, j - 1);
-                      fact.insert_sym (u, j - 1);
+                  FloatCHOL fact;
+                  fact.set (R);
+                  err = 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 ())
-                    {
-                      // real case
-                      FloatMatrix R = argr.float_matrix_value ();
-                      FloatColumnVector u = argu.float_column_vector_value ();
-
-                      FloatCHOL fact;
-                      fact.set (R);
-                      err = fact.insert_sym (u, j-1);
-
-                      retval(0) = get_chol_r (fact);
-                    }
-                  else
-                    {
-                      // complex case
-                      FloatComplexMatrix R = argr.float_complex_matrix_value ();
-                      FloatComplexColumnVector u =
-                        argu.float_complex_column_vector_value ();
-
-                      FloatComplexCHOL fact;
-                      fact.set (R);
-                      err = fact.insert_sym (u, j-1);
-
-                      retval(0) = get_chol_r (fact);
-                    }
+                  retval(0) = get_chol_r (fact);
                 }
               else
                 {
-                  if (argr.is_real_type () && argu.is_real_type ())
-                    {
-                      // real case
-                      Matrix R = argr.matrix_value ();
-                      ColumnVector u = argu.column_vector_value ();
+                  // complex case
+                  FloatComplexMatrix R = argr.float_complex_matrix_value ();
+                  FloatComplexColumnVector u =
+                    argu.float_complex_column_vector_value ();
 
-                      CHOL fact;
-                      fact.set (R);
-                      err = fact.insert_sym (u, j-1);
+                  FloatComplexCHOL fact;
+                  fact.set (R);
+                  err = fact.insert_sym (u, j-1);
 
-                      retval(0) = get_chol_r (fact);
-                    }
-                  else
-                    {
-                      // complex case
-                      ComplexMatrix R = argr.complex_matrix_value ();
-                      ComplexColumnVector u =
-                        argu.complex_column_vector_value ();
-
-                      ComplexCHOL fact;
-                      fact.set (R);
-                      err = fact.insert_sym (u, j-1);
-
-                      retval(0) = get_chol_r (fact);
-                    }
+                  retval(0) = get_chol_r (fact);
                 }
-
-              if (nargout > 1)
-                retval(1) = err;
-              else if (err == 1)
-                error ("cholinsert: insertion violates positiveness");
-              else if (err == 2)
-                error ("cholinsert: singular matrix");
-              else if (err == 3)
-                error ("cholinsert: diagonal element must be real");
             }
           else
-            error ("cholinsert: index J out of range");
+            {
+              if (argr.is_real_type () && argu.is_real_type ())
+                {
+                  // real case
+                  Matrix R = argr.matrix_value ();
+                  ColumnVector u = argu.column_vector_value ();
+
+                  CHOL fact;
+                  fact.set (R);
+                  err = fact.insert_sym (u, j-1);
+
+                  retval(0) = get_chol_r (fact);
+                }
+              else
+                {
+                  // complex case
+                  ComplexMatrix R = argr.complex_matrix_value ();
+                  ComplexColumnVector u =
+                    argu.complex_column_vector_value ();
+
+                  ComplexCHOL fact;
+                  fact.set (R);
+                  err = fact.insert_sym (u, j-1);
+
+                  retval(0) = get_chol_r (fact);
+                }
+            }
+
+          if (nargout > 1)
+            retval(1) = err;
+          else if (err == 1)
+            error ("cholinsert: insertion violates positiveness");
+          else if (err == 2)
+            error ("cholinsert: singular matrix");
+          else if (err == 3)
+            error ("cholinsert: diagonal element must be real");
         }
       else
-        error ("cholinsert: dimension mismatch between R and U");
+        error ("cholinsert: index J out of range");
     }
   else
     print_usage ();
@@ -1122,13 +1170,19 @@
 
 %!testif HAVE_CHOLMOD
 %! a = sparse ([4, 2, 0; 2, 3, 0; 0, 0, 1]);
+%! r = chol (a);
+%!
 %! 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);
+%!
+%! ap = a(2, 2);
+%! u = a(:, [1, 3]);
+%! rp = chol (ap);
+%! rp = cholinsert (rp, [1, 3], u);
+%! assert (rp, r, 16 * eps);
 
 %!error <sparse complex matrix not supported>
 %! cholinsert (sparse (1i), 2, [1; 1]);
@@ -1136,6 +1190,12 @@
 %! cholinsert (sparse (1), 2, [1i; 1]);
 %!error <Invalid Cholesky factor>
 %! cholinsert (sparse ([1, 0; 1, 1]), 3, [1; 1; 1]);
+%!error <dimension mismatch for U>
+%! cholinsert (eye (2), [1, 1], ones (4, 2));
+%!error <insert positions should be increasing>
+%! cholinsert (sparse (eye (2)), [1, 1], ones (4, 2));
+%!error <dimension mismatch between J and U>
+%! cholinsert (sparse (eye (2)), 1, ones (4, 2));
 */
 
 DEFUN_DLD (choldelete, args, ,
diff -r b92bf24ae3b3 -r f1cad7943333 liboctave/numeric/sparse-base-chol.cc
--- a/liboctave/numeric/sparse-base-chol.cc	Fri Sep 12 08:48:47 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.cc	Fri Sep 12 11:51:08 2014 +0200
@@ -105,7 +105,7 @@
 template <class chol_type, class chol_elt, class p_type>
 void
 sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::set
-  (const chol_type& matL, octave_idx_type k)
+  (const chol_type& matL, const RowVector& ks)
 {
   const octave_idx_type nr = matL.rows ();
   const octave_idx_type nc = matL.cols ();
@@ -114,10 +114,25 @@
       (*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);
+  const octave_idx_type nk = ks.length ();
+  for (octave_idx_type k = 0; k < nk; ++k)
+    {
+      if (k + 1 < nk && ks(k) >= ks(k + 1))
+        {
+          (*current_liboctave_error_handler)
+            ("The insert positions should be increasing");
+          return;
+        }
+      if (ks(k) < 0 || ks(k) >= nr + nk)
+        {
+          (*current_liboctave_error_handler) ("Insert position out-of-bounds");
+          return;
+        }
+    }
+
+  /* Size of the result may be larger if we add rows.  */
+  const octave_idx_type n = nr + nk;
   
   clear_factor ();
 
@@ -138,7 +153,7 @@
 #endif
 
   const octave_idx_type nnzIn = matL.nnz ();
-  const octave_idx_type nnz = (k == -1 ? nnzIn : nnzIn + 1);
+  const octave_idx_type nnz = nnzIn + nk;
 
   Lfactor->nzmax = nnz;
   octave_idx_type *p = new octave_idx_type[n + 1];
@@ -148,11 +163,21 @@
   octave_idx_type *next = new octave_idx_type[n + 2];
   octave_idx_type *prev = new octave_idx_type[n + 2];
 
+  Lfactor->p = p;
+  Lfactor->i = i;
+  Lfactor->x = x;
+  Lfactor->nz = nz;
+  Lfactor->next = next;
+  Lfactor->prev = prev;
+
+  octave_idx_type inserted = 0;
   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;
+    {
+      p[c] = matL.cidx (c - inserted) + inserted;
+      if (inserted < nk && c == ks(inserted))
+        ++inserted;
+    }
+  assert (inserted == nk);
   for (octave_idx_type c = 0; c < n; ++c)
     nz[c] = p[c + 1] - p[c];
 
@@ -160,24 +185,31 @@
      (and column) if applicable.  */
   octave_idx_type in = 0;
   octave_idx_type out = 0;
+  inserted = 0;
   for (octave_idx_type c = 0; c < n; ++c)
-    if (c == k)
+    if (inserted < nk && c == ks(inserted))
       {
         i[out] = c;
         x[out] = 1.;
         ++out;
+        ++inserted;
       }
     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);
+      {
+        octave_idx_type offs = 0;
+        for (octave_idx_type el = 0; el < nz[c]; ++el)
+          {
+            i[out] = matL.ridx (in);
+            while (offs < nk && i[out] >= ks(offs))
+              ++offs;
+            i[out] += offs;
+
+            x[out] = matL.data (in);
+            ++in;
+            ++out;
+          }
+      }
+  assert (in == nnzIn && out == nnz && inserted == nk);
 
   /* Columns are monotonic.  Reflect this in the next/prev lists.  */
   for (octave_idx_type j = 0; j < n; ++j)
@@ -189,13 +221,6 @@
   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;
-
   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
   if (! CHOLMOD_NAME (check_factor) (Lfactor, &Common))
     (*current_liboctave_error_handler) ("Invalid Cholesky factor");
@@ -365,10 +390,10 @@
 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)
+  (const chol_type& matR, const RowVector& ks)
 {
 #ifdef HAVE_CHOLMOD
-  rep->set (matR.hermitian (), k);
+  rep->set (matR.hermitian (), ks);
 #else
   (*current_liboctave_error_handler)
     ("Missing CHOLMOD. Sparse cholesky factorization disabled");
diff -r b92bf24ae3b3 -r f1cad7943333 liboctave/numeric/sparse-base-chol.h
--- a/liboctave/numeric/sparse-base-chol.h	Fri Sep 12 08:48:47 2014 +0200
+++ b/liboctave/numeric/sparse-base-chol.h	Fri Sep 12 11:51:08 2014 +0200
@@ -72,9 +72,10 @@
 
     /* 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);
+       objects.  It should be in the L L' factorisation type.  Additionally,
+       insert identity rows/columns at the given indices.  They
+       should be in increasing order.  */
+    void set (const chol_type& matL, const RowVector& ks);
 
     octave_idx_type insert_sym (const chol_type& u, octave_idx_type k);
 
@@ -167,7 +168,7 @@
 
     ~sparse_base_chol_rep (void) { }
 
-    void set (const chol_type& mat, octave_idx_type k = -1) { }
+    void set (const chol_type&, const RowVector&) { }
 
     octave_idx_type insert_sym (const chol_type& u, octave_idx_type k) { }
 
@@ -255,7 +256,13 @@
   /* 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);
+  void
+  set (const chol_type& matL)
+  {
+    RowVector dummy;
+    set (matL, dummy);
+  }
+  void set (const chol_type& matL, const RowVector& ks);
 
   octave_idx_type
   insert_sym (const chol_type& u, octave_idx_type k)
