# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1453958133 18000
# Node ID ea9c050148093d0ad5841f84e642b47c053c0f93
# Parent  307096fb67e14016af0ee0f9d856588633afd843
revamp sparse LU factorization classes

* sparse-lu.h, sparse-lu.cc: Rename from sparse-base-lu.h and
sparse-base-lu.cc, respectively.
(class sparse_lu): Rename from sparse_base_lu.  Incorporate code
from SparseCmplxLU and SparsedbleLU classes into the sparse_lu
template.

* sparse-lu-inst.cc: New file.

* SparseCmplxLU.cc, SparseCmplxLU.h, SparsedbleLU.cc, SparsedbleLU.h:
Delete.

* lu.cc, luinc.cc, CSparse.cc, dSparse.cc, eigs-base.cc: Change
all uses of SparsedbleLU and SparseCmplxLU to use new
sparse_lu template class.

* liboctave/numeric/module.mk: Update.

diff -r 307096fb67e1 -r ea9c05014809 libinterp/corefcn/lu.cc
--- a/libinterp/corefcn/lu.cc	Tue Jan 26 16:30:54 2016 -0500
+++ b/libinterp/corefcn/lu.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -28,8 +28,7 @@
 #include "dbleLU.h"
 #include "fCmplxLU.h"
 #include "floatLU.h"
-#include "SparseCmplxLU.h"
-#include "SparsedbleLU.h"
+#include "sparse-lu.h"
 
 #include "defun.h"
 #include "error.h"
@@ -208,7 +207,7 @@
               ColumnVector Qinit (nc);
               for (octave_idx_type i = 0; i < nc; i++)
                 Qinit(i) = i;
-              SparseLU fact (m, Qinit, thres, false, true);
+              sparse_lu<SparseMatrix> fact (m, Qinit, thres, false, true);
 
               if (nargout < 2)
                 retval(0) = fact.Y ();
@@ -242,7 +241,7 @@
           else
             {
               retval.resize (scale ? 5 : 4);
-              SparseLU fact (m, thres, scale);
+              sparse_lu<SparseMatrix> fact (m, thres, scale);
 
               retval(0) = octave_value (fact.L (),
                                         MatrixType (MatrixType::Lower));
@@ -273,7 +272,7 @@
               ColumnVector Qinit (nc);
               for (octave_idx_type i = 0; i < nc; i++)
                 Qinit(i) = i;
-              SparseComplexLU fact (m, Qinit, thres, false, true);
+              sparse_lu<SparseComplexMatrix> fact (m, Qinit, thres, false, true);
 
               if (nargout < 2)
                 retval(0) = fact.Y ();
@@ -306,7 +305,7 @@
           else
             {
               retval.resize (scale ? 5 : 4);
-              SparseComplexLU fact (m, thres, scale);
+              sparse_lu<SparseComplexMatrix> fact (m, thres, scale);
 
               retval(0) = octave_value (fact.L (),
                                         MatrixType (MatrixType::Lower));
diff -r 307096fb67e1 -r ea9c05014809 libinterp/corefcn/luinc.cc
--- a/libinterp/corefcn/luinc.cc	Tue Jan 26 16:30:54 2016 -0500
+++ b/libinterp/corefcn/luinc.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -32,8 +32,7 @@
 #include "oct-map.h"
 
 #include "MatrixType.h"
-#include "SparseCmplxLU.h"
-#include "SparsedbleLU.h"
+#include "sparse-lu.h"
 #include "ov-re-sparse.h"
 #include "ov-cx-sparse.h"
 
@@ -195,7 +194,7 @@
         case 1:
         case 2:
           {
-            SparseLU fact (sm, Qinit, thresh, false, true, droptol,
+            sparse_lu<SparseMatrix> fact (sm, Qinit, thresh, false, true, droptol,
                            milu, udiag);
 
             SparseMatrix P = fact.Pr ();
@@ -212,7 +211,7 @@
 
         case 3:
           {
-            SparseLU fact (sm, Qinit, thresh, false, true, droptol,
+            sparse_lu<SparseMatrix> fact (sm, Qinit, thresh, false, true, droptol,
                            milu, udiag);
 
             if (vecout)
@@ -231,7 +230,7 @@
         case 4:
         default:
           {
-            SparseLU fact (sm, Qinit, thresh, false, false, droptol,
+            sparse_lu<SparseMatrix> fact (sm, Qinit, thresh, false, false, droptol,
                            milu, udiag);
 
             if (vecout)
@@ -270,7 +269,7 @@
         case 1:
         case 2:
           {
-            SparseComplexLU fact (sm, Qinit, thresh, false, true,
+            sparse_lu<SparseComplexMatrix> fact (sm, Qinit, thresh, false, true,
                                   droptol, milu, udiag);
 
             SparseMatrix P = fact.Pr ();
@@ -287,7 +286,7 @@
 
         case 3:
           {
-            SparseComplexLU fact (sm, Qinit, thresh, false, true,
+            sparse_lu<SparseComplexMatrix> fact (sm, Qinit, thresh, false, true,
                                   droptol, milu, udiag);
 
             if (vecout)
@@ -306,7 +305,7 @@
         case 4:
         default:
           {
-            SparseComplexLU fact (sm, Qinit, thresh, false, false,
+            sparse_lu<SparseComplexMatrix> fact (sm, Qinit, thresh, false, false,
                                   droptol, milu, udiag);
 
             if (vecout)
diff -r 307096fb67e1 -r ea9c05014809 liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc	Tue Jan 26 16:30:54 2016 -0500
+++ b/liboctave/array/CSparse.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -51,7 +51,7 @@
 #include "dSparse.h"
 #include "functor.h"
 #include "oct-spparms.h"
-#include "SparseCmplxLU.h"
+#include "sparse-lu.h"
 #include "oct-sparse.h"
 #include "sparse-util.h"
 #include "sparse-chol.h"
@@ -1105,7 +1105,7 @@
             Qinit(i) = i;
 
           MatrixType tmp_typ (MatrixType::Upper);
-          SparseComplexLU fact (*this, Qinit, Matrix (), false, false);
+          sparse_lu<SparseComplexMatrix> fact (*this, Qinit, Matrix (), false, false);
           rcond = fact.rcond ();
           double rcond2;
           SparseComplexMatrix InvL = fact.L ().transpose ().
diff -r 307096fb67e1 -r ea9c05014809 liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc	Tue Jan 26 16:30:54 2016 -0500
+++ b/liboctave/array/dSparse.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -45,7 +45,7 @@
 #include "dSparse.h"
 #include "functor.h"
 #include "oct-spparms.h"
-#include "SparsedbleLU.h"
+#include "sparse-lu.h"
 #include "MatrixType.h"
 #include "oct-sparse.h"
 #include "sparse-util.h"
@@ -1194,7 +1194,7 @@
             Qinit(i) = i;
 
           MatrixType tmp_typ (MatrixType::Upper);
-          SparseLU fact (*this, Qinit, Matrix (), false, false);
+          sparse_lu<SparseMatrix> fact (*this, Qinit, Matrix (), false, false);
           rcond = fact.rcond ();
           double rcond2;
           SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/SparseCmplxLU.cc
--- a/liboctave/numeric/SparseCmplxLU.cc	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,485 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <vector>
-
-#include "lo-error.h"
-#include "oct-locbuf.h"
-
-#include "SparseCmplxLU.h"
-#include "oct-spparms.h"
-
-// Instantiate the base LU class for the types we need.
-
-#include "sparse-base-lu.h"
-#include "sparse-base-lu.cc"
-
-template class sparse_base_lu <SparseComplexMatrix, Complex,
-                               SparseMatrix, double>;
-
-#include "oct-sparse.h"
-
-SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
-                                  const Matrix& piv_thres, bool scale)
-{
-#ifdef HAVE_UMFPACK
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  UMFPACK_ZNAME (defaults) (control);
-
-  double tmp = octave_sparse_params::get_key ("spumoni");
-  if (! xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-  if (piv_thres.numel () == 2)
-    {
-      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-  else
-    {
-      tmp = octave_sparse_params::get_key ("piv_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-
-      tmp = octave_sparse_params::get_key ("sym_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-
-  // Set whether we are allowed to modify Q or not
-  tmp = octave_sparse_params::get_key ("autoamd");
-  if (! xisnan (tmp))
-    Control (UMFPACK_FIXQ) = tmp;
-
-  // Turn-off UMFPACK scaling for LU
-  if (scale)
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
-  else
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  UMFPACK_ZNAME (report_control) (control);
-
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const Complex *Ax = a.data ();
-
-  UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
-                                 reinterpret_cast<const double *> (Ax),
-                                 0, 1, control);
-
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
-                                          reinterpret_cast<const double *> (Ax),
-                                          0, 0,
-                                          &Symbolic, control, info);
-
-  if (status < 0)
-    {
-      UMFPACK_ZNAME (report_status) (control, status);
-      UMFPACK_ZNAME (report_info) (control, info);
-
-      UMFPACK_ZNAME (free_symbolic) (&Symbolic);
-
-      (*current_liboctave_error_handler)
-        ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
-    }
-  else
-    {
-      UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
-
-      void *Numeric;
-      status = UMFPACK_ZNAME (numeric) (Ap, Ai,
-                                        reinterpret_cast<const double *> (Ax),
-                                        0, Symbolic, &Numeric, control,
-                                        info);
-      UMFPACK_ZNAME (free_symbolic) (&Symbolic);
-
-      cond = Info (UMFPACK_RCOND);
-
-      if (status < 0)
-        {
-          UMFPACK_ZNAME (report_status) (control, status);
-          UMFPACK_ZNAME (report_info) (control, info);
-
-          UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-          (*current_liboctave_error_handler)
-            ("SparseComplexLU::SparseComplexLU numeric factorization failed");
-        }
-      else
-        {
-          UMFPACK_ZNAME (report_numeric) (Numeric, control);
-
-          octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
-          status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1,
-                                             &ignore2, &ignore3, Numeric);
-
-          if (status < 0)
-            {
-              UMFPACK_ZNAME (report_status) (control, status);
-              UMFPACK_ZNAME (report_info) (control, info);
-
-              UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-              (*current_liboctave_error_handler)
-                ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
-            }
-          else
-            {
-              octave_idx_type n_inner = (nr < nc ? nr : nc);
-
-              if (lnz < 1)
-                Lfact = SparseComplexMatrix (n_inner, nr,
-                                             static_cast<octave_idx_type> (1));
-              else
-                Lfact = SparseComplexMatrix (n_inner, nr, lnz);
-
-              octave_idx_type *Ltp = Lfact.cidx ();
-              octave_idx_type *Ltj = Lfact.ridx ();
-              Complex *Ltx = Lfact.data ();
-
-              if (unz < 1)
-                Ufact = SparseComplexMatrix (n_inner, nc,
-                                             static_cast<octave_idx_type> (1));
-              else
-                Ufact = SparseComplexMatrix (n_inner, nc, unz);
-
-              octave_idx_type *Up = Ufact.cidx ();
-              octave_idx_type *Uj = Ufact.ridx ();
-              Complex *Ux = Ufact.data ();
-
-              Rfact = SparseMatrix (nr, nr, nr);
-              for (octave_idx_type i = 0; i < nr; i++)
-                {
-                  Rfact.xridx (i) = i;
-                  Rfact.xcidx (i) = i;
-                }
-              Rfact.xcidx (nr) = nr;
-              double *Rx = Rfact.data ();
-
-              P.resize (dim_vector (nr, 1));
-              octave_idx_type *p = P.fortran_vec ();
-
-              Q.resize (dim_vector (nc, 1));
-              octave_idx_type *q = Q.fortran_vec ();
-
-              octave_idx_type do_recip;
-              status = UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
-                                                    reinterpret_cast<double *> (Ltx),
-                                                    0, Up, Uj,
-                                                    reinterpret_cast <double *> (Ux),
-                                                    0, p, q, 0, 0,
-                                                    &do_recip, Rx, Numeric);
-
-              UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-              if (status < 0)
-                {
-                  UMFPACK_ZNAME (report_status) (control, status);
-
-                  (*current_liboctave_error_handler)
-                    ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
-                }
-              else
-                {
-                  Lfact = Lfact.transpose ();
-
-                  if (do_recip)
-                    for (octave_idx_type i = 0; i < nr; i++)
-                      Rx[i] = 1.0 / Rx[i];
-
-                  UMFPACK_ZNAME (report_matrix) (nr, n_inner,
-                                                 Lfact.cidx (), Lfact.ridx (),
-                                                 reinterpret_cast<double *> (Lfact.data ()),
-                                                 0, 1, control);
-
-                  UMFPACK_ZNAME (report_matrix) (n_inner, nc,
-                                                 Ufact.cidx (), Ufact.ridx (),
-                                                 reinterpret_cast<double *> (Ufact.data ()),
-                                                 0, 1, control);
-                  UMFPACK_ZNAME (report_perm) (nr, p, control);
-                  UMFPACK_ZNAME (report_perm) (nc, q, control);
-                }
-
-              UMFPACK_ZNAME (report_info) (control, info);
-            }
-        }
-    }
-
-#else
-  (*current_liboctave_error_handler)
-    ("support for UMFPACK was unavailable or disabled when liboctave was built");
-#endif
-}
-
-SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a,
-                                  const ColumnVector& Qinit,
-                                  const Matrix& piv_thres, bool scale,
-                                  bool FixedQ, double droptol,
-                                  bool milu, bool udiag)
-{
-#ifdef HAVE_UMFPACK
-  if (milu)
-    (*current_liboctave_error_handler)
-      ("Modified incomplete LU not implemented");
-
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  UMFPACK_ZNAME (defaults) (control);
-
-  double tmp = octave_sparse_params::get_key ("spumoni");
-  if (! xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-  if (piv_thres.numel () == 2)
-    {
-      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-  else
-    {
-      tmp = octave_sparse_params::get_key ("piv_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-
-      tmp = octave_sparse_params::get_key ("sym_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-
-  if (droptol >= 0.)
-    Control (UMFPACK_DROPTOL) = droptol;
-
-  // Set whether we are allowed to modify Q or not
-  if (FixedQ)
-    Control (UMFPACK_FIXQ) = 1.0;
-  else
-    {
-      tmp = octave_sparse_params::get_key ("autoamd");
-      if (! xisnan (tmp))
-        Control (UMFPACK_FIXQ) = tmp;
-    }
-
-  // Turn-off UMFPACK scaling for LU
-  if (scale)
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
-  else
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  UMFPACK_ZNAME (report_control) (control);
-
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const Complex *Ax = a.data ();
-
-  UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
-                                 reinterpret_cast<const double *> (Ax), 0,
-                                 1, control);
-
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status;
-
-  // Null loop so that qinit is imediately deallocated when not
-  // needed
-  do
-    {
-      OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
-
-      for (octave_idx_type i = 0; i < nc; i++)
-        qinit[i] = static_cast<octave_idx_type> (Qinit (i));
-
-      status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
-                                          reinterpret_cast<const double *> (Ax),
-                                          0, qinit, &Symbolic, control,
-                                          info);
-    }
-  while (0);
-
-  if (status < 0)
-    {
-      UMFPACK_ZNAME (report_status) (control, status);
-      UMFPACK_ZNAME (report_info) (control, info);
-
-      UMFPACK_ZNAME (free_symbolic) (&Symbolic);
-
-      (*current_liboctave_error_handler)
-        ("SparseComplexLU::SparseComplexLU symbolic factorization failed");
-    }
-  else
-    {
-      UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
-
-      void *Numeric;
-      status = UMFPACK_ZNAME (numeric) (Ap, Ai,
-                                        reinterpret_cast<const double *> (Ax), 0,
-                                        Symbolic, &Numeric, control, info);
-      UMFPACK_ZNAME (free_symbolic) (&Symbolic);
-
-      cond = Info (UMFPACK_RCOND);
-
-      if (status < 0)
-        {
-          UMFPACK_ZNAME (report_status) (control, status);
-          UMFPACK_ZNAME (report_info) (control, info);
-
-          UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-          (*current_liboctave_error_handler)
-            ("SparseComplexLU::SparseComplexLU numeric factorization failed");
-        }
-      else
-        {
-          UMFPACK_ZNAME (report_numeric) (Numeric, control);
-
-          octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
-          status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz,
-                                             &ignore1, &ignore2, &ignore3,
-                                             Numeric);
-
-          if (status < 0)
-            {
-              UMFPACK_ZNAME (report_status) (control, status);
-              UMFPACK_ZNAME (report_info) (control, info);
-
-              UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-              (*current_liboctave_error_handler)
-                ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
-            }
-          else
-            {
-              octave_idx_type n_inner = (nr < nc ? nr : nc);
-
-              if (lnz < 1)
-                Lfact = SparseComplexMatrix (n_inner, nr,
-                                             static_cast<octave_idx_type> (1));
-              else
-                Lfact = SparseComplexMatrix (n_inner, nr, lnz);
-
-              octave_idx_type *Ltp = Lfact.cidx ();
-              octave_idx_type *Ltj = Lfact.ridx ();
-              Complex *Ltx = Lfact.data ();
-
-              if (unz < 1)
-                Ufact = SparseComplexMatrix (n_inner, nc,
-                                             static_cast<octave_idx_type> (1));
-              else
-                Ufact = SparseComplexMatrix  (n_inner, nc, unz);
-
-              octave_idx_type *Up = Ufact.cidx ();
-              octave_idx_type *Uj = Ufact.ridx ();
-              Complex *Ux = Ufact.data ();
-
-              Rfact = SparseMatrix (nr, nr, nr);
-              for (octave_idx_type i = 0; i < nr; i++)
-                {
-                  Rfact.xridx (i) = i;
-                  Rfact.xcidx (i) = i;
-                }
-              Rfact.xcidx (nr) = nr;
-              double *Rx = Rfact.data ();
-
-              P.resize (dim_vector (nr, 1));
-              octave_idx_type *p = P.fortran_vec ();
-
-              Q.resize (dim_vector (nc, 1));
-              octave_idx_type *q = Q.fortran_vec ();
-
-              octave_idx_type do_recip;
-              status =
-                UMFPACK_ZNAME (get_numeric) (Ltp, Ltj,
-                                             reinterpret_cast<double *> (Ltx),
-                                             0, Up, Uj,
-                                             reinterpret_cast<double *> (Ux),
-                                             0, p, q, 0, 0,
-                                             &do_recip, Rx, Numeric);
-
-              UMFPACK_ZNAME (free_numeric) (&Numeric);
-
-              if (status < 0)
-                {
-                  UMFPACK_ZNAME (report_status) (control, status);
-
-                  (*current_liboctave_error_handler)
-                    ("SparseComplexLU::SparseComplexLU extracting LU factors failed");
-                }
-              else
-                {
-                  Lfact = Lfact.transpose ();
-
-                  if (do_recip)
-                    for (octave_idx_type i = 0; i < nr; i++)
-                      Rx[i] = 1.0 / Rx[i];
-
-                  UMFPACK_ZNAME (report_matrix) (nr, n_inner,
-                                                 Lfact.cidx (),
-                                                 Lfact.ridx (),
-                                                 reinterpret_cast<double *> (Lfact.data ()),
-                                                 0, 1, control);
-
-                  UMFPACK_ZNAME (report_matrix) (n_inner, nc,
-                                                 Ufact.cidx (),
-                                                 Ufact.ridx (),
-                                                 reinterpret_cast<double *> (Ufact.data ()),
-                                                 0, 1, control);
-                  UMFPACK_ZNAME (report_perm) (nr, p, control);
-                  UMFPACK_ZNAME (report_perm) (nc, q, control);
-                }
-
-              UMFPACK_ZNAME (report_info) (control, info);
-            }
-        }
-    }
-
-  if (udiag)
-    (*current_liboctave_error_handler)
-      ("Option udiag of incomplete LU not implemented");
-
-#else
-  (*current_liboctave_error_handler)
-    ("support for UMFPACK was unavailable or disabled when liboctave was built");
-#endif
-}
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/SparseCmplxLU.h
--- a/liboctave/numeric/SparseCmplxLU.h	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,67 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_SparseCmplxLU_h)
-#define octave_SparseCmplxLU_h 1
-
-#include "sparse-base-lu.h"
-#include "dSparse.h"
-#include "CSparse.h"
-
-class
-OCTAVE_API
-SparseComplexLU
-  : public sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>
-{
-public:
-
-  SparseComplexLU (void)
-    : sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double> () { }
-
-  SparseComplexLU (const SparseComplexMatrix& a,
-                   const Matrix& piv_thres = Matrix (),
-                   bool scale = false);
-
-  SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit,
-                   const Matrix& piv_thres = Matrix (),
-                   bool scale = false, bool FixedQ = false,
-                   double droptol = -1., bool milu = false,
-                   bool udiag = false);
-
-  SparseComplexLU (const SparseComplexLU& a)
-    : sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double> (a)
-  { }
-
-  SparseComplexLU& operator = (const SparseComplexLU& a)
-  {
-    if (this != &a)
-      sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double>
-                     :: operator = (a);
-
-    return *this;
-  }
-
-  ~SparseComplexLU (void) { }
-};
-
-#endif
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/SparsedbleLU.cc
--- a/liboctave/numeric/SparsedbleLU.cc	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,462 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <vector>
-
-#include "lo-error.h"
-#include "oct-locbuf.h"
-
-#include "SparsedbleLU.h"
-#include "oct-spparms.h"
-
-// Instantiate the base LU class for the types we need.
-
-#include "sparse-base-lu.h"
-#include "sparse-base-lu.cc"
-
-template class sparse_base_lu <SparseMatrix, double, SparseMatrix, double>;
-
-#include "oct-sparse.h"
-
-SparseLU::SparseLU (const SparseMatrix& a, const Matrix& piv_thres, bool scale)
-{
-#ifdef HAVE_UMFPACK
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  UMFPACK_DNAME (defaults) (control);
-
-  double tmp = octave_sparse_params::get_key ("spumoni");
-  if (! xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-
-  if (piv_thres.numel () == 2)
-    {
-      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-  else
-    {
-      tmp = octave_sparse_params::get_key ("piv_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-
-      tmp = octave_sparse_params::get_key ("sym_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-
-  // Set whether we are allowed to modify Q or not
-  tmp = octave_sparse_params::get_key ("autoamd");
-  if (! xisnan (tmp))
-    Control (UMFPACK_FIXQ) = tmp;
-
-  if (scale)
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
-  else
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  UMFPACK_DNAME (report_control) (control);
-
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const double *Ax = a.data ();
-
-  UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
-
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
-                                          &Symbolic, control, info);
-
-  if (status < 0)
-    {
-      UMFPACK_DNAME (report_status) (control, status);
-      UMFPACK_DNAME (report_info) (control, info);
-
-      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
-
-      (*current_liboctave_error_handler)
-        ("SparseLU::SparseLU symbolic factorization failed");
-    }
-  else
-    {
-      UMFPACK_DNAME (report_symbolic) (Symbolic, control);
-
-      void *Numeric;
-      status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
-                                        &Numeric, control, info) ;
-      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
-
-      cond = Info (UMFPACK_RCOND);
-
-      if (status < 0)
-        {
-          UMFPACK_DNAME (report_status) (control, status);
-          UMFPACK_DNAME (report_info) (control, info);
-
-          UMFPACK_DNAME (free_numeric) (&Numeric);
-
-          (*current_liboctave_error_handler)
-            ("SparseLU::SparseLU numeric factorization failed");
-        }
-      else
-        {
-          UMFPACK_DNAME (report_numeric) (Numeric, control);
-
-          octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
-          status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1,
-                                             &ignore2, &ignore3, Numeric) ;
-
-          if (status < 0)
-            {
-              UMFPACK_DNAME (report_status) (control, status);
-              UMFPACK_DNAME (report_info) (control, info);
-
-              UMFPACK_DNAME (free_numeric) (&Numeric);
-
-              (*current_liboctave_error_handler)
-                ("SparseLU::SparseLU extracting LU factors failed");
-            }
-          else
-            {
-              octave_idx_type n_inner = (nr < nc ? nr : nc);
-
-              if (lnz < 1)
-                Lfact = SparseMatrix (n_inner, nr,
-                                      static_cast<octave_idx_type> (1));
-              else
-                Lfact = SparseMatrix (n_inner, nr, lnz);
-
-              octave_idx_type *Ltp = Lfact.cidx ();
-              octave_idx_type *Ltj = Lfact.ridx ();
-              double *Ltx = Lfact.data ();
-
-              if (unz < 1)
-                Ufact = SparseMatrix (n_inner, nc,
-                                      static_cast<octave_idx_type> (1));
-              else
-                Ufact = SparseMatrix (n_inner, nc, unz);
-
-              octave_idx_type *Up = Ufact.cidx ();
-              octave_idx_type *Uj = Ufact.ridx ();
-              double *Ux = Ufact.data ();
-
-              Rfact = SparseMatrix (nr, nr, nr);
-              for (octave_idx_type i = 0; i < nr; i++)
-                {
-                  Rfact.xridx (i) = i;
-                  Rfact.xcidx (i) = i;
-                }
-              Rfact.xcidx (nr) = nr;
-              double *Rx = Rfact.data ();
-
-              P.resize (dim_vector (nr, 1));
-              octave_idx_type *p = P.fortran_vec ();
-
-              Q.resize (dim_vector (nc, 1));
-              octave_idx_type *q = Q.fortran_vec ();
-
-              octave_idx_type do_recip;
-              status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx,
-                                                    Up, Uj, Ux, p, q, 0,
-                                                    &do_recip, Rx,
-                                                    Numeric) ;
-
-              UMFPACK_DNAME (free_numeric) (&Numeric) ;
-
-              if (status < 0)
-                {
-                  UMFPACK_DNAME (report_status) (control, status);
-
-                  (*current_liboctave_error_handler)
-                    ("SparseLU::SparseLU extracting LU factors failed");
-                }
-              else
-                {
-                  Lfact = Lfact.transpose ();
-
-                  if (do_recip)
-                    for (octave_idx_type i = 0; i < nr; i++)
-                      Rx[i] = 1.0 / Rx[i];
-
-                  UMFPACK_DNAME (report_matrix) (nr, n_inner,
-                                                 Lfact.cidx (), Lfact.ridx (),
-                                                 Lfact.data (), 1, control);
-                  UMFPACK_DNAME (report_matrix) (n_inner, nc,
-                                                 Ufact.cidx (), Ufact.ridx (),
-                                                 Ufact.data (), 1, control);
-                  UMFPACK_DNAME (report_perm) (nr, p, control);
-                  UMFPACK_DNAME (report_perm) (nc, q, control);
-                }
-
-              UMFPACK_DNAME (report_info) (control, info);
-            }
-        }
-    }
-
-#else
-  (*current_liboctave_error_handler)
-    ("support for UMFPACK was unavailable or disabled when liboctave was built");
-#endif
-}
-
-SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
-                    const Matrix& piv_thres, bool scale, bool FixedQ,
-                    double droptol, bool milu, bool udiag)
-{
-#ifdef HAVE_UMFPACK
-  if (milu)
-    (*current_liboctave_error_handler)
-      ("Modified incomplete LU not implemented");
-
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-
-  // Setup the control parameters
-  Matrix Control (UMFPACK_CONTROL, 1);
-  double *control = Control.fortran_vec ();
-  UMFPACK_DNAME (defaults) (control);
-
-  double tmp = octave_sparse_params::get_key ("spumoni");
-  if (! xisnan (tmp))
-    Control (UMFPACK_PRL) = tmp;
-
-  if (piv_thres.numel () == 2)
-    {
-      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-  else
-    {
-      tmp = octave_sparse_params::get_key ("piv_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
-
-      tmp = octave_sparse_params::get_key ("sym_tol");
-      if (! xisnan (tmp))
-        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
-    }
-
-  if (droptol >= 0.)
-    Control (UMFPACK_DROPTOL) = droptol;
-
-
-  // Set whether we are allowed to modify Q or not
-  if (FixedQ)
-    Control (UMFPACK_FIXQ) = 1.0;
-  else
-    {
-      tmp = octave_sparse_params::get_key ("autoamd");
-      if (! xisnan (tmp))
-        Control (UMFPACK_FIXQ) = tmp;
-    }
-
-  if (scale)
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
-  else
-    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
-
-  UMFPACK_DNAME (report_control) (control);
-
-  const octave_idx_type *Ap = a.cidx ();
-  const octave_idx_type *Ai = a.ridx ();
-  const double *Ax = a.data ();
-
-  UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1,
-                                 control);
-
-  void *Symbolic;
-  Matrix Info (1, UMFPACK_INFO);
-  double *info = Info.fortran_vec ();
-  int status;
-
-  // Null loop so that qinit is imediately deallocated when not needed
-  do
-    {
-      OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
-
-      for (octave_idx_type i = 0; i < nc; i++)
-        qinit[i] = static_cast<octave_idx_type> (Qinit (i));
-
-      status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax,
-                                          qinit, &Symbolic, control, info);
-    }
-  while (0);
-
-  if (status < 0)
-    {
-      UMFPACK_DNAME (report_status) (control, status);
-      UMFPACK_DNAME (report_info) (control, info);
-
-      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
-
-      (*current_liboctave_error_handler)
-        ("SparseLU::SparseLU symbolic factorization failed");
-    }
-  else
-    {
-      UMFPACK_DNAME (report_symbolic) (Symbolic, control);
-
-      void *Numeric;
-      status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
-                                        &Numeric, control, info) ;
-      UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
-
-      cond = Info (UMFPACK_RCOND);
-
-      if (status < 0)
-        {
-          UMFPACK_DNAME (report_status) (control, status);
-          UMFPACK_DNAME (report_info) (control, info);
-
-          UMFPACK_DNAME (free_numeric) (&Numeric);
-
-          (*current_liboctave_error_handler)
-            ("SparseLU::SparseLU numeric factorization failed");
-        }
-      else
-        {
-          UMFPACK_DNAME (report_numeric) (Numeric, control);
-
-          octave_idx_type lnz, unz, ignore1, ignore2, ignore3;
-          status = UMFPACK_DNAME (get_lunz) (&lnz, &unz,
-                                             &ignore1, &ignore2, &ignore3,
-                                             Numeric) ;
-
-          if (status < 0)
-            {
-              UMFPACK_DNAME (report_status) (control, status);
-              UMFPACK_DNAME (report_info) (control, info);
-
-              UMFPACK_DNAME (free_numeric) (&Numeric);
-
-              (*current_liboctave_error_handler)
-                ("SparseLU::SparseLU extracting LU factors failed");
-            }
-          else
-            {
-              octave_idx_type n_inner = (nr < nc ? nr : nc);
-
-              if (lnz < 1)
-                Lfact = SparseMatrix (n_inner, nr,
-                                      static_cast<octave_idx_type> (1));
-              else
-                Lfact = SparseMatrix (n_inner, nr, lnz);
-
-              octave_idx_type *Ltp = Lfact.cidx ();
-              octave_idx_type *Ltj = Lfact.ridx ();
-              double *Ltx = Lfact.data ();
-
-              if (unz < 1)
-                Ufact = SparseMatrix (n_inner, nc,
-                                      static_cast<octave_idx_type> (1));
-              else
-                Ufact = SparseMatrix (n_inner, nc, unz);
-
-              octave_idx_type *Up = Ufact.cidx ();
-              octave_idx_type *Uj = Ufact.ridx ();
-              double *Ux = Ufact.data ();
-
-              Rfact = SparseMatrix (nr, nr, nr);
-              for (octave_idx_type i = 0; i < nr; i++)
-                {
-                  Rfact.xridx (i) = i;
-                  Rfact.xcidx (i) = i;
-                }
-              Rfact.xcidx (nr) = nr;
-              double *Rx = Rfact.data ();
-
-              P.resize (dim_vector (nr, 1));
-              octave_idx_type *p = P.fortran_vec ();
-
-              Q.resize (dim_vector (nc, 1));
-              octave_idx_type *q = Q.fortran_vec ();
-
-              octave_idx_type do_recip;
-              status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj,
-                                                    Ltx, Up, Uj, Ux, p, q,
-                                                    0, &do_recip,
-                                                    Rx, Numeric) ;
-
-              UMFPACK_DNAME (free_numeric) (&Numeric) ;
-
-              if (status < 0)
-                {
-                  UMFPACK_DNAME (report_status) (control, status);
-
-                  (*current_liboctave_error_handler)
-                    ("SparseLU::SparseLU extracting LU factors failed");
-                }
-              else
-                {
-                  Lfact = Lfact.transpose ();
-
-                  if (do_recip)
-                    for (octave_idx_type i = 0; i < nr; i++)
-                      Rx[i] = 1.0 / Rx[i];
-
-                  UMFPACK_DNAME (report_matrix) (nr, n_inner,
-                                                 Lfact.cidx (),
-                                                 Lfact.ridx (),
-                                                 Lfact.data (),
-                                                 1, control);
-                  UMFPACK_DNAME (report_matrix) (n_inner, nc,
-                                                 Ufact.cidx (),
-                                                 Ufact.ridx (),
-                                                 Ufact.data (),
-                                                 1, control);
-                  UMFPACK_DNAME (report_perm) (nr, p, control);
-                  UMFPACK_DNAME (report_perm) (nc, q, control);
-                }
-
-              UMFPACK_DNAME (report_info) (control, info);
-            }
-        }
-    }
-
-  if (udiag)
-    (*current_liboctave_error_handler)
-      ("Option udiag of incomplete LU not implemented");
-
-#else
-  (*current_liboctave_error_handler)
-    ("support for UMFPACK was unavailable or disabled when liboctave was built");
-#endif
-}
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/SparsedbleLU.h
--- a/liboctave/numeric/SparsedbleLU.h	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,62 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_SparsedbleLU_h)
-#define octave_SparsedbleLU_h 1
-
-#include "sparse-base-lu.h"
-#include "dSparse.h"
-
-class
-OCTAVE_API
-SparseLU : public sparse_base_lu <SparseMatrix, double, SparseMatrix, double>
-{
-public:
-
-  SparseLU (void)
-    : sparse_base_lu <SparseMatrix, double, SparseMatrix, double> () { }
-
-  SparseLU (const SparseMatrix& a, const Matrix& piv_thres = Matrix (),
-            bool scale = false);
-
-  SparseLU (const SparseMatrix& a, const ColumnVector& Qinit,
-            const Matrix& piv_thres = Matrix (), bool scale = false,
-            bool FixedQ = false, double droptol = -1.,
-            bool milu = false, bool udiag = false);
-
-  SparseLU (const SparseLU& a)
-    : sparse_base_lu <SparseMatrix, double, SparseMatrix, double> (a) { }
-
-  SparseLU& operator = (const SparseLU& a)
-  {
-    if (this != &a)
-      sparse_base_lu <SparseMatrix, double, SparseMatrix, double>
-      :: operator = (a);
-
-    return *this;
-  }
-
-  ~SparseLU (void) { }
-};
-
-#endif
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/eigs-base.cc
--- a/liboctave/numeric/eigs-base.cc	Tue Jan 26 16:30:54 2016 -0500
+++ b/liboctave/numeric/eigs-base.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -32,8 +32,7 @@
 #include "f77-fcn.h"
 #include "oct-locbuf.h"
 #include "quit.h"
-#include "SparsedbleLU.h"
-#include "SparseCmplxLU.h"
+#include "sparse-lu.h"
 #include "dSparse.h"
 #include "CSparse.h"
 #include "MatrixType.h"
@@ -473,7 +472,7 @@
       AminusSigmaB -= sigmat;
     }
 
-  SparseLU fact (AminusSigmaB);
+  sparse_lu<SparseMatrix> fact (AminusSigmaB);
 
   L = fact.L ();
   U = fact.U ();
@@ -637,7 +636,7 @@
       AminusSigmaB -= sigmat;
     }
 
-  SparseComplexLU fact (AminusSigmaB);
+  sparse_lu<SparseComplexMatrix> fact (AminusSigmaB);
 
   L = fact.L ();
   U = fact.U ();
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Tue Jan 26 16:30:54 2016 -0500
+++ b/liboctave/numeric/module.mk	Thu Jan 28 00:15:33 2016 -0500
@@ -80,10 +80,8 @@
   liboctave/numeric/randmtzig.h \
   liboctave/numeric/randpoisson.h \
   liboctave/numeric/sparse-chol.h \
-  liboctave/numeric/sparse-base-lu.h \
-  liboctave/numeric/SparseCmplxLU.h \
+  liboctave/numeric/sparse-lu.h \
   liboctave/numeric/SparseCmplxQR.h \
-  liboctave/numeric/SparsedbleLU.h \
   liboctave/numeric/SparseQR.h
 
 NUMERIC_C_SRC = \
@@ -144,12 +142,11 @@
   liboctave/numeric/oct-spparms.cc \
   liboctave/numeric/ODES.cc \
   liboctave/numeric/Quad.cc \
-  liboctave/numeric/SparseCmplxLU.cc \
   liboctave/numeric/SparseCmplxQR.cc \
-  liboctave/numeric/SparsedbleLU.cc \
   liboctave/numeric/SparseQR.cc \
   liboctave/numeric/sparse-chol-inst.cc \
-  $(NUMERIC_C_SRC)
+  liboctave/numeric/sparse-lu-inst.cc \
+$(NUMERIC_C_SRC)
 
 LIBOCTAVE_TEMPLATE_SRC += \
   liboctave/numeric/base-lu.cc \
@@ -157,7 +154,7 @@
   liboctave/numeric/bsxfun-defs.cc \
   liboctave/numeric/eigs-base.cc \
   liboctave/numeric/sparse-chol.cc \
-  liboctave/numeric/sparse-base-lu.cc \
+  liboctave/numeric/sparse-lu.cc \
   liboctave/numeric/sparse-dmsolve.cc
 
 ## Special rules for sources which must be built before rest of compilation.
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/sparse-base-lu.cc
--- a/liboctave/numeric/sparse-base-lu.cc	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,148 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "sparse-base-lu.h"
-
-#include "PermMatrix.h"
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-lu_type
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Y (void) const
-{
-  octave_idx_type nr = Lfact.rows ();
-  octave_idx_type nz = Lfact.cols ();
-  octave_idx_type nc = Ufact.cols ();
-
-  lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
-  octave_idx_type ii = 0;
-  Yout.xcidx (0) = 0;
-
-  for (octave_idx_type j = 0; j < nc; j++)
-    {
-      for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
-        {
-          Yout.xridx (ii) = Ufact.ridx (i);
-          Yout.xdata (ii++) = Ufact.data (i);
-        }
-      if (j < nz)
-        {
-          // Note the +1 skips the 1.0 on the diagonal
-          for (octave_idx_type i = Lfact.cidx (j) + 1;
-               i < Lfact.cidx (j +1); i++)
-            {
-              Yout.xridx (ii) = Lfact.ridx (i);
-              Yout.xdata (ii++) = Lfact.data (i);
-            }
-        }
-      Yout.xcidx (j + 1) = ii;
-    }
-
-  return Yout;
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-p_type
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr (void) const
-{
-
-  octave_idx_type nr = Lfact.rows ();
-
-  p_type Pout (nr, nr, nr);
-
-  for (octave_idx_type i = 0; i < nr; i++)
-    {
-      Pout.cidx (i) = i;
-      Pout.ridx (P (i)) = i;
-      Pout.data (i) = 1;
-    }
-  Pout.cidx (nr) = nr;
-
-  return Pout;
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-ColumnVector
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr_vec (void) const
-{
-
-  octave_idx_type nr = Lfact.rows ();
-
-  ColumnVector Pout (nr);
-
-  for (octave_idx_type i = 0; i < nr; i++)
-    Pout.xelem (i) = static_cast<double> (P(i) + 1);
-
-  return Pout;
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-PermMatrix
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr_mat (void) const
-{
-  return PermMatrix (P, false);
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-p_type
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc (void) const
-{
-  octave_idx_type nc = Ufact.cols ();
-
-  p_type Pout (nc, nc, nc);
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    {
-      Pout.cidx (i) = i;
-      Pout.ridx (i) = Q (i);
-      Pout.data (i) = 1;
-    }
-  Pout.cidx (nc) = nc;
-
-  return Pout;
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-ColumnVector
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc_vec (void) const
-{
-
-  octave_idx_type nc = Ufact.cols ();
-
-  ColumnVector Pout (nc);
-
-  for (octave_idx_type i = 0; i < nc; i++)
-    Pout.xelem (i) = static_cast<double> (Q(i) + 1);
-
-  return Pout;
-}
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-PermMatrix
-sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc_mat (void) const
-{
-  return PermMatrix (Q, true);
-}
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/sparse-base-lu.h
--- a/liboctave/numeric/sparse-base-lu.h	Tue Jan 26 16:30:54 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-/*
-
-Copyright (C) 2004-2015 David Bateman
-Copyright (C) 1998-2004 Andy Adler
-
-This file is part of Octave.
-
-Octave is free software; you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
-option) any later version.
-
-Octave is distributed in the hope that it will be useful, but WITHOUT
-ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-for more details.
-
-You should have received a copy of the GNU General Public License
-along with Octave; see the file COPYING.  If not, see
-<http://www.gnu.org/licenses/>.
-
-*/
-
-
-#if ! defined (octave_sparse_base_lu_h)
-#define octave_sparse_base_lu_h 1
-
-#include "MArray.h"
-#include "dSparse.h"
-
-template <typename lu_type, typename lu_elt_type, typename p_type, typename p_elt_type>
-class
-sparse_base_lu
-{
-public:
-
-  sparse_base_lu (void)
-    : Lfact (), Ufact (), Rfact (), cond (0), P (), Q () { }
-
-  sparse_base_lu (const sparse_base_lu& a)
-    : Lfact (a.Lfact), Ufact (a.Ufact), Rfact (), cond (a.cond),
-      P (a.P), Q (a.Q)
-  { }
-
-  sparse_base_lu& operator = (const sparse_base_lu& a)
-  {
-    if (this != &a)
-      {
-        Lfact = a.Lfact;
-        Ufact = a.Ufact;
-        cond = a.cond;
-        P = a.P;
-        Q = a.Q;
-      }
-    return *this;
-  }
-
-  virtual ~sparse_base_lu (void) { }
-
-  lu_type L (void) const { return Lfact; }
-
-  lu_type U (void) const { return Ufact; }
-
-  SparseMatrix R (void) const { return Rfact; }
-
-  lu_type Y (void) const;
-
-  p_type Pc (void) const;
-
-  p_type Pr (void) const;
-
-  ColumnVector Pc_vec (void) const;
-
-  ColumnVector Pr_vec (void) const;
-
-  PermMatrix Pc_mat (void) const;
-
-  PermMatrix Pr_mat (void) const;
-
-  const octave_idx_type * row_perm (void) const { return P.fortran_vec (); }
-
-  const octave_idx_type * col_perm (void) const { return Q.fortran_vec (); }
-
-  double rcond (void) const { return cond; }
-
-protected:
-
-  lu_type Lfact;
-  lu_type Ufact;
-  SparseMatrix Rfact;
-
-  double cond;
-
-  MArray<octave_idx_type> P;
-  MArray<octave_idx_type> Q;
-};
-
-#endif
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/sparse-lu-inst.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-lu-inst.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -0,0 +1,32 @@
+/*
+
+Copyright (C) 2016 John W. Eaton
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "sparse-lu.h"
+#include "sparse-lu.cc"
+
+template class sparse_lu<SparseMatrix>;
+
+template class sparse_lu<SparseComplexMatrix>;
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/sparse-lu.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-lu.cc	Thu Jan 28 00:15:33 2016 -0500
@@ -0,0 +1,884 @@
+/*
+
+Copyright (C) 2004-2015 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CSparse.h"
+#include "PermMatrix.h"
+#include "dSparse.h"
+#include "lo-error.h"
+#include "oct-locbuf.h"
+#include "oct-sparse.h"
+#include "oct-spparms.h"
+#include "sparse-lu.h"
+
+// Wrappers for SuiteSparse (formerly UMFPACK) functions that have
+// different names depending on the sparse matrix data type.
+//
+// All of these functions must be specialized to forward to the correct
+// SuiteSparse functions.
+
+template <typename T>
+void
+umfpack_defaults (double *Control);
+
+template <typename T>
+void
+umfpack_free_numeric (void **Numeric);
+
+template <typename T>
+void
+umfpack_free_symbolic (void **Symbolic);
+
+template <typename T>
+octave_idx_type
+umfpack_get_lunz (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric);
+
+template <typename T>
+octave_idx_type
+umfpack_get_numeric (octave_idx_type *Lp, octave_idx_type *Lj,
+                     T *Lx, // Or Lz_packed
+                     octave_idx_type *Up, octave_idx_type *Ui,
+                     T *Ux, // Or Uz_packed
+                     octave_idx_type *p, octave_idx_type *q,
+                     double *Dz_packed, octave_idx_type *do_recip,
+                     double *Rs, void *Numeric);
+
+template <typename T>
+octave_idx_type
+umfpack_numeric (const octave_idx_type *Ap, const octave_idx_type *Ai,
+                 const T *Ax, // Or Az_packed
+                 void *Symbolic, void **Numeric,
+                 const double *Control, double *Info);
+
+template <typename T>
+octave_idx_type
+umfpack_qsymbolic (octave_idx_type n_row, octave_idx_type n_col,
+                   const octave_idx_type *Ap, const octave_idx_type *Ai,
+                   const T *Ax, // Or Az_packed
+                   const octave_idx_type *Qinit, void **Symbolic,
+                   const double *Control, double *Info);
+
+template <typename T>
+void
+umfpack_report_control (const double *Control);
+
+template <typename T>
+void
+umfpack_report_info (const double *Control, const double *Info);
+
+template <typename T>
+void
+umfpack_report_matrix (octave_idx_type n_row, octave_idx_type n_col,
+                       const octave_idx_type *Ap, const octave_idx_type *Ai,
+                       const T *Ax, // Or Az_packed
+                       octave_idx_type col_form, const double *Control);
+
+template <typename T>
+void
+umfpack_report_numeric (void *Numeric, const double *Control);
+
+template <typename T>
+void
+umfpack_report_perm (octave_idx_type np, const octave_idx_type *Perm,
+                     const double *Control);
+
+template <typename T>
+void
+umfpack_report_status (double *Control, octave_idx_type status);
+
+template <typename T>
+void
+umfpack_report_symbolic (void *Symbolic, const double *Control);
+
+// SparseMatrix Specialization.
+
+template <>
+inline void
+umfpack_defaults<double> (double *Control)
+{
+  UMFPACK_DNAME (defaults) (Control);
+}
+
+template <>
+inline void
+umfpack_free_numeric<double> (void **Numeric)
+{
+  UMFPACK_DNAME (free_numeric) (Numeric);
+}
+
+template <>
+inline void
+umfpack_free_symbolic<double> (void **Symbolic)
+{
+  UMFPACK_DNAME (free_symbolic) (Symbolic);
+}
+
+template <>
+inline octave_idx_type
+umfpack_get_lunz<double>
+  (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
+{
+  octave_idx_type ignore1, ignore2, ignore3;
+
+  return UMFPACK_DNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
+                                   &ignore3, Numeric);
+}
+
+template <>
+inline octave_idx_type
+umfpack_get_numeric<double>
+  (octave_idx_type *Lp, octave_idx_type *Lj, double *Lx,
+   octave_idx_type *Up, octave_idx_type *Ui, double *Ux,
+   octave_idx_type *p, octave_idx_type *q, double *Dx,
+   octave_idx_type *do_recip, double *Rs, void *Numeric)
+{
+  return UMFPACK_DNAME (get_numeric) (Lp, Lj, Lx, Up, Ui, Ux, p, q, Dx,
+                                      do_recip, Rs, Numeric);
+}
+
+template <>
+inline octave_idx_type
+umfpack_numeric<double>
+  (const octave_idx_type *Ap, const octave_idx_type *Ai,
+   const double *Ax, void *Symbolic, void **Numeric,
+   const double *Control, double *Info)
+{
+  return UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, Numeric, Control,
+                                  Info);
+}
+
+template <>
+inline octave_idx_type
+umfpack_qsymbolic<double>
+  (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
+   const octave_idx_type *Ai, const double *Ax,
+   const octave_idx_type *Qinit, void **Symbolic,
+   const double *Control, double *Info)
+{
+  return UMFPACK_DNAME (qsymbolic) (n_row, n_col, Ap, Ai, Ax, Qinit,
+                                    Symbolic, Control, Info);
+}
+
+template <>
+inline void
+umfpack_report_control<double> (const double *Control)
+{
+  UMFPACK_DNAME (report_control) (Control);
+}
+
+template <>
+inline void
+umfpack_report_info<double> (const double *Control, const double *Info)
+{
+  UMFPACK_DNAME (report_info) (Control, Info);
+}
+
+template <>
+inline void
+umfpack_report_matrix<double>
+  (octave_idx_type n_row, octave_idx_type n_col, const octave_idx_type *Ap,
+   const octave_idx_type *Ai, const double *Ax, octave_idx_type col_form,
+   const double *Control)
+{
+  UMFPACK_DNAME (report_matrix) (n_row, n_col, Ap, Ai, Ax, col_form, Control);
+}
+
+template <>
+inline void
+umfpack_report_numeric<double> (void *Numeric, const double *Control)
+{
+  UMFPACK_DNAME (report_numeric) (Numeric, Control);
+}
+
+template <>
+inline void
+umfpack_report_perm<double>
+  (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
+{
+  UMFPACK_DNAME (report_perm) (np, Perm, Control);
+}
+
+template <>
+inline void
+umfpack_report_status<double> (double *Control, octave_idx_type status)
+{
+  UMFPACK_DNAME (report_status) (Control, status);
+}
+
+template <>
+inline void
+umfpack_report_symbolic<double> (void *Symbolic, const double *Control)
+{
+  UMFPACK_DNAME (report_symbolic) (Symbolic, Control);
+}
+
+// SparseComplexMatrix specialization.
+
+template <>
+inline void
+umfpack_defaults<Complex> (double *Control)
+{
+  UMFPACK_ZNAME (defaults) (Control);
+}
+
+template <>
+inline void
+umfpack_free_numeric<Complex> (void **Numeric)
+{
+  UMFPACK_ZNAME (free_numeric) (Numeric);
+}
+
+template <>
+inline void
+umfpack_free_symbolic<Complex> (void **Symbolic)
+{
+  UMFPACK_ZNAME (free_symbolic) (Symbolic);
+}
+
+template <>
+inline octave_idx_type
+umfpack_get_lunz<Complex>
+  (octave_idx_type *lnz, octave_idx_type *unz, void *Numeric)
+{
+  octave_idx_type ignore1, ignore2, ignore3;
+
+  return UMFPACK_ZNAME (get_lunz) (lnz, unz, &ignore1, &ignore2,
+                                   &ignore3, Numeric);
+}
+
+template <>
+inline octave_idx_type
+umfpack_get_numeric<Complex>
+  (octave_idx_type *Lp, octave_idx_type *Lj, Complex *Lz,
+   octave_idx_type *Up, octave_idx_type *Ui, Complex *Uz,
+   octave_idx_type *p, octave_idx_type *q, double *Dz,
+   octave_idx_type *do_recip, double *Rs, void *Numeric)
+{
+  return UMFPACK_ZNAME (get_numeric) (Lp, Lj,
+                                      reinterpret_cast<double *> (Lz),
+                                      0, Up, Ui,
+                                      reinterpret_cast<double *> (Uz),
+                                      0, p, q,
+                                      reinterpret_cast<double *> (Dz),
+                                      0, do_recip, Rs, Numeric);
+}
+
+template <>
+inline octave_idx_type
+umfpack_numeric<Complex>
+  (const octave_idx_type *Ap, const octave_idx_type *Ai,
+   const Complex *Az, void *Symbolic, void **Numeric,
+   const double *Control, double *Info)
+{
+  return UMFPACK_ZNAME (numeric) (Ap, Ai,
+                                  reinterpret_cast<const double *> (Az),
+                                  0, Symbolic, Numeric, Control, Info);
+}
+
+
+template <>
+inline octave_idx_type
+umfpack_qsymbolic<Complex>
+  (octave_idx_type n_row, octave_idx_type n_col,
+   const octave_idx_type *Ap, const octave_idx_type *Ai,
+   const Complex *Az, const octave_idx_type *Qinit,
+   void **Symbolic, const double *Control, double *Info)
+{
+  return UMFPACK_ZNAME (qsymbolic) (n_row, n_col, Ap, Ai,
+                                    reinterpret_cast<const double *> (Az),
+                                    0, Qinit, Symbolic, Control, Info);
+}
+
+
+template <>
+inline void
+umfpack_report_control<Complex> (const double *Control)
+{
+  UMFPACK_ZNAME (report_control) (Control);
+}
+
+template <>
+inline void
+umfpack_report_info<Complex> (const double *Control, const double *Info)
+{
+  UMFPACK_ZNAME (report_info) (Control, Info);
+}
+
+template <>
+inline void
+umfpack_report_matrix<Complex>
+  (octave_idx_type n_row, octave_idx_type n_col,
+   const octave_idx_type *Ap, const octave_idx_type *Ai,
+   const Complex *Az, octave_idx_type col_form, const double *Control)
+{
+  UMFPACK_ZNAME (report_matrix) (n_row, n_col, Ap, Ai,
+                                 reinterpret_cast<const double *> (Az),
+                                 0, col_form, Control);
+}
+
+template <>
+inline void
+umfpack_report_numeric<Complex> (void *Numeric, const double *Control)
+{
+  UMFPACK_ZNAME (report_numeric) (Numeric, Control);
+}
+
+template <>
+inline void
+umfpack_report_perm<Complex>
+  (octave_idx_type np, const octave_idx_type *Perm, const double *Control)
+{
+  UMFPACK_ZNAME (report_perm) (np, Perm, Control);
+}
+
+template <>
+inline void
+umfpack_report_status<Complex> (double *Control, octave_idx_type status)
+{
+  UMFPACK_ZNAME (report_status) (Control, status);
+}
+
+template <>
+inline void
+umfpack_report_symbolic <Complex> (void *Symbolic, const double *Control)
+{
+  UMFPACK_ZNAME (report_symbolic) (Symbolic, Control);
+}
+
+template <typename lu_type>
+sparse_lu<lu_type>::sparse_lu (const lu_type& a, const Matrix& piv_thres,
+                               bool scale)
+  : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
+{
+#ifdef HAVE_UMFPACK
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  // Setup the control parameters
+  Matrix Control (UMFPACK_CONTROL, 1);
+  double *control = Control.fortran_vec ();
+  umfpack_defaults<lu_elt_type> (control);
+
+  double tmp = octave_sparse_params::get_key ("spumoni");
+  if (! xisnan (tmp))
+    Control (UMFPACK_PRL) = tmp;
+
+  if (piv_thres.numel () == 2)
+    {
+      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
+      if (! xisnan (tmp))
+        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
+      if (! xisnan (tmp))
+        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+    }
+  else
+    {
+      tmp = octave_sparse_params::get_key ("piv_tol");
+      if (! xisnan (tmp))
+        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+
+      tmp = octave_sparse_params::get_key ("sym_tol");
+      if (! xisnan (tmp))
+        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+    }
+
+  // Set whether we are allowed to modify Q or not
+  tmp = octave_sparse_params::get_key ("autoamd");
+  if (! xisnan (tmp))
+    Control (UMFPACK_FIXQ) = tmp;
+
+  // Turn-off UMFPACK scaling for LU
+  if (scale)
+    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
+  else
+    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
+
+  umfpack_report_control<lu_elt_type> (control);
+
+  const octave_idx_type *Ap = a.cidx ();
+  const octave_idx_type *Ai = a.ridx ();
+  const lu_elt_type *Ax = a.data ();
+
+  umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);
+
+  void *Symbolic;
+  Matrix Info (1, UMFPACK_INFO);
+  double *info = Info.fortran_vec ();
+  int status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, 0, &Symbolic, control, info);
+
+  if (status < 0)
+    {
+      umfpack_report_status<lu_elt_type> (control, status);
+      umfpack_report_info<lu_elt_type> (control, info);
+
+      umfpack_free_symbolic<lu_elt_type> (&Symbolic);
+
+      (*current_liboctave_error_handler)
+        ("sparse_lu: symbolic factorization failed");
+    }
+  else
+    {
+      umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
+
+      void *Numeric;
+      status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
+      umfpack_free_symbolic<lu_elt_type> (&Symbolic);
+
+      cond = Info (UMFPACK_RCOND);
+
+      if (status < 0)
+        {
+          umfpack_report_status<lu_elt_type> (control, status);
+          umfpack_report_info<lu_elt_type> (control, info);
+
+          umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+          (*current_liboctave_error_handler)
+            ("sparse_lu: numeric factorization failed");
+        }
+      else
+        {
+          umfpack_report_numeric<lu_elt_type> (Numeric, control);
+
+          octave_idx_type lnz, unz;
+          status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
+
+          if (status < 0)
+            {
+              umfpack_report_status<lu_elt_type> (control, status);
+              umfpack_report_info<lu_elt_type> (control, info);
+
+              umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+              (*current_liboctave_error_handler)
+                ("sparse_lu: extracting LU factors failed");
+            }
+          else
+            {
+              octave_idx_type n_inner = (nr < nc ? nr : nc);
+
+              if (lnz < 1)
+                Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
+              else
+                Lfact = lu_type (n_inner, nr, lnz);
+
+              octave_idx_type *Ltp = Lfact.cidx ();
+              octave_idx_type *Ltj = Lfact.ridx ();
+              lu_elt_type *Ltx = Lfact.data ();
+
+              if (unz < 1)
+                Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
+              else
+                Ufact = lu_type (n_inner, nc, unz);
+
+              octave_idx_type *Up = Ufact.cidx ();
+              octave_idx_type *Uj = Ufact.ridx ();
+              lu_elt_type *Ux = Ufact.data ();
+
+              Rfact = SparseMatrix (nr, nr, nr);
+              for (octave_idx_type i = 0; i < nr; i++)
+                {
+                  Rfact.xridx (i) = i;
+                  Rfact.xcidx (i) = i;
+                }
+              Rfact.xcidx (nr) = nr;
+              double *Rx = Rfact.data ();
+
+              P.resize (dim_vector (nr, 1));
+              octave_idx_type *p = P.fortran_vec ();
+
+              Q.resize (dim_vector (nc, 1));
+              octave_idx_type *q = Q.fortran_vec ();
+
+              octave_idx_type do_recip;
+              status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);
+
+              umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+              if (status < 0)
+                {
+                  umfpack_report_status<lu_elt_type> (control, status);
+
+                  (*current_liboctave_error_handler)
+                    ("sparse_lu: extracting LU factors failed");
+                }
+              else
+                {
+                  Lfact = Lfact.transpose ();
+
+                  if (do_recip)
+                    for (octave_idx_type i = 0; i < nr; i++)
+                      Rx[i] = 1.0 / Rx[i];
+
+                  umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
+                  umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
+                  umfpack_report_perm<lu_elt_type> (nr, p, control);
+                  umfpack_report_perm<lu_elt_type> (nc, q, control);
+                }
+
+              umfpack_report_info<lu_elt_type> (control, info);
+            }
+        }
+    }
+
+#else
+  (*current_liboctave_error_handler)
+    ("support for UMFPACK was unavailable or disabled when liboctave was built");
+#endif
+}
+
+template <typename lu_type>
+sparse_lu<lu_type>::sparse_lu (const lu_type& a,
+                               const ColumnVector& Qinit,
+                               const Matrix& piv_thres, bool scale,
+                               bool FixedQ, double droptol,
+                               bool milu, bool udiag)
+  : Lfact (), Ufact (), Rfact (), cond (0), P (), Q ()
+{
+#ifdef HAVE_UMFPACK
+  if (milu)
+    (*current_liboctave_error_handler)
+      ("Modified incomplete LU not implemented");
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  // Setup the control parameters
+  Matrix Control (UMFPACK_CONTROL, 1);
+  double *control = Control.fortran_vec ();
+  umfpack_defaults<lu_elt_type> (control);
+
+  double tmp = octave_sparse_params::get_key ("spumoni");
+  if (! xisnan (tmp))
+    Control (UMFPACK_PRL) = tmp;
+
+  if (piv_thres.numel () == 2)
+    {
+      tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0));
+      if (! xisnan (tmp))
+        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+      tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1));
+      if (! xisnan (tmp))
+        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+    }
+  else
+    {
+      tmp = octave_sparse_params::get_key ("piv_tol");
+      if (! xisnan (tmp))
+        Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
+
+      tmp = octave_sparse_params::get_key ("sym_tol");
+      if (! xisnan (tmp))
+        Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
+    }
+
+  if (droptol >= 0.)
+    Control (UMFPACK_DROPTOL) = droptol;
+
+  // Set whether we are allowed to modify Q or not
+  if (FixedQ)
+    Control (UMFPACK_FIXQ) = 1.0;
+  else
+    {
+      tmp = octave_sparse_params::get_key ("autoamd");
+      if (! xisnan (tmp))
+        Control (UMFPACK_FIXQ) = tmp;
+    }
+
+  // Turn-off UMFPACK scaling for LU
+  if (scale)
+    Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM;
+  else
+    Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
+
+  umfpack_report_control<lu_elt_type> (control);
+
+  const octave_idx_type *Ap = a.cidx ();
+  const octave_idx_type *Ai = a.ridx ();
+  const lu_elt_type *Ax = a.data ();
+
+  umfpack_report_matrix<lu_elt_type> (nr, nc, Ap, Ai, Ax, static_cast<octave_idx_type> (1), control);
+
+  void *Symbolic;
+  Matrix Info (1, UMFPACK_INFO);
+  double *info = Info.fortran_vec ();
+  int status;
+
+  // Null loop so that qinit is imediately deallocated when not needed
+  do
+    {
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc);
+
+      for (octave_idx_type i = 0; i < nc; i++)
+        qinit[i] = static_cast<octave_idx_type> (Qinit (i));
+
+      status = umfpack_qsymbolic<lu_elt_type> (nr, nc, Ap, Ai, Ax, qinit, &Symbolic, control, info);
+    }
+  while (0);
+
+  if (status < 0)
+    {
+      umfpack_report_status<lu_elt_type> (control, status);
+      umfpack_report_info<lu_elt_type> (control, info);
+
+      umfpack_free_symbolic<lu_elt_type> (&Symbolic);
+
+      (*current_liboctave_error_handler)
+        ("sparse_lu: symbolic factorization failed");
+    }
+  else
+    {
+      umfpack_report_symbolic<lu_elt_type> (Symbolic, control);
+
+      void *Numeric;
+      status = umfpack_numeric<lu_elt_type> (Ap, Ai, Ax, Symbolic, &Numeric, control, info);
+      umfpack_free_symbolic<lu_elt_type> (&Symbolic);
+
+      cond = Info (UMFPACK_RCOND);
+
+      if (status < 0)
+        {
+          umfpack_report_status<lu_elt_type> (control, status);
+          umfpack_report_info<lu_elt_type> (control, info);
+
+          umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+          (*current_liboctave_error_handler)
+            ("sparse_lu: numeric factorization failed");
+        }
+      else
+        {
+          umfpack_report_numeric<lu_elt_type> (Numeric, control);
+
+          octave_idx_type lnz, unz;
+          status = umfpack_get_lunz<lu_elt_type> (&lnz, &unz, Numeric);
+
+          if (status < 0)
+            {
+              umfpack_report_status<lu_elt_type> (control, status);
+              umfpack_report_info<lu_elt_type> (control, info);
+
+              umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+              (*current_liboctave_error_handler)
+                ("sparse_lu: extracting LU factors failed");
+            }
+          else
+            {
+              octave_idx_type n_inner = (nr < nc ? nr : nc);
+
+              if (lnz < 1)
+                Lfact = lu_type (n_inner, nr, static_cast<octave_idx_type> (1));
+              else
+                Lfact = lu_type (n_inner, nr, lnz);
+
+              octave_idx_type *Ltp = Lfact.cidx ();
+              octave_idx_type *Ltj = Lfact.ridx ();
+              lu_elt_type *Ltx = Lfact.data ();
+
+              if (unz < 1)
+                Ufact = lu_type (n_inner, nc, static_cast<octave_idx_type> (1));
+              else
+                Ufact = lu_type (n_inner, nc, unz);
+
+              octave_idx_type *Up = Ufact.cidx ();
+              octave_idx_type *Uj = Ufact.ridx ();
+              lu_elt_type *Ux = Ufact.data ();
+
+              Rfact = SparseMatrix (nr, nr, nr);
+              for (octave_idx_type i = 0; i < nr; i++)
+                {
+                  Rfact.xridx (i) = i;
+                  Rfact.xcidx (i) = i;
+                }
+              Rfact.xcidx (nr) = nr;
+              double *Rx = Rfact.data ();
+
+              P.resize (dim_vector (nr, 1));
+              octave_idx_type *p = P.fortran_vec ();
+
+              Q.resize (dim_vector (nc, 1));
+              octave_idx_type *q = Q.fortran_vec ();
+
+              octave_idx_type do_recip;
+              status = umfpack_get_numeric<lu_elt_type> (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, 0, &do_recip, Rx, Numeric);
+
+              umfpack_free_numeric<lu_elt_type> (&Numeric);
+
+              if (status < 0)
+                {
+                  umfpack_report_status<lu_elt_type> (control, status);
+
+                  (*current_liboctave_error_handler)
+                    ("sparse_lu: extracting LU factors failed");
+                }
+              else
+                {
+                  Lfact = Lfact.transpose ();
+
+                  if (do_recip)
+                    for (octave_idx_type i = 0; i < nr; i++)
+                      Rx[i] = 1.0 / Rx[i];
+
+                  umfpack_report_matrix<lu_elt_type> (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), static_cast<octave_idx_type> (1), control);
+                  umfpack_report_matrix<lu_elt_type> (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), static_cast<octave_idx_type> (1), control);
+                  umfpack_report_perm<lu_elt_type> (nr, p, control);
+                  umfpack_report_perm<lu_elt_type> (nc, q, control);
+                }
+
+              umfpack_report_info<lu_elt_type> (control, info);
+            }
+        }
+    }
+
+  if (udiag)
+    (*current_liboctave_error_handler)
+      ("Option udiag of incomplete LU not implemented");
+
+#else
+  (*current_liboctave_error_handler)
+    ("support for UMFPACK was unavailable or disabled when liboctave was built");
+#endif
+}
+
+template <typename lu_type>
+lu_type
+sparse_lu<lu_type>::Y (void) const
+{
+  octave_idx_type nr = Lfact.rows ();
+  octave_idx_type nz = Lfact.cols ();
+  octave_idx_type nc = Ufact.cols ();
+
+  lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz () - (nr<nz?nr:nz));
+  octave_idx_type ii = 0;
+  Yout.xcidx (0) = 0;
+
+  for (octave_idx_type j = 0; j < nc; j++)
+    {
+      for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
+        {
+          Yout.xridx (ii) = Ufact.ridx (i);
+          Yout.xdata (ii++) = Ufact.data (i);
+        }
+      if (j < nz)
+        {
+          // Note the +1 skips the 1.0 on the diagonal
+          for (octave_idx_type i = Lfact.cidx (j) + 1;
+               i < Lfact.cidx (j +1); i++)
+            {
+              Yout.xridx (ii) = Lfact.ridx (i);
+              Yout.xdata (ii++) = Lfact.data (i);
+            }
+        }
+      Yout.xcidx (j + 1) = ii;
+    }
+
+  return Yout;
+}
+
+template <typename lu_type>
+SparseMatrix
+sparse_lu<lu_type>::Pr (void) const
+{
+
+  octave_idx_type nr = Lfact.rows ();
+
+  SparseMatrix Pout (nr, nr, nr);
+
+  for (octave_idx_type i = 0; i < nr; i++)
+    {
+      Pout.cidx (i) = i;
+      Pout.ridx (P (i)) = i;
+      Pout.data (i) = 1;
+    }
+  Pout.cidx (nr) = nr;
+
+  return Pout;
+}
+
+template <typename lu_type>
+ColumnVector
+sparse_lu<lu_type>::Pr_vec (void) const
+{
+
+  octave_idx_type nr = Lfact.rows ();
+
+  ColumnVector Pout (nr);
+
+  for (octave_idx_type i = 0; i < nr; i++)
+    Pout.xelem (i) = static_cast<double> (P(i) + 1);
+
+  return Pout;
+}
+
+template <typename lu_type>
+PermMatrix
+sparse_lu<lu_type>::Pr_mat (void) const
+{
+  return PermMatrix (P, false);
+}
+
+template <typename lu_type>
+SparseMatrix
+sparse_lu<lu_type>::Pc (void) const
+{
+  octave_idx_type nc = Ufact.cols ();
+
+  SparseMatrix Pout (nc, nc, nc);
+
+  for (octave_idx_type i = 0; i < nc; i++)
+    {
+      Pout.cidx (i) = i;
+      Pout.ridx (i) = Q (i);
+      Pout.data (i) = 1;
+    }
+  Pout.cidx (nc) = nc;
+
+  return Pout;
+}
+
+template <typename lu_type>
+ColumnVector
+sparse_lu<lu_type>::Pc_vec (void) const
+{
+
+  octave_idx_type nc = Ufact.cols ();
+
+  ColumnVector Pout (nc);
+
+  for (octave_idx_type i = 0; i < nc; i++)
+    Pout.xelem (i) = static_cast<double> (Q(i) + 1);
+
+  return Pout;
+}
+
+template <typename lu_type>
+PermMatrix
+sparse_lu<lu_type>::Pc_mat (void) const
+{
+  return PermMatrix (Q, true);
+}
diff -r 307096fb67e1 -r ea9c05014809 liboctave/numeric/sparse-lu.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-lu.h	Thu Jan 28 00:15:33 2016 -0500
@@ -0,0 +1,109 @@
+/*
+
+Copyright (C) 2004-2015 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+
+This file is part of Octave.
+
+Octave is free software; you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+Octave is distributed in the hope that it will be useful, but WITHOUT
+ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with Octave; see the file COPYING.  If not, see
+<http://www.gnu.org/licenses/>.
+
+*/
+
+
+#if ! defined (octave_sparse_lu_h)
+#define octave_sparse_lu_h 1
+
+#include "MArray.h"
+#include "dSparse.h"
+
+template <typename lu_type>
+class
+sparse_lu
+{
+public:
+
+  typedef typename lu_type::element_type lu_elt_type;
+  
+  sparse_lu (void)
+    : Lfact (), Ufact (), Rfact (), cond (0), P (), Q () { }
+
+  sparse_lu (const lu_type& a, const Matrix& piv_thres = Matrix (),
+             bool scale = false);
+
+  sparse_lu (const lu_type& a, const ColumnVector& Qinit,
+             const Matrix& piv_thres, bool scale = false,
+             bool FixedQ = false, double droptol = -1.0,
+             bool milu = false, bool udiag = false);
+
+  sparse_lu (const sparse_lu& a)
+    : Lfact (a.Lfact), Ufact (a.Ufact), Rfact (), cond (a.cond),
+      P (a.P), Q (a.Q)
+  { }
+
+  sparse_lu& operator = (const sparse_lu& a)
+  {
+    if (this != &a)
+      {
+        Lfact = a.Lfact;
+        Ufact = a.Ufact;
+        cond = a.cond;
+        P = a.P;
+        Q = a.Q;
+      }
+
+    return *this;
+  }
+
+  virtual ~sparse_lu (void) { }
+
+  lu_type L (void) const { return Lfact; }
+
+  lu_type U (void) const { return Ufact; }
+
+  SparseMatrix R (void) const { return Rfact; }
+
+  lu_type Y (void) const;
+
+  SparseMatrix Pc (void) const;
+
+  SparseMatrix Pr (void) const;
+
+  ColumnVector Pc_vec (void) const;
+
+  ColumnVector Pr_vec (void) const;
+
+  PermMatrix Pc_mat (void) const;
+
+  PermMatrix Pr_mat (void) const;
+
+  const octave_idx_type * row_perm (void) const { return P.fortran_vec (); }
+
+  const octave_idx_type * col_perm (void) const { return Q.fortran_vec (); }
+
+  double rcond (void) const { return cond; }
+
+protected:
+
+  lu_type Lfact;
+  lu_type Ufact;
+  SparseMatrix Rfact;
+
+  double cond;
+
+  MArray<octave_idx_type> P;
+  MArray<octave_idx_type> Q;
+};
+
+#endif