Mercurial > octave

changeset 21145:307096fb67e1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
revamp sparse Cholesky factorization classes * sparse-chol.h, sparse-chol.cc: Rename from sparse-base-chol.h and sparse-base-chol.cc, respectively. (class sparse_chol): Rename from sparse_base_chol. Incorporate code from SparseCmplxCHOL and SparsedbleCHOL classes into the sparse_chol template. Hide representation and HAVE_CHOLMOD macro from public interface. * sparse-chol-inst.cc: New file. * SparseCmplxCHOL.cc, SparseCmplxCHOL.h, SparsedbleCHOL.cc, SparsedbleCHOL.h: Delete. * chol.cc, symbfact.cc, CSparse.cc, dSparse.cc, eigs-base.cc: Change all uses of SparsedbleCHOL and SparseCmplxCHOL to use new sparse_chol template class. * liboctave/numeric/module.mk: Update.
author John W. Eaton <jwe@octave.org>
date Tue, 26 Jan 2016 16:30:54 -0500
parents 76e0ef020dae
children ea9c05014809
files libinterp/dldfcn/chol.cc libinterp/dldfcn/symbfact.cc liboctave/array/CSparse.cc liboctave/array/dSparse.cc liboctave/numeric/SparseCmplxCHOL.cc liboctave/numeric/SparseCmplxCHOL.h liboctave/numeric/SparsedbleCHOL.cc liboctave/numeric/SparsedbleCHOL.h liboctave/numeric/eigs-base.cc liboctave/numeric/module.mk liboctave/numeric/sparse-base-chol.cc liboctave/numeric/sparse-base-chol.h liboctave/numeric/sparse-chol-inst.cc liboctave/numeric/sparse-chol.cc liboctave/numeric/sparse-chol.h
diffstat 15 files changed, 688 insertions(+), 870 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/dldfcn/chol.cc	Wed Jan 27 14:15:17 2016 +0100
+++ b/libinterp/dldfcn/chol.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -31,8 +31,7 @@
 #include "dbleCHOL.h"
 #include "fCmplxCHOL.h"
 #include "floatCHOL.h"
-#include "SparseCmplxCHOL.h"
-#include "SparsedbleCHOL.h"
+#include "sparse-chol.h"
 #include "oct-spparms.h"
 #include "sparse-util.h"
 
@@ -194,7 +193,7 @@
         {
           SparseMatrix m = arg.sparse_matrix_value ();
 
-          SparseCHOL fact (m, info, natural, force);
+          sparse_chol<SparseMatrix> fact (m, info, natural, force);
 
           if (nargout == 3)
             {
@@ -219,7 +218,7 @@
         {
           SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
 
-          SparseComplexCHOL fact (m, info, natural, force);
+          sparse_chol<SparseComplexMatrix> fact (m, info, natural, force);
 
           if (nargout == 3)
             {
@@ -358,7 +357,7 @@
             {
               SparseMatrix m = arg.sparse_matrix_value ();
 
-              SparseCHOL chol (m, info);
+              sparse_chol<SparseMatrix> chol (m, info);
 
               if (info == 0)
                 retval = chol.inverse ();
@@ -369,7 +368,7 @@
             {
               SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
 
-              SparseComplexCHOL chol (m, info);
+              sparse_chol<SparseComplexMatrix> chol (m, info);
 
               if (info == 0)
                 retval = chol.inverse ();
--- a/libinterp/dldfcn/symbfact.cc	Wed Jan 27 14:15:17 2016 +0100
+++ b/libinterp/dldfcn/symbfact.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -25,8 +25,7 @@
 #include <config.h>
 #endif
 
-#include "SparseCmplxCHOL.h"
-#include "SparsedbleCHOL.h"
+#include "sparse-chol.h"
 #include "oct-spparms.h"
 #include "sparse-util.h"
 #include "oct-locbuf.h"
--- a/liboctave/array/CSparse.cc	Wed Jan 27 14:15:17 2016 +0100
+++ b/liboctave/array/CSparse.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -54,7 +54,7 @@
 #include "SparseCmplxLU.h"
 #include "oct-sparse.h"
 #include "sparse-util.h"
-#include "SparseCmplxCHOL.h"
+#include "sparse-chol.h"
 #include "SparseCmplxQR.h"
 
 #include "Sparse-op-defs.h"
@@ -1079,7 +1079,7 @@
       if (mattype.is_hermitian ())
         {
           MatrixType tmp_typ (MatrixType::Upper);
-          SparseComplexCHOL fact (*this, info, false);
+          sparse_chol<SparseComplexMatrix> fact (*this, info, false);
           rcond = fact.rcond ();
           if (info == 0)
             {
--- a/liboctave/array/dSparse.cc	Wed Jan 27 14:15:17 2016 +0100
+++ b/liboctave/array/dSparse.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -49,7 +49,7 @@
 #include "MatrixType.h"
 #include "oct-sparse.h"
 #include "sparse-util.h"
-#include "SparsedbleCHOL.h"
+#include "sparse-chol.h"
 #include "SparseQR.h"
 
 #include "Sparse-op-defs.h"
@@ -1169,7 +1169,7 @@
       if (mattype.is_hermitian ())
         {
           MatrixType tmp_typ (MatrixType::Upper);
-          SparseCHOL fact (*this, info, false);
+          sparse_chol<SparseMatrix> fact (*this, info, false);
           rcond = fact.rcond ();
           if (info == 0)
             {
--- a/liboctave/numeric/SparseCmplxCHOL.cc	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,67 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "SparseCmplxCHOL.h"
-
-// Instantiate the base CHOL class for the type we need
-#define OCTAVE_CHOLMOD_TYPE CHOLMOD_COMPLEX
-#include "sparse-base-chol.h"
-#include "sparse-base-chol.cc"
-template class sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix>;
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-SparseComplexMatrix
-chol2inv (const SparseComplexMatrix& r)
-{
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("U must be a square matrix");
-
-  SparseComplexMatrix retval;
-  MatrixType mattype (r);
-  int typ = mattype.type (false);
-  double rcond;
-  octave_idx_type info;
-  SparseComplexMatrix rinv;
-
-  if (typ == MatrixType::Upper)
-    {
-      rinv = r.inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
-    }
-  else if (typ == MatrixType::Lower)
-    {
-      rinv = r.transpose ().inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
-    }
-  else
-    (*current_liboctave_error_handler) ("U must be a triangular matrix");
-
-  return retval;
-}
--- a/liboctave/numeric/SparseCmplxCHOL.h	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_SparseCmplxCHOL_h)
-#define octave_SparseCmplxCHOL_h 1
-
-#include "sparse-base-chol.h"
-#include "dSparse.h"
-#include "CSparse.h"
-
-class
-OCTAVE_API
-SparseComplexCHOL
-  : public sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix>
-{
-public:
-
-  SparseComplexCHOL (void)
-    : sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> () { }
-
-  SparseComplexCHOL (const SparseComplexMatrix& a, bool natural = true,
-                     bool force = false)
-    : sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>
-       (a, natural, force) { }
-
-  SparseComplexCHOL (const SparseComplexMatrix& a, octave_idx_type& info,
-                     bool natural = true, bool force = false)
-    : sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>
-       (a, info, natural, force) { }
-
-  SparseComplexCHOL (const SparseComplexCHOL& a)
-    : sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix> (a) { }
-
-  ~SparseComplexCHOL (void) { }
-
-  SparseComplexCHOL& operator = (const SparseComplexCHOL& a)
-  {
-    if (this != &a)
-      sparse_base_chol <SparseComplexMatrix, Complex, SparseMatrix> ::
-      operator = (a);
-
-    return *this;
-  }
-
-  SparseComplexMatrix chol_matrix (void) const { return R (); }
-
-  SparseComplexMatrix L (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>:: L ();
-  }
-
-  SparseComplexMatrix R (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>:: R ();
-  }
-
-  octave_idx_type P (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>:: P ();
-  }
-
-  ColumnVector perm (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex,
-                            SparseMatrix>:: perm ();
-  }
-
-  SparseMatrix Q (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex, SparseMatrix>:: Q ();
-  }
-
-  double rcond (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex,
-                            SparseMatrix>:: rcond ();
-  }
-
-  // Compute the inverse of a matrix using the Cholesky factorization.
-  SparseComplexMatrix inverse (void) const
-  {
-    return sparse_base_chol<SparseComplexMatrix, Complex,
-                            SparseMatrix>:: inverse ();
-  }
-};
-
-SparseComplexMatrix OCTAVE_API chol2inv (const SparseComplexMatrix& r);
-
-#endif
--- a/liboctave/numeric/SparsedbleCHOL.cc	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,67 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "SparsedbleCHOL.h"
-
-// Instantiate the base CHOL class for the type we need
-#define OCTAVE_CHOLMOD_TYPE CHOLMOD_REAL
-#include "sparse-base-chol.h"
-#include "sparse-base-chol.cc"
-template class sparse_base_chol <SparseMatrix, double, SparseMatrix>;
-
-// Compute the inverse of a matrix using the Cholesky factorization.
-SparseMatrix
-chol2inv (const SparseMatrix& r)
-{
-  octave_idx_type r_nr = r.rows ();
-  octave_idx_type r_nc = r.cols ();
-  SparseMatrix retval;
-
-  if (r_nr != r_nc)
-    (*current_liboctave_error_handler) ("U must be a square matrix");
-
-  MatrixType mattype (r);
-  int typ = mattype.type (false);
-  double rcond;
-  octave_idx_type info;
-  SparseMatrix rinv;
-
-  if (typ == MatrixType::Upper)
-    {
-      rinv = r.inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
-    }
-  else if (typ == MatrixType::Lower)
-    {
-      rinv = r.transpose ().inverse (mattype, info, rcond, true, false);
-      retval = rinv.transpose () * rinv;
-    }
-  else
-    (*current_liboctave_error_handler) ("U must be a triangular matrix");
-
-  return retval;
-}
--- a/liboctave/numeric/SparsedbleCHOL.h	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,91 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_SparsedbleCHOL_h)
-#define octave_SparsedbleCHOL_h 1
-
-#include "sparse-base-chol.h"
-#include "dSparse.h"
-
-class
-OCTAVE_API
-SparseCHOL : public sparse_base_chol <SparseMatrix, double, SparseMatrix>
-{
-public:
-
-  SparseCHOL (void) : sparse_base_chol<SparseMatrix, double, SparseMatrix> ()
-  { }
-
-  SparseCHOL (const SparseMatrix& a, bool natural = true, bool force = false)
-    : sparse_base_chol<SparseMatrix, double, SparseMatrix> (a, natural, force)
-  { }
-
-  SparseCHOL (const SparseMatrix& a, octave_idx_type& info,
-              bool natural = false, bool force = false)
-    : sparse_base_chol<SparseMatrix, double, SparseMatrix> (a, info, natural,
-                                                            force)
-  { }
-
-  SparseCHOL (const SparseCHOL& a) :
-    sparse_base_chol<SparseMatrix, double, SparseMatrix> (a) { }
-
-  ~SparseCHOL (void) { }
-
-  SparseCHOL& operator = (const SparseCHOL& a)
-  {
-    if (this != &a)
-      sparse_base_chol <SparseMatrix, double, SparseMatrix> :: operator = (a);
-
-    return *this;
-  }
-
-  SparseMatrix chol_matrix (void) const { return R (); }
-
-  SparseMatrix L (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: L (); }
-
-  SparseMatrix R (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: R (); }
-
-  octave_idx_type P (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: P (); }
-
-  ColumnVector perm (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: perm (); }
-
-  SparseMatrix Q (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: Q (); }
-
-  double rcond (void) const
-  { return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: rcond (); }
-
-  // Compute the inverse of a matrix using the Cholesky factorization.
-  SparseMatrix inverse (void) const
-  {
-    return sparse_base_chol<SparseMatrix, double, SparseMatrix>:: inverse ();
-  }
-};
-
-SparseMatrix OCTAVE_API chol2inv (const SparseMatrix& r);
-
-#endif
--- a/liboctave/numeric/eigs-base.cc	Wed Jan 27 14:15:17 2016 +0100
+++ b/liboctave/numeric/eigs-base.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -37,8 +37,7 @@
 #include "dSparse.h"
 #include "CSparse.h"
 #include "MatrixType.h"
-#include "SparsedbleCHOL.h"
-#include "SparseCmplxCHOL.h"
+#include "sparse-chol.h"
 #include "oct-rand.h"
 #include "dbleCHOL.h"
 #include "CmplxCHOL.h"
@@ -370,7 +369,7 @@
 make_cholb (SparseMatrix& b, SparseMatrix& bt, ColumnVector& permB)
 {
   octave_idx_type info;
-  SparseCHOL fact (b, info, false);
+  sparse_chol<SparseMatrix> fact (b, info, false);
 
   if (fact.P () != 0)
     return false;
@@ -408,7 +407,7 @@
             ColumnVector& permB)
 {
   octave_idx_type info;
-  SparseComplexCHOL fact (b, info, false);
+  sparse_chol<SparseComplexMatrix> fact (b, info, false);
 
   if (fact.P () != 0)
     return false;
--- a/liboctave/numeric/module.mk	Wed Jan 27 14:15:17 2016 +0100
+++ b/liboctave/numeric/module.mk	Tue Jan 26 16:30:54 2016 -0500
@@ -79,12 +79,10 @@
   liboctave/numeric/randgamma.h \
   liboctave/numeric/randmtzig.h \
   liboctave/numeric/randpoisson.h \
-  liboctave/numeric/sparse-base-chol.h \
+  liboctave/numeric/sparse-chol.h \
   liboctave/numeric/sparse-base-lu.h \
-  liboctave/numeric/SparseCmplxCHOL.h \
   liboctave/numeric/SparseCmplxLU.h \
   liboctave/numeric/SparseCmplxQR.h \
-  liboctave/numeric/SparsedbleCHOL.h \
   liboctave/numeric/SparsedbleLU.h \
   liboctave/numeric/SparseQR.h
 
@@ -146,12 +144,11 @@
   liboctave/numeric/oct-spparms.cc \
   liboctave/numeric/ODES.cc \
   liboctave/numeric/Quad.cc \
-  liboctave/numeric/SparseCmplxCHOL.cc \
   liboctave/numeric/SparseCmplxLU.cc \
   liboctave/numeric/SparseCmplxQR.cc \
-  liboctave/numeric/SparsedbleCHOL.cc \
   liboctave/numeric/SparsedbleLU.cc \
   liboctave/numeric/SparseQR.cc \
+  liboctave/numeric/sparse-chol-inst.cc \
   $(NUMERIC_C_SRC)
 
 LIBOCTAVE_TEMPLATE_SRC += \
@@ -159,7 +156,7 @@
   liboctave/numeric/base-qr.cc \
   liboctave/numeric/bsxfun-defs.cc \
   liboctave/numeric/eigs-base.cc \
-  liboctave/numeric/sparse-base-chol.cc \
+  liboctave/numeric/sparse-chol.cc \
   liboctave/numeric/sparse-base-lu.cc \
   liboctave/numeric/sparse-dmsolve.cc
 
--- a/liboctave/numeric/sparse-base-chol.cc	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,289 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "sparse-base-chol.h"
-#include "sparse-util.h"
-#include "lo-error.h"
-#include "oct-sparse.h"
-#include "oct-spparms.h"
-#include "quit.h"
-#include "MatrixType.h"
-
-#ifdef HAVE_CHOLMOD
-// Can't use CHOLMOD_NAME(drop)(0.0, S, cm). It doesn't treat complex matrices
-template <typename chol_type, typename chol_elt, typename p_type>
-void
-sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::drop_zeros
-  (const cholmod_sparse* S)
-{
-  chol_elt sik;
-  octave_idx_type *Sp, *Si;
-  chol_elt *Sx;
-  octave_idx_type pdest, k, ncol, p, pend;
-
-  if (! S)
-    return;
-
-  Sp = static_cast<octave_idx_type *>(S->p);
-  Si = static_cast<octave_idx_type *>(S->i);
-  Sx = static_cast<chol_elt *>(S->x);
-  pdest = 0;
-  ncol = S->ncol;
-
-  for (k = 0; k < ncol; k++)
-    {
-      p = Sp[k];
-      pend = Sp[k+1];
-      Sp[k] = pdest;
-      for (; p < pend; p++)
-        {
-          sik = Sx[p];
-          if (CHOLMOD_IS_NONZERO (sik))
-            {
-              if (p != pdest)
-                {
-                  Si[pdest] = Si[p];
-                  Sx[pdest] = sik;
-                }
-              pdest++;
-            }
-        }
-    }
-  Sp[ncol] = pdest;
-}
-#endif
-
-template <typename chol_type, typename chol_elt, typename p_type>
-octave_idx_type
-sparse_base_chol<chol_type, chol_elt, p_type>::sparse_base_chol_rep::init
-  (const chol_type& a, bool natural, bool force)
-{
-  volatile octave_idx_type info = 0;
-
-#ifdef HAVE_CHOLMOD
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("SparseCHOL requires square matrix");
-
-  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;
-      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
-    }
-  else
-    {
-      cm->print = static_cast<int> (spu) + 2;
-      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
-    }
-
-  cm->error_handler = &SparseCholError;
-  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
-  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, 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;
-  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;
-#if defined (ENABLE_64)
-  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 ();
-
-  // use natural ordering if no q output parameter
-  if (natural)
-    {
-      cm->nmethods = 1 ;
-      cm->method[0].ordering = CHOLMOD_NATURAL ;
-      cm->postorder = false ;
-    }
-
-  cholmod_factor *Lfactor;
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
-  CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-  is_pd = cm->status == CHOLMOD_OK;
-  info = (is_pd ? 0 : cm->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;
-    }
-
-  return info;
-
-#else
-  (*current_liboctave_error_handler)
-    ("support for CHOLMOD was unavailable or disabled when liboctave was built");
-#endif
-}
-
-template <typename chol_type, typename chol_elt, typename p_type>
-chol_type
-sparse_base_chol<chol_type, chol_elt, p_type>::L (void) const
-{
-#ifdef HAVE_CHOLMOD
-  cholmod_sparse *m = rep->L ();
-  octave_idx_type nc = m->ncol;
-  octave_idx_type nnz = m->nzmax;
-  chol_type ret (m->nrow, nc, nnz);
-  for (octave_idx_type j = 0; j < nc+1; j++)
-    ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j];
-  for (octave_idx_type i = 0; i < nnz; i++)
-    {
-      ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i];
-      ret.xdata (i) = static_cast<chol_elt *>(m->x)[i];
-    }
-  return ret;
-#else
-  return chol_type ();
-#endif
-}
-
-template <typename chol_type, typename chol_elt, typename p_type>
-p_type
-sparse_base_chol<chol_type, chol_elt, p_type>::
-sparse_base_chol_rep::Q (void) const
-{
-#ifdef HAVE_CHOLMOD
-  octave_idx_type n = Lsparse->nrow;
-  p_type p (n, n, n);
-
-  for (octave_idx_type i = 0; i < n; i++)
-    {
-      p.xcidx (i) = i;
-      p.xridx (i) = static_cast<octave_idx_type>(perms (i));
-      p.xdata (i) = 1;
-    }
-  p.xcidx (n) = n;
-
-  return p;
-#else
-  return p_type ();
-#endif
-}
-
-template <typename chol_type, typename chol_elt, typename p_type>
-chol_type
-sparse_base_chol<chol_type, chol_elt, p_type>::inverse (void) const
-{
-  chol_type retval;
-#ifdef HAVE_CHOLMOD
-  cholmod_sparse *m = rep->L ();
-  octave_idx_type n = m->ncol;
-  ColumnVector perms = rep->perm ();
-  double rcond2;
-  octave_idx_type info;
-  MatrixType mattype (MatrixType::Upper);
-  chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
-
-  if (perms.numel () == n)
-    {
-      p_type Qc = Q ();
-      retval = Qc * linv * linv.hermitian () * Qc.transpose ();
-    }
-  else
-    retval = linv * linv.hermitian ();
-#endif
-  return retval;
-}
--- a/liboctave/numeric/sparse-base-chol.h	Wed Jan 27 14:15:17 2016 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,225 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-Copyright (C) 1998-2005 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_sparse_base_chol_h)
-#define octave_sparse_base_chol_h 1
-
-#include "oct-sparse.h"
-#include "dColVector.h"
-
-template <typename chol_type, typename chol_elt, typename p_type>
-class
-sparse_base_chol
-{
-protected:
-#ifdef HAVE_CHOLMOD
-  class sparse_base_chol_rep
-  {
-  public:
-    sparse_base_chol_rep (void)
-      : count (1), Lsparse (0), Common (), is_pd (false), minor_p (0),
-        perms (), cond (0)
-    { }
-
-    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)
-    {
-      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)
-    {
-      info = init (a, natural, force);
-    }
-
-    ~sparse_base_chol_rep (void)
-    {
-      if (is_pd)
-        CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
-    }
-
-    cholmod_sparse * L (void) const { return Lsparse; }
-
-    octave_idx_type P (void) const
-    {
-      return (minor_p == static_cast<octave_idx_type>(Lsparse->ncol) ?
-              0 : minor_p + 1);
-    }
-
-    ColumnVector perm (void) const { return perms + 1; }
-
-    p_type Q (void) const;
-
-    bool is_positive_definite (void) const { return is_pd; }
-
-    double rcond (void) const { return cond; }
-
-    octave_refcount<int> count;
-
-  private:
-    cholmod_sparse *Lsparse;
-
-    cholmod_common Common;
-
-    bool is_pd;
-
-    octave_idx_type minor_p;
-
-    ColumnVector perms;
-
-    double cond;
-
-    octave_idx_type init (const chol_type& a, bool natural, bool force);
-
-    void drop_zeros (const cholmod_sparse* S);
-
-    // No copying!
-
-    sparse_base_chol_rep (const sparse_base_chol_rep&);
-
-    sparse_base_chol_rep& operator = (const sparse_base_chol_rep&);
-  };
-#else
-  class sparse_base_chol_rep
-  {
-  public:
-    sparse_base_chol_rep (void)
-      : count (1), is_pd (false), minor_p (0), perms (), cond (0) { }
-
-    sparse_base_chol_rep (const chol_type& a, bool natural, bool force)
-      : count (1), is_pd (false), minor_p (0), perms (), cond (0)
-    {
-      init (a, natural, force);
-    }
-
-    sparse_base_chol_rep (const chol_type& a, octave_idx_type& info,
-                          bool natural, bool force)
-      : count (1), is_pd (false), minor_p (0), perms (), cond (0)
-    {
-      info = init (a, natural, force);
-    }
-
-    ~sparse_base_chol_rep (void) { }
-
-    octave_idx_type P (void) const { return 0; }
-
-    ColumnVector perm (void) const { return perms + 1; }
-
-    p_type Q (void) const;
-
-    bool is_positive_definite (void) const { return is_pd; }
-
-    double rcond (void) const { return cond; }
-
-    octave_refcount<int> count;
-
-  private:
-    bool is_pd;
-
-    octave_idx_type minor_p;
-
-    ColumnVector perms;
-
-    double cond;
-
-    octave_idx_type init (const chol_type& a, bool natural, bool force);
-
-    // No copying!
-
-    sparse_base_chol_rep (const sparse_base_chol_rep&);
-
-    sparse_base_chol_rep& operator = (const sparse_base_chol_rep&);
-  };
-#endif
-
-private:
-  sparse_base_chol_rep *rep;
-
-public:
-
-  sparse_base_chol (void)
-    : rep (new typename
-           sparse_base_chol<chol_type, chol_elt, p_type>
-           ::sparse_base_chol_rep ())
-  { }
-
-  sparse_base_chol (const chol_type& a, bool natural, bool force)
-    : rep (new typename
-           sparse_base_chol<chol_type, chol_elt, p_type>
-           ::sparse_base_chol_rep (a, natural, force))
-  { }
-
-  sparse_base_chol (const chol_type& a, octave_idx_type& info,
-                    bool natural, bool force)
-    : rep (new typename
-           sparse_base_chol<chol_type, chol_elt, p_type>
-           ::sparse_base_chol_rep (a, info, natural, force))
-  { }
-
-  sparse_base_chol (const sparse_base_chol<chol_type, chol_elt, p_type>& a)
-    : rep (a.rep)
-  { rep->count++; }
-
-  virtual ~sparse_base_chol (void)
-  {
-    if (--rep->count == 0)
-      delete rep;
-  }
-
-  sparse_base_chol& operator = (const sparse_base_chol& a)
-  {
-    if (this != &a)
-      {
-        if (--rep->count == 0)
-          delete rep;
-
-        rep = a.rep;
-        rep->count++;
-      }
-
-    return *this;
-  }
-
-  chol_type L (void) const;
-
-  chol_type R (void) const { return L ().hermitian (); }
-
-  octave_idx_type P (void) const { return rep->P (); }
-
-  ColumnVector perm (void) const { return rep->perm (); }
-
-  p_type Q (void) const { return rep->Q (); }
-
-  bool is_positive_definite (void) const
-  { return rep->is_positive_definite (); }
-
-  double rcond (void) const { return rep->rcond (); }
-
-  chol_type inverse (void) const;
-};
-
-#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-chol-inst.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -0,0 +1,38 @@
+/*
+
+Copyright (C) 2016 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "sparse-chol.h"
+#include "sparse-chol.cc"
+
+template class sparse_chol<SparseMatrix>;
+
+template class sparse_chol<SparseComplexMatrix>;
+
+template SparseMatrix
+chol2inv<SparseMatrix> (const SparseMatrix& r);
+
+template SparseComplexMatrix
+chol2inv<SparseComplexMatrix> (const SparseComplexMatrix& r);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-chol.cc	Tue Jan 26 16:30:54 2016 -0500
@@ -0,0 +1,542 @@
+/*
+
+Copyright (C) 2005-2015 David Bateman
+Copyright (C) 1998-2005 Andy Adler
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "sparse-chol.h"
+#include "sparse-util.h"
+#include "lo-error.h"
+#include "oct-sparse.h"
+#include "oct-spparms.h"
+#include "quit.h"
+#include "MatrixType.h"
+
+template <typename chol_type>
+class sparse_chol<chol_type>::sparse_chol_rep
+{
+public:
+
+  sparse_chol_rep (void)
+    : count (1), is_pd (false), minor_p (0), perms (), cond (0),
+#ifdef HAVE_CHOLMOD
+      Lsparse (0), Common ()
+#endif
+  { }
+
+  sparse_chol_rep (const chol_type& a, bool natural, bool force)
+    : count (1), is_pd (false), minor_p (0), perms (), cond (0),
+#ifdef HAVE_CHOLMOD
+      Lsparse (0), Common ()
+#endif
+  {
+    init (a, natural, force);
+  }
+
+  sparse_chol_rep (const chol_type& a, octave_idx_type& info,
+                   bool natural, bool force)
+    : count (1), is_pd (false), minor_p (0), perms (), cond (0),
+#ifdef HAVE_CHOLMOD
+      Lsparse (0), Common ()
+#endif
+  {
+    info = init (a, natural, force);
+  }
+
+  ~sparse_chol_rep (void)
+  {
+#ifdef HAVE_CHOLMOD
+    if (is_pd)
+      CHOLMOD_NAME (free_sparse) (&Lsparse, &Common);
+#endif
+  }
+
+  cholmod_sparse *L (void) const
+  {
+#ifdef HAVE_CHOLMOD
+    return Lsparse;
+#else
+    return 0;
+#endif
+  }
+
+  octave_idx_type P (void) const
+  {
+#ifdef HAVE_CHOLMOD
+    return (minor_p == static_cast<octave_idx_type>(Lsparse->ncol) ?
+            0 : minor_p + 1);
+#else
+    return 0;
+#endif
+  }
+
+  ColumnVector perm (void) const { return perms + 1; }
+
+  SparseMatrix Q (void) const;
+
+  bool is_positive_definite (void) const { return is_pd; }
+
+  double rcond (void) const { return cond; }
+
+  octave_refcount<int> count;
+
+private:
+
+  bool is_pd;
+
+  octave_idx_type minor_p;
+
+  ColumnVector perms;
+
+  double cond;
+
+#ifdef HAVE_CHOLMOD
+  cholmod_sparse *Lsparse;
+
+  cholmod_common Common;
+
+  void drop_zeros (const cholmod_sparse *S);
+#endif
+
+  octave_idx_type init (const chol_type& a, bool natural, bool force);
+
+  // No copying!
+
+  sparse_chol_rep (const sparse_chol_rep&);
+
+  sparse_chol_rep& operator = (const sparse_chol_rep&);
+};
+
+#ifdef HAVE_CHOLMOD
+
+// Can't use CHOLMOD_NAME(drop)(0.0, S, cm) because it doesn't treat
+// complex matrices.
+
+template <typename chol_type>
+void
+sparse_chol<chol_type>::sparse_chol_rep::drop_zeros (const cholmod_sparse *S)
+{
+  chol_elt sik;
+  octave_idx_type *Sp, *Si;
+  chol_elt *Sx;
+  octave_idx_type pdest, k, ncol, p, pend;
+
+  if (! S)
+    return;
+
+  Sp = static_cast<octave_idx_type *>(S->p);
+  Si = static_cast<octave_idx_type *>(S->i);
+  Sx = static_cast<chol_elt *>(S->x);
+  pdest = 0;
+  ncol = S->ncol;
+
+  for (k = 0; k < ncol; k++)
+    {
+      p = Sp[k];
+      pend = Sp[k+1];
+      Sp[k] = pdest;
+      for (; p < pend; p++)
+        {
+          sik = Sx[p];
+          if (CHOLMOD_IS_NONZERO (sik))
+            {
+              if (p != pdest)
+                {
+                  Si[pdest] = Si[p];
+                  Sx[pdest] = sik;
+                }
+              pdest++;
+            }
+        }
+    }
+  Sp[ncol] = pdest;
+}
+
+// Must provide a specialization for this function.
+template <typename T>
+int
+get_xtype (void);
+
+template <>
+inline int
+get_xtype<double> (void)
+{
+  return CHOLMOD_REAL;
+}
+
+template <>
+inline int
+get_xtype<Complex> (void)
+{
+  return CHOLMOD_COMPLEX;
+}
+
+#endif
+
+template <typename chol_type>
+octave_idx_type
+sparse_chol<chol_type>::sparse_chol_rep::init (const chol_type& a,
+                                               bool natural, bool force)
+{
+  volatile octave_idx_type info = 0;
+
+#ifdef HAVE_CHOLMOD
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("sparse_chol requires square matrix");
+
+  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;
+      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
+    }
+  else
+    {
+      cm->print = static_cast<int> (spu) + 2;
+      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
+    }
+
+  cm->error_handler = &SparseCholError;
+  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
+  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, 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;
+  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;
+#if defined (ENABLE_64)
+  ac->itype = CHOLMOD_LONG;
+#else
+  ac->itype = CHOLMOD_INT;
+#endif
+  ac->dtype = CHOLMOD_DOUBLE;
+  ac->stype = 1;
+  ac->xtype = get_xtype<chol_elt> ();
+
+  if (a_nr < 1)
+    ac->x = &dummy;
+  else
+    ac->x = a.data ();
+
+  // use natural ordering if no q output parameter
+  if (natural)
+    {
+      cm->nmethods = 1 ;
+      cm->method[0].ordering = CHOLMOD_NATURAL ;
+      cm->postorder = false ;
+    }
+
+  cholmod_factor *Lfactor;
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  Lfactor = CHOLMOD_NAME(analyze) (ac, cm);
+  CHOLMOD_NAME(factorize) (ac, Lfactor, cm);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  is_pd = cm->status == CHOLMOD_OK;
+  info = (is_pd ? 0 : cm->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;
+    }
+
+  return info;
+
+#else
+  (*current_liboctave_error_handler)
+    ("support for CHOLMOD was unavailable or disabled when liboctave was built");
+#endif
+}
+
+template <typename chol_type>
+SparseMatrix
+sparse_chol<chol_type>::sparse_chol_rep::Q (void) const
+{
+#ifdef HAVE_CHOLMOD
+  octave_idx_type n = Lsparse->nrow;
+  SparseMatrix p (n, n, n);
+
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      p.xcidx (i) = i;
+      p.xridx (i) = static_cast<octave_idx_type>(perms (i));
+      p.xdata (i) = 1;
+    }
+  p.xcidx (n) = n;
+
+  return p;
+#else
+  return SparseMatrix ();
+#endif
+}
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (void)
+  : rep (new typename sparse_chol<chol_type>::sparse_chol_rep ())
+{ }
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (const chol_type& a, bool natural,
+                                     bool force)
+  : rep (new typename
+         sparse_chol<chol_type>::sparse_chol_rep (a, natural, force))
+{ }
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (const chol_type& a, octave_idx_type& info,
+                                     bool natural, bool force)
+  : rep (new typename
+         sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, force))
+{ }
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (const chol_type& a, octave_idx_type& info,
+                                     bool natural)
+  : rep (new typename
+         sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, false))
+{ }
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (const chol_type& a, octave_idx_type& info)
+  : rep (new typename
+         sparse_chol<chol_type>::sparse_chol_rep (a, info, false, false))
+{ }
+
+template <typename chol_type>
+sparse_chol<chol_type>::sparse_chol (const sparse_chol<chol_type>& a)
+  : rep (a.rep)
+{
+  rep->count++;
+}
+
+template <typename chol_type>
+sparse_chol<chol_type>::~sparse_chol (void)
+{
+  if (--rep->count == 0)
+    delete rep;
+}
+
+template <typename chol_type>
+sparse_chol<chol_type>&
+sparse_chol<chol_type>::operator = (const sparse_chol& a)
+{
+  if (this != &a)
+    {
+      if (--rep->count == 0)
+        delete rep;
+
+      rep = a.rep;
+      rep->count++;
+    }
+
+  return *this;
+}
+
+template <typename chol_type>
+chol_type
+sparse_chol<chol_type>::L (void) const
+{
+#ifdef HAVE_CHOLMOD
+  cholmod_sparse *m = rep->L ();
+  octave_idx_type nc = m->ncol;
+  octave_idx_type nnz = m->nzmax;
+  chol_type ret (m->nrow, nc, nnz);
+  for (octave_idx_type j = 0; j < nc+1; j++)
+    ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j];
+  for (octave_idx_type i = 0; i < nnz; i++)
+    {
+      ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i];
+      ret.xdata (i) = static_cast<chol_elt *>(m->x)[i];
+    }
+  return ret;
+#else
+  return chol_type ();
+#endif
+}
+
+template <typename chol_type>
+octave_idx_type
+sparse_chol<chol_type>::P (void) const
+{
+  return rep->P ();
+}
+
+template <typename chol_type>
+ColumnVector
+sparse_chol<chol_type>::perm (void) const
+{
+  return rep->perm ();
+}
+
+template <typename chol_type>
+SparseMatrix
+sparse_chol<chol_type>::Q (void) const
+{
+  return rep->Q ();
+}
+
+template <typename chol_type>
+bool
+sparse_chol<chol_type>::is_positive_definite (void) const
+{
+  return rep->is_positive_definite ();
+}
+
+template <typename chol_type>
+double
+sparse_chol<chol_type>::rcond (void) const
+{
+  return rep->rcond ();
+}
+
+template <typename chol_type>
+chol_type
+sparse_chol<chol_type>::inverse (void) const
+{
+  chol_type retval;
+#ifdef HAVE_CHOLMOD
+  cholmod_sparse *m = rep->L ();
+  octave_idx_type n = m->ncol;
+  ColumnVector perms = rep->perm ();
+  double rcond2;
+  octave_idx_type info;
+  MatrixType mattype (MatrixType::Upper);
+  chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0);
+
+  if (perms.numel () == n)
+    {
+      SparseMatrix Qc = Q ();
+      retval = Qc * linv * linv.hermitian () * Qc.transpose ();
+    }
+  else
+    retval = linv * linv.hermitian ();
+#endif
+  return retval;
+}
+
+template <typename chol_type>
+chol_type
+chol2inv (const chol_type& r)
+{
+  octave_idx_type r_nr = r.rows ();
+  octave_idx_type r_nc = r.cols ();
+  chol_type retval;
+
+  if (r_nr != r_nc)
+    (*current_liboctave_error_handler) ("U must be a square matrix");
+
+  MatrixType mattype (r);
+  int typ = mattype.type (false);
+  double rcond;
+  octave_idx_type info;
+  chol_type rinv;
+
+  if (typ == MatrixType::Upper)
+    {
+      rinv = r.inverse (mattype, info, rcond, true, false);
+      retval = rinv.transpose () * rinv;
+    }
+  else if (typ == MatrixType::Lower)
+    {
+      rinv = r.transpose ().inverse (mattype, info, rcond, true, false);
+      retval = rinv.transpose () * rinv;
+    }
+  else
+    (*current_liboctave_error_handler) ("U must be a triangular matrix");
+
+  return retval;
+}
+
+// SparseComplexMatrix specialization (the value for the NATURAL
+// parameter in the sparse_chol<T>::sparse_chol_rep constructor is
+// different from the default).
+
+template <>
+sparse_chol<SparseComplexMatrix>::sparse_chol (const SparseComplexMatrix& a,
+                                               octave_idx_type& info)
+  : rep (new typename
+         sparse_chol<SparseComplexMatrix>::sparse_chol_rep (a, info, true, false))
+{ }
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-chol.h	Tue Jan 26 16:30:54 2016 -0500
@@ -0,0 +1,92 @@
+/*
+
+Copyright (C) 2005-2015 David Bateman
+Copyright (C) 1998-2005 Andy Adler
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#if ! defined (octave_sparse_chol_h)
+#define octave_sparse_chol_h 1
+
+#include "CSparse.h"
+#include "dColVector.h"
+#include "dSparse.h"
+#include "oct-sparse.h"
+
+template <typename chol_type>
+class
+sparse_chol
+{
+public:
+
+  sparse_chol (void);
+
+  sparse_chol (const chol_type& a, bool natural, bool force);
+
+  sparse_chol (const chol_type& a, octave_idx_type& info,
+               bool natural, bool force);
+
+  sparse_chol (const chol_type& a, octave_idx_type& info, bool natural);
+
+  sparse_chol (const chol_type& a, octave_idx_type& info);
+
+  sparse_chol (const sparse_chol<chol_type>& a);
+
+  virtual ~sparse_chol (void);
+
+  sparse_chol& operator = (const sparse_chol& a);
+
+  chol_type L (void) const;
+
+  chol_type R (void) const { return L ().hermitian (); }
+
+  octave_idx_type P (void) const;
+
+  ColumnVector perm (void) const;
+
+  SparseMatrix Q (void) const;
+
+  bool is_positive_definite (void) const;
+
+  double rcond (void) const;
+
+  chol_type inverse (void) const;
+
+protected:
+
+  typedef typename chol_type::element_type chol_elt;
+
+  class sparse_chol_rep;
+
+private:
+
+  sparse_chol_rep *rep;
+};
+
+template <typename chol_type>
+chol_type
+chol2inv (const chol_type& r);
+
+// SparseComplexMatrix specialization.
+
+template <>
+sparse_chol<SparseComplexMatrix>::sparse_chol (const SparseComplexMatrix& a,
+                                               octave_idx_type& info);
+
+#endif