# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1454434446 18000
# Node ID 791dcb32b657722923a0fa9d0466ca6de8e3626e
# Parent  f45f4f888db5339654c7b965eafa79d2702d16c1
revamp sparse QR factorizatino classes

* sparse-qr-inst.cc, sparse-qr.cc, sparse-qr.h: New files, adapted
from SparseCmplxQR.cc, SparseCmplxQR.h, SparseQR.cc, and SparseQR.h.
(sparse_qr): New template class.
* SparseCmplxQR.cc, SparseCmplxQR.h, SparseQR.cc, SparseQR.h: Delete.
* liboctave/numeric/module.mk: Update.

* dmperm.cc, qr.cc, CSparse.cc, dSparse.cc, sparse-dmsolve.cc:
Use new sparse_qr class.

* oct-sparse.h (CXSPARSE_DNAME, CXSPARSE_ZNAME):
Move macros here from SparseQR.h and SparseCmplxQR.h.

diff -r f45f4f888db5 -r 791dcb32b657 libinterp/dldfcn/dmperm.cc
--- a/libinterp/dldfcn/dmperm.cc	Tue Feb 02 12:33:38 2016 -0500
+++ b/libinterp/dldfcn/dmperm.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -34,8 +34,7 @@
 #include "oct-sparse.h"
 #include "ov-re-sparse.h"
 #include "ov-cx-sparse.h"
-#include "SparseQR.h"
-#include "SparseCmplxQR.h"
+#include "sparse-qr.h"
 
 #if defined (ENABLE_64)
 # define CXSPARSE_NAME(name) cs_dl ## name
diff -r f45f4f888db5 -r 791dcb32b657 libinterp/dldfcn/qr.cc
--- a/libinterp/dldfcn/qr.cc	Tue Feb 02 12:33:38 2016 -0500
+++ b/libinterp/dldfcn/qr.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -34,8 +34,7 @@
 #include "fCmplxQRP.h"
 #include "floatQR.h"
 #include "floatQRP.h"
-#include "SparseQR.h"
-#include "SparseCmplxQR.h"
+#include "sparse-qr.h"
 
 
 #include "defun-dld.h"
@@ -239,7 +238,7 @@
 
       if (is_cmplx)
         {
-          SparseComplexQR q (arg.sparse_complex_matrix_value ());
+          sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ());
 
           if (have_b > 0)
             {
@@ -255,7 +254,7 @@
         }
       else
         {
-          SparseQR q (arg.sparse_matrix_value ());
+          sparse_qr<SparseMatrix> q (arg.sparse_matrix_value ());
 
           if (have_b > 0)
             {
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/array/CSparse.cc
--- a/liboctave/array/CSparse.cc	Tue Feb 02 12:33:38 2016 -0500
+++ b/liboctave/array/CSparse.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -55,7 +55,7 @@
 #include "oct-sparse.h"
 #include "sparse-util.h"
 #include "sparse-chol.h"
-#include "SparseCmplxQR.h"
+#include "sparse-qr.h"
 
 #include "Sparse-op-defs.h"
 
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/array/dSparse.cc
--- a/liboctave/array/dSparse.cc	Tue Feb 02 12:33:38 2016 -0500
+++ b/liboctave/array/dSparse.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -50,7 +50,7 @@
 #include "oct-sparse.h"
 #include "sparse-util.h"
 #include "sparse-chol.h"
-#include "SparseQR.h"
+#include "sparse-qr.h"
 
 #include "Sparse-op-defs.h"
 
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/SparseCmplxQR.cc
--- a/liboctave/numeric/SparseCmplxQR.cc	Tue Feb 02 12:33:38 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,736 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-
-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 "SparseCmplxQR.h"
-#include "oct-locbuf.h"
-
-SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep
-  (OCTAVE_UNUSED const SparseComplexMatrix& a, OCTAVE_UNUSED int order)
-  : count (1), nrows (0)
-#ifdef HAVE_CXSPARSE
-    , S (0), N (0)
-#endif
-{
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_ZNAME () A;
-  A.nzmax = a.nnz ();
-  A.m = a.rows ();
-  A.n = a.cols ();
-  nrows = A.m;
-  // Cast away const on A, with full knowledge that CSparse won't touch it
-  // Prevents the methods below making a copy of the data.
-  A.p = const_cast<octave_idx_type *>(a.cidx ());
-  A.i = const_cast<octave_idx_type *>(a.ridx ());
-  A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *>
-                                   (a.data ()));
-  A.nz = -1;
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
-  N = CXSPARSE_ZNAME (_qr) (&A, S);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  if (! N)
-    (*current_liboctave_error_handler)
-      ("SparseComplexQR: sparse matrix QR factorization filled");
-
-  count = 1;
-
-#else
-  (*current_liboctave_error_handler)
-    ("SparseComplexQR: support for CXSparse was unavailable or disabled when liboctave was built");
-#endif
-}
-
-SparseComplexQR::SparseComplexQR_rep::~SparseComplexQR_rep (void)
-{
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_ZNAME (_sfree) (S);
-  CXSPARSE_ZNAME (_nfree) (N);
-#endif
-}
-
-SparseComplexMatrix
-SparseComplexQR::SparseComplexQR_rep::V (void) const
-{
-#ifdef HAVE_CXSPARSE
-  // Drop zeros from V and sort
-  // FIXME: Is the double transpose to sort necessary?
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  CXSPARSE_ZNAME (_dropzeros) (N->L);
-  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
-  CXSPARSE_ZNAME (_spfree) (N->L);
-  N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
-  CXSPARSE_ZNAME (_spfree) (D);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nz = N->L->nzmax;
-  SparseComplexMatrix ret (N->L->m, nc, nz);
-  for (octave_idx_type j = 0; j < nc+1; j++)
-    ret.xcidx (j) = N->L->p[j];
-  for (octave_idx_type j = 0; j < nz; j++)
-    {
-      ret.xridx (j) = N->L->i[j];
-      ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j];
-    }
-  return ret;
-#else
-  return SparseComplexMatrix ();
-#endif
-}
-
-ColumnVector
-SparseComplexQR::SparseComplexQR_rep::Pinv (void) const
-{
-#ifdef HAVE_CXSPARSE
-  ColumnVector ret(N->L->m);
-  for (octave_idx_type i = 0; i < N->L->m; i++)
-    ret.xelem (i) = S->pinv[i];
-  return ret;
-#else
-  return ColumnVector ();
-#endif
-}
-
-ColumnVector
-SparseComplexQR::SparseComplexQR_rep::P (void) const
-{
-#ifdef HAVE_CXSPARSE
-  ColumnVector ret(N->L->m);
-  for (octave_idx_type i = 0; i < N->L->m; i++)
-    ret.xelem (S->pinv[i]) = i;
-  return ret;
-#else
-  return ColumnVector ();
-#endif
-}
-
-SparseComplexMatrix
-SparseComplexQR::SparseComplexQR_rep::R (const bool econ) const
-{
-#ifdef HAVE_CXSPARSE
-  // Drop zeros from R and sort
-  // FIXME: Is the double transpose to sort necessary?
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  CXSPARSE_ZNAME (_dropzeros) (N->U);
-  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
-  CXSPARSE_ZNAME (_spfree) (N->U);
-  N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
-  CXSPARSE_ZNAME (_spfree) (D);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-  octave_idx_type nc = N->U->n;
-  octave_idx_type nz = N->U->nzmax;
-  SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
-  for (octave_idx_type j = 0; j < nc+1; j++)
-    ret.xcidx (j) = N->U->p[j];
-  for (octave_idx_type j = 0; j < nz; j++)
-    {
-      ret.xridx (j) = N->U->i[j];
-      ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j];
-    }
-  return ret;
-#else
-  return SparseComplexMatrix ();
-#endif
-}
-
-ComplexMatrix
-SparseComplexQR::SparseComplexQR_rep::C (const ComplexMatrix &b) const
-{
-#ifdef HAVE_CXSPARSE
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nr = nrows;
-  const cs_complex_t *bvec =
-    reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
-  ComplexMatrix ret(b_nr, b_nc);
-  Complex *vec = ret.fortran_vec ();
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler) ("matrix dimension mismatch");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
-      for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
-        {
-          octave_quit ();
-          volatile octave_idx_type nm = (nr < nc ? nr : nc);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec)
-            (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type i = 0; i < nm; i++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly)
-                (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          for (octave_idx_type i = 0; i < b_nr; i++)
-            vec[i+idx] = buf[i];
-        }
-    }
-  return ret;
-#else
-  return ComplexMatrix ();
-#endif
-}
-
-ComplexMatrix
-SparseComplexQR::SparseComplexQR_rep::Q (void) const
-{
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nr = nrows;
-  ComplexMatrix ret(nr, nr);
-  Complex *vec = ret.fortran_vec ();
-  if (nr < 0 || nc < 0)
-    (*current_liboctave_error_handler) ("matrix dimension mismatch");
-
-  if (nr == 0 || nc == 0)
-    ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr);
-      for (octave_idx_type i = 0; i < nr; i++)
-        bvec[i] = cs_complex_t (0.0, 0.0);
-      OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
-      for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
-        {
-          octave_quit ();
-          bvec[j] = cs_complex_t (1.0, 0.0);
-          volatile octave_idx_type nm = (nr < nc ? nr : nc);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec)
-            (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type i = 0; i < nm; i++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly)
-                (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          for (octave_idx_type i = 0; i < nr; i++)
-            vec[i+idx] = buf[i];
-          bvec[j] = cs_complex_t (0.0, 0.0);
-        }
-    }
-  return ret.hermitian ();
-#else
-  return ComplexMatrix ();
-#endif
-}
-
-ComplexMatrix
-qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
-         octave_idx_type& info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  ComplexMatrix x;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
-  else if (nr >= nc)
-    {
-      SparseComplexQR q (a, 2);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec)
-            (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseComplexMatrix at = a.hermitian ();
-      SparseComplexQR q (at, 2);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
-      OCTAVE_LOCAL_BUFFER (double, B, nr);
-      for (octave_idx_type i = 0; i < nr; i++)
-        B[i] = q.N ()->B[i];
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec)
-            (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-
-  return x;
-#else
-  return ComplexMatrix ();
-#endif
-}
-
-SparseComplexMatrix
-qrsolve (const SparseComplexMatrix&a, const SparseMatrix &b,
-         octave_idx_type &info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  SparseComplexMatrix x;
-  volatile octave_idx_type ii, x_nz;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = SparseComplexMatrix (nc, b_nc);
-  else if (nr >= nc)
-    {
-      SparseComplexQR q (a, 2);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec)
-            (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_ZNAME (_ipvec)
-            (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseComplexMatrix at = a.hermitian ();
-      SparseComplexQR q (at, 2);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
-
-      OCTAVE_LOCAL_BUFFER (double, B, nr);
-      for (octave_idx_type i = 0; i < nr; i++)
-        B[i] = q.N ()->B[i];
-
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec)
-            (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec)
-            (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-
-  x.maybe_compress ();
-  return x;
-#else
-  return SparseComplexMatrix ();
-#endif
-}
-
-ComplexMatrix
-qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
-         octave_idx_type& info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  const cs_complex_t *bvec =
-    reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
-  ComplexMatrix x;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
-  else if (nr >= nc)
-    {
-      SparseComplexQR q (a, 2);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
-                          (x.fortran_vec ());
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
-           i++, idx+=nc, bidx+=b_nr)
-        {
-          octave_quit ();
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseComplexMatrix at = a.hermitian ();
-      SparseComplexQR q (at, 2);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
-      OCTAVE_LOCAL_BUFFER (double, B, nr);
-      for (octave_idx_type i = 0; i < nr; i++)
-        B[i] = q.N ()->B[i];
-      for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
-           i++, idx+=nc, bidx+=b_nr)
-        {
-          octave_quit ();
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
-          CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-
-  return x;
-#else
-  return ComplexMatrix ();
-#endif
-}
-
-SparseComplexMatrix
-qrsolve (const SparseComplexMatrix&a, const SparseComplexMatrix &b,
-         octave_idx_type &info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  SparseComplexMatrix x;
-  volatile octave_idx_type ii, x_nz;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = SparseComplexMatrix (nc, b_nc);
-  else if (nr >= nc)
-    {
-      SparseComplexQR q (a, 2);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_ipvec)
-            (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_ZNAME (_ipvec)
-            (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseComplexMatrix at = a.hermitian ();
-      SparseComplexQR q (at, 2);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
-      OCTAVE_LOCAL_BUFFER (double, B, nr);
-      for (octave_idx_type i = 0; i < nr; i++)
-        B[i] = q.N ()->B[i];
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = cs_complex_t (0.0, 0.0);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec)
-            (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
-          CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_ZNAME (_pvec)
-            (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-
-  x.maybe_compress ();
-  return x;
-#else
-  return SparseComplexMatrix ();
-#endif
-}
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/SparseCmplxQR.h
--- a/liboctave/numeric/SparseCmplxQR.h	Tue Feb 02 12:33:38 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,173 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-
-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_SparseCmplxQR_h)
-#define octave_SparseCmplxQR_h 1
-
-#include <iosfwd>
-
-#include "dMatrix.h"
-#include "CMatrix.h"
-#include "dSparse.h"
-#include "CSparse.h"
-#include "oct-sparse.h"
-
-#if defined (ENABLE_64)
-#define CXSPARSE_ZNAME(name) cs_cl ## name
-#else
-#define CXSPARSE_ZNAME(name) cs_ci ## name
-#endif
-
-class
-OCTAVE_API
-SparseComplexQR
-{
-protected:
-  class SparseComplexQR_rep
-  {
-  public:
-    SparseComplexQR_rep (const SparseComplexMatrix& a, int order);
-
-    ~SparseComplexQR_rep (void);
-#ifdef HAVE_CXSPARSE
-    bool ok (void) const { return (N && S); }
-#else
-    bool ok (void) const { return false; }
-#endif
-    SparseComplexMatrix V (void) const;
-
-    ColumnVector Pinv (void) const;
-
-    ColumnVector P (void) const;
-
-    SparseComplexMatrix R (const bool econ) const;
-
-    ComplexMatrix C (const ComplexMatrix &b) const;
-
-    ComplexMatrix Q (void) const;
-
-    octave_refcount<int> count;
-
-    octave_idx_type nrows;
-#ifdef HAVE_CXSPARSE
-    CXSPARSE_ZNAME (s) *S;
-
-    CXSPARSE_ZNAME (n) *N;
-#endif
-  private:
-
-    // No copying!
-
-    SparseComplexQR_rep (const SparseComplexQR_rep&);
-
-    SparseComplexQR_rep operator = (const SparseComplexQR_rep&);
-
-  };
-private:
-  SparseComplexQR_rep *rep;
-
-public:
-  SparseComplexQR (void) :
-    rep (new SparseComplexQR_rep (SparseComplexMatrix (), 0)) { }
-
-  SparseComplexQR (const SparseComplexMatrix& a, int order = 0) :
-    rep (new SparseComplexQR_rep (a, order)) { }
-
-  SparseComplexQR (const SparseComplexQR& a) : rep (a.rep) { rep->count++; }
-
-  ~SparseComplexQR (void)
-  {
-    if (--rep->count == 0)
-      delete rep;
-  }
-
-  SparseComplexQR& operator = (const SparseComplexQR& a)
-  {
-    if (this != &a)
-      {
-        if (--rep->count == 0)
-          delete rep;
-
-        rep = a.rep;
-        rep->count++;
-      }
-    return *this;
-  }
-
-  bool ok (void) const { return rep->ok (); }
-
-  SparseComplexMatrix V (void) const { return rep->V (); }
-
-  ColumnVector Pinv (void) const { return rep->P (); }
-
-  ColumnVector P (void) const { return rep->P (); }
-
-  SparseComplexMatrix R (const bool econ = false) const
-  { return rep->R(econ); }
-
-  ComplexMatrix C (const ComplexMatrix &b) const { return rep->C(b); }
-
-  ComplexMatrix Q (void) const { return rep->Q (); }
-
-  friend ComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                const MArray<double> &b,
-                                octave_idx_type &info);
-
-  friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                      const SparseMatrix &b,
-                                      octave_idx_type &info);
-
-  friend ComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                const MArray<Complex> &b,
-                                octave_idx_type &info);
-
-  friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                      const SparseComplexMatrix &b,
-                                      octave_idx_type &info);
-
-protected:
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_ZNAME (s) * S (void) { return rep->S; }
-
-  CXSPARSE_ZNAME (n) * N (void) { return rep->N; }
-#endif
-};
-
-
-// Publish externally used friend functions.
-
-extern ComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                              const MArray<double> &b,
-                              octave_idx_type &info);
-
-extern SparseComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                    const SparseMatrix &b,
-                                    octave_idx_type &info);
-
-extern ComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                              const MArray<Complex> &b,
-                              octave_idx_type &info);
-
-extern SparseComplexMatrix qrsolve (const SparseComplexMatrix &a,
-                                    const SparseComplexMatrix &b,
-                                    octave_idx_type &info);
-#endif
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/SparseQR.cc
--- a/liboctave/numeric/SparseQR.cc	Tue Feb 02 12:33:38 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,795 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-
-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 "SparseQR.h"
-#include "oct-locbuf.h"
-
-SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order)
-  : count (1), nrows (0)
-#ifdef HAVE_CXSPARSE
-    , S (0), N (0)
-#endif
-{
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_DNAME () A;
-  A.nzmax = a.nnz ();
-  A.m = a.rows ();
-  A.n = a.cols ();
-  nrows = A.m;
-  // Cast away const on A, with full knowledge that CSparse won't touch it
-  // Prevents the methods below making a copy of the data.
-  A.p = const_cast<octave_idx_type *>(a.cidx ());
-  A.i = const_cast<octave_idx_type *>(a.ridx ());
-  A.x = const_cast<double *>(a.data ());
-  A.nz = -1;
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
-
-  N = CXSPARSE_DNAME (_qr) (&A, S);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  if (! N)
-    (*current_liboctave_error_handler)
-      ("SparseQR: sparse matrix QR factorization filled");
-
-  count = 1;
-
-#else
-  (*current_liboctave_error_handler)
-    ("SparseQR: support for CXSparse was unavailable or disabled when liboctave was built");
-#endif
-}
-
-SparseQR::SparseQR_rep::~SparseQR_rep (void)
-{
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_DNAME (_sfree) (S);
-  CXSPARSE_DNAME (_nfree) (N);
-#endif
-}
-
-SparseMatrix
-SparseQR::SparseQR_rep::V (void) const
-{
-#ifdef HAVE_CXSPARSE
-  // Drop zeros from V and sort
-  // FIXME: Is the double transpose to sort necessary?
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  CXSPARSE_DNAME (_dropzeros) (N->L);
-  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
-  CXSPARSE_DNAME (_spfree) (N->L);
-  N->L = CXSPARSE_DNAME (_transpose) (D, 1);
-  CXSPARSE_DNAME (_spfree) (D);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nz = N->L->nzmax;
-  SparseMatrix ret (N->L->m, nc, nz);
-  for (octave_idx_type j = 0; j < nc+1; j++)
-    ret.xcidx (j) = N->L->p[j];
-  for (octave_idx_type j = 0; j < nz; j++)
-    {
-      ret.xridx (j) = N->L->i[j];
-      ret.xdata (j) = N->L->x[j];
-    }
-  return ret;
-#else
-  return SparseMatrix ();
-#endif
-}
-
-ColumnVector
-SparseQR::SparseQR_rep::Pinv (void) const
-{
-#ifdef HAVE_CXSPARSE
-  ColumnVector ret(N->L->m);
-  for (octave_idx_type i = 0; i < N->L->m; i++)
-    ret.xelem (i) = S->pinv[i];
-  return ret;
-#else
-  return ColumnVector ();
-#endif
-}
-
-ColumnVector
-SparseQR::SparseQR_rep::P (void) const
-{
-#ifdef HAVE_CXSPARSE
-  ColumnVector ret(N->L->m);
-  for (octave_idx_type i = 0; i < N->L->m; i++)
-    ret.xelem (S->pinv[i]) = i;
-  return ret;
-#else
-  return ColumnVector ();
-#endif
-}
-
-SparseMatrix
-SparseQR::SparseQR_rep::R (const bool econ) const
-{
-#ifdef HAVE_CXSPARSE
-  // Drop zeros from R and sort
-  // FIXME: Is the double transpose to sort necessary?
-  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-  CXSPARSE_DNAME (_dropzeros) (N->U);
-  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
-  CXSPARSE_DNAME (_spfree) (N->U);
-  N->U = CXSPARSE_DNAME (_transpose) (D, 1);
-  CXSPARSE_DNAME (_spfree) (D);
-  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-  octave_idx_type nc = N->U->n;
-  octave_idx_type nz = N->U->nzmax;
-
-  SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
-
-  for (octave_idx_type j = 0; j < nc+1; j++)
-    ret.xcidx (j) = N->U->p[j];
-  for (octave_idx_type j = 0; j < nz; j++)
-    {
-      ret.xridx (j) = N->U->i[j];
-      ret.xdata (j) = N->U->x[j];
-    }
-  return ret;
-#else
-  return SparseMatrix ();
-#endif
-}
-
-Matrix
-SparseQR::SparseQR_rep::C (const Matrix &b) const
-{
-#ifdef HAVE_CXSPARSE
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nr = nrows;
-  const double *bvec = b.fortran_vec ();
-  Matrix ret (b_nr, b_nc);
-  double *vec = ret.fortran_vec ();
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler) ("matrix dimension mismatch");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    ret = Matrix (nc, b_nc, 0.0);
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
-      for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
-        {
-          octave_quit ();
-          for (octave_idx_type i = nr; i < S->m2; i++)
-            buf[i] = 0.;
-          volatile octave_idx_type nm = (nr < nc ? nr : nc);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (volatile octave_idx_type i = 0; i < nm; i++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          for (octave_idx_type i = 0; i < b_nr; i++)
-            vec[i+idx] = buf[i];
-        }
-    }
-  return ret;
-#else
-  return Matrix ();
-#endif
-}
-
-Matrix
-SparseQR::SparseQR_rep::Q (void) const
-{
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nc = N->L->n;
-  octave_idx_type nr = nrows;
-  Matrix ret (nr, nr);
-  double *vec = ret.fortran_vec ();
-  if (nr < 0 || nc < 0)
-    (*current_liboctave_error_handler) ("matrix dimension mismatch");
-
-  if (nr == 0 || nc == 0)
-    ret = Matrix (nc, nr, 0.0);
-  else
-    {
-      OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
-      for (octave_idx_type i = 0; i < nr; i++)
-        bvec[i] = 0.;
-      OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
-      for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
-        {
-          octave_quit ();
-          bvec[j] = 1.0;
-          for (octave_idx_type i = nr; i < S->m2; i++)
-            buf[i] = 0.;
-          volatile octave_idx_type nm = (nr < nc ? nr : nc);
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (volatile octave_idx_type i = 0; i < nm; i++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          for (octave_idx_type i = 0; i < nr; i++)
-            vec[i+idx] = buf[i];
-          bvec[j] = 0.0;
-        }
-    }
-  return ret.transpose ();
-#else
-  return Matrix ();
-#endif
-}
-
-Matrix
-qrsolve (const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  const double *bvec = b.fortran_vec ();
-  Matrix x;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = Matrix (nc, b_nc, 0.0);
-  else if (nr >= nc)
-    {
-      SparseQR q (a, 3);
-      if (! q.ok ())
-        return Matrix ();
-      x.resize (nc, b_nc);
-      double *vec = x.fortran_vec ();
-      OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
-           i++, idx+=nc, bidx+=b_nr)
-        {
-          octave_quit ();
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseMatrix at = a.hermitian ();
-      SparseQR q (at, 3);
-      if (! q.ok ())
-        return Matrix ();
-      x.resize (nc, b_nc);
-      double *vec = x.fortran_vec ();
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
-      for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
-           i++, idx+=nc, bidx+=b_nr)
-        {
-          octave_quit ();
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-        }
-      info = 0;
-    }
-
-  return x;
-#else
-  return Matrix ();
-#endif
-}
-
-SparseMatrix
-qrsolve (const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-  SparseMatrix x;
-  volatile octave_idx_type ii, x_nz;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = SparseMatrix (nc, b_nc);
-  else if (nr >= nc)
-    {
-      SparseQR q (a, 3);
-      if (! q.ok ())
-        return SparseMatrix ();
-      x = SparseMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              double tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseMatrix at = a.hermitian ();
-      SparseQR q (at, 3);
-      if (! q.ok ())
-        return SparseMatrix ();
-      x = SparseMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            Xx[j] = b.xelem (j,i);
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              double tmp = Xx[j];
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-
-  x.maybe_compress ();
-  return x;
-#else
-  return SparseMatrix ();
-#endif
-}
-
-ComplexMatrix
-qrsolve (const SparseMatrix& a, const MArray<Complex>& b, octave_idx_type& info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nc = b.cols ();
-  octave_idx_type b_nr = b.rows ();
-  ComplexMatrix x;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
-  else if (nr >= nc)
-    {
-      SparseQR q (a, 3);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      Complex *vec = x.fortran_vec ();
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            {
-              Complex c = b.xelem (j,i);
-              Xx[j] = std::real (c);
-              Xz[j] = std::imag (c);
-            }
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (octave_idx_type j = 0; j < nc; j++)
-            vec[j+idx] = Complex (Xx[j], Xz[j]);
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseMatrix at = a.hermitian ();
-      SparseQR q (at, 3);
-      if (! q.ok ())
-        return ComplexMatrix ();
-      x.resize (nc, b_nc);
-      Complex *vec = x.fortran_vec ();
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            {
-              Complex c = b.xelem (j,i);
-              Xx[j] = std::real (c);
-              Xz[j] = std::imag (c);
-            }
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (octave_idx_type j = 0; j < nc; j++)
-            vec[j+idx] = Complex (Xx[j], Xz[j]);
-        }
-      info = 0;
-    }
-
-  return x;
-#else
-  return ComplexMatrix ();
-#endif
-}
-
-SparseComplexMatrix
-qrsolve (const SparseMatrix&a, const SparseComplexMatrix &b,
-         octave_idx_type &info)
-{
-  info = -1;
-#ifdef HAVE_CXSPARSE
-  octave_idx_type nr = a.rows ();
-  octave_idx_type nc = a.cols ();
-  octave_idx_type b_nr = b.rows ();
-  octave_idx_type b_nc = b.cols ();
-  SparseComplexMatrix x;
-  volatile octave_idx_type ii, x_nz;
-
-  if (nr < 0 || nc < 0 || nr != b_nr)
-    (*current_liboctave_error_handler)
-      ("matrix dimension mismatch in solution of minimum norm problem");
-
-  if (nr == 0 || nc == 0 || b_nc == 0)
-    x = SparseComplexMatrix (nc, b_nc);
-  else if (nr >= nc)
-    {
-      SparseQR q (a, 3);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            {
-              Complex c = b.xelem (j,i);
-              Xx[j] = std::real (c);
-              Xz[j] = std::imag (c);
-            }
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (octave_idx_type j = nr; j < q.S ()->m2; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = 0; j < nc; j++)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_usolve) (q.N ()->U, buf);
-          CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Complex (Xx[j], Xz[j]);
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-  else
-    {
-      SparseMatrix at = a.hermitian ();
-      SparseQR q (at, 3);
-      if (! q.ok ())
-        return SparseComplexMatrix ();
-      x = SparseComplexMatrix (nc, b_nc, b.nnz ());
-      x.xcidx (0) = 0;
-      x_nz = b.nnz ();
-      ii = 0;
-      volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2);
-      OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
-      OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
-      for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
-        {
-          octave_quit ();
-          for (octave_idx_type j = 0; j < b_nr; j++)
-            {
-              Complex c = b.xelem (j,i);
-              Xx[j] = std::real (c);
-              Xz[j] = std::imag (c);
-            }
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (octave_idx_type j = nr; j < nbuf; j++)
-            buf[j] = 0.;
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr);
-          CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          for (volatile octave_idx_type j = nr-1; j >= 0; j--)
-            {
-              octave_quit ();
-              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-              CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf);
-              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-            }
-          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-          CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc);
-          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
-
-          for (octave_idx_type j = 0; j < nc; j++)
-            {
-              Complex tmp = Complex (Xx[j], Xz[j]);
-              if (tmp != 0.0)
-                {
-                  if (ii == x_nz)
-                    {
-                      // Resize the sparse matrix
-                      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
-                      sz = (sz > 10 ? sz : 10) + x_nz;
-                      x.change_capacity (sz);
-                      x_nz = sz;
-                    }
-                  x.xdata (ii) = tmp;
-                  x.xridx (ii++) = j;
-                }
-            }
-          x.xcidx (i+1) = ii;
-        }
-      info = 0;
-    }
-
-  x.maybe_compress ();
-  return x;
-#else
-  return SparseComplexMatrix ();
-#endif
-}
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/SparseQR.h
--- a/liboctave/numeric/SparseQR.h	Tue Feb 02 12:33:38 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,169 +0,0 @@
-/*
-
-Copyright (C) 2005-2015 David Bateman
-
-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_SparseQR_h)
-#define octave_SparseQR_h 1
-
-#include <iosfwd>
-
-#include "dMatrix.h"
-#include "CMatrix.h"
-#include "dSparse.h"
-#include "CSparse.h"
-#include "oct-sparse.h"
-
-#if defined (ENABLE_64)
-#define CXSPARSE_DNAME(name) cs_dl ## name
-#else
-#define CXSPARSE_DNAME(name) cs_di ## name
-#endif
-
-class
-OCTAVE_API
-SparseQR
-{
-protected:
-  class SparseQR_rep
-  {
-  public:
-    SparseQR_rep (const SparseMatrix& a, int order);
-
-    ~SparseQR_rep (void);
-#ifdef HAVE_CXSPARSE
-    bool ok (void) const { return (N && S); }
-#else
-    bool ok (void) const { return false; }
-#endif
-    SparseMatrix V (void) const;
-
-    ColumnVector Pinv (void) const;
-
-    ColumnVector P (void) const;
-
-    SparseMatrix R (const bool econ) const;
-
-    Matrix C (const Matrix &b) const;
-
-    Matrix Q (void) const;
-
-    octave_refcount<int> count;
-
-    octave_idx_type nrows;
-#ifdef HAVE_CXSPARSE
-    CXSPARSE_DNAME (s) *S;
-
-    CXSPARSE_DNAME (n) *N;
-#endif
-
-  private:
-
-    // No copying!
-
-    SparseQR_rep (const SparseQR_rep&);
-
-    SparseQR_rep& operator = (const SparseQR_rep&);
-  };
-
-private:
-
-  SparseQR_rep *rep;
-
-public:
-
-  SparseQR (void) : rep (new SparseQR_rep (SparseMatrix (), 0)) { }
-
-  SparseQR (const SparseMatrix& a, int order = 0) :
-    rep (new SparseQR_rep (a, order)) { }
-
-  SparseQR (const SparseQR& a) : rep (a.rep) { rep->count++; }
-
-  ~SparseQR (void)
-  {
-    if (--rep->count == 0)
-      delete rep;
-  }
-
-  SparseQR& operator = (const SparseQR& a)
-  {
-    if (this != &a)
-      {
-        if (--rep->count == 0)
-          delete rep;
-
-        rep = a.rep;
-        rep->count++;
-      }
-    return *this;
-  }
-
-  bool ok (void) const { return rep->ok (); }
-
-  SparseMatrix V (void) const { return rep->V (); }
-
-  ColumnVector Pinv (void) const { return rep->P (); }
-
-  ColumnVector P (void) const { return rep->P (); }
-
-  SparseMatrix R (const bool econ = false) const { return rep->R(econ); }
-
-  Matrix C (const Matrix &b) const { return rep->C(b); }
-
-  Matrix Q (void) const { return rep->Q (); }
-
-  friend Matrix qrsolve (const SparseMatrix &a, const MArray<double> &b,
-                         octave_idx_type &info);
-
-  friend SparseMatrix qrsolve (const SparseMatrix &a, const SparseMatrix &b,
-                               octave_idx_type &info);
-
-  friend ComplexMatrix qrsolve (const SparseMatrix &a, const MArray<Complex> &b,
-                                octave_idx_type &info);
-
-  friend SparseComplexMatrix qrsolve (const SparseMatrix &a,
-                                      const SparseComplexMatrix &b,
-                                      octave_idx_type &info);
-
-protected:
-#ifdef HAVE_CXSPARSE
-  CXSPARSE_DNAME (s) * S (void) { return rep->S; }
-
-  CXSPARSE_DNAME (n) * N (void) { return rep->N; }
-#endif
-};
-
-
-// Publish externally used friend functions.
-
-extern Matrix qrsolve (const SparseMatrix &a, const MArray<double> &b,
-                       octave_idx_type &info);
-
-extern SparseMatrix qrsolve (const SparseMatrix &a, const SparseMatrix &b,
-                             octave_idx_type &info);
-
-extern ComplexMatrix qrsolve (const SparseMatrix &a, const MArray<Complex> &b,
-                              octave_idx_type &info);
-
-extern SparseComplexMatrix qrsolve (const SparseMatrix &a,
-                                    const SparseComplexMatrix &b,
-                                    octave_idx_type &info);
-
-#endif
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Tue Feb 02 12:33:38 2016 -0500
+++ b/liboctave/numeric/module.mk	Tue Feb 02 12:34:06 2016 -0500
@@ -81,8 +81,7 @@
   liboctave/numeric/randpoisson.h \
   liboctave/numeric/sparse-chol.h \
   liboctave/numeric/sparse-lu.h \
-  liboctave/numeric/SparseCmplxQR.h \
-  liboctave/numeric/SparseQR.h
+  liboctave/numeric/sparse-qr.h
 
 NUMERIC_C_SRC = \
   liboctave/numeric/randgamma.c \
@@ -142,11 +141,10 @@
   liboctave/numeric/oct-spparms.cc \
   liboctave/numeric/ODES.cc \
   liboctave/numeric/Quad.cc \
-  liboctave/numeric/SparseCmplxQR.cc \
-  liboctave/numeric/SparseQR.cc \
   liboctave/numeric/sparse-chol-inst.cc \
   liboctave/numeric/sparse-lu-inst.cc \
-$(NUMERIC_C_SRC)
+  liboctave/numeric/sparse-qr-inst.cc \
+  $(NUMERIC_C_SRC)
 
 LIBOCTAVE_TEMPLATE_SRC += \
   liboctave/numeric/base-lu.cc \
@@ -155,6 +153,7 @@
   liboctave/numeric/eigs-base.cc \
   liboctave/numeric/sparse-chol.cc \
   liboctave/numeric/sparse-lu.cc \
+  liboctave/numeric/sparse-qr.cc \
   liboctave/numeric/sparse-dmsolve.cc
 
 ## Special rules for sources which must be built before rest of compilation.
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/sparse-dmsolve.cc
--- a/liboctave/numeric/sparse-dmsolve.cc	Tue Feb 02 12:33:38 2016 -0500
+++ b/liboctave/numeric/sparse-dmsolve.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -28,8 +28,7 @@
 
 #include "MArray.h"
 #include "MSparse.h"
-#include "SparseQR.h"
-#include "SparseCmplxQR.h"
+#include "sparse-qr.h"
 #include "MatrixType.h"
 #include "oct-sort.h"
 #include "oct-locbuf.h"
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/sparse-qr-inst.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-qr-inst.cc	Tue Feb 02 12:34:06 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-qr.h"
+#include "sparse-qr.cc"
+
+template class sparse_qr<SparseMatrix>;
+
+template class sparse_qr<SparseComplexMatrix>;
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/sparse-qr.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-qr.cc	Tue Feb 02 12:34:06 2016 -0500
@@ -0,0 +1,1975 @@
+/*
+
+Copyright (C) 2005-2015 David Bateman
+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 "lo-error.h"
+#include "sparse-qr.h"
+#include "oct-locbuf.h"
+
+template <typename SPARSE_T>
+class
+cxsparse_types
+{
+public:
+  typedef void symbolic_type;
+  typedef void numeric_type;
+};
+
+template <>
+class
+cxsparse_types<SparseMatrix>
+{
+public:
+#ifdef HAVE_CXSPARSE
+  typedef CXSPARSE_DNAME (s) symbolic_type;
+  typedef CXSPARSE_DNAME (n) numeric_type;
+#endif
+};
+
+template <>
+class
+cxsparse_types<SparseComplexMatrix>
+{
+public:
+#ifdef HAVE_CXSPARSE
+  typedef CXSPARSE_ZNAME (s) symbolic_type;
+  typedef CXSPARSE_ZNAME (n) numeric_type;
+#endif
+};
+
+template <typename SPARSE_T>
+class sparse_qr<SPARSE_T>::sparse_qr_rep
+{
+public:
+
+  sparse_qr_rep (const SPARSE_T& a, int order);
+
+  ~sparse_qr_rep (void);
+
+  bool ok (void) const
+  {
+#ifdef HAVE_CXSPARSE
+    return (N && S);
+#else
+    return false;
+#endif
+  }
+
+  SPARSE_T V (void) const;
+
+  ColumnVector Pinv (void) const;
+
+  ColumnVector P (void) const;
+
+  SPARSE_T R (bool econ) const;
+
+  typename SPARSE_T::dense_matrix_type
+  C (const typename SPARSE_T::dense_matrix_type& b) const;
+
+  typename SPARSE_T::dense_matrix_type
+  Q (void) const;
+
+  octave_refcount<int> count;
+
+  octave_idx_type nrows;
+  octave_idx_type ncols;
+
+  typename cxsparse_types<SPARSE_T>::symbolic_type *S;
+  typename cxsparse_types<SPARSE_T>::numeric_type *N;
+
+  template <typename RHS_T, typename RET_T>
+  RET_T
+  tall_solve (const RHS_T& b, octave_idx_type& info) const;
+
+  template <typename RHS_T, typename RET_T>
+  RET_T
+  wide_solve (const RHS_T& b, octave_idx_type& info) const;
+
+private:
+
+  // No copying!
+
+  sparse_qr_rep (const sparse_qr_rep&);
+
+  sparse_qr_rep& operator = (const sparse_qr_rep&);
+};
+
+template <typename SPARSE_T>
+ColumnVector
+sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const
+{
+#ifdef HAVE_CXSPARSE
+  ColumnVector ret (N->L->m);
+  for (octave_idx_type i = 0; i < N->L->m; i++)
+    ret.xelem (i) = S->pinv[i];
+  return ret;
+#else
+  return ColumnVector ();
+#endif
+}
+
+template <typename SPARSE_T>
+ColumnVector
+sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const
+{
+#ifdef HAVE_CXSPARSE
+  ColumnVector ret (N->L->m);
+  for (octave_idx_type i = 0; i < N->L->m; i++)
+    ret.xelem (S->pinv[i]) = i;
+  return ret;
+#else
+  return ColumnVector ();
+#endif
+}
+
+// Specializations.
+
+// Real-valued matrices.
+
+template <>
+sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep
+  (const SparseMatrix& a, int order)
+    : count (1), nrows (a.rows ()), ncols (a.columns ())
+#ifdef HAVE_CXSPARSE
+    , S (0), N (0)
+#endif
+{
+#ifdef HAVE_CXSPARSE
+  CXSPARSE_DNAME () A;
+  A.nzmax = a.nnz ();
+  A.m = nrows;
+  A.n = ncols;
+  // Cast away const on A, with full knowledge that CSparse won't touch it
+  // Prevents the methods below making a copy of the data.
+  A.p = const_cast<octave_idx_type *>(a.cidx ());
+  A.i = const_cast<octave_idx_type *>(a.ridx ());
+  A.x = const_cast<double *>(a.data ());
+  A.nz = -1;
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
+  N = CXSPARSE_DNAME (_qr) (&A, S);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  if (! N)
+    (*current_liboctave_error_handler)
+      ("sparse_qr: sparse matrix QR factorization filled");
+
+  count = 1;
+
+#else
+  (*current_liboctave_error_handler)
+    ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
+#endif
+}
+
+template <>
+sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
+{
+#ifdef HAVE_CXSPARSE
+  CXSPARSE_DNAME (_sfree) (S);
+  CXSPARSE_DNAME (_nfree) (N);
+#endif
+}
+
+template <>
+SparseMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::V (void) const
+{
+#ifdef HAVE_CXSPARSE
+  // Drop zeros from V and sort
+  // FIXME: Is the double transpose to sort necessary?
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  CXSPARSE_DNAME (_dropzeros) (N->L);
+  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
+  CXSPARSE_DNAME (_spfree) (N->L);
+  N->L = CXSPARSE_DNAME (_transpose) (D, 1);
+  CXSPARSE_DNAME (_spfree) (D);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nz = N->L->nzmax;
+  SparseMatrix ret (N->L->m, nc, nz);
+  for (octave_idx_type j = 0; j < nc+1; j++)
+    ret.xcidx (j) = N->L->p[j];
+  for (octave_idx_type j = 0; j < nz; j++)
+    {
+      ret.xridx (j) = N->L->i[j];
+      ret.xdata (j) = N->L->x[j];
+    }
+  return ret;
+#else
+  return SparseMatrix ();
+#endif
+}
+
+template <>
+SparseMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const
+{
+#ifdef HAVE_CXSPARSE
+  // Drop zeros from R and sort
+  // FIXME: Is the double transpose to sort necessary?
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  CXSPARSE_DNAME (_dropzeros) (N->U);
+  CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
+  CXSPARSE_DNAME (_spfree) (N->U);
+  N->U = CXSPARSE_DNAME (_transpose) (D, 1);
+  CXSPARSE_DNAME (_spfree) (D);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  octave_idx_type nc = N->U->n;
+  octave_idx_type nz = N->U->nzmax;
+
+  SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
+
+  for (octave_idx_type j = 0; j < nc+1; j++)
+    ret.xcidx (j) = N->U->p[j];
+  for (octave_idx_type j = 0; j < nz; j++)
+    {
+      ret.xridx (j) = N->U->i[j];
+      ret.xdata (j) = N->U->x[j];
+    }
+  return ret;
+#else
+  return SparseMatrix ();
+#endif
+}
+
+template <>
+Matrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const
+{
+#ifdef HAVE_CXSPARSE
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nr = nrows;
+  const double *bvec = b.fortran_vec ();
+  Matrix ret (b_nr, b_nc);
+  double *vec = ret.fortran_vec ();
+  if (nr < 0 || nc < 0 || nr != b_nr)
+    (*current_liboctave_error_handler) ("matrix dimension mismatch");
+
+  if (nr == 0 || nc == 0 || b_nc == 0)
+    ret = Matrix (nc, b_nc, 0.0);
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+      for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
+        {
+          octave_quit ();
+          for (octave_idx_type i = nr; i < S->m2; i++)
+            buf[i] = 0.;
+          volatile octave_idx_type nm = (nr < nc ? nr : nc);
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          for (volatile octave_idx_type i = 0; i < nm; i++)
+            {
+              octave_quit ();
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+          for (octave_idx_type i = 0; i < b_nr; i++)
+            vec[i+idx] = buf[i];
+        }
+    }
+  return ret;
+#else
+  return Matrix ();
+#endif
+}
+
+template <>
+Matrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::Q (void) const
+{
+#ifdef HAVE_CXSPARSE
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nr = nrows;
+  Matrix ret (nr, nr);
+  double *vec = ret.fortran_vec ();
+  if (nr < 0 || nc < 0)
+    (*current_liboctave_error_handler) ("matrix dimension mismatch");
+
+  if (nr == 0 || nc == 0)
+    ret = Matrix (nc, nr, 0.0);
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
+      for (octave_idx_type i = 0; i < nr; i++)
+        bvec[i] = 0.;
+      OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+      for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
+        {
+          octave_quit ();
+          bvec[j] = 1.0;
+          for (octave_idx_type i = nr; i < S->m2; i++)
+            buf[i] = 0.;
+          volatile octave_idx_type nm = (nr < nc ? nr : nc);
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          for (volatile octave_idx_type i = 0; i < nm; i++)
+            {
+              octave_quit ();
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+          for (octave_idx_type i = 0; i < nr; i++)
+            vec[i+idx] = buf[i];
+          bvec[j] = 0.0;
+        }
+    }
+  return ret.transpose ();
+#else
+  return Matrix ();
+#endif
+}
+
+template <>
+template <>
+Matrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix>
+  (const MArray<double>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  const double *bvec = b.data ();
+
+  Matrix x (nc, b_nc);
+  double *vec = x.fortran_vec ();
+  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
+       i++, idx+=nc, bidx+=b_nr)
+    {
+      octave_quit ();
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return Matrix ();
+
+#endif
+}
+
+template <>
+template <>
+Matrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix>
+  (const MArray<double>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  const double *bvec = b.data ();
+
+  Matrix x (nc, b_nc);
+  double *vec = x.fortran_vec ();
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
+  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
+       i++, idx+=nc, bidx+=b_nr)
+    {
+      octave_quit ();
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, bvec + bidx, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return Matrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix>
+  (const SparseMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return SparseMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseMatrix>
+  (const SparseMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          double tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
+  (const MArray<Complex>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  ComplexMatrix x (nc, b_nc);
+  Complex *vec = x.fortran_vec ();
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        {
+          Complex c = b.xelem (j,i);
+          Xx[j] = std::real (c);
+          Xz[j] = std::imag (c);
+        }
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (octave_idx_type j = 0; j < nc; j++)
+        vec[j+idx] = Complex (Xx[j], Xz[j]);
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix>
+  (const MArray<Complex>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  ComplexMatrix x (nc, b_nc);
+  Complex *vec = x.fortran_vec ();
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        {
+          Complex c = b.xelem (j,i);
+          Xx[j] = std::real (c);
+          Xz[j] = std::imag (c);
+        }
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (octave_idx_type j = 0; j < nc; j++)
+        vec[j+idx] = Complex (Xx[j], Xz[j]);
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+// Complex-valued matrices.
+
+template <>
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep
+  (const SparseComplexMatrix& a, int order)
+    : count (1), nrows (a.rows ()), ncols (a.columns ())
+#ifdef HAVE_CXSPARSE
+    , S (0), N (0)
+#endif
+{
+#ifdef HAVE_CXSPARSE
+  CXSPARSE_ZNAME () A;
+  A.nzmax = a.nnz ();
+  A.m = nrows;
+  A.n = ncols;
+  // Cast away const on A, with full knowledge that CSparse won't touch it
+  // Prevents the methods below making a copy of the data.
+  A.p = const_cast<octave_idx_type *>(a.cidx ());
+  A.i = const_cast<octave_idx_type *>(a.ridx ());
+  A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *> (a.data ()));
+  A.nz = -1;
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
+  N = CXSPARSE_ZNAME (_qr) (&A, S);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  if (! N)
+    (*current_liboctave_error_handler)
+      ("sparse_qr: sparse matrix QR factorization filled");
+
+  count = 1;
+
+#else
+  (*current_liboctave_error_handler)
+    ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
+#endif
+}
+
+template <>
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
+{
+#ifdef HAVE_CXSPARSE
+  CXSPARSE_ZNAME (_sfree) (S);
+  CXSPARSE_ZNAME (_nfree) (N);
+#endif
+}
+
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V (void) const
+{
+#ifdef HAVE_CXSPARSE
+  // Drop zeros from V and sort
+  // FIXME: Is the double transpose to sort necessary?
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  CXSPARSE_ZNAME (_dropzeros) (N->L);
+  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
+  CXSPARSE_ZNAME (_spfree) (N->L);
+  N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
+  CXSPARSE_ZNAME (_spfree) (D);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nz = N->L->nzmax;
+  SparseComplexMatrix ret (N->L->m, nc, nz);
+  for (octave_idx_type j = 0; j < nc+1; j++)
+    ret.xcidx (j) = N->L->p[j];
+  for (octave_idx_type j = 0; j < nz; j++)
+    {
+      ret.xridx (j) = N->L->i[j];
+      ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j];
+    }
+  return ret;
+#else
+  return SparseComplexMatrix ();
+#endif
+}
+
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
+{
+#ifdef HAVE_CXSPARSE
+  // Drop zeros from R and sort
+  // FIXME: Is the double transpose to sort necessary?
+  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+  CXSPARSE_ZNAME (_dropzeros) (N->U);
+  CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
+  CXSPARSE_ZNAME (_spfree) (N->U);
+  N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
+  CXSPARSE_ZNAME (_spfree) (D);
+  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+  octave_idx_type nc = N->U->n;
+  octave_idx_type nz = N->U->nzmax;
+
+  SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
+
+  for (octave_idx_type j = 0; j < nc+1; j++)
+    ret.xcidx (j) = N->U->p[j];
+  for (octave_idx_type j = 0; j < nz; j++)
+    {
+      ret.xridx (j) = N->U->i[j];
+      ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j];
+    }
+  return ret;
+#else
+  return SparseComplexMatrix ();
+#endif
+}
+
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const
+{
+#ifdef HAVE_CXSPARSE
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nr = nrows;
+  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
+  ComplexMatrix ret (b_nr, b_nc);
+  Complex *vec = ret.fortran_vec ();
+  if (nr < 0 || nc < 0 || nr != b_nr)
+    (*current_liboctave_error_handler) ("matrix dimension mismatch");
+
+  if (nr == 0 || nc == 0 || b_nc == 0)
+    ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
+      for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr)
+        {
+          octave_quit ();
+          volatile octave_idx_type nm = (nr < nc ? nr : nc);
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          for (volatile octave_idx_type i = 0; i < nm; i++)
+            {
+              octave_quit ();
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+          for (octave_idx_type i = 0; i < b_nr; i++)
+            vec[i+idx] = buf[i];
+        }
+    }
+  return ret;
+#else
+  return ComplexMatrix ();
+#endif
+}
+
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const
+{
+#ifdef HAVE_CXSPARSE
+  octave_idx_type nc = N->L->n;
+  octave_idx_type nr = nrows;
+  ComplexMatrix ret (nr, nr);
+  Complex *vec = ret.fortran_vec ();
+  if (nr < 0 || nc < 0)
+    (*current_liboctave_error_handler) ("matrix dimension mismatch");
+
+  if (nr == 0 || nc == 0)
+    ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
+  else
+    {
+      OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr);
+      for (octave_idx_type i = 0; i < nr; i++)
+        bvec[i] = cs_complex_t (0.0, 0.0);
+      OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
+      for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
+        {
+          octave_quit ();
+          bvec[j] = cs_complex_t (1.0, 0.0);
+          volatile octave_idx_type nm = (nr < nc ? nr : nc);
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          for (volatile octave_idx_type i = 0; i < nm; i++)
+            {
+              octave_quit ();
+              BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+              CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf));
+              END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+            }
+          for (octave_idx_type i = 0; i < nr; i++)
+            vec[i+idx] = buf[i];
+          bvec[j] = cs_complex_t (0.0, 0.0);
+        }
+    }
+  return ret.hermitian ();
+#else
+  return ComplexMatrix ();
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix>
+  (const SparseComplexMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        {
+          Complex c = b.xelem (j,i);
+          Xx[j] = std::real (c);
+          Xz[j] = std::imag (c);
+        }
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_usolve) (N->U, buf);
+      CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Complex (Xx[j], Xz[j]);
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix>
+  (const SparseComplexMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nr = b.rows ();
+  octave_idx_type b_nc = b.cols ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        {
+          Complex c = b.xelem (j,i);
+          Xx[j] = std::real (c);
+          Xz[j] = std::imag (c);
+        }
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = 0.;
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
+      CXSPARSE_DNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Complex (Xx[j], Xz[j]);
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>, ComplexMatrix>
+  (const MArray<double>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  ComplexMatrix x (nc, b_nc);
+  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_usolve) (N->U, buf);
+      CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>, ComplexMatrix>
+  (const MArray<double>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseComplexMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  ComplexMatrix x (nc, b_nc);
+  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
+
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
+  OCTAVE_LOCAL_BUFFER (double, B, nr);
+  for (octave_idx_type i = 0; i < nr; i++)
+    B[i] = N->B[i];
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      CXSPARSE_ZNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseComplexMatrix>
+  (const SparseMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_usolve) (N->U, buf);
+      CXSPARSE_ZNAME (_ipvec) (S->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseComplexMatrix>
+  (const SparseMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseComplexMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
+
+  OCTAVE_LOCAL_BUFFER (double, B, nr);
+  for (octave_idx_type i = 0; i < nr; i++)
+    B[i] = N->B[i];
+
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      CXSPARSE_ZNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
+  (const MArray<Complex>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
+
+  ComplexMatrix x (nc, b_nc);
+  cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
+    (x.fortran_vec ());
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
+       i++, idx+=nc, bidx+=b_nr)
+    {
+      octave_quit ();
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_usolve) (N->U, buf);
+      CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+ComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix>
+  (const MArray<Complex>& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseComplexMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ());
+
+  ComplexMatrix x (nc, b_nc);
+  cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ());
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
+  OCTAVE_LOCAL_BUFFER (double, B, nr);
+  for (octave_idx_type i = 0; i < nr; i++)
+    B[i] = N->B[i];
+  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
+       i++, idx+=nc, bidx+=b_nr)
+    {
+      octave_quit ();
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr);
+      CXSPARSE_ZNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+    }
+  info = 0;
+
+  return x;
+
+#else
+
+  return ComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix>
+  (const SparseComplexMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  octave_idx_type nr = nrows;
+  octave_idx_type nc = ncols;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < S->m2; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = 0; j < nc; j++)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_usolve) (N->U, buf);
+      CXSPARSE_ZNAME (_ipvec) (S->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <>
+template <>
+SparseComplexMatrix
+sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix>
+  (const SparseComplexMatrix& b, octave_idx_type& info) const
+{
+  info = -1;
+
+#ifdef HAVE_CXSPARSE
+
+  // These are swapped because the original matrix was transposed in
+  // sparse_qr<SparseComplexMatrix>::solve.
+
+  octave_idx_type nr = ncols;
+  octave_idx_type nc = nrows;
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  volatile octave_idx_type ii, x_nz;
+
+  SparseComplexMatrix x (nc, b_nc, b.nnz ());
+  x.xcidx (0) = 0;
+  x_nz = b.nnz ();
+  ii = 0;
+  volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
+  OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
+  OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
+  OCTAVE_LOCAL_BUFFER (double, B, nr);
+  for (octave_idx_type i = 0; i < nr; i++)
+    B[i] = N->B[i];
+  for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
+    {
+      octave_quit ();
+      for (octave_idx_type j = 0; j < b_nr; j++)
+        Xx[j] = b.xelem (j,i);
+      for (octave_idx_type j = nr; j < nbuf; j++)
+        buf[j] = cs_complex_t (0.0, 0.0);
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr);
+      CXSPARSE_ZNAME (_utsolve) (N->U, buf);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      for (volatile octave_idx_type j = nr-1; j >= 0; j--)
+        {
+          octave_quit ();
+          BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+          CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
+          END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+        }
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      CXSPARSE_ZNAME (_pvec) (S->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      for (octave_idx_type j = 0; j < nc; j++)
+        {
+          Complex tmp = Xx[j];
+          if (tmp != 0.0)
+            {
+              if (ii == x_nz)
+                {
+                  // Resize the sparse matrix
+                  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+                  sz = (sz > 10 ? sz : 10) + x_nz;
+                  x.change_capacity (sz);
+                  x_nz = sz;
+                }
+              x.xdata (ii) = tmp;
+              x.xridx (ii++) = j;
+            }
+        }
+      x.xcidx (i+1) = ii;
+    }
+
+  info = 0;
+
+  x.maybe_compress ();
+
+  return x;
+
+#else
+
+  return SparseComplexMatrix ();
+
+#endif
+}
+
+template <typename SPARSE_T>
+sparse_qr<SPARSE_T>::sparse_qr (void)
+  : rep (new sparse_qr_rep (SPARSE_T (), 0))
+{ }
+
+template <typename SPARSE_T>
+sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order)
+  : rep (new sparse_qr_rep (a, order))
+{ }
+
+template <typename SPARSE_T>
+sparse_qr<SPARSE_T>::sparse_qr (const sparse_qr<SPARSE_T>& a)
+  : rep (a.rep)
+{
+  rep->count++;
+}
+
+template <typename SPARSE_T>
+sparse_qr<SPARSE_T>::~sparse_qr (void)
+{
+  if (--rep->count == 0)
+    delete rep;
+}
+
+template <typename SPARSE_T>
+sparse_qr<SPARSE_T>&
+sparse_qr<SPARSE_T>::operator = (const sparse_qr<SPARSE_T>& a)
+{
+  if (this != &a)
+    {
+      if (--rep->count == 0)
+        delete rep;
+
+      rep = a.rep;
+      rep->count++;
+    }
+
+  return *this;
+}
+
+template <typename SPARSE_T>
+bool
+sparse_qr<SPARSE_T>::ok (void) const
+{
+  return rep->ok ();
+}
+
+template <typename SPARSE_T>
+SPARSE_T
+sparse_qr<SPARSE_T>::V (void) const
+{
+  return rep->V ();
+}
+
+template <typename SPARSE_T>
+ColumnVector
+sparse_qr<SPARSE_T>::Pinv (void) const
+{
+  return rep->P ();
+}
+
+template <typename SPARSE_T>
+ColumnVector
+sparse_qr<SPARSE_T>::P (void) const
+{
+  return rep->P ();
+}
+
+template <typename SPARSE_T>
+SPARSE_T
+sparse_qr<SPARSE_T>::R (bool econ) const
+{
+  return rep->R (econ);
+}
+
+template <typename SPARSE_T>
+typename SPARSE_T::dense_matrix_type
+sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b) const
+{
+  return rep->C (b);
+}
+
+template <typename SPARSE_T>
+typename SPARSE_T::dense_matrix_type
+sparse_qr<SPARSE_T>::Q (void) const
+{
+  return rep->Q ();
+}
+
+// FIXME: Why is the "order" of the QR calculation as used in the
+// CXSparse function sqr 3 for real matrices and 2 for complex?  These
+// values seem to be required but there was no explanation in David
+// Bateman's original code.
+
+template <typename SPARSE_T>
+class
+cxsparse_defaults
+{
+public:
+  enum { order = -1 };
+};
+
+template <>
+class
+cxsparse_defaults<SparseMatrix>
+{
+public:
+  enum { order = 3 };
+};
+
+template <>
+class
+cxsparse_defaults<SparseComplexMatrix>
+{
+public:
+  enum { order = 2 };
+};
+
+template <typename SPARSE_T>
+template <typename RHS_T, typename RET_T>
+RET_T
+sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
+                            octave_idx_type& info)
+{
+  info = -1;
+
+  octave_idx_type nr = a.rows ();
+  octave_idx_type nc = a.cols ();
+
+  octave_idx_type b_nc = b.cols ();
+  octave_idx_type b_nr = b.rows ();
+
+  int order = cxsparse_defaults<SPARSE_T>::order;
+
+  if (nr < 0 || nc < 0 || nr != b_nr)
+    (*current_liboctave_error_handler)
+      ("matrix dimension mismatch in solution of minimum norm problem");
+
+  if (nr == 0 || nc == 0 || b_nc == 0)
+    {
+      info = 0;
+
+      return RET_T (nc, b_nc, 0.0);
+    }
+  else if (nr >= nc)
+    {
+      sparse_qr<SPARSE_T> q (a, order);
+
+      return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T ();
+    }
+  else
+    {
+      sparse_qr<SPARSE_T> q (a.hermitian (), order);
+
+      return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
+    }
+}
+
+template <typename SPARSE_T>
+template <typename RHS_T, typename RET_T>
+RET_T
+sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const
+{
+  return rep->tall_solve<RHS_T, RET_T> (b, info);
+}
+
+template <typename SPARSE_T>
+template <typename RHS_T, typename RET_T>
+RET_T
+sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const
+{
+  return rep->wide_solve<RHS_T, RET_T> (b, info);
+}
+
+Matrix
+qrsolve (const SparseMatrix& a, const MArray<double>& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b, info);
+}
+
+SparseMatrix
+qrsolve (const SparseMatrix& a, const SparseMatrix& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b, info);
+}
+
+ComplexMatrix
+qrsolve (const SparseMatrix& a, const MArray<Complex>& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseMatrix>::solve<MArray<Complex>, ComplexMatrix> (a, b, info);
+}
+
+SparseComplexMatrix
+qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix> (a, b, info);
+}
+
+ComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseComplexMatrix>::solve<MArray<double>, ComplexMatrix> (a, b, info);
+}
+
+SparseComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix, SparseComplexMatrix> (a, b, info);
+}
+
+ComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>, ComplexMatrix> (a, b, info);
+}
+
+SparseComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+         octave_idx_type& info)
+{
+  return sparse_qr<SparseComplexMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix> (a, b, info);
+}
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/numeric/sparse-qr.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/sparse-qr.h	Tue Feb 02 12:34:06 2016 -0500
@@ -0,0 +1,122 @@
+/*
+
+Copyright (C) 2005-2015 David Bateman
+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/>.
+
+*/
+
+#if ! defined (octave_sparse_qr_h)
+#define octave_sparse_qr_h 1
+
+#include "dMatrix.h"
+#include "CMatrix.h"
+#include "dSparse.h"
+#include "CSparse.h"
+#include "oct-sparse.h"
+
+template <typename SPARSE_T>
+class
+sparse_qr
+{
+public:
+
+  sparse_qr (void);
+
+  sparse_qr (const SPARSE_T& a, int order = 0);
+
+  sparse_qr (const sparse_qr& a);
+
+  ~sparse_qr (void);
+
+  sparse_qr& operator = (const sparse_qr& a);
+
+  bool ok (void) const;
+
+  SPARSE_T V (void) const;
+
+  ColumnVector Pinv (void) const;
+
+  ColumnVector P (void) const;
+
+  SPARSE_T R (bool econ = false) const;
+
+  typename SPARSE_T::dense_matrix_type
+  C (const typename SPARSE_T::dense_matrix_type& b) const;
+
+  typename SPARSE_T::dense_matrix_type
+  Q (void) const;
+
+  template <typename RHS_T, typename RET_T>
+  static RET_T
+  solve (const SPARSE_T& a, const RHS_T& b,
+         octave_idx_type& info);
+
+private:
+
+  class sparse_qr_rep;
+
+  sparse_qr_rep *rep;
+
+  template <typename RHS_T, typename RET_T>
+  RET_T
+  tall_solve (const RHS_T& b, octave_idx_type& info) const;
+
+  template <typename RHS_T, typename RET_T>
+  RET_T
+  wide_solve (const RHS_T& b, octave_idx_type& info) const;
+};
+
+// Provide qrsolve for backward compatibility.
+
+extern Matrix
+qrsolve (const SparseMatrix& a, const MArray<double>& b,
+         octave_idx_type& info);
+
+extern SparseMatrix
+qrsolve (const SparseMatrix& a, const SparseMatrix& b,
+         octave_idx_type& info);
+
+extern ComplexMatrix
+qrsolve (const SparseMatrix& a, const MArray<Complex>& b,
+         octave_idx_type& info);
+
+extern SparseComplexMatrix
+qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b,
+         octave_idx_type& info);
+
+extern ComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
+         octave_idx_type& info);
+
+extern SparseComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b,
+         octave_idx_type& info);
+
+extern ComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
+         octave_idx_type& info);
+
+extern SparseComplexMatrix
+qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+         octave_idx_type& info);
+
+typedef sparse_qr<SparseMatrix> SparseQR;
+typedef sparse_qr<SparseComplexMatrix> SparseComplexQR;
+
+#endif
diff -r f45f4f888db5 -r 791dcb32b657 liboctave/util/oct-sparse.h
--- a/liboctave/util/oct-sparse.h	Tue Feb 02 12:33:38 2016 -0500
+++ b/liboctave/util/oct-sparse.h	Tue Feb 02 12:34:06 2016 -0500
@@ -107,4 +107,14 @@
 # endif
 #endif
 
+#if defined (HAVE_CXSPARSE)
+# if defined (ENABLE_64)
+#  define CXSPARSE_DNAME(name) cs_dl ## name
+#  define CXSPARSE_ZNAME(name) cs_cl ## name
+# else
+#  define CXSPARSE_DNAME(name) cs_di ## name
+#  define CXSPARSE_ZNAME(name) cs_ci ## name
+# endif
 #endif
+
+#endif