# HG changeset patch
# User Daniel Kraft
# 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
+%! cholinsert (sparse (1i), 2, [1; 1]);
+%!error
+%! cholinsert (sparse (1), 2, [1i; 1]);
+%!error
+%! 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
+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
void
sparse_base_chol::sparse_base_chol_rep::init_common
@@ -69,11 +104,136 @@
template
void
-sparse_base_chol::sparse_base_chol_rep
- ::factor_to_sparse (bool natural)
+sparse_base_chol::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
+octave_idx_type
+sparse_base_chol::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
+void
+sparse_base_chol::
+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
+void
+sparse_base_chol::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
chol_type
sparse_base_chol::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;