--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/__eigs__.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,675 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 2005-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <limits>
+#include <string>
+
+#include "Matrix.h"
+#include "eigs-base.h"
+#include "unwind-prot.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "interpreter-private.h"
+#include "oct-map.h"
+#include "ov.h"
+#include "ovl.h"
+#include "pager.h"
+#include "parse.h"
+#include "variables.h"
+
+#if defined (HAVE_ARPACK)
+
+// Global pointer for user defined function.
+static octave_value eigs_fcn;
+
+// Have we warned about imaginary values returned from user function?
+static bool warned_imaginary = false;
+
+// Is this a recursive call?
+static int call_depth = 0;
+
+ColumnVector
+eigs_func (const ColumnVector& x, int& eigs_error)
+{
+  ColumnVector retval;
+  octave_value_list args;
+  args(0) = x;
+
+  if (eigs_fcn.is_defined ())
+    {
+      octave_value_list tmp;
+
+      try
+        {
+          tmp = octave::feval (eigs_fcn, args, 1);
+        }
+      catch (octave::execution_exception& e)
+        {
+          err_user_supplied_eval (e, "eigs");
+        }
+
+      if (tmp.length () && tmp(0).is_defined ())
+        {
+          if (! warned_imaginary && tmp(0).iscomplex ())
+            {
+              warning ("eigs: ignoring imaginary part returned from user-supplied function");
+              warned_imaginary = true;
+            }
+
+          retval = tmp(0).xvector_value ("eigs: evaluation of user-supplied function failed");
+        }
+      else
+        {
+          eigs_error = 1;
+          err_user_supplied_eval ("eigs");
+        }
+    }
+
+  return retval;
+}
+
+ComplexColumnVector
+eigs_complex_func (const ComplexColumnVector& x, int& eigs_error)
+{
+  ComplexColumnVector retval;
+  octave_value_list args;
+  args(0) = x;
+
+  if (eigs_fcn.is_defined ())
+    {
+      octave_value_list tmp;
+
+      try
+        {
+          tmp = octave::feval (eigs_fcn, args, 1);
+        }
+      catch (octave::execution_exception& e)
+        {
+          err_user_supplied_eval (e, "eigs");
+        }
+
+      if (tmp.length () && tmp(0).is_defined ())
+        {
+          retval = tmp(0).xcomplex_vector_value ("eigs: evaluation of user-supplied function failed");
+        }
+      else
+        {
+          eigs_error = 1;
+          err_user_supplied_eval ("eigs");
+        }
+    }
+
+  return retval;
+}
+
+#endif
+
+DEFMETHOD (__eigs__, interp, args, nargout,
+           doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{d} =} __eigs__ (@var{A})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{k})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{k}, @var{sigma}, @var{opts})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{B})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{A}, @var{B}, @var{k}, @var{sigma}, @var{opts})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{k}, @var{sigma}, @var{opts})
+@deftypefnx {} {@var{d} =} __eigs__ (@var{af}, @var{n}, @var{B}, @var{k}, @var{sigma}, @var{opts})
+@deftypefnx {} {[@var{V}, @var{d}] =} __eigs__ (@var{A}, @dots{})
+@deftypefnx {} {[@var{V}, @var{d}] =} __eigs__ (@var{af}, @var{n}, @dots{})
+@deftypefnx {} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{A}, @dots{})
+@deftypefnx {} {[@var{V}, @var{d}, @var{flag}] =} __eigs__ (@var{af}, @var{n}, @dots{})
+Undocumented internal function.
+@end deftypefn */)
+{
+#if defined (HAVE_ARPACK)
+
+  int nargin = args.length ();
+
+  if (nargin == 0)
+    print_usage ();
+
+  octave_value_list retval;
+
+  std::string fcn_name;
+  octave_idx_type n = 0;
+  octave_idx_type k = 6;
+  Complex sigma = 0.0;
+  double sigmar, sigmai;
+  bool have_sigma = false;
+  std::string typ = "LM";
+  Matrix amm, bmm, bmt;
+  ComplexMatrix acm, bcm, bct;
+  SparseMatrix asmm, bsmm, bsmt;
+  SparseComplexMatrix ascm, bscm, bsct;
+  int b_arg = 0;
+  bool have_b = false;
+  bool have_a_fun = false;
+  bool a_is_complex = false;
+  bool b_is_complex = false;
+  bool symmetric = false;
+  bool sym_tested = false;
+  bool cholB = false;
+  bool a_is_sparse = false;
+  bool b_is_sparse = false;
+  ColumnVector permB;
+  int arg_offset = 0;
+  double tol = std::numeric_limits<double>::epsilon ();
+  int maxit = 300;
+  int disp = 0;
+  octave_idx_type p = -1;
+  ColumnVector resid;
+  ComplexColumnVector cresid;
+  octave_idx_type info = 1;
+
+  warned_imaginary = false;
+
+  octave::unwind_protect frame;
+
+  frame.protect_var (call_depth);
+  call_depth++;
+
+  if (call_depth > 1)
+    error ("eigs: invalid recursive call");
+
+  if (args(0).is_function_handle () || args(0).is_inline_function ()
+      || args(0).is_string ())
+    {
+      eigs_fcn = octave::get_function_handle (interp, args(0), "x");
+
+      if (eigs_fcn.is_undefined ())
+        error ("eigs: unknown function");
+
+      if (nargin < 2)
+        error ("eigs: incorrect number of arguments");
+
+      n = args(1).nint_value ();
+      arg_offset = 1;
+      have_a_fun = true;
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        {
+          if (args(0).issparse ())
+            {
+              ascm = (args(0).sparse_complex_matrix_value ());
+              a_is_sparse = true;
+            }
+          else
+            acm = (args(0).complex_matrix_value ());
+          a_is_complex = true;
+        }
+      else
+        {
+          if (args(0).issparse ())
+            {
+              asmm = (args(0).sparse_matrix_value ());
+              a_is_sparse = true;
+            }
+          else
+            {
+              amm = (args(0).matrix_value ());
+            }
+        }
+    }
+
+  // Note hold off reading B until later to avoid issues of double
+  // copies of the matrix if B is full/real while A is complex.
+  if (nargin > 1 + arg_offset
+      && ! (args(1 + arg_offset).is_real_scalar ()))
+    {
+      if (args(1+arg_offset).iscomplex ())
+        {
+          b_arg = 1+arg_offset;
+          if (args(b_arg).issparse ())
+            {
+              bscm = (args(b_arg).sparse_complex_matrix_value ());
+              b_is_sparse = true;
+            }
+          else
+            bcm = (args(b_arg).complex_matrix_value ());
+          have_b = true;
+          b_is_complex = true;
+          arg_offset++;
+        }
+      else
+        {
+          b_arg = 1+arg_offset;
+          if (args(b_arg).issparse ())
+            {
+              bsmm = (args(b_arg).sparse_matrix_value ());
+              b_is_sparse = true;
+            }
+          else
+            bmm = (args(b_arg).matrix_value ());
+          have_b = true;
+          arg_offset++;
+        }
+    }
+
+  if (nargin > (1+arg_offset))
+    k = args(1+arg_offset).nint_value ();
+
+  if (nargin > (2+arg_offset))
+    {
+      if (args(2+arg_offset).is_string ())
+        {
+          typ = args(2+arg_offset).string_value ();
+
+          // Use STL function to convert to upper case
+          transform (typ.begin (), typ.end (), typ.begin (), toupper);
+
+          sigma = 0.0;
+        }
+      else
+        {
+          sigma = args(2+arg_offset).xcomplex_value ("eigs: SIGMA must be a scalar or a string");
+
+          have_sigma = true;
+        }
+    }
+
+  sigmar = sigma.real ();
+  sigmai = sigma.imag ();
+
+  if (nargin > (3+arg_offset))
+    {
+      if (! args(3+arg_offset).isstruct ())
+        error ("eigs: OPTS argument must be a structure");
+
+      octave_scalar_map map = args(3+arg_offset).xscalar_map_value ("eigs: OPTS argument must be a scalar structure");
+
+      octave_value tmp;
+
+      // issym is ignored for complex matrix inputs
+      tmp = map.getfield ("issym");
+      if (tmp.is_defined ())
+        {
+          if (tmp.numel () != 1)
+            error ("eigs: OPTS.issym must be a scalar value");
+
+          symmetric = tmp.xbool_value ("eigs: OPTS.issym must be a logical value");
+          sym_tested = true;
+        }
+
+      // isreal is ignored if A is not a function
+      if (have_a_fun)
+        {
+          tmp = map.getfield ("isreal");
+          if (tmp.is_defined ())
+            {
+              if (tmp.numel () != 1)
+                error ("eigs: OPTS.isreal must be a scalar value");
+
+              a_is_complex = ! tmp.xbool_value ("eigs: OPTS.isreal must be a logical value");
+            }
+        }
+
+      tmp = map.getfield ("tol");
+      if (tmp.is_defined ())
+        tol = tmp.double_value ();
+
+      tmp = map.getfield ("maxit");
+      if (tmp.is_defined ())
+        maxit = tmp.nint_value ();
+
+      tmp = map.getfield ("p");
+      if (tmp.is_defined ())
+        p = tmp.nint_value ();
+
+      tmp = map.getfield ("v0");
+      if (tmp.is_defined ())
+        {
+          if (a_is_complex || b_is_complex)
+            cresid = ComplexColumnVector (tmp.complex_vector_value ());
+          else
+            resid = ColumnVector (tmp.vector_value ());
+        }
+
+      tmp = map.getfield ("disp");
+      if (tmp.is_defined ())
+        disp = tmp.nint_value ();
+
+      tmp = map.getfield ("cholB");
+      if (tmp.is_defined ())
+        {
+          if (tmp.numel () != 1)
+            error ("eigs: OPTS.cholB must be a scalar value");
+
+          cholB = tmp.xbool_value ("eigs: OPTS.cholB must be a logical value");
+        }
+
+      tmp = map.getfield ("permB");
+      if (tmp.is_defined ())
+        permB = ColumnVector (tmp.vector_value ()) - 1.0;
+    }
+
+  if (nargin > (4+arg_offset))
+    error ("eigs: incorrect number of arguments");
+
+  // Test undeclared (no issym) matrix inputs for symmetry
+  if (! sym_tested && ! have_a_fun)
+    {
+      if (a_is_complex)
+        {
+          if (a_is_sparse)
+            symmetric = ascm.ishermitian ();
+          else
+            symmetric = acm.ishermitian ();
+        }
+      else
+        {
+          if (a_is_sparse)
+            symmetric = asmm.issymmetric ();
+          else
+            symmetric = amm.issymmetric ();
+        }
+    }
+
+  if (have_b)
+    {
+      if (a_is_complex || b_is_complex)
+        {
+          if (b_is_sparse)
+            bscm = args(b_arg).sparse_complex_matrix_value ();
+          else
+            bcm = args(b_arg).complex_matrix_value ();
+        }
+      else
+        {
+          if (b_is_sparse)
+            bsmm = args(b_arg).sparse_matrix_value ();
+          else
+            bmm = args(b_arg).matrix_value ();
+        }
+    }
+
+  // Mode 1 for SM mode seems unstable for some reason.
+  // Use Mode 3 instead, with sigma = 0.
+  if (! have_sigma && typ == "SM")
+    have_sigma = true;
+
+  octave_idx_type nconv;
+  if (a_is_complex || b_is_complex)
+    {
+      ComplexMatrix eig_vec;
+      ComplexColumnVector eig_val;
+
+      if (have_a_fun)
+        {
+          if (b_is_sparse)
+            nconv = EigsComplexNonSymmetricFunc
+              (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
+               eig_val, bscm, permB, cresid, octave_stdout, tol,
+               (nargout > 1), cholB, disp, maxit);
+          else
+            nconv = EigsComplexNonSymmetricFunc
+              (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
+               eig_val, bcm, permB, cresid, octave_stdout, tol,
+               (nargout > 1), cholB, disp, maxit);
+        }
+      else if (have_sigma)
+        {
+          if (a_is_sparse)
+            nconv = EigsComplexNonSymmetricMatrixShift
+                    (ascm, sigma, k, p, info, eig_vec, eig_val, bscm, permB,
+                     cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+                     maxit);
+          else
+            nconv = EigsComplexNonSymmetricMatrixShift
+                    (acm, sigma, k, p, info, eig_vec, eig_val, bcm, permB,
+                     cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+                     maxit);
+        }
+      else
+        {
+          if (a_is_sparse)
+            nconv = EigsComplexNonSymmetricMatrix
+                    (ascm, typ, k, p, info, eig_vec, eig_val, bscm, permB,
+                     cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+                     maxit);
+          else
+            nconv = EigsComplexNonSymmetricMatrix
+                    (acm, typ, k, p, info, eig_vec, eig_val, bcm, permB,
+                     cresid, octave_stdout, tol, (nargout > 1), cholB, disp,
+                     maxit);
+        }
+
+      if (nargout < 2)
+        {
+          if (symmetric)
+            retval(0) = real (eig_val);
+          else
+            retval(0) = eig_val;
+        }
+      else
+        {
+          if (symmetric)
+            retval = ovl (eig_vec, DiagMatrix (real (eig_val)), double (info));
+          else
+            retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
+        }
+    }
+  else if (sigmai != 0.0)
+    {
+      // Promote real problem to a complex one.
+      ComplexMatrix eig_vec;
+      ComplexColumnVector eig_val;
+
+      if (have_a_fun)
+        {
+          if (b_is_sparse)
+            nconv = EigsComplexNonSymmetricFunc
+              (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
+               eig_val, bscm, permB, cresid, octave_stdout, tol,
+               (nargout > 1), cholB, disp, maxit);
+          else
+            nconv = EigsComplexNonSymmetricFunc
+              (eigs_complex_func, n, typ, sigma, k, p, info, eig_vec,
+               eig_val, bcm, permB, cresid, octave_stdout, tol,
+               (nargout > 1), cholB, disp, maxit);
+        }
+      else
+        {
+          if (a_is_sparse)
+            nconv = EigsComplexNonSymmetricMatrixShift
+                    (SparseComplexMatrix (asmm), sigma, k, p, info, eig_vec,
+                     eig_val, SparseComplexMatrix (bsmm), permB, cresid,
+                     octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+          else
+            nconv = EigsComplexNonSymmetricMatrixShift
+                    (ComplexMatrix (amm), sigma, k, p, info, eig_vec,
+                     eig_val, ComplexMatrix (bmm), permB, cresid,
+                     octave_stdout, tol, (nargout > 1), cholB, disp, maxit);
+        }
+
+      if (nargout < 2)
+        {
+          if (symmetric)
+            retval(0) = real (eig_val);
+          else
+            retval(0) = eig_val;
+        }
+      else
+        {
+          if (symmetric)
+            retval = ovl (eig_vec, DiagMatrix (real (eig_val)), double (info));
+          else
+            retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
+        }
+    }
+  else
+    {
+      if (symmetric)
+        {
+          Matrix eig_vec;
+          ColumnVector eig_val;
+
+          if (have_a_fun)
+            {
+              if (b_is_sparse)
+                nconv = EigsRealSymmetricFunc
+                       (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
+                        eig_val, bsmm, permB, resid, octave_stdout, tol,
+                        (nargout > 1), cholB, disp, maxit);
+              else
+                nconv = EigsRealSymmetricFunc
+                       (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
+                        eig_val, bmm, permB, resid, octave_stdout, tol,
+                        (nargout > 1), cholB, disp, maxit);
+            }
+          else if (have_sigma)
+            {
+              if (a_is_sparse)
+                nconv = EigsRealSymmetricMatrixShift
+                        (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+              else
+                nconv = EigsRealSymmetricMatrixShift
+                        (amm, sigmar, k, p, info, eig_vec, eig_val, bmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+            }
+          else
+            {
+              if (a_is_sparse)
+                nconv = EigsRealSymmetricMatrix
+                        (asmm, typ, k, p, info, eig_vec, eig_val, bsmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+              else
+                nconv = EigsRealSymmetricMatrix
+                        (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+                         resid, octave_stdout, tol, (nargout > 1), cholB,
+                         disp, maxit);
+            }
+
+          if (nargout < 2)
+            retval(0) = eig_val;
+          else
+            retval = ovl (eig_vec, DiagMatrix (eig_val), double (info));
+        }
+      else
+        {
+          ComplexMatrix eig_vec;
+          ComplexColumnVector eig_val;
+
+          if (have_a_fun)
+            {
+              if (b_is_sparse)
+                nconv = EigsRealNonSymmetricFunc
+                        (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
+                         eig_val, bsmm, permB, resid, octave_stdout, tol,
+                         (nargout > 1), cholB, disp, maxit);
+              else
+                nconv = EigsRealNonSymmetricFunc
+                        (eigs_func, n, typ, sigmar, k, p, info, eig_vec,
+                         eig_val, bmm, permB, resid, octave_stdout, tol,
+                         (nargout > 1), cholB, disp, maxit);
+            }
+          else if (have_sigma)
+            {
+              if (a_is_sparse)
+                nconv = EigsRealNonSymmetricMatrixShift
+                        (asmm, sigmar, k, p, info, eig_vec, eig_val, bsmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+              else
+                nconv = EigsRealNonSymmetricMatrixShift
+                        (amm, sigmar, k, p, info, eig_vec, eig_val, bmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+            }
+          else
+            {
+              if (a_is_sparse)
+                nconv = EigsRealNonSymmetricMatrix
+                        (asmm, typ, k, p, info, eig_vec, eig_val, bsmm,
+                         permB, resid, octave_stdout, tol, (nargout > 1),
+                         cholB, disp, maxit);
+              else
+                nconv = EigsRealNonSymmetricMatrix
+                        (amm, typ, k, p, info, eig_vec, eig_val, bmm, permB,
+                         resid, octave_stdout, tol, (nargout > 1), cholB,
+                         disp, maxit);
+            }
+
+          if (nargout < 2)
+            retval(0) = eig_val;
+          else
+            retval = ovl (eig_vec, ComplexDiagMatrix (eig_val), double (info));
+        }
+    }
+
+  if (nconv <= 0)
+    warning_with_id ("Octave:eigs:UnconvergedEigenvalues",
+                     "eigs: None of the %" OCTAVE_IDX_TYPE_FORMAT
+                     " requested eigenvalues converged", k);
+  else if (nconv < k)
+    warning_with_id ("Octave:eigs:UnconvergedEigenvalues",
+                     "eigs: Only %" OCTAVE_IDX_TYPE_FORMAT
+                     " of the %" OCTAVE_IDX_TYPE_FORMAT
+                     " requested eigenvalues converged",
+                     nconv, k);
+
+  if (! fcn_name.empty ())
+    {
+      octave::symbol_table& symtab = interp.get_symbol_table ();
+
+      symtab.clear_function (fcn_name);
+    }
+
+  return retval;
+
+#else
+
+  octave_unused_parameter (interp);
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("eigs", "ARPACK");
+
+#endif
+}
+
+/*
+## No test needed for internal helper function.
+%!assert (1)
+*/

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/amd.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,202 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 2008-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+// This is the octave interface to amd, which bore the copyright given
+// in the help of the functions.
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <cstdlib>
+
+#include "CSparse.h"
+#include "Sparse.h"
+#include "dMatrix.h"
+#include "oct-locbuf.h"
+#include "oct-sparse.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "oct-map.h"
+#include "ov.h"
+#include "ovl.h"
+#include "parse.h"
+
+DEFUN (amd, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} amd (@var{S})
+@deftypefnx {} {@var{p} =} amd (@var{S}, @var{opts})
+
+Return the approximate minimum degree permutation of a matrix.
+
+This is a permutation such that the Cholesky@tie{}factorization of
+@code{@var{S} (@var{p}, @var{p})} tends to be sparser than the
+Cholesky@tie{}factorization of @var{S} itself.  @code{amd} is typically
+faster than @code{symamd} but serves a similar purpose.
+
+The optional parameter @var{opts} is a structure that controls the behavior
+of @code{amd}.  The fields of the structure are
+
+@table @asis
+@item @var{opts}.dense
+Determines what @code{amd} considers to be a dense row or column of the
+input matrix.  Rows or columns with more than @code{max (16, (dense *
+sqrt (@var{n})))} entries, where @var{n} is the order of the matrix @var{S},
+are ignored by @code{amd} during the calculation of the permutation.
+The value of dense must be a positive scalar and the default value is 10.0
+
+@item @var{opts}.aggressive
+If this value is a nonzero scalar, then @code{amd} performs aggressive
+absorption.  The default is not to perform aggressive absorption.
+@end table
+
+The author of the code itself is Timothy A. Davis
+(see @url{http://faculty.cse.tamu.edu/davis/suitesparse.html}).
+@seealso{symamd, colamd}
+@end deftypefn */)
+{
+#if defined (HAVE_AMD)
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+
+  octave_idx_type n_row, n_col;
+  const octave::suitesparse_integer *ridx, *cidx;
+  SparseMatrix sm;
+  SparseComplexMatrix scm;
+
+  if (args(0).issparse ())
+    {
+      if (args(0).iscomplex ())
+        {
+          scm = args(0).sparse_complex_matrix_value ();
+          n_row = scm.rows ();
+          n_col = scm.cols ();
+          ridx = octave::to_suitesparse_intptr (scm.xridx ());
+          cidx = octave::to_suitesparse_intptr (scm.xcidx ());
+        }
+      else
+        {
+          sm = args(0).sparse_matrix_value ();
+          n_row = sm.rows ();
+          n_col = sm.cols ();
+          ridx = octave::to_suitesparse_intptr (sm.xridx ());
+          cidx = octave::to_suitesparse_intptr (sm.xcidx ());
+        }
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+      else
+        sm = SparseMatrix (args(0).matrix_value ());
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      ridx = octave::to_suitesparse_intptr (sm.xridx ());
+      cidx = octave::to_suitesparse_intptr (sm.xcidx ());
+    }
+
+  if (n_row != n_col)
+    err_square_matrix_required ("amd", "S");
+
+  OCTAVE_LOCAL_BUFFER (double, Control, AMD_CONTROL);
+  AMD_NAME (_defaults) (Control);
+  if (nargin > 1)
+    {
+      octave_scalar_map arg1 = args(1).xscalar_map_value ("amd: OPTS argument must be a scalar structure");
+
+      octave_value tmp;
+
+      tmp = arg1.getfield ("dense");
+      if (tmp.is_defined ())
+        Control[AMD_DENSE] = tmp.double_value ();
+
+      tmp = arg1.getfield ("aggressive");
+      if (tmp.is_defined ())
+        Control[AMD_AGGRESSIVE] = tmp.double_value ();
+    }
+
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, P, n_col);
+  Matrix xinfo (AMD_INFO, 1);
+  double *Info = xinfo.fortran_vec ();
+
+  // FIXME: how can we manage the memory allocation of amd
+  //        in a cleaner manner?
+  SUITESPARSE_ASSIGN_FPTR (malloc_func, amd_malloc, malloc);
+  SUITESPARSE_ASSIGN_FPTR (free_func, amd_free, free);
+  SUITESPARSE_ASSIGN_FPTR (calloc_func, amd_calloc, calloc);
+  SUITESPARSE_ASSIGN_FPTR (realloc_func, amd_realloc, realloc);
+  SUITESPARSE_ASSIGN_FPTR (printf_func, amd_printf, printf);
+
+  octave_idx_type result = AMD_NAME (_order) (n_col, cidx, ridx, P, Control,
+                                              Info);
+
+  if (result == AMD_OUT_OF_MEMORY)
+    error ("amd: out of memory");
+  else if (result == AMD_INVALID)
+    error ("amd: matrix S is corrupted");
+
+  Matrix Pout (1, n_col);
+  for (octave_idx_type i = 0; i < n_col; i++)
+    Pout.xelem (i) = P[i] + 1;
+
+  if (nargout > 1)
+    return ovl (Pout, xinfo);
+  else
+    return ovl (Pout);
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("amd", "AMD");
+
+#endif
+}
+
+/*
+%!shared A, A2, opts
+%! A = ones (20, 30);
+%! A2 = ones (30, 30);
+
+%!testif HAVE_AMD
+%! assert(amd (A2), [1:30]);
+%! opts.dense = 25;
+%! assert(amd (A2, opts), [1:30]);
+%! opts.aggressive = 1;
+%! assert(amd (A2, opts), [1:30]);
+
+%!testif HAVE_AMD
+%! assert (amd ([]), zeros (1,0))
+
+%!error <S must be a square matrix|was unavailable or disabled> amd (A)
+%!error amd (A2, 2)
+*/

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/ccolamd.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,583 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 2005-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+// This is the octave interface to ccolamd, which bore the copyright given
+// in the help of the functions.
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <cstdlib>
+
+#include "CSparse.h"
+#include "Sparse.h"
+#include "dNDArray.h"
+#include "oct-locbuf.h"
+#include "oct-sparse.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "ov.h"
+#include "pager.h"
+
+DEFUN (ccolamd, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} ccolamd (@var{S})
+@deftypefnx {} {@var{p} =} ccolamd (@var{S}, @var{knobs})
+@deftypefnx {} {@var{p} =} ccolamd (@var{S}, @var{knobs}, @var{cmember})
+@deftypefnx {} {[@var{p}, @var{stats}] =} ccolamd (@dots{})
+
+Constrained column approximate minimum degree permutation.
+
+@code{@var{p} = ccolamd (@var{S})} returns the column approximate minimum
+degree permutation vector for the sparse matrix @var{S}.  For a
+non-symmetric matrix @var{S}, @code{@var{S}(:, @var{p})} tends to have
+sparser LU@tie{}factors than @var{S}.
+@code{chol (@var{S}(:, @var{p})' * @var{S}(:, @var{p}))} also tends to be
+sparser than @code{chol (@var{S}' * @var{S})}.
+@code{@var{p} = ccolamd (@var{S}, 1)} optimizes the ordering for
+@code{lu (@var{S}(:, @var{p}))}.  The ordering is followed by a column
+elimination tree post-ordering.
+
+@var{knobs} is an optional 1-element to 5-element input vector, with a
+default value of @code{[0 10 10 1 0]} if not present or empty.  Entries not
+present are set to their defaults.
+
+@table @code
+@item @var{knobs}(1)
+if nonzero, the ordering is optimized for @code{lu (S(:, p))}.  It will be a
+poor ordering for @code{chol (@var{S}(:, @var{p})' * @var{S}(:, @var{p}))}.
+This is the most important knob for ccolamd.
+
+@item @var{knobs}(2)
+if @var{S} is m-by-n, rows with more than
+@code{max (16, @var{knobs}(2) * sqrt (n))} entries are ignored.
+
+@item @var{knobs}(3)
+columns with more than
+@code{max (16, @var{knobs}(3) * sqrt (min (@var{m}, @var{n})))} entries are
+ignored and ordered last in the output permutation
+(subject to the cmember constraints).
+
+@item @var{knobs}(4)
+if nonzero, aggressive absorption is performed.
+
+@item @var{knobs}(5)
+if nonzero, statistics and knobs are printed.
+
+@end table
+
+@var{cmember} is an optional vector of length @math{n}.  It defines the
+constraints on the column ordering.  If @code{@var{cmember}(j) = @var{c}},
+then column @var{j} is in constraint set @var{c} (@var{c} must be in the
+range 1 to n).  In the output permutation @var{p}, all columns in set 1
+appear first, followed by all columns in set 2, and so on.
+@code{@var{cmember} = ones (1,n)} if not present or empty.
+@code{ccolamd (@var{S}, [], 1 : n)} returns @code{1 : n}
+
+@code{@var{p} = ccolamd (@var{S})} is about the same as
+@code{@var{p} = colamd (@var{S})}.  @var{knobs} and its default values
+differ.  @code{colamd} always does aggressive absorption, and it finds an
+ordering suitable for both @code{lu (@var{S}(:, @var{p}))} and @code{chol
+(@var{S}(:, @var{p})' * @var{S}(:, @var{p}))}; it cannot optimize its
+ordering for @code{lu (@var{S}(:, @var{p}))} to the extent that
+@code{ccolamd (@var{S}, 1)} can.
+
+@var{stats} is an optional 20-element output vector that provides data
+about the ordering and the validity of the input matrix @var{S}.  Ordering
+statistics are in @code{@var{stats}(1 : 3)}.  @code{@var{stats}(1)} and
+@code{@var{stats}(2)} are the number of dense or empty rows and columns
+ignored by @sc{ccolamd} and @code{@var{stats}(3)} is the number of garbage
+collections performed on the internal data structure used by @sc{ccolamd}
+(roughly of size @code{2.2 * nnz (@var{S}) + 4 * @var{m} + 7 * @var{n}}
+integers).
+
+@code{@var{stats}(4 : 7)} provide information if CCOLAMD was able to
+continue.  The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if
+invalid.  @code{@var{stats}(5)} is the rightmost column index that is
+unsorted or contains duplicate entries, or zero if no such column exists.
+@code{@var{stats}(6)} is the last seen duplicate or out-of-order row
+index in the column index given by @code{@var{stats}(5)}, or zero if no
+such row index exists.  @code{@var{stats}(7)} is the number of duplicate
+or out-of-order row indices.  @code{@var{stats}(8 : 20)} is always zero in
+the current version of @sc{ccolamd} (reserved for future use).
+
+The authors of the code itself are @nospell{S. Larimore, T. Davis} and
+@nospell{S. Rajamanickam} in collaboration with @nospell{J. Bilbert and E. Ng}.
+Supported by the National Science Foundation
+@nospell{(DMS-9504974, DMS-9803599, CCR-0203270)}, and a grant from
+@nospell{Sandia} National Lab.
+See @url{http://faculty.cse.tamu.edu/davis/suitesparse.html} for ccolamd,
+csymamd, amd, colamd, symamd, and other related orderings.
+@seealso{colamd, csymamd}
+@end deftypefn */)
+{
+#if defined (HAVE_CCOLAMD)
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+  int spumoni = 0;
+
+  // Get knobs
+  static_assert (CCOLAMD_KNOBS <= 40,
+                 "ccolamd: # of CCOLAMD_KNOBS exceeded.  Please report this to bugs.octave.org");
+  double knob_storage[CCOLAMD_KNOBS];
+  double *knobs = &knob_storage[0];
+  CCOLAMD_NAME (_set_defaults) (knobs);
+
+  // Check for user-passed knobs
+  if (nargin > 1)
+    {
+      NDArray User_knobs = args(1).array_value ();
+      int nel_User_knobs = User_knobs.numel ();
+
+      if (nel_User_knobs > 0)
+        knobs[CCOLAMD_LU] = (User_knobs(0) != 0);
+      if (nel_User_knobs > 1)
+        knobs[CCOLAMD_DENSE_ROW] = User_knobs(1);
+      if (nel_User_knobs > 2)
+        knobs[CCOLAMD_DENSE_COL] = User_knobs(2);
+      if (nel_User_knobs > 3)
+        knobs[CCOLAMD_AGGRESSIVE] = (User_knobs(3) != 0);
+      if (nel_User_knobs > 4)
+        spumoni = (User_knobs(4) != 0);
+
+      // print knob settings if spumoni is set
+      if (spumoni)
+        {
+          octave_stdout << "\nccolamd version " << CCOLAMD_MAIN_VERSION << '.'
+                        <<  CCOLAMD_SUB_VERSION << ", " << CCOLAMD_DATE
+                        << ":\nknobs(1): " << User_knobs(0) << ", order for ";
+          if (knobs[CCOLAMD_LU] != 0)
+            octave_stdout << "lu (A)\n";
+          else
+            octave_stdout << "chol (A'*A)\n";
+
+          if (knobs[CCOLAMD_DENSE_ROW] >= 0)
+            octave_stdout << "knobs(2): " << User_knobs(1)
+                          << ", rows with > max (16,"
+                          << knobs[CCOLAMD_DENSE_ROW]
+                          << "*sqrt (size(A,2)))"
+                          << " entries removed\n";
+          else
+            octave_stdout << "knobs(2): " << User_knobs(1)
+                          << ", no dense rows removed\n";
+
+          if (knobs[CCOLAMD_DENSE_COL] >= 0)
+            octave_stdout << "knobs(3): " << User_knobs(2)
+                          << ", cols with > max (16,"
+                          << knobs[CCOLAMD_DENSE_COL] << "*sqrt (size(A)))"
+                          << " entries removed\n";
+          else
+            octave_stdout << "knobs(3): " << User_knobs(2)
+                          << ", no dense columns removed\n";
+
+          if (knobs[CCOLAMD_AGGRESSIVE] != 0)
+            octave_stdout << "knobs(4): " << User_knobs(3)
+                          << ", aggressive absorption: yes";
+          else
+            octave_stdout << "knobs(4): " << User_knobs(3)
+                          << ", aggressive absorption: no";
+
+          octave_stdout << "knobs(5): " << User_knobs(4)
+                        << ", statistics and knobs printed\n";
+        }
+    }
+
+  octave_idx_type n_row, n_col, nnz;
+  octave_idx_type *ridx, *cidx;
+  SparseComplexMatrix scm;
+  SparseMatrix sm;
+
+  if (args(0).issparse ())
+    {
+      if (args(0).iscomplex ())
+        {
+          scm = args(0).sparse_complex_matrix_value ();
+          n_row = scm.rows ();
+          n_col = scm.cols ();
+          nnz = scm.nnz ();
+          ridx = scm.xridx ();
+          cidx = scm.xcidx ();
+        }
+      else
+        {
+          sm = args(0).sparse_matrix_value ();
+
+          n_row = sm.rows ();
+          n_col = sm.cols ();
+          nnz = sm.nnz ();
+          ridx = sm.xridx ();
+          cidx = sm.xcidx ();
+        }
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+      else
+        sm = SparseMatrix (args(0).matrix_value ());
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      nnz = sm.nnz ();
+      ridx = sm.xridx ();
+      cidx = sm.xcidx ();
+    }
+
+  // Allocate workspace for ccolamd
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, p, n_col+1);
+  for (octave_idx_type i = 0; i < n_col+1; i++)
+    p[i] = cidx[i];
+
+  octave_idx_type Alen = CCOLAMD_NAME (_recommended) (nnz, n_row, n_col);
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, A, Alen);
+  for (octave_idx_type i = 0; i < nnz; i++)
+    A[i] = ridx[i];
+
+  static_assert (CCOLAMD_STATS <= 40,
+                 "ccolamd: # of CCOLAMD_STATS exceeded.  Please report this to bugs.octave.org");
+  octave::suitesparse_integer stats_storage[CCOLAMD_STATS];
+  octave::suitesparse_integer *stats = &stats_storage[0];
+
+  if (nargin > 2)
+    {
+      NDArray in_cmember = args(2).array_value ();
+      octave_idx_type cslen = in_cmember.numel ();
+      OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, cmember, cslen);
+      for (octave_idx_type i = 0; i < cslen; i++)
+        // convert cmember from 1-based to 0-based
+        cmember[i] = static_cast<octave::suitesparse_integer>(in_cmember(i) - 1);
+
+      if (cslen != n_col)
+        error ("ccolamd: CMEMBER must be of length equal to #cols of A");
+
+      // Order the columns (destroys A)
+      if (! CCOLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats,cmember))
+        {
+          CCOLAMD_NAME (_report) (stats);
+
+          error ("ccolamd: internal error!");
+        }
+    }
+  else
+    {
+      // Order the columns (destroys A)
+      if (! CCOLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats, nullptr))
+        {
+          CCOLAMD_NAME (_report) (stats);
+
+          error ("ccolamd: internal error!");
+        }
+    }
+
+  // return the permutation vector
+  NDArray out_perm (dim_vector (1, n_col));
+  for (octave_idx_type i = 0; i < n_col; i++)
+    out_perm(i) = p[i] + 1;
+
+  retval(0) = out_perm;
+
+  // print stats if spumoni > 0
+  if (spumoni > 0)
+    CCOLAMD_NAME (_report) (stats);
+
+  // Return the stats vector
+  if (nargout == 2)
+    {
+      NDArray out_stats (dim_vector (1, CCOLAMD_STATS));
+      for (octave_idx_type i = 0 ; i < CCOLAMD_STATS ; i++)
+        out_stats(i) = stats[i];
+      retval(1) = out_stats;
+
+      // fix stats (5) and (6), for 1-based information on
+      // jumbled matrix.  note that this correction doesn't
+      // occur if symamd returns FALSE
+      out_stats(CCOLAMD_INFO1)++;
+      out_stats(CCOLAMD_INFO2)++;
+    }
+
+  return retval;
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("ccolamd", "CCOLAMD");
+
+#endif
+}
+
+DEFUN (csymamd, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} csymamd (@var{S})
+@deftypefnx {} {@var{p} =} csymamd (@var{S}, @var{knobs})
+@deftypefnx {} {@var{p} =} csymamd (@var{S}, @var{knobs}, @var{cmember})
+@deftypefnx {} {[@var{p}, @var{stats}] =} csymamd (@dots{})
+
+For a symmetric positive definite matrix @var{S}, return the permutation
+vector @var{p} such that @code{@var{S}(@var{p},@var{p})} tends to have a
+sparser Cholesky@tie{}factor than @var{S}.
+
+Sometimes @code{csymamd} works well for symmetric indefinite matrices too.
+The matrix @var{S} is assumed to be symmetric; only the strictly lower
+triangular part is referenced.  @var{S} must be square.  The ordering is
+followed by an elimination tree post-ordering.
+
+@var{knobs} is an optional 1-element to 3-element input vector, with a
+default value of @code{[10 1 0]}.  Entries not present are set to their
+defaults.
+
+@table @code
+@item @var{knobs}(1)
+If @var{S} is n-by-n, then rows and columns with more than
+@code{max(16,@var{knobs}(1)*sqrt(n))} entries are ignored, and ordered
+last in the output permutation (subject to the cmember constraints).
+
+@item @var{knobs}(2)
+If nonzero, aggressive absorption is performed.
+
+@item @var{knobs}(3)
+If nonzero, statistics and knobs are printed.
+
+@end table
+
+@var{cmember} is an optional vector of length n.  It defines the constraints
+on the ordering.  If @code{@var{cmember}(j) = @var{S}}, then row/column j is
+in constraint set @var{c} (@var{c} must be in the range 1 to n).  In the
+output permutation @var{p}, rows/columns in set 1 appear first, followed
+by all rows/columns in set 2, and so on.  @code{@var{cmember} = ones (1,n)}
+if not present or empty.  @code{csymamd (@var{S},[],1:n)} returns
+@code{1:n}.
+
+@code{@var{p} = csymamd (@var{S})} is about the same as
+@code{@var{p} = symamd (@var{S})}.  @var{knobs} and its default values
+differ.
+
+@code{@var{stats}(4:7)} provide information if CCOLAMD was able to
+continue.  The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if
+invalid.  @code{@var{stats}(5)} is the rightmost column index that is
+unsorted or contains duplicate entries, or zero if no such column exists.
+@code{@var{stats}(6)} is the last seen duplicate or out-of-order row
+index in the column index given by @code{@var{stats}(5)}, or zero if no
+such row index exists.  @code{@var{stats}(7)} is the number of duplicate
+or out-of-order row indices.  @code{@var{stats}(8:20)} is always zero in
+the current version of @sc{ccolamd} (reserved for future use).
+
+The authors of the code itself are @nospell{S. Larimore, T. Davis} and
+@nospell{S. Rajamanickam} in collaboration with @nospell{J. Bilbert and E. Ng}.
+Supported by the National Science Foundation
+@nospell{(DMS-9504974, DMS-9803599, CCR-0203270)}, and a grant from
+@nospell{Sandia} National Lab.
+See @url{http://faculty.cse.tamu.edu/davis/suitesparse.html} for ccolamd,
+colamd, csymamd, amd, colamd, symamd, and other related orderings.
+@seealso{symamd, ccolamd}
+@end deftypefn */)
+{
+#if defined (HAVE_CCOLAMD)
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 3)
+    print_usage ();
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+  int spumoni = 0;
+
+  // Get knobs
+  static_assert (CCOLAMD_KNOBS <= 40,
+                 "csymamd: # of CCOLAMD_KNOBS exceeded.  Please report this to bugs.octave.org");
+  double knob_storage[CCOLAMD_KNOBS];
+  double *knobs = &knob_storage[0];
+  CCOLAMD_NAME (_set_defaults) (knobs);
+
+  // Check for user-passed knobs
+  if (nargin > 1)
+    {
+      NDArray User_knobs = args(1).array_value ();
+      int nel_User_knobs = User_knobs.numel ();
+
+      if (nel_User_knobs > 0)
+        knobs[CCOLAMD_DENSE_ROW] = User_knobs(0);
+      if (nel_User_knobs > 1)
+        knobs[CCOLAMD_AGGRESSIVE] = User_knobs(1);
+      if (nel_User_knobs > 2)
+        spumoni = static_cast<int> (User_knobs(2));
+
+      // print knob settings if spumoni is set
+      if (spumoni)
+        {
+          octave_stdout << "\ncsymamd version " << CCOLAMD_MAIN_VERSION
+                        << '.' << CCOLAMD_SUB_VERSION
+                        << ", " << CCOLAMD_DATE << "\n";
+
+          if (knobs[CCOLAMD_DENSE_ROW] >= 0)
+            octave_stdout << "knobs(1): " << User_knobs(0)
+                          << ", rows/cols with > max (16,"
+                          << knobs[CCOLAMD_DENSE_ROW]
+                          << "*sqrt (size(A,2)))"
+                          << " entries removed\n";
+          else
+            octave_stdout << "knobs(1): " << User_knobs(0)
+                          << ", no dense rows/cols removed\n";
+
+          if (knobs[CCOLAMD_AGGRESSIVE] != 0)
+            octave_stdout << "knobs(2): " << User_knobs(1)
+                          << ", aggressive absorption: yes";
+          else
+            octave_stdout << "knobs(2): " << User_knobs(1)
+                          << ", aggressive absorption: no";
+
+          octave_stdout << "knobs(3): " << User_knobs(2)
+                        << ", statistics and knobs printed\n";
+        }
+    }
+
+  octave_idx_type n_row, n_col;
+  octave_idx_type *ridx, *cidx;
+  SparseMatrix sm;
+  SparseComplexMatrix scm;
+
+  if (args(0).issparse ())
+    {
+      if (args(0).iscomplex ())
+        {
+          scm = args(0).sparse_complex_matrix_value ();
+          n_row = scm.rows ();
+          n_col = scm.cols ();
+          ridx = scm.xridx ();
+          cidx = scm.xcidx ();
+        }
+      else
+        {
+          sm = args(0).sparse_matrix_value ();
+          n_row = sm.rows ();
+          n_col = sm.cols ();
+          ridx = sm.xridx ();
+          cidx = sm.xcidx ();
+        }
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+      else
+        sm = SparseMatrix (args(0).matrix_value ());
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      ridx = sm.xridx ();
+      cidx = sm.xcidx ();
+    }
+
+  if (n_row != n_col)
+    err_square_matrix_required ("csymamd", "S");
+
+  // Allocate workspace for symamd
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, perm, n_col+1);
+  static_assert (CCOLAMD_STATS <= 40,
+                 "csymamd: # of CCOLAMD_STATS exceeded.  Please report this to bugs.octave.org");
+  octave::suitesparse_integer stats_storage[CCOLAMD_STATS];
+  octave::suitesparse_integer *stats = &stats_storage[0];
+
+  if (nargin > 2)
+    {
+      NDArray in_cmember = args(2).array_value ();
+      octave_idx_type cslen = in_cmember.numel ();
+      OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, cmember, cslen);
+      for (octave_idx_type i = 0; i < cslen; i++)
+        // convert cmember from 1-based to 0-based
+        cmember[i] = static_cast<octave_idx_type> (in_cmember(i) - 1);
+
+      if (cslen != n_col)
+        error ("csymamd: CMEMBER must be of length equal to #cols of A");
+
+      if (! CSYMAMD_NAME () (n_col,
+                             octave::to_suitesparse_intptr (ridx),
+                             octave::to_suitesparse_intptr (cidx),
+                             perm, knobs, stats, &calloc, &free, cmember, -1))
+        {
+          CSYMAMD_NAME (_report)(stats);
+
+          error ("csymamd: internal error!");
+        }
+    }
+  else
+    {
+      if (! CSYMAMD_NAME () (n_col,
+                             octave::to_suitesparse_intptr (ridx),
+                             octave::to_suitesparse_intptr (cidx),
+                             perm, knobs, stats, &calloc, &free, nullptr, -1))
+        {
+          CSYMAMD_NAME (_report)(stats);
+
+          error ("csymamd: internal error!");
+        }
+    }
+
+  // return the permutation vector
+  NDArray out_perm (dim_vector (1, n_col));
+  for (octave_idx_type i = 0; i < n_col; i++)
+    out_perm(i) = perm[i] + 1;
+
+  retval(0) = out_perm;
+
+  // print stats if spumoni > 0
+  if (spumoni > 0)
+    CSYMAMD_NAME (_report)(stats);
+
+  // Return the stats vector
+  if (nargout == 2)
+    {
+      NDArray out_stats (dim_vector (1, CCOLAMD_STATS));
+      for (octave_idx_type i = 0 ; i < CCOLAMD_STATS ; i++)
+        out_stats(i) = stats[i];
+      retval(1) = out_stats;
+
+      // fix stats (5) and (6), for 1-based information on
+      // jumbled matrix.  note that this correction doesn't
+      // occur if symamd returns FALSE
+      out_stats(CCOLAMD_INFO1)++;
+      out_stats(CCOLAMD_INFO2)++;
+    }
+
+  return retval;
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("csymamd", "CCOLAMD");
+
+#endif
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/chol.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,1331 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 1996-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <string>
+
+#include "Matrix.h"
+#include "chol.h"
+#include "oct-string.h"
+#include "sparse-chol.h"
+#include "sparse-util.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "ov.h"
+#include "ovl.h"
+
+template <typename CHOLT>
+static octave_value
+get_chol (const CHOLT& fact)
+{
+  return octave_value (fact.chol_matrix());
+}
+
+template <typename CHOLT>
+static octave_value
+get_chol_r (const CHOLT& fact)
+{
+  return octave_value (fact.chol_matrix (),
+                       MatrixType (MatrixType::Upper));
+}
+
+template <typename CHOLT>
+static octave_value
+get_chol_l (const CHOLT& fact)
+{
+  return octave_value (fact.chol_matrix ().transpose (),
+                       MatrixType (MatrixType::Lower));
+}
+
+DEFUN (chol, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{R} =} chol (@var{A})
+@deftypefnx {} {[@var{R}, @var{p}] =} chol (@var{A})
+@deftypefnx {} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{A})
+@deftypefnx {} {[@var{R}, @var{p}, @var{Q}] =} chol (@var{A}, "vector")
+@deftypefnx {} {[@var{L}, @dots{}] =} chol (@dots{}, "lower")
+@deftypefnx {} {[@var{R}, @dots{}] =} chol (@dots{}, "upper")
+@cindex Cholesky factorization
+Compute the upper Cholesky@tie{}factor, @var{R}, of the real symmetric
+or complex Hermitian positive definite matrix @var{A}.
+
+The upper Cholesky@tie{}factor @var{R} is computed by using the upper
+triangular part of matrix @var{A} and is defined by
+@tex
+$ R^T R = A $.
+@end tex
+@ifnottex
+
+@example
+@var{R}' * @var{R} = @var{A}.
+@end example
+
+@end ifnottex
+
+Calling @code{chol} using the optional @qcode{"upper"} flag has the
+same behavior.  In contrast, using the optional @qcode{"lower"} flag,
+@code{chol} returns the lower triangular factorization, computed by using
+the lower triangular part of matrix @var{A}, such that
+@tex
+$ L L^T = A $.
+@end tex
+@ifnottex
+
+@example
+@var{L} * @var{L}' = @var{A}.
+@end example
+
+@end ifnottex
+
+Called with one output argument @code{chol} fails if matrix @var{A} is
+not positive definite.  Note that if matrix @var{A} is not real symmetric
+or complex Hermitian then the lower triangular part is considered to be
+the (complex conjugate) transpose of the upper triangular part, or vice
+versa, given the @qcode{"lower"} flag.
+
+Called with two or more output arguments @var{p} flags whether the matrix
+@var{A} was positive definite and @code{chol} does not fail.  A zero value
+of @var{p} indicates that matrix @var{A} is positive definite and @var{R}
+gives the factorization.  Otherwise, @var{p} will have a positive value.
+
+If called with three output arguments matrix @var{A} must be sparse and
+a sparsity preserving row/column permutation is applied to matrix @var{A}
+prior to the factorization.  That is @var{R} is the factorization of
+@code{@var{A}(@var{Q},@var{Q})} such that
+@tex
+$ R^T R = Q^T A Q$.
+@end tex
+@ifnottex
+
+@example
+@var{R}' * @var{R} = @var{Q}' * @var{A} * @var{Q}.
+@end example
+
+@end ifnottex
+
+The sparsity preserving permutation is generally returned as a matrix.
+However, given the optional flag @qcode{"vector"}, @var{Q} will be
+returned as a vector such that
+@tex
+$ R^T R = A (Q, Q)$.
+@end tex
+@ifnottex
+
+@example
+@var{R}' * @var{R} = @var{A}(@var{Q}, @var{Q}).
+@end example
+
+@end ifnottex
+
+In general the lower triangular factorization is significantly faster for
+sparse matrices.
+@seealso{hess, lu, qr, qz, schur, svd, ichol, cholinv, chol2inv, cholupdate, cholinsert, choldelete, cholshift}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 3 || nargout > 3)
+    print_usage ();
+  if (nargout > 2 && ! args(0).issparse ())
+    error ("chol: using three output arguments, matrix A must be sparse");
+
+  bool LLt = false;
+  bool vecout = false;
+
+  int n = 1;
+  while (n < nargin)
+    {
+      std::string tmp = args(n++).xstring_value ("chol: optional arguments must be strings");
+
+      if (octave::string::strcmpi (tmp, "vector"))
+        vecout = true;
+      else if (octave::string::strcmpi (tmp, "lower"))
+        LLt = true;
+      else if (octave::string::strcmpi (tmp, "upper"))
+        LLt = false;
+      else
+        error (R"(chol: optional argument must be one of "vector", "lower", or "upper")");
+    }
+
+  octave_value_list retval;
+  octave_value arg = args(0);
+
+  if (arg.isempty ())
+    return ovl (Matrix ());
+
+  if (arg.issparse ())
+    {
+      octave_idx_type info;
+      bool natural = (nargout != 3);
+      bool force = nargout > 1;
+
+      if (arg.isreal ())
+        {
+          SparseMatrix m = arg.sparse_matrix_value ();
+
+          octave::math::sparse_chol<SparseMatrix> fact (m, info, natural, force);
+
+          if (nargout == 3)
+            {
+              if (vecout)
+                retval(2) = fact.perm ();
+              else
+                retval(2) = fact.Q ();
+            }
+
+          if (nargout >= 2 || info == 0)
+            {
+              retval(1) = info;
+              if (LLt)
+                retval(0) = fact.L ();
+              else
+                retval(0) = fact.R ();
+            }
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else if (arg.iscomplex ())
+        {
+          SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
+
+          octave::math::sparse_chol<SparseComplexMatrix> fact (m, info, natural, force);
+
+          if (nargout == 3)
+            {
+              if (vecout)
+                retval(2) = fact.perm ();
+              else
+                retval(2) = fact.Q ();
+            }
+
+          if (nargout >= 2 || info == 0)
+            {
+              retval(1) = info;
+              if (LLt)
+                retval(0) = fact.L ();
+              else
+                retval(0) = fact.R ();
+            }
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else
+        err_wrong_type_arg ("chol", arg);
+    }
+  else if (arg.is_single_type ())
+    {
+      if (vecout)
+        error (R"(chol: A must be sparse for the "vector" option)");
+      if (arg.isreal ())
+        {
+          FloatMatrix m = arg.float_matrix_value ();
+
+          octave_idx_type info;
+
+          octave::math::chol<FloatMatrix> fact (m, info, LLt != true);
+
+          if (nargout == 2 || info == 0)
+            retval = ovl (get_chol (fact), info);
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else if (arg.iscomplex ())
+        {
+          FloatComplexMatrix m = arg.float_complex_matrix_value ();
+
+          octave_idx_type info;
+
+          octave::math::chol<FloatComplexMatrix> fact (m, info, LLt != true);
+
+          if (nargout == 2 || info == 0)
+            retval = ovl (get_chol (fact), info);
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else
+        err_wrong_type_arg ("chol", arg);
+    }
+  else
+    {
+      if (vecout)
+        error (R"(chol: A must be sparse for the "vector" option)");
+      if (arg.isreal ())
+        {
+          Matrix m = arg.matrix_value ();
+
+          octave_idx_type info;
+
+          octave::math::chol<Matrix> fact (m, info, LLt != true);
+
+          if (nargout == 2 || info == 0)
+            retval = ovl (get_chol (fact), info);
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else if (arg.iscomplex ())
+        {
+          ComplexMatrix m = arg.complex_matrix_value ();
+
+          octave_idx_type info;
+
+          octave::math::chol<ComplexMatrix> fact (m, info, LLt != true);
+
+          if (nargout == 2 || info == 0)
+            retval = ovl (get_chol (fact), info);
+          else
+            error ("chol: input matrix must be positive definite");
+        }
+      else
+        err_wrong_type_arg ("chol", arg);
+    }
+
+  return retval;
+}
+
+/*
+%!assert (chol ([2, 1; 1, 1]), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)], sqrt (eps))
+%!assert (chol (single ([2, 1; 1, 1])), single ([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]), sqrt (eps ("single")))
+
+%!assert (chol ([2, 1; 1, 1], "upper"), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)],
+%!        sqrt (eps))
+%!assert (chol ([2, 1; 1, 1], "lower"), [sqrt(2), 0; 1/sqrt(2), 1/sqrt(2)],
+%!        sqrt (eps))
+
+%!assert (chol ([2, 1; 1, 1], "lower"), chol ([2, 1; 1, 1], "LoweR"))
+%!assert (chol ([2, 1; 1, 1], "upper"), chol ([2, 1; 1, 1], "Upper"))
+
+## Check the "vector" option which only affects the 3rd argument and
+## is only valid for sparse input.
+%!testif HAVE_CHOLMOD
+%! a = sparse ([2 1; 1 1]);
+%! r = sparse ([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]);
+%! [rd, pd, qd] = chol (a);
+%! [rv, pv, qv] = chol (a, "vector");
+%! assert (r, rd, eps)
+%! assert (r, rv, eps)
+%! assert (pd, 0)
+%! assert (pd, pv)
+%! assert (qd, sparse (eye (2)))
+%! assert (qv, [1 2])
+%!
+%! [rv, pv, qv] = chol (a, "Vector"); # check case sensitivity
+%! assert (r, rv, eps)
+%! assert (pd, pv)
+%! assert (qv, [1 2])
+
+%!testif HAVE_CHOLMOD <*42587>
+%! A = sparse ([1 0 8;0 1 8;8 8 1]);
+%! [Q, p] = chol (A);
+%! assert (p != 0);
+
+%!error chol ()
+%!error <matrix must be positive definite> chol ([1, 2; 3, 4])
+%!error <requires square matrix> chol ([1, 2; 3, 4; 5, 6])
+%!error <optional arguments must be strings> chol (1, 2)
+%!error <optional argument must be one of "vector", "lower"> chol (1, "foobar")
+%!error <matrix A must be sparse> [L,p,Q] = chol ([1, 2; 3, 4])
+%!error <A must be sparse> [L, p] = chol ([1, 2; 3, 4], "vector")
+*/
+
+DEFUN (cholinv, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {} cholinv (@var{A})
+Compute the inverse of the symmetric positive definite matrix @var{A} using
+the Cholesky@tie{}factorization.
+@seealso{chol, chol2inv, inv}
+@end deftypefn */)
+{
+  if (args.length () != 1)
+    print_usage ();
+
+  octave_value retval;
+  octave_value arg = args(0);
+
+  octave_idx_type nr = arg.rows ();
+  octave_idx_type nc = arg.columns ();
+
+  if (nr == 0 || nc == 0)
+    retval = Matrix ();
+  else
+    {
+      if (arg.issparse ())
+        {
+          octave_idx_type info;
+
+          if (arg.isreal ())
+            {
+              SparseMatrix m = arg.sparse_matrix_value ();
+
+              octave::math::sparse_chol<SparseMatrix> chol (m, info);
+
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else if (arg.iscomplex ())
+            {
+              SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
+
+              octave::math::sparse_chol<SparseComplexMatrix> chol (m, info);
+
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else
+            err_wrong_type_arg ("cholinv", arg);
+        }
+      else if (arg.is_single_type ())
+        {
+          if (arg.isreal ())
+            {
+              FloatMatrix m = arg.float_matrix_value ();
+
+              octave_idx_type info;
+              octave::math::chol<FloatMatrix> chol (m, info);
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else if (arg.iscomplex ())
+            {
+              FloatComplexMatrix m = arg.float_complex_matrix_value ();
+
+              octave_idx_type info;
+              octave::math::chol<FloatComplexMatrix> chol (m, info);
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else
+            err_wrong_type_arg ("chol", arg);
+        }
+      else
+        {
+          if (arg.isreal ())
+            {
+              Matrix m = arg.matrix_value ();
+
+              octave_idx_type info;
+              octave::math::chol<Matrix> chol (m, info);
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else if (arg.iscomplex ())
+            {
+              ComplexMatrix m = arg.complex_matrix_value ();
+
+              octave_idx_type info;
+              octave::math::chol<ComplexMatrix> chol (m, info);
+              if (info == 0)
+                retval = chol.inverse ();
+              else
+                error ("cholinv: A must be positive definite");
+            }
+          else
+            err_wrong_type_arg ("chol", arg);
+        }
+    }
+
+  return retval;
+}
+
+/*
+%!shared A, Ainv
+%! A = [2,0.2;0.2,1];
+%! Ainv = inv (A);
+%!test
+%! Ainv1 = cholinv (A);
+%! assert (norm (Ainv-Ainv1), 0, 1e-10);
+%!testif HAVE_CHOLMOD
+%! Ainv2 = inv (sparse (A));
+%! assert (norm (Ainv-Ainv2), 0, 1e-10);
+%!testif HAVE_CHOLMOD
+%! Ainv3 = cholinv (sparse (A));
+%! assert (norm (Ainv-Ainv3), 0, 1e-10);
+*/
+
+DEFUN (chol2inv, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {} chol2inv (@var{U})
+Invert a symmetric, positive definite square matrix from its Cholesky
+decomposition, @var{U}.
+
+Note that @var{U} should be an upper-triangular matrix with positive
+diagonal elements.  @code{chol2inv (@var{U})} provides
+@code{inv (@var{U}'*@var{U})} but it is much faster than using @code{inv}.
+@seealso{chol, cholinv, inv}
+@end deftypefn */)
+{
+  if (args.length () != 1)
+    print_usage ();
+
+  octave_value retval;
+
+  octave_value arg = args(0);
+
+  octave_idx_type nr = arg.rows ();
+  octave_idx_type nc = arg.columns ();
+
+  if (nr == 0 || nc == 0)
+    retval = Matrix ();
+  else
+    {
+      if (arg.issparse ())
+        {
+          if (arg.isreal ())
+            {
+              SparseMatrix r = arg.sparse_matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else if (arg.iscomplex ())
+            {
+              SparseComplexMatrix r = arg.sparse_complex_matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else
+            err_wrong_type_arg ("chol2inv", arg);
+        }
+      else if (arg.is_single_type ())
+        {
+          if (arg.isreal ())
+            {
+              FloatMatrix r = arg.float_matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else if (arg.iscomplex ())
+            {
+              FloatComplexMatrix r = arg.float_complex_matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else
+            err_wrong_type_arg ("chol2inv", arg);
+
+        }
+      else
+        {
+          if (arg.isreal ())
+            {
+              Matrix r = arg.matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else if (arg.iscomplex ())
+            {
+              ComplexMatrix r = arg.complex_matrix_value ();
+
+              retval = octave::math::chol2inv (r);
+            }
+          else
+            err_wrong_type_arg ("chol2inv", arg);
+        }
+    }
+
+  return retval;
+}
+
+/*
+
+## Test for bug #36437
+%!function sparse_chol2inv (T, tol)
+%!  iT = inv (T);
+%!  ciT = chol2inv (chol (T));
+%!  assert (ciT, iT, tol);
+%!  assert (chol2inv (chol ( full (T))), ciT, tol*2);
+%!endfunction
+
+%!testif HAVE_CHOLMOD
+%! A = gallery ("poisson", 3);
+%! sparse_chol2inv (A, eps);
+
+%!testif HAVE_CHOLMOD
+%! n = 10;
+%! B = spdiags (ones (n, 1) * [1 2 1], [-1 0 1], n, n);
+%! sparse_chol2inv (B, eps*100);
+
+%!testif HAVE_CHOLMOD
+%! C = gallery("tridiag", 5);
+%! sparse_chol2inv (C, eps*10);
+
+%!testif HAVE_CHOLMOD
+%! D = gallery("wathen", 1, 1);
+%! sparse_chol2inv (D, eps*10^4);
+
+*/
+
+DEFUN (cholupdate, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})
+Update or downdate a Cholesky@tie{}factorization.
+
+Given an upper triangular matrix @var{R} and a column vector @var{u},
+attempt to determine another upper triangular matrix @var{R1} such that
+
+@itemize @bullet
+@item
+@var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'
+if @var{op} is @qcode{"+"}
+
+@item
+@var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'
+if @var{op} is @qcode{"-"}
+@end itemize
+
+If @var{op} is @qcode{"-"}, @var{info} is set to
+
+@itemize
+@item 0 if the downdate was successful,
+
+@item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,
+
+@item 2 if @var{R} is singular.
+@end itemize
+
+If @var{info} is not present, an error message is printed in cases 1 and 2.
+@seealso{chol, cholinsert, choldelete, cholshift}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  if (nargin < 2 || nargin > 3)
+    print_usage ();
+
+  octave_value argr = args(0);
+  octave_value argu = args(1);
+
+  if (! argr.isnumeric () || ! argu.isnumeric ()
+      || (nargin > 2 && ! args(2).is_string ()))
+    print_usage ();
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+
+  octave_idx_type n = argr.rows ();
+
+  std::string op = (nargin < 3) ? "+" : args(2).string_value ();
+
+  bool down = (op == "-");
+
+  if (! down && op != "+")
+    error (R"(cholupdate: OP must be "+" or "-")");
+
+  if (argr.columns () != n || argu.rows () != n || argu.columns () != 1)
+    error ("cholupdate: dimension mismatch between R and U");
+
+  int err = 0;
+  if (argr.is_single_type () || argu.is_single_type ())
+    {
+      if (argr.isreal () && argu.isreal ())
+        {
+          // real case
+          FloatMatrix R = argr.float_matrix_value ();
+          FloatColumnVector u = argu.float_column_vector_value ();
+
+          octave::math::chol<FloatMatrix> fact;
+          fact.set (R);
+
+          if (down)
+            err = fact.downdate (u);
+          else
+            fact.update (u);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          FloatComplexMatrix R = argr.float_complex_matrix_value ();
+          FloatComplexColumnVector u
+            = argu.float_complex_column_vector_value ();
+
+          octave::math::chol<FloatComplexMatrix> fact;
+          fact.set (R);
+
+          if (down)
+            err = fact.downdate (u);
+          else
+            fact.update (u);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+  else
+    {
+      if (argr.isreal () && argu.isreal ())
+        {
+          // real case
+          Matrix R = argr.matrix_value ();
+          ColumnVector u = argu.column_vector_value ();
+
+          octave::math::chol<Matrix> fact;
+          fact.set (R);
+
+          if (down)
+            err = fact.downdate (u);
+          else
+            fact.update (u);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          ComplexMatrix R = argr.complex_matrix_value ();
+          ComplexColumnVector u = argu.complex_column_vector_value ();
+
+          octave::math::chol<ComplexMatrix> fact;
+          fact.set (R);
+
+          if (down)
+            err = fact.downdate (u);
+          else
+            fact.update (u);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+
+  if (nargout > 1)
+    retval(1) = err;
+  else if (err == 1)
+    error ("cholupdate: downdate violates positiveness");
+  else if (err == 2)
+    error ("cholupdate: singular matrix");
+
+  return retval;
+}
+
+/*
+%!shared A, u, Ac, uc
+%! A = [  0.436997  -0.131721   0.124120  -0.061673 ;
+%!       -0.131721   0.738529   0.019851  -0.140295 ;
+%!        0.124120   0.019851   0.354879  -0.059472 ;
+%!       -0.061673  -0.140295  -0.059472   0.600939 ];
+%!
+%! u = [  0.98950 ;
+%!        0.39844 ;
+%!        0.63484 ;
+%!        0.13351 ];
+%! Ac = [  0.5585528 + 0.0000000i  -0.1662088 - 0.0315341i   0.0107873 + 0.0236411i  -0.0276775 - 0.0186073i ;
+%!        -0.1662088 + 0.0315341i   0.6760061 + 0.0000000i   0.0011452 - 0.0475528i   0.0145967 + 0.0247641i ;
+%!         0.0107873 - 0.0236411i   0.0011452 + 0.0475528i   0.6263149 - 0.0000000i  -0.1585837 - 0.0719763i ;
+%!        -0.0276775 + 0.0186073i   0.0145967 - 0.0247641i  -0.1585837 + 0.0719763i   0.6034234 - 0.0000000i ];
+%!
+%! uc = [ 0.54267 + 0.91519i ;
+%!        0.99647 + 0.43141i ;
+%!        0.83760 + 0.68977i ;
+%!        0.39160 + 0.90378i ];
+
+%!test
+%! R = chol (A);
+%! R1 = cholupdate (R, u);
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - R'*R - u*u', Inf) < 1e1*eps);
+%!
+%! R1 = cholupdate (R1, u, "-");
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1 - R, Inf) < 1e1*eps);
+
+%!test
+%! R = chol (Ac);
+%! R1 = cholupdate (R, uc);
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - R'*R - uc*uc', Inf) < 1e1*eps);
+%!
+%! R1 = cholupdate (R1, uc, "-");
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1 - R, Inf) < 1e1*eps);
+
+%!test
+%! R = chol (single (A));
+%! R1 = cholupdate (R, single (u));
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1'*R1 - R'*R - single (u*u'), Inf) < 1e1*eps ("single"));
+%!
+%! R1 = cholupdate (R1, single (u), "-");
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1 - R, Inf) < 2e1*eps ("single"));
+
+%!test
+%! R = chol (single (Ac));
+%! R1 = cholupdate (R, single (uc));
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1'*R1 - R'*R - single (uc*uc'), Inf) < 1e1*eps ("single"));
+%!
+%! R1 = cholupdate (R1, single (uc), "-");
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1 - R, Inf) < 2e1*eps ("single"));
+*/
+
+DEFUN (cholinsert, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{R1} =} cholinsert (@var{R}, @var{j}, @var{u})
+@deftypefnx {} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})
+Update a Cholesky factorization given a row or column to insert in the
+original factored matrix.
+
+Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian
+positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper
+triangular, return the Cholesky@tie{}factorization of
+@var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and
+@w{p = [1:j-1,j+1:n+1]}.  @w{u(j)} should be positive.
+
+On return, @var{info} is set to
+
+@itemize
+@item 0 if the insertion was successful,
+
+@item 1 if @var{A1} is not positive definite,
+
+@item 2 if @var{R} is singular.
+@end itemize
+
+If @var{info} is not present, an error message is printed in cases 1 and 2.
+@seealso{chol, cholupdate, choldelete, cholshift}
+@end deftypefn */)
+{
+  if (args.length () != 3)
+    print_usage ();
+
+  octave_value argr = args(0);
+  octave_value argj = args(1);
+  octave_value argu = args(2);
+
+  if (! argr.isnumeric () || ! argu.isnumeric ()
+      || ! argj.is_real_scalar ())
+    print_usage ();
+
+  octave_idx_type n = argr.rows ();
+  octave_idx_type j = argj.scalar_value ();
+
+  if (argr.columns () != n || argu.rows () != n+1 || argu.columns () != 1)
+    error ("cholinsert: dimension mismatch between R and U");
+
+  if (j < 1 || j > n+1)
+    error ("cholinsert: index J out of range");
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+
+  int err = 0;
+  if (argr.is_single_type () || argu.is_single_type ())
+    {
+      if (argr.isreal () && argu.isreal ())
+        {
+          // real case
+          FloatMatrix R = argr.float_matrix_value ();
+          FloatColumnVector u = argu.float_column_vector_value ();
+
+          octave::math::chol<FloatMatrix> fact;
+          fact.set (R);
+          err = fact.insert_sym (u, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          FloatComplexMatrix R = argr.float_complex_matrix_value ();
+          FloatComplexColumnVector u
+            = argu.float_complex_column_vector_value ();
+
+          octave::math::chol<FloatComplexMatrix> fact;
+          fact.set (R);
+          err = fact.insert_sym (u, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+  else
+    {
+      if (argr.isreal () && argu.isreal ())
+        {
+          // real case
+          Matrix R = argr.matrix_value ();
+          ColumnVector u = argu.column_vector_value ();
+
+          octave::math::chol<Matrix> fact;
+          fact.set (R);
+          err = fact.insert_sym (u, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          ComplexMatrix R = argr.complex_matrix_value ();
+          ComplexColumnVector u = argu.complex_column_vector_value ();
+
+          octave::math::chol<ComplexMatrix> fact;
+          fact.set (R);
+          err = fact.insert_sym (u, j-1);
+
+          retval = ovl (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");
+
+  return retval;
+}
+
+/*
+%!test
+%! u2 = [  0.35080 ;
+%!         0.63930 ;
+%!         3.31057 ;
+%!        -0.13825 ;
+%!         0.45266 ];
+%!
+%! R = chol (A);
+%!
+%! j = 3;  p = [1:j-1, j+1:5];
+%! R1 = cholinsert (R, j, u2);
+%! A1 = R1'*R1;
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (A1(p,p) - A, Inf) < 1e1*eps);
+
+%!test
+%! u2 = [  0.35080  + 0.04298i;
+%!         0.63930  + 0.23778i;
+%!         3.31057  + 0.00000i;
+%!        -0.13825  + 0.19879i;
+%!         0.45266  + 0.50020i];
+%!
+%! R = chol (Ac);
+%!
+%! j = 3;  p = [1:j-1, j+1:5];
+%! R1 = cholinsert (R, j, u2);
+%! A1 = R1'*R1;
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (A1(p,p) - Ac, Inf) < 1e1*eps);
+
+%!test
+%! u2 = single ([  0.35080 ;
+%!                 0.63930 ;
+%!                 3.31057 ;
+%!                -0.13825 ;
+%!                 0.45266 ]);
+%!
+%! R = chol (single (A));
+%!
+%! j = 3;  p = [1:j-1, j+1:5];
+%! R1 = cholinsert (R, j, u2);
+%! A1 = R1'*R1;
+%!
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (A1(p,p) - A, Inf) < 1e1*eps ("single"));
+
+%!test
+%! u2 = single ([  0.35080  + 0.04298i;
+%!                 0.63930  + 0.23778i;
+%!                 3.31057  + 0.00000i;
+%!                -0.13825  + 0.19879i;
+%!                 0.45266  + 0.50020i]);
+%!
+%! R = chol (single (Ac));
+%!
+%! j = 3;  p = [1:j-1, j+1:5];
+%! R1 = cholinsert (R, j, u2);
+%! A1 = R1'*R1;
+%!
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (A1(p,p) - single (Ac), Inf) < 2e1*eps ("single"));
+
+%!test
+%! cu = chol (triu (A), "upper");
+%! cl = chol (tril (A), "lower");
+%! assert (cu, cl', eps);
+
+%!test
+%! cca   = chol (Ac);
+%!
+%! ccal  = chol (Ac, "lower");
+%! ccal2 = chol (tril (Ac), "lower");
+%!
+%! ccau  = chol (Ac, "upper");
+%! ccau2 = chol (triu (Ac), "upper");
+%!
+%! assert (cca'*cca,     Ac, eps);
+%! assert (ccau'*ccau,   Ac, eps);
+%! assert (ccau2'*ccau2, Ac, eps);
+%!
+%! assert (cca, ccal',  eps);
+%! assert (cca, ccau,   eps);
+%! assert (cca, ccal2', eps);
+%! assert (cca, ccau2,  eps);
+
+%!test
+%! cca   = chol (single (Ac));
+%!
+%! ccal  = chol (single (Ac), "lower");
+%! ccal2 = chol (tril (single (Ac)), "lower");
+%!
+%! ccau  = chol (single (Ac), "upper");
+%! ccau2 = chol (triu (single (Ac)), "upper");
+%!
+%! assert (cca'*cca,     single (Ac), eps ("single"));
+%! assert (ccau'*ccau,   single (Ac), eps ("single"));
+%! assert (ccau2'*ccau2, single (Ac), eps ("single"));
+%!
+%! assert (cca, ccal',  eps ("single"));
+%! assert (cca, ccau,   eps ("single"));
+%! assert (cca, ccal2', eps ("single"));
+%! assert (cca, ccau2,  eps ("single"));
+
+%!test
+%! a = [12,  2,  3,  4;
+%!       2, 14,  5,  3;
+%!       3,  5, 16,  6;
+%!       4,  3,  6, 16];
+%!
+%! b = [0,  1,  2,  3;
+%!     -1,  0,  1,  2;
+%!     -2, -1,  0,  1;
+%!     -3, -2, -1,  0];
+%!
+%! ca = a + i*b;
+%!
+%! cca   = chol (ca);
+%!
+%! ccal  = chol (ca, "lower");
+%! ccal2 = chol (tril (ca), "lower");
+%!
+%! ccau  = chol (ca, "upper");
+%! ccau2 = chol (triu (ca), "upper");
+%!
+%! assert (cca'*cca,     ca, 16*eps);
+%! assert (ccau'*ccau,   ca, 16*eps);
+%! assert (ccau2'*ccau2, ca, 16*eps);
+%!
+%! assert (cca, ccal',  16*eps);
+%! assert (cca, ccau,   16*eps);
+%! assert (cca, ccal2', 16*eps);
+%! assert (cca, ccau2,  16*eps);
+*/
+
+DEFUN (choldelete, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {@var{R1} =} choldelete (@var{R}, @var{j})
+Update a Cholesky factorization given a row or column to delete from the
+original factored matrix.
+
+Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian
+positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper
+triangular, return the Cholesky@tie{}factorization of @w{A(p,p)}, where
+@w{p = [1:j-1,j+1:n+1]}.
+@seealso{chol, cholupdate, cholinsert, cholshift}
+@end deftypefn */)
+{
+  if (args.length () != 2)
+    print_usage ();
+
+  octave_value argr = args(0);
+  octave_value argj = args(1);
+
+  if (! argr.isnumeric () || ! argj.is_real_scalar ())
+    print_usage ();
+
+  octave_idx_type n = argr.rows ();
+  octave_idx_type j = argj.scalar_value ();
+
+  if (argr.columns () != n)
+    err_square_matrix_required ("choldelete", "R");
+
+  if (j < 0 && j > n)
+    error ("choldelete: index J out of range");
+
+  octave_value_list retval;
+
+  if (argr.is_single_type ())
+    {
+      if (argr.isreal ())
+        {
+          // real case
+          FloatMatrix R = argr.float_matrix_value ();
+
+          octave::math::chol<FloatMatrix> fact;
+          fact.set (R);
+          fact.delete_sym (j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          FloatComplexMatrix R = argr.float_complex_matrix_value ();
+
+          octave::math::chol<FloatComplexMatrix> fact;
+          fact.set (R);
+          fact.delete_sym (j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+  else
+    {
+      if (argr.isreal ())
+        {
+          // real case
+          Matrix R = argr.matrix_value ();
+
+          octave::math::chol<Matrix> fact;
+          fact.set (R);
+          fact.delete_sym (j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          ComplexMatrix R = argr.complex_matrix_value ();
+
+          octave::math::chol<ComplexMatrix> fact;
+          fact.set (R);
+          fact.delete_sym (j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+
+  return retval;
+}
+
+/*
+%!test
+%! R = chol (A);
+%!
+%! j = 3;  p = [1:j-1,j+1:4];
+%! R1 = choldelete (R, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
+
+%!test
+%! R = chol (Ac);
+%!
+%! j = 3;  p = [1:j-1,j+1:4];
+%! R1 = choldelete (R, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
+
+%!test
+%! R = chol (single (A));
+%!
+%! j = 3;  p = [1:j-1,j+1:4];
+%! R1 = choldelete (R, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
+
+%!test
+%! R = chol (single (Ac));
+%!
+%! j = 3;  p = [1:j-1,j+1:4];
+%! R1 = choldelete (R,j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), single (0));
+%! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
+*/
+
+DEFUN (cholshift, args, ,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})
+Update a Cholesky factorization given a range of columns to shift in the
+original factored matrix.
+
+Given a Cholesky@tie{}factorization of a real symmetric or complex Hermitian
+positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper
+triangular, return the Cholesky@tie{}factorization of
+@w{@var{A}(p,p)}, where @w{p} is the permutation @*
+@code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*
+ or @*
+@code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}.  @*
+
+@seealso{chol, cholupdate, cholinsert, choldelete}
+@end deftypefn */)
+{
+  if (args.length () != 3)
+    print_usage ();
+
+  octave_value argr = args(0);
+  octave_value argi = args(1);
+  octave_value argj = args(2);
+
+  if (! argr.isnumeric () || ! argi.is_real_scalar ()
+      || ! argj.is_real_scalar ())
+    print_usage ();
+
+  octave_idx_type n = argr.rows ();
+  octave_idx_type i = argi.scalar_value ();
+  octave_idx_type j = argj.scalar_value ();
+
+  if (argr.columns () != n)
+    err_square_matrix_required ("cholshift", "R");
+
+  if (j < 0 || j > n+1 || i < 0 || i > n+1)
+    error ("cholshift: index I or J is out of range");
+
+  octave_value_list retval;
+
+  if (argr.is_single_type () && argi.is_single_type ()
+      && argj.is_single_type ())
+    {
+      if (argr.isreal ())
+        {
+          // real case
+          FloatMatrix R = argr.float_matrix_value ();
+
+          octave::math::chol<FloatMatrix> fact;
+          fact.set (R);
+          fact.shift_sym (i-1, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          FloatComplexMatrix R = argr.float_complex_matrix_value ();
+
+          octave::math::chol<FloatComplexMatrix> fact;
+          fact.set (R);
+          fact.shift_sym (i-1, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+  else
+    {
+      if (argr.isreal ())
+        {
+          // real case
+          Matrix R = argr.matrix_value ();
+
+          octave::math::chol<Matrix> fact;
+          fact.set (R);
+          fact.shift_sym (i-1, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+      else
+        {
+          // complex case
+          ComplexMatrix R = argr.complex_matrix_value ();
+
+          octave::math::chol<ComplexMatrix> fact;
+          fact.set (R);
+          fact.shift_sym (i-1, j-1);
+
+          retval = ovl (get_chol_r (fact));
+        }
+    }
+
+  return retval;
+}
+
+/*
+%!test
+%! R = chol (A);
+%!
+%! i = 1;  j = 3;  p = [1:i-1, shift(i:j,-1), j+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
+%!
+%! j = 1;  i = 3;  p = [1:j-1, shift(j:i,+1), i+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1) - R1, Inf), 0);
+%! assert (norm (R1'*R1 - A(p,p), Inf) < 1e1*eps);
+
+%!test
+%! R = chol (Ac);
+%!
+%! i = 1;  j = 3;  p = [1:i-1, shift(i:j,-1), j+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
+%!
+%! j = 1;  i = 3;  p = [1:j-1, shift(j:i,+1), i+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - Ac(p,p), Inf) < 1e1*eps);
+
+%!test
+%! R = chol (single (A));
+%!
+%! i = 1;  j = 3;  p = [1:i-1, shift(i:j,-1), j+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
+%!
+%! j = 1;  i = 3;  p = [1:j-1, shift(j:i,+1), i+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - single (A(p,p)), Inf) < 1e1*eps ("single"));
+
+%!test
+%! R = chol (single (Ac));
+%!
+%! i = 1;  j = 3;  p = [1:i-1, shift(i:j,-1), j+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
+%!
+%! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
+%! R1 = cholshift (R, i, j);
+%!
+%! assert (norm (triu (R1)-R1, Inf), 0);
+%! assert (norm (R1'*R1 - single (Ac(p,p)), Inf) < 1e1*eps ("single"));
+*/

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/colamd.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,774 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 1998-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+// This is the octave interface to colamd, which bore the copyright given
+// in the help of the functions.
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <cstdlib>
+
+#include <string>
+
+#include "CSparse.h"
+#include "dNDArray.h"
+#include "dSparse.h"
+#include "oct-locbuf.h"
+#include "oct-sparse.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "ovl.h"
+#include "pager.h"
+
+// The symmetric column elimination tree code take from the Davis LDL code.
+// Copyright given elsewhere in this file.
+static void
+symetree (const octave_idx_type *ridx, const octave_idx_type *cidx,
+          octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
+{
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
+  if (P)
+    // If P is present then compute Pinv, the inverse of P
+    for (octave_idx_type k = 0 ; k < n ; k++)
+      Pinv[P[k]] = k;
+
+  for (octave_idx_type k = 0 ; k < n ; k++)
+    {
+      // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k)
+      Parent[k] = n ;                // parent of k is not yet known
+      Flag[k] = k ;                  // mark node k as visited
+      octave_idx_type kk = (P ? P[k] : k); // kth original, or permuted, column
+      octave_idx_type p2 = cidx[kk+1];
+      for (octave_idx_type p = cidx[kk] ; p < p2 ; p++)
+        {
+          // A (i,k) is nonzero (original or permuted A)
+          octave_idx_type i = (P ? Pinv[ridx[p]] : ridx[p]);
+          if (i < k)
+            {
+              // follow path from i to root of etree, stop at flagged node
+              for ( ; Flag[i] != k ; i = Parent[i])
+                {
+                  // find parent of i if not yet determined
+                  if (Parent[i] == n)
+                    Parent[i] = k;
+                  Flag[i] = k ;        // mark i as visited
+                }
+            }
+        }
+    }
+}
+
+// The elimination tree post-ordering code below is taken from SuperLU
+static inline octave_idx_type
+make_set (octave_idx_type i, octave_idx_type *pp)
+{
+  pp[i] = i;
+  return i;
+}
+
+static inline octave_idx_type
+link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
+{
+  pp[s] = t;
+  return t;
+}
+
+static inline octave_idx_type
+find (octave_idx_type i, octave_idx_type *pp)
+{
+  octave_idx_type p = pp[i];
+  octave_idx_type gp = pp[p];
+
+  while (gp != p)
+    {
+      pp[i] = gp;
+      i = gp;
+      p = pp[i];
+      gp = pp[p];
+    }
+
+  return p;
+}
+
+static octave_idx_type
+etdfs (octave_idx_type v, octave_idx_type *first_kid,
+       octave_idx_type *next_kid, octave_idx_type *post,
+       octave_idx_type postnum)
+{
+  for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w])
+    postnum = etdfs (w, first_kid, next_kid, post, postnum);
+
+  post[postnum++] = v;
+
+  return postnum;
+}
+
+static void
+tree_postorder (octave_idx_type n, octave_idx_type *parent,
+                octave_idx_type *post)
+{
+  // Allocate storage for working arrays and results
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
+
+  // Set up structure describing children
+  for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1)
+    ; // do nothing
+
+  for (octave_idx_type v = n-1; v >= 0; v--)
+    {
+      octave_idx_type dad = parent[v];
+      next_kid[v] = first_kid[dad];
+      first_kid[dad] = v;
+    }
+
+  // Depth-first search from dummy root vertex #n
+  etdfs (n, first_kid, next_kid, post, 0);
+}
+
+static void
+coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
+          octave_idx_type *colend, octave_idx_type *parent,
+          octave_idx_type nr, octave_idx_type nc)
+{
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
+
+  // Compute firstcol[row] = first nonzero column in row
+  for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc)
+    ; // do nothing
+
+  for (octave_idx_type col = 0; col < nc; col++)
+    for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
+      {
+        octave_idx_type row = ridx[p];
+        if (firstcol[row] > col)
+          firstcol[row] = col;
+      }
+
+  // Compute etree by Liu's algorithm for symmetric matrices,
+  // except use (firstcol[r],c) in place of an edge (r,c) of A.
+  // Thus each row clique in A'*A is replaced by a star
+  // centered at its first vertex, which has the same fill.
+  for (octave_idx_type col = 0; col < nc; col++)
+    {
+      octave_idx_type cset = make_set (col, pp);
+      root[cset] = col;
+      parent[col] = nc;
+      for (octave_idx_type p = colbeg[col]; p < colend[col]; p++)
+        {
+          octave_idx_type row = firstcol[ridx[p]];
+          if (row >= col)
+            continue;
+          octave_idx_type rset = find (row, pp);
+          octave_idx_type rroot = root[rset];
+          if (rroot != col)
+            {
+              parent[rroot] = col;
+              cset = link (cset, rset, pp);
+              root[cset] = col;
+            }
+        }
+    }
+}
+
+DEFUN (colamd, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} colamd (@var{S})
+@deftypefnx {} {@var{p} =} colamd (@var{S}, @var{knobs})
+@deftypefnx {} {[@var{p}, @var{stats}] =} colamd (@var{S})
+@deftypefnx {} {[@var{p}, @var{stats}] =} colamd (@var{S}, @var{knobs})
+
+Compute the column approximate minimum degree permutation.
+
+@code{@var{p} = colamd (@var{S})} returns the column approximate minimum
+degree permutation vector for the sparse matrix @var{S}.  For a
+non-symmetric matrix @var{S}, @code{@var{S}(:,@var{p})} tends to have
+sparser LU@tie{}factors than @var{S}.  The Cholesky@tie{}factorization of
+@code{@var{S}(:,@var{p})' * @var{S}(:,@var{p})} also tends to be sparser
+than that of @code{@var{S}' * @var{S}}.
+
+@var{knobs} is an optional one- to three-element input vector.  If @var{S}
+is m-by-n, then rows with more than @code{max(16,@var{knobs}(1)*sqrt(n))}
+entries are ignored.  Columns with more than
+@code{max (16,@var{knobs}(2)*sqrt(min(m,n)))} entries are removed prior to
+ordering, and ordered last in the output permutation @var{p}.  Only
+completely dense rows or columns are removed if @code{@var{knobs}(1)} and
+@code{@var{knobs}(2)} are < 0, respectively.  If @code{@var{knobs}(3)} is
+nonzero, @var{stats} and @var{knobs} are printed.  The default is
+@code{@var{knobs} = [10 10 0]}.  Note that @var{knobs} differs from earlier
+versions of colamd.
+
+@var{stats} is an optional 20-element output vector that provides data
+about the ordering and the validity of the input matrix @var{S}.  Ordering
+statistics are in @code{@var{stats}(1:3)}.  @code{@var{stats}(1)} and
+@code{@var{stats}(2)} are the number of dense or empty rows and columns
+ignored by @sc{colamd} and @code{@var{stats}(3)} is the number of garbage
+collections performed on the internal data structure used by @sc{colamd}
+(roughly of size @code{2.2 * nnz(@var{S}) + 4 * @var{m} + 7 * @var{n}}
+integers).
+
+Octave built-in functions are intended to generate valid sparse matrices,
+with no duplicate entries, with ascending row indices of the nonzeros
+in each column, with a non-negative number of entries in each column (!)
+and so on.  If a matrix is invalid, then @sc{colamd} may or may not be able
+to continue.  If there are duplicate entries (a row index appears two or
+more times in the same column) or if the row indices in a column are out
+of order, then @sc{colamd} can correct these errors by ignoring the
+duplicate entries and sorting each column of its internal copy of the
+matrix @var{S} (the input matrix @var{S} is not repaired, however).  If a
+matrix is invalid in other ways then @sc{colamd} cannot continue, an error
+message is printed, and no output arguments (@var{p} or @var{stats}) are
+returned.
+@sc{colamd} is thus a simple way to check a sparse matrix to see if it's
+valid.
+
+@code{@var{stats}(4:7)} provide information if @sc{colamd} was able to
+continue.  The matrix is OK if @code{@var{stats}(4)} is zero, or 1 if
+invalid.  @code{@var{stats}(5)} is the rightmost column index that is
+unsorted or contains duplicate entries, or zero if no such column exists.
+@code{@var{stats}(6)} is the last seen duplicate or out-of-order row
+index in the column index given by @code{@var{stats}(5)}, or zero if no
+such row index exists.  @code{@var{stats}(7)} is the number of duplicate
+or out-of-order row indices.  @code{@var{stats}(8:20)} is always zero in
+the current version of @sc{colamd} (reserved for future use).
+
+The ordering is followed by a column elimination tree post-ordering.
+
+The authors of the code itself are @nospell{Stefan I. Larimore} and
+@nospell{Timothy A. Davis}.  The algorithm was developed in collaboration with
+@nospell{John Gilbert}, Xerox PARC, and @nospell{Esmond Ng}, Oak Ridge National
+Laboratory.  (see @url{http://faculty.cse.tamu.edu/davis/suitesparse.html})
+@seealso{colperm, symamd, ccolamd}
+@end deftypefn */)
+{
+#if defined (HAVE_COLAMD)
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+  int spumoni = 0;
+
+  // Get knobs
+  static_assert (COLAMD_KNOBS <= 40,
+                 "colamd: # of COLAMD_KNOBS exceeded.  Please report this to bugs.octave.org");
+  double knob_storage[COLAMD_KNOBS];
+  double *knobs = &knob_storage[0];
+  COLAMD_NAME (_set_defaults) (knobs);
+
+  // Check for user-passed knobs
+  if (nargin == 2)
+    {
+      NDArray User_knobs = args(1).array_value ();
+      int nel_User_knobs = User_knobs.numel ();
+
+      if (nel_User_knobs > 0)
+        knobs[COLAMD_DENSE_ROW] = User_knobs(0);
+      if (nel_User_knobs > 1)
+        knobs[COLAMD_DENSE_COL] = User_knobs(1);
+      if (nel_User_knobs > 2)
+        spumoni = static_cast<int> (User_knobs(2));
+
+      // print knob settings if spumoni is set
+      if (spumoni)
+        {
+
+          octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION
+                        << '.' <<  COLAMD_SUB_VERSION
+                        << ", " << COLAMD_DATE << ":\n";
+
+          if (knobs[COLAMD_DENSE_ROW] >= 0)
+            octave_stdout << "knobs(1): " << User_knobs (0)
+                          << ", rows with > max (16,"
+                          << knobs[COLAMD_DENSE_ROW] << "*sqrt (size(A,2)))"
+                          << " entries removed\n";
+          else
+            octave_stdout << "knobs(1): " << User_knobs (0)
+                          << ", only completely dense rows removed\n";
+
+          if (knobs[COLAMD_DENSE_COL] >= 0)
+            octave_stdout << "knobs(2): " << User_knobs (1)
+                          << ", cols with > max (16,"
+                          << knobs[COLAMD_DENSE_COL] << "*sqrt (size(A)))"
+                          << " entries removed\n";
+          else
+            octave_stdout << "knobs(2): " << User_knobs (1)
+                          << ", only completely dense columns removed\n";
+
+          octave_stdout << "knobs(3): " << User_knobs (2)
+                        << ", statistics and knobs printed\n";
+
+        }
+    }
+
+  octave_idx_type n_row, n_col, nnz;
+  octave_idx_type *ridx, *cidx;
+  SparseComplexMatrix scm;
+  SparseMatrix sm;
+
+  if (args(0).issparse ())
+    {
+      if (args(0).iscomplex ())
+        {
+          scm = args(0).sparse_complex_matrix_value ();
+          n_row = scm.rows ();
+          n_col = scm.cols ();
+          nnz = scm.nnz ();
+          ridx = scm.xridx ();
+          cidx = scm.xcidx ();
+        }
+      else
+        {
+          sm = args(0).sparse_matrix_value ();
+
+          n_row = sm.rows ();
+          n_col = sm.cols ();
+          nnz = sm.nnz ();
+          ridx = sm.xridx ();
+          cidx = sm.xcidx ();
+        }
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+      else
+        sm = SparseMatrix (args(0).matrix_value ());
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      nnz = sm.nnz ();
+      ridx = sm.xridx ();
+      cidx = sm.xcidx ();
+    }
+
+  // Allocate workspace for colamd
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, p, n_col+1);
+  for (octave_idx_type i = 0; i < n_col+1; i++)
+    p[i] = cidx[i];
+
+  octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
+  OCTAVE_LOCAL_BUFFER (octave::suitesparse_integer, A, Alen);
+  for (octave_idx_type i = 0; i < nnz; i++)
+    A[i] = ridx[i];
+
+  // Order the columns (destroys A)
+  static_assert (COLAMD_STATS <= 40,
+                 "colamd: # of COLAMD_STATS exceeded.  Please report this to bugs.octave.org");
+  octave::suitesparse_integer stats_storage[COLAMD_STATS];
+  octave::suitesparse_integer *stats = &stats_storage[0];
+  if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
+    {
+      COLAMD_NAME (_report)(stats);
+
+      error ("colamd: internal error!");
+    }
+
+  // column elimination tree post-ordering (reuse variables)
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
+
+  for (octave_idx_type i = 0; i < n_col; i++)
+    {
+      colbeg[i] = cidx[p[i]];
+      colend[i] = cidx[p[i]+1];
+    }
+
+  coletree (ridx, colbeg, colend, etree, n_row, n_col);
+
+  // Calculate the tree post-ordering
+  tree_postorder (n_col, etree, colbeg);
+
+  // return the permutation vector
+  NDArray out_perm (dim_vector (1, n_col));
+  for (octave_idx_type i = 0; i < n_col; i++)
+    out_perm(i) = p[colbeg[i]] + 1;
+
+  retval(0) = out_perm;
+
+  // print stats if spumoni > 0
+  if (spumoni > 0)
+    COLAMD_NAME (_report)(stats);
+
+  // Return the stats vector
+  if (nargout == 2)
+    {
+      NDArray out_stats (dim_vector (1, COLAMD_STATS));
+      for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
+        out_stats(i) = stats[i];
+      retval(1) = out_stats;
+
+      // fix stats (5) and (6), for 1-based information on
+      // jumbled matrix.  note that this correction doesn't
+      // occur if symamd returns FALSE
+      out_stats(COLAMD_INFO1)++;
+      out_stats(COLAMD_INFO2)++;
+    }
+
+  return retval;
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("colamd", "COLAMD");
+
+#endif
+}
+
+DEFUN (symamd, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} symamd (@var{S})
+@deftypefnx {} {@var{p} =} symamd (@var{S}, @var{knobs})
+@deftypefnx {} {[@var{p}, @var{stats}] =} symamd (@var{S})
+@deftypefnx {} {[@var{p}, @var{stats}] =} symamd (@var{S}, @var{knobs})
+
+For a symmetric positive definite matrix @var{S}, returns the permutation
+vector p such that @code{@var{S}(@var{p}, @var{p})} tends to have a
+sparser Cholesky@tie{}factor than @var{S}.
+
+Sometimes @code{symamd} works well for symmetric indefinite matrices too.
+The matrix @var{S} is assumed to be symmetric; only the strictly lower
+triangular part is referenced.  @var{S} must be square.
+
+@var{knobs} is an optional one- to two-element input vector.  If @var{S} is
+n-by-n, then rows and columns with more than
+@code{max (16,@var{knobs}(1)*sqrt(n))} entries are removed prior to
+ordering, and ordered last in the output permutation @var{p}.  No
+rows/columns are removed if @code{@var{knobs}(1) < 0}.  If
+@code{@var{knobs}(2)} is nonzero, @var{stats} and @var{knobs} are
+printed.  The default is @code{@var{knobs} = [10 0]}.  Note that
+@var{knobs} differs from earlier versions of @code{symamd}.
+
+@var{stats} is an optional 20-element output vector that provides data
+about the ordering and the validity of the input matrix @var{S}.  Ordering
+statistics are in @code{@var{stats}(1:3)}.
+@code{@var{stats}(1) = @var{stats}(2)} is the number of dense or empty rows
+and columns ignored by SYMAMD and @code{@var{stats}(3)} is the number of
+garbage collections performed on the internal data structure used by SYMAMD
+(roughly of size @code{8.4 * nnz (tril (@var{S}, -1)) + 9 * @var{n}}
+integers).
+
+Octave built-in functions are intended to generate valid sparse matrices,
+with no duplicate entries, with ascending row indices of the nonzeros
+in each column, with a non-negative number of entries in each column (!)
+and so on.  If a matrix is invalid, then SYMAMD may or may not be able
+to continue.  If there are duplicate entries (a row index appears two or
+more times in the same column) or if the row indices in a column are out
+of order, then SYMAMD can correct these errors by ignoring the duplicate
+entries and sorting each column of its internal copy of the matrix S (the
+input matrix S is not repaired, however).  If a matrix is invalid in
+other ways then SYMAMD cannot continue, an error message is printed, and
+no output arguments (@var{p} or @var{stats}) are returned.  SYMAMD is
+thus a simple way to check a sparse matrix to see if it's valid.
+
+@code{@var{stats}(4:7)} provide information if SYMAMD was able to
+continue.  The matrix is OK if @code{@var{stats} (4)} is zero, or 1
+if invalid.  @code{@var{stats}(5)} is the rightmost column index that
+is unsorted or contains duplicate entries, or zero if no such column
+exists.  @code{@var{stats}(6)} is the last seen duplicate or out-of-order
+row index in the column index given by @code{@var{stats}(5)}, or zero
+if no such row index exists.  @code{@var{stats}(7)} is the number of
+duplicate or out-of-order row indices.  @code{@var{stats}(8:20)} is
+always zero in the current version of SYMAMD (reserved for future use).
+
+The ordering is followed by a column elimination tree post-ordering.
+
+The authors of the code itself are @nospell{Stefan I. Larimore} and
+@nospell{Timothy A. Davis}.  The algorithm was developed in collaboration with
+@nospell{John Gilbert}, Xerox PARC, and @nospell{Esmond Ng}, Oak Ridge National
+Laboratory.  (see @url{http://faculty.cse.tamu.edu/davis/suitesparse.html})
+@seealso{colperm, colamd}
+@end deftypefn */)
+{
+#if defined (HAVE_COLAMD)
+
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+
+  octave_value_list retval (nargin == 2 ? 2 : 1);
+  int spumoni = 0;
+
+  // Get knobs
+  static_assert (COLAMD_KNOBS <= 40,
+                 "symamd: # of COLAMD_KNOBS exceeded.  Please report this to bugs.octave.org");
+  double knob_storage[COLAMD_KNOBS];
+  double *knobs = &knob_storage[0];
+  COLAMD_NAME (_set_defaults) (knobs);
+
+  // Check for user-passed knobs
+  if (nargin == 2)
+    {
+      NDArray User_knobs = args(1).array_value ();
+      int nel_User_knobs = User_knobs.numel ();
+
+      if (nel_User_knobs > 0)
+        knobs[COLAMD_DENSE_ROW] = User_knobs(COLAMD_DENSE_ROW);
+      if (nel_User_knobs > 1)
+        spumoni = static_cast<int> (User_knobs (1));
+    }
+
+  // print knob settings if spumoni is set
+  if (spumoni > 0)
+    octave_stdout << "symamd: dense row/col fraction: "
+                  << knobs[COLAMD_DENSE_ROW] << std::endl;
+
+  octave_idx_type n_row, n_col;
+  octave_idx_type *ridx, *cidx;
+  SparseMatrix sm;
+  SparseComplexMatrix scm;
+
+  if (args(0).issparse ())
+    {
+      if (args(0).iscomplex ())
+        {
+          scm = args(0).sparse_complex_matrix_value ();
+          n_row = scm.rows ();
+          n_col = scm.cols ();
+          ridx = scm.xridx ();
+          cidx = scm.xcidx ();
+        }
+      else
+        {
+          sm = args(0).sparse_matrix_value ();
+          n_row = sm.rows ();
+          n_col = sm.cols ();
+          ridx = sm.xridx ();
+          cidx = sm.xcidx ();
+        }
+    }
+  else
+    {
+      if (args(0).iscomplex ())
+        sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+      else
+        sm = SparseMatrix (args(0).matrix_value ());
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      ridx = sm.xridx ();
+      cidx = sm.xcidx ();
+    }
+
+  if (n_row != n_col)
+    err_square_matrix_required ("symamd", "S");
+
+  // Allocate workspace for symamd
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1);
+  static_assert (COLAMD_STATS <= 40,
+                 "symamd: # of COLAMD_STATS exceeded.  Please report this to bugs.octave.org");
+  octave::suitesparse_integer stats_storage[COLAMD_STATS];
+  octave::suitesparse_integer *stats = &stats_storage[0];
+  if (! SYMAMD_NAME () (n_col, octave::to_suitesparse_intptr (ridx),
+                        octave::to_suitesparse_intptr (cidx),
+                        octave::to_suitesparse_intptr (perm),
+                        knobs, stats, &calloc, &free))
+    {
+      SYMAMD_NAME (_report)(stats);
+
+      error ("symamd: internal error!");
+    }
+
+  // column elimination tree post-ordering
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
+  symetree (ridx, cidx, etree, perm, n_col);
+
+  // Calculate the tree post-ordering
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
+  tree_postorder (n_col, etree, post);
+
+  // return the permutation vector
+  NDArray out_perm (dim_vector (1, n_col));
+  for (octave_idx_type i = 0; i < n_col; i++)
+    out_perm(i) = perm[post[i]] + 1;
+
+  retval(0) = out_perm;
+
+  // print stats if spumoni > 0
+  if (spumoni > 0)
+    SYMAMD_NAME (_report)(stats);
+
+  // Return the stats vector
+  if (nargout == 2)
+    {
+      NDArray out_stats (dim_vector (1, COLAMD_STATS));
+      for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
+        out_stats(i) = stats[i];
+      retval(1) = out_stats;
+
+      // fix stats (5) and (6), for 1-based information on
+      // jumbled matrix.  note that this correction doesn't
+      // occur if symamd returns FALSE
+      out_stats(COLAMD_INFO1)++;
+      out_stats(COLAMD_INFO2)++;
+    }
+
+  return retval;
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("symamd", "COLAMD");
+
+#endif
+}
+
+DEFUN (etree, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} etree (@var{S})
+@deftypefnx {} {@var{p} =} etree (@var{S}, @var{typ})
+@deftypefnx {} {[@var{p}, @var{q}] =} etree (@var{S}, @var{typ})
+
+Return the elimination tree for the matrix @var{S}.
+
+By default @var{S} is assumed to be symmetric and the symmetric elimination
+tree is returned.  The argument @var{typ} controls whether a symmetric or
+column elimination tree is returned.  Valid values of @var{typ} are
+@qcode{"sym"} or @qcode{"col"}, for symmetric or column elimination tree
+respectively.
+
+Called with a second argument, @code{etree} also returns the postorder
+permutations on the tree.
+@end deftypefn */)
+{
+  int nargin = args.length ();
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+
+  octave_value_list retval (nargout == 2 ? 2 : 1);
+
+  octave_idx_type n_row = 0;
+  octave_idx_type n_col = 0;
+  octave_idx_type *ridx = nullptr;
+  octave_idx_type *cidx = nullptr;
+
+  if (! args(0).issparse ())
+    error ("etree: S must be a sparse matrix");
+
+  if (args(0).iscomplex ())
+    {
+      SparseComplexMatrix scm = args(0).sparse_complex_matrix_value ();
+
+      n_row = scm.rows ();
+      n_col = scm.cols ();
+      ridx = scm.xridx ();
+      cidx = scm.xcidx ();
+    }
+  else
+    {
+      SparseMatrix sm = args(0).sparse_matrix_value ();
+
+      n_row = sm.rows ();
+      n_col = sm.cols ();
+      ridx = sm.xridx ();
+      cidx = sm.xcidx ();
+    }
+
+  bool is_sym = true;
+
+  if (nargin == 2)
+    {
+      std::string str = args(1).xstring_value ("etree: TYP must be a string");
+      if (str.find ('C') == 0 || str.find ('c') == 0)
+        is_sym = false;
+    }
+
+  // column elimination tree post-ordering (reuse variables)
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
+
+  if (is_sym)
+    {
+      if (n_row != n_col)
+        error ("etree: S is marked as symmetric, but is not square");
+
+      symetree (ridx, cidx, etree, nullptr, n_col);
+    }
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
+
+      for (octave_idx_type i = 0; i < n_col; i++)
+        {
+          colbeg[i] = cidx[i];
+          colend[i] = cidx[i+1];
+        }
+
+      coletree (ridx, colbeg, colend, etree, n_row, n_col);
+    }
+
+  NDArray tree (dim_vector (1, n_col));
+  for (octave_idx_type i = 0; i < n_col; i++)
+    // We flag a root with n_col while Matlab does it with zero
+    // Convert for matlab compatible output
+    if (etree[i] == n_col)
+      tree(i) = 0;
+    else
+      tree(i) = etree[i] + 1;
+
+  retval(0) = tree;
+
+  if (nargout == 2)
+    {
+      // Calculate the tree post-ordering
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
+      tree_postorder (n_col, etree, post);
+
+      NDArray postorder (dim_vector (1, n_col));
+      for (octave_idx_type i = 0; i < n_col; i++)
+        postorder(i) = post[i] + 1;
+
+      retval(1) = postorder;
+    }
+
+  return retval;
+}
+
+/*
+%!assert (etree (speye (2)), [0, 0]);
+%!assert (etree (gallery ("poisson", 16)), [2:256, 0]);
+
+%!error etree ()
+%!error etree (1, 2, 3)
+%!error <S must be a sparse matrix> etree ([1, 2; 3, 4])
+%!error <TYP must be a string> etree (speye (2), 3)
+%!error <is not square> etree (sprand (2, 4, .25))
+*/

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/dmperm.cc	Wed Jan 29 06:30:40 2020 -0500
@@ -0,0 +1,215 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 1998-2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// 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
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include "CSparse.h"
+#include "dRowVector.h"
+#include "dSparse.h"
+#include "oct-sparse.h"
+
+#include "defun.h"
+#include "errwarn.h"
+#include "ov.h"
+#include "ovl.h"
+#include "utils.h"
+
+#if defined (OCTAVE_ENABLE_64)
+#  define CXSPARSE_NAME(name) cs_dl ## name
+#else
+#  define CXSPARSE_NAME(name) cs_di ## name
+#endif
+
+#if defined (HAVE_CXSPARSE)
+
+static RowVector
+put_int (octave::suitesparse_integer *p, octave_idx_type n)
+{
+  RowVector ret (n);
+  for (octave_idx_type i = 0; i < n; i++)
+    ret.xelem (i) = p[i] + 1;
+  return ret;
+}
+
+static octave_value_list
+dmperm_internal (bool rank, const octave_value arg, int nargout)
+{
+  octave_value_list retval;
+  octave_idx_type nr = arg.rows ();
+  octave_idx_type nc = arg.columns ();
+  SparseMatrix m;
+  SparseComplexMatrix cm;
+  CXSPARSE_NAME () csm;
+  csm.m = nr;
+  csm.n = nc;
+  csm.x = nullptr;
+  csm.nz = -1;
+
+  if (arg.isreal ())
+    {
+      m = arg.sparse_matrix_value ();
+      csm.nzmax = m.nnz ();
+      csm.p = octave::to_suitesparse_intptr (m.xcidx ());
+      csm.i = octave::to_suitesparse_intptr (m.xridx ());
+    }
+  else
+    {
+      cm = arg.sparse_complex_matrix_value ();
+      csm.nzmax = cm.nnz ();
+      csm.p = octave::to_suitesparse_intptr (cm.xcidx ());
+      csm.i = octave::to_suitesparse_intptr (cm.xridx ());
+    }
+
+  if (nargout <= 1 || rank)
+    {
+      octave::suitesparse_integer *jmatch = CXSPARSE_NAME (_maxtrans) (&csm, 0);
+      if (rank)
+        {
+          octave_idx_type r = 0;
+          for (octave_idx_type i = 0; i < nc; i++)
+            if (jmatch[nr+i] >= 0)
+              r++;
+          retval(0) = static_cast<double> (r);
+        }
+      else
+        retval(0) = put_int (jmatch + nr, nc);
+      CXSPARSE_NAME (_free) (jmatch);
+    }
+  else
+    {
+      CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm, 0);
+
+      //retval(5) = put_int (dm->rr, 5);
+      //retval(4) = put_int (dm->cc, 5);
+      retval = ovl (put_int (dm->p, nr), put_int (dm->q, nc),
+                    put_int (dm->r, dm->nb+1), put_int (dm->s, dm->nb+1));
+
+      CXSPARSE_NAME (_dfree) (dm);
+    }
+
+  return retval;
+}
+
+#endif
+
+DEFUN (dmperm, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn  {} {@var{p} =} dmperm (@var{S})
+@deftypefnx {} {[@var{p}, @var{q}, @var{r}, @var{S}] =} dmperm (@var{S})
+
+@cindex @nospell{Dulmage-Mendelsohn} decomposition
+Perform a @nospell{Dulmage-Mendelsohn} permutation of the sparse matrix
+@var{S}.
+
+With a single output argument @code{dmperm} performs the row permutations
+@var{p} such that @code{@var{S}(@var{p},:)} has no zero elements on the
+diagonal.
+
+Called with two or more output arguments, returns the row and column
+permutations, such that @code{@var{S}(@var{p}, @var{q})} is in block
+triangular form.  The values of @var{r} and @var{S} define the boundaries
+of the blocks.  If @var{S} is square then @code{@var{r} == @var{S}}.
+
+The method used is described in: @nospell{A. Pothen & C.-J. Fan.}
+@cite{Computing the Block Triangular Form of a Sparse Matrix}.
+@nospell{ACM} Trans.@: Math.@: Software, 16(4):303--324, 1990.
+@seealso{colamd, ccolamd}
+@end deftypefn */)
+{
+#if defined (HAVE_CXSPARSE)
+
+  if (args.length () != 1)
+    print_usage ();
+
+  return dmperm_internal (false, args(0), nargout);
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("dmperm", "CXSparse");
+
+#endif
+}
+
+/*
+%!testif HAVE_CXSPARSE
+%! n = 20;
+%! a = speye (n,n);
+%! a = a(randperm (n),:);
+%! assert (a(dmperm (a),:), speye (n));
+
+%!testif HAVE_CXSPARSE
+%! n = 20;
+%! d = 0.2;
+%! a = tril (sprandn (n,n,d), -1) + speye (n,n);
+%! a = a(randperm (n), randperm (n));
+%! [p,q,r,s] = dmperm (a);
+%! assert (tril (a(p,q), -1), sparse (n, n));
+*/
+
+DEFUN (sprank, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn {} {@var{p} =} sprank (@var{S})
+@cindex structural rank
+
+Calculate the structural rank of the sparse matrix @var{S}.
+
+Note that only the structure of the matrix is used in this calculation based
+on a @nospell{Dulmage-Mendelsohn} permutation to block triangular form.  As
+such the numerical rank of the matrix @var{S} is bounded by
+@code{sprank (@var{S}) >= rank (@var{S})}.  Ignoring floating point errors
+@code{sprank (@var{S}) == rank (@var{S})}.
+@seealso{dmperm}
+@end deftypefn */)
+{
+#if defined (HAVE_CXSPARSE)
+
+  if (args.length () != 1)
+    print_usage ();
+
+  return dmperm_internal (true, args(0), nargout);
+
+#else
+
+  octave_unused_parameter (args);
+  octave_unused_parameter (nargout);
+
+  err_disabled_feature ("sprank", "CXSparse");
+
+#endif
+}
+
+/*
+%!testif HAVE_CXSPARSE
+%! assert (sprank (speye (20)), 20);
+%!testif HAVE_CXSPARSE
+%! assert (sprank ([1,0,2,0;2,0,4,0]), 2);
+
+%!error sprank (1,2)
+*/

--- a/libinterp/corefcn/module.mk	Sat Jan 25 15:26:07 2020 +0100
+++ b/libinterp/corefcn/module.mk	Wed Jan 29 06:30:40 2020 -0500
@@ -113,6 +113,7 @@
   %reldir%/__betainc__.cc \
   %reldir%/__contourc__.cc \
   %reldir%/__dsearchn__.cc \
+  %reldir%/__eigs__.cc \
   %reldir%/__expint__.cc \
   %reldir%/__ftp__.cc \
   %reldir%/__gammainc__.cc \
@@ -122,16 +123,20 @@
   %reldir%/__magick_read__.cc \
   %reldir%/__pchip_deriv__.cc \
   %reldir%/__qp__.cc \
+  %reldir%/amd.cc \
   %reldir%/balance.cc \
   %reldir%/besselj.cc \
   %reldir%/bitfcns.cc \
   %reldir%/bsxfun.cc \
   %reldir%/c-file-ptr-stream.cc \
   %reldir%/call-stack.cc \
+  %reldir%/ccolamd.cc \
   %reldir%/cdisplay.c \
   %reldir%/cellfun.cc \
+  %reldir%/chol.cc \
+  %reldir%/coct-hdf5-types.c \
+  %reldir%/colamd.cc \
   %reldir%/colloc.cc \
-  %reldir%/coct-hdf5-types.c \
   %reldir%/conv2.cc \
   %reldir%/daspk.cc \
   %reldir%/dasrt.cc \
@@ -144,6 +149,7 @@
   %reldir%/dirfns.cc \
   %reldir%/display.cc \
   %reldir%/dlmread.cc \
+  %reldir%/dmperm.cc \
   %reldir%/dot.cc \
   %reldir%/dynamic-ld.cc \
   %reldir%/eig.cc \
@@ -178,9 +184,9 @@
   %reldir%/hex2num.cc \
   %reldir%/hook-fcn.cc \
   %reldir%/input.cc \
-  %reldir%/inv.cc \
   %reldir%/interpreter-private.cc \
   %reldir%/interpreter.cc \
+  %reldir%/inv.cc \
   %reldir%/kron.cc \
   %reldir%/load-path.cc \
   %reldir%/load-save.cc \
@@ -190,8 +196,8 @@
   %reldir%/ls-mat-ascii.cc \
   %reldir%/ls-mat4.cc \
   %reldir%/ls-mat5.cc \
+  %reldir%/ls-oct-binary.cc \
   %reldir%/ls-oct-text.cc \
-  %reldir%/ls-oct-binary.cc \
   %reldir%/ls-utils.cc \
   %reldir%/lsode.cc \
   %reldir%/lu.cc \
@@ -221,6 +227,7 @@
   %reldir%/pr-output.cc \
   %reldir%/procstream.cc \
   %reldir%/psi.cc \
+  %reldir%/qr.cc \
   %reldir%/quad.cc \
   %reldir%/quadcc.cc \
   %reldir%/qz.cc \
@@ -242,21 +249,23 @@
   %reldir%/sub2ind.cc \
   %reldir%/svd.cc \
   %reldir%/sylvester.cc \
+  %reldir%/symbfact.cc \
   %reldir%/syminfo.cc \
+  %reldir%/symrcm.cc \
   %reldir%/symrec.cc \
   %reldir%/symscope.cc \
   %reldir%/symtab.cc \
   %reldir%/syscalls.cc \
   %reldir%/sysdep.cc \
-  %reldir%/time.cc \
   %reldir%/text-engine.cc \
   %reldir%/text-renderer.cc \
+  %reldir%/time.cc \
   %reldir%/toplev.cc \
   %reldir%/tril.cc \
   %reldir%/tsearch.cc \
   %reldir%/typecast.cc \
+  %reldir%/url-handle-manager.cc \
   %reldir%/urlwrite.cc \
-  %reldir%/url-handle-manager.cc \
   %reldir%/utils.cc \
   %reldir%/variables.cc \
   %reldir%/xdiv.cc \